text:
const CKSP = Ptr{Cvoid}
const CKSPType = Cstring

mutable struct KSP{T} <: Factorization{T}
    ptr::CKSP
    comm::MPI.Comm
    # keep around so that they don't get gc'ed
    A  # Operator
    P  # preconditioning operator
    opts::Options{T}
end

scalartype(::KSP{T}) where {T} = T

# allows us to pass XXMat objects directly into CMat ccall signatures
Base.cconvert(::Type{CKSP}, obj::KSP) = obj.ptr
# allows us to pass XXMat objects directly into Ptr{CMat} ccall signatures
Base.unsafe_convert(::Type{Ptr{CKSP}}, obj::KSP) =
    convert(Ptr{CKSP}, pointer_from_objref(obj))

Base.eltype(::KSP{T}) where {T} = T
LinearAlgebra.transpose(ksp) = LinearAlgebra.Transpose(ksp)
LinearAlgebra.adjoint(ksp) = LinearAlgebra.Adjoint(ksp)

@for_libpetsc begin

    function KSP{$PetscScalar}(comm::MPI.Comm; kwargs...)
        initialize($PetscScalar)
        opts = Options{$PetscScalar}(kwargs...)
        ksp = KSP{$PetscScalar}(C_NULL, comm, nothing, nothing, opts)
        @chk ccall((:KSPCreate, $libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CKSP}), comm, ksp)
        if comm == MPI.COMM_SELF
            finalizer(destroy, ksp)
        end
        return ksp
    end

    function destroy(ksp::KSP{$PetscScalar})
        finalized($PetscScalar) ||
        @chk ccall((:KSPDestroy, $libpetsc), PetscErrorCode, (Ptr{CKSP},), ksp)
        return nothing
    end

    function setoperators!(ksp::KSP{$PetscScalar}, A::AbstractMat{$PetscScalar}, P::AbstractMat{$PetscScalar})
        @chk ccall((:KSPSetOperators, $libpetsc), PetscErrorCode, (CKSP, CMat, CMat), ksp, A, P)
        ksp.A = A
        ksp.P = P
        return nothing
    end

    function settolerances!(ksp::KSP{$PetscScalar}; rtol=PETSC_DEFAULT, atol=PETSC_DEFAULT, divtol=PETSC_DEFAULT, max_it=PETSC_DEFAULT)
        @chk ccall((:KSPSetTolerances, $libpetsc), PetscErrorCode,
                    (CKSP, $PetscReal, $PetscReal, $PetscReal, $PetscInt),
                    ksp, rtol, atol, divtol, max_it)
        return nothing
    end

    function setfromoptions!(ksp::KSP{$PetscScalar})
        @chk ccall((:KSPSetFromOptions, $libpetsc), PetscErrorCode, (CKSP,), ksp)
    end

    function gettype(ksp::KSP{$PetscScalar})
        t_r = Ref{CKSPType}()
        @chk ccall((:KSPGetType, $libpetsc), PetscErrorCode, (CKSP, Ptr{CKSPType}), ksp, t_r)
        return unsafe_string(t_r[])
    end

    function iters(ksp::KSP{$PetscScalar})
        r_its = Ref{$PetscInt}()
        @chk ccall((:KSPGetIterationNumber, $libpetsc), PetscErrorCode,
        (KSP, Ptr{$PetscInt}), ksp, r_its)
        return r_its[]
    end

    function view(ksp::KSP{$PetscScalar}, viewer::Viewer{$PetscScalar}=ViewerStdout{$PetscScalar}(ksp.comm))
        @chk ccall((:KSPView, $libpetsc), PetscErrorCode,
                    (CKSP, CPetscViewer),
                ksp, viewer);
        return nothing
    end

    function resnorm(ksp::KSP{$PetscScalar})
        r_rnorm = Ref{$PetscReal}()
        @chk ccall((:KSPGetResidualNorm, $libpetsc), PetscErrorCode,
        (KSP, Ptr{$PetscReal}), ksp, r_rnorm)
        return r_rnorm[]
    end

    function solve!(x::AbstractVec{$PetscScalar}, ksp::KSP{$PetscScalar}, b::AbstractVec{$PetscScalar})
        with(ksp.opts) do
            @chk ccall((:KSPSolve, $libpetsc), PetscErrorCode,
            (CKSP, CVec, CVec), ksp, b, x)
        end
        return x
    end

    function solve!(x::AbstractVec{$PetscScalar}, tksp::Transpose{T,K}, b::AbstractVec{$PetscScalar}) where {T,K <: KSP{$PetscScalar}}
        ksp = parent(tksp)
        with(ksp.opts) do
            @chk ccall((:KSPSolveTranspose, $libpetsc), PetscErrorCode,
            (CKSP, CVec, CVec), ksp, b, x)
        end
        return x
    end

end

# no generic Adjoint solve defined, but for Real we can use Adjoint
solve!(x::AbstractVec{T}, aksp::Adjoint{T,K}, b::AbstractVec{T}) where {K <: KSP{T}} where {T<:Real} =
    solve!(x, transpose(parent(aksp)), b)

const KSPAT{T} = Union{KSP{T}, Transpose{T, KSP{T}}, Adjoint{T, KSP{T}}}

LinearAlgebra.ldiv!(x::AbstractVec{T}, ksp::KSPAT{T}, b::AbstractVec{T}) where {T} = solve!(x, ksp, b)
function LinearAlgebra.ldiv!(x::AbstractVector{T}, ksp::KSPAT{T}, b::AbstractVector{T}) where {T}
    parent(solve!(AbstractVec(x), ksp, AbstractVec(b)))
end
Base.:\(ksp::KSPAT{T}, b::AbstractVector{T}) where {T} = ldiv!(similar(b), ksp, b)

"""
    KSP(A, P; options...)

Construct a PETSc Krylov subspace solver.

Any PETSc options prefixed with `ksp_` and `pc_` can be passed as keywords.
"""
function KSP(A::AbstractMat{T}, P::AbstractMat{T}=A; kwargs...) where {T}
    ksp = KSP{T}(A.comm; kwargs...)
    setoperators!(ksp, A, P)
    with(ksp.opts) do
        setfromoptions!(ksp)
    end
    return ksp
end

Base.show(io::IO, ksp::KSP) = _show(io, ksp)

"""
    iters(ksp::KSP)

Gets the current iteration number; if the `solve!` is complete, returns the number of iterations used.

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPGetIterationNumber.html
"""
iters

"""
    resnorm(ksp::KSP)

Gets the last (approximate preconditioned) residual norm that has been computed.

https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPGetResidualNorm.html
"""
resnorm


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo

# ## Mechanism
mech = get_mechanism(:hopper,
    timestep=0.1,
    gravity=-9.81,
    # friction_coefficient=0.5,
    # damper=25.0,
    # spring=1.0,
    contact_foot=true,
    contact_body=true,
   )

# ## Simulate
Dojo.initialize!(mech, :hopper)

storage = simulate!(mech, 2.5,
    record=true,
    opts=SolverOptions(rtol=1.0e-6, btol=1.0e-6, verbose=true))

# ## Visualize
vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis, show_contact=true)

# ## Contact interpenetration
res = get_sdf(mech, storage) # distance from floor to each contact
minimum(minimum([min.(0.0, r) for r in res]))

ext: jl


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

ArgParse.parse_item(::Type{Vector{String}}, x::AbstractString) = split(x, ",")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--min_nv"
            help = "Minimum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemin(Int)
        "--max_nv"
            help = "Maximum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemax(Int)
        "--serialize_path"
            help = "Where to store the serialize networks."
            default = "data/networks_serialize/"
        "--graphs_path"
            help = "Where to store the graphs."
            default = "data/graphs/"
        "--recompute"
            help = "Wether to restore the serialize `PowerFlowNetwork` object if the path already exists in the database."
            default = false
            action = :store_true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    serialize_instances!(db;
                         serialize_path=parsed_args["serialize_path"],
                         graphs_path=parsed_args["graphs_path"],
                         min_nv=parsed_args["min_nv"],
                         max_nv=parsed_args["max_nv"],
                         recompute=parsed_args["recompute"])
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
   IO functions for the NEID spectrograph
   https://neid.psu.edu/
Author: Eric Ford and collaborators
Created: August 2020
"""

"""Create Dataframe containing filenames and key data for all files neid*.fits in directory"""
function make_manifest(data_path::String ; max_spectra_to_use::Int = 1000 )
    df_filenames = EchelleInstruments.make_manifest(data_path,r"^neidL[1,2]_\d+[T\.]\d+\.fits$")
    #=
    dir_filelist = readdir(data_path,join=true)
    idx_spectra = map(fn->occursin(r"^neidL1_\d+[T\.]\d+\.fits$", last(split(fn,'/')) ),dir_filelist)
    spectra_filelist = dir_filelist[idx_spectra]
    =#
    #=
    df_files = DataFrame(Filename = String[], target = String[], bjd = Float64[], ssbz=Float64[] )
    map(fn->add_metadata_from_fits!(df_files,fn),spectra_filelist)
    df_files
    =#
    #@assert length(spectra_filelist) >= 1
    #df_files = DataFrame(read_metadata(spectra_filelist[1]))
    df_files = DataFrame(read_metadata(df_filenames.Filename[1]))
    keys = propertynames(df_files)
    allowmissing!(df_files, keys[map(k->k∉[:Filename, :bjd, :target, :airmass],keys)] )
    #println("df_files = ", df_files)

    #if length(spectra_filelist) >= 2
    if length(df_filenames.Filename) >= 2
        max_idx = min(max_spectra_to_use,length(df_filenames.Filename))
        map(fn->add_metadata_from_fits!(df_files,fn),df_filenames.Filename[2:max_idx])
    end
    #=
    for i in 2:length(spectra_filelist)
        try
            add_metadata_from_fits!(df_files,spectra_filelist[i])

        catch
            println("# Problem reading ", spectra_filelist[i])
            fields_to_save = metadata_symbols_default(NEID2D())
            fields_str_to_save = metadata_strings_default(NEID2D())

            metadata_tmp = read_metadata_from_fits(spectra_filelist[i],fields=fields_to_save,fields_str=fields_str_to_save)
            println("Before: ", metadata_tmp)
            #println("After : ", metadata_tmp)
            #push!(df_files,metadata_tmp)
        end
    end
    =#
    df_files

end

"""Create Dict containing filename and default metadata from file."""
function read_metadata(fn::String)
    f = FITS(fn)
    read_metadata(f)
end

function read_metadata(f::FITS)
    fields_to_save = metadata_symbols_default(NEID2D())
    fields_str_to_save = metadata_strings_default(NEID2D())
    dict = read_metadata_from_fits(f,fields=fields_to_save,fields_str=fields_str_to_save)
    check_metadata_fields_expected_present(dict,fields_to_save)
    dict[:Filename] = f.filename
    if haskey(dict,:airmass) && typeof(dict[:airmass]) == String
       dict[:airmass] = parse(Float64,dict[:airmass])
    end
    hdr1 = read_header(f[1])
    if hdr1["DATALVL"] == 2
        dict_ccf = read_metadata_ccf(f)
        dict_expmeter = get_exposure_meter_summary(f)
        dict_activity = read_activity_table(f)
        dict = merge(dict,dict_ccf,dict_expmeter,dict_activity)
    end
    return dict
end

function read_metadata_ccf(f::FITS)
    fields_to_save = [ :drp_extsnr, :drp_ccfjdmod, :drp_ccfrvmod, :drp_dvrmsmod]
    fields_str_to_save = [ "EXTSNR", "CCFJDMOD", "CCFRVMOD", "DVRMSMOD"]
    dict = read_metadata_from_fits(f,fields=fields_to_save,fields_str=fields_str_to_save, hdu="CCFS")
end

function read_metadata_ccf(fn::String)
    f = FITS(fn)
    read_metadata_ccf(f)
end

function read_activity_table(fn::String)
    f = FITS(fn)
    read_activity_table(f)
end

"""Create Dict containing filename and default metadata from file."""
function read_activity_table(f::FITS)
    activity_index = read(f["ACTIVITY"],"INDEX")
    activity_value = read(f["ACTIVITY"],"VALUE")
    activity_uncert = read(f["ACTIVITY"],"UNCERTAINTY")
    telluric_zenith = read_header(f["TELLURIC"])["ZENITH"]
    telluric_wvapor = read_header(f["TELLURIC"])["WVAPOR"]
    Dict(vcat(Symbol.(activity_index) .=> read(f["ACTIVITY"],"VALUE"),
              Symbol.("σ_" .* activity_index) .=> read(f["ACTIVITY"],"UNCERTAINTY"),
              Symbol("telluric_zenith") => telluric_zenith,
              Symbol("telluric_wvapor") => telluric_wvapor ) )
end

"""Create Dict containing filename and default metadata from file."""
function read_exposure_meter(f::FITS)
    hdr1 = read_header(f[1])
    exptime = hdr1["EXPTIME"]
    λ_range = 38:97
    fiber = 5
    #expmeter_data = read(f["EXPMETER"],t_range,λ_range,fiber)
    expmeter_data = read(f["EXPMETER"]) # ,:,λ_range,fiber)
    t_len = size(expmeter_data,1)
    t_range = 2:(min(round(Int64,exptime)-1,t_len))
    if t_len < floor(Int64,exptime)
       @warn("*** Exposure meter contains fewer readings than floor(EXPTIME) in " * string(f.filename) )
    end
    view(expmeter_data, t_range,λ_range,fiber)
end

function read_exposure_meter(fn::String)
    f = FITS(fn)
    read_exposure_meter(f)
end

function get_exposure_meter_summary(f::Union{FITS,String})
    try
       expmeter_data = read_exposure_meter(f)
       mean_expmeter = mean(sum(expmeter_data,dims=2),dims=1)[1,1]
       mean_expmeter_blue = mean(sum(view(expmeter_data,:,1:20),dims=2),dims=1)[1,1]
       mean_expmeter_green = mean(sum(view(expmeter_data,:,21:40),dims=2),dims=1)[1,1]
       mean_expmeter_red = mean(sum(view(expmeter_data,:,41:60),dims=2),dims=1)[1,1]
       rms_expmeter = sqrt(var(sum(expmeter_data,dims=2),mean=mean_expmeter,corrected=false))

       expmeter_data_winsor = similar(expmeter_data)
       for i in 1:size(expmeter_data_winsor,2)
           expmeter_data_winsor[:,i] .= winsor(view(expmeter_data,:,i),count=2)
       end
       expmeter_mean_winsor = mean(sum(expmeter_data_winsor,dims=2),dims=1)[1,1]
       expmeter_rms_winsor = sqrt(var(sum(expmeter_data_winsor,dims=2),mean=expmeter_mean_winsor,corrected=false))
       return Dict(:expmeter_mean => mean_expmeter, :expmeter_rms => rms_expmeter, :expmeter_mean_winsor=>expmeter_mean_winsor, :expmeter_rms_winsor=>expmeter_rms_winsor, :expmeter_mean_red=>mean_expmeter_red, :expmeter_mean_green=>mean_expmeter_green, :expmeter_mean_blue=>mean_expmeter_blue )
    catch ex
       if typeof(f) == String
          @warn("*** Error extracting exposure meter data for " * string(f))
       else
          @warn("*** Error extracting exposure meter data for " * string(f.filename))
       end
       return Dict(:expmeter_mean => missing, :expmeter_rms => missing, :expmeter_mean_winsor=>expmeter_mean_winsor, :expmeter_rms_winsor=>expmeter_rms_winsor, :expmeter_mean_red=>missing, :expmeter_mean_green=>missing, :expmeter_mean_blue=>missing )
    end
end

function add_metadata_from_fits!(df::DataFrame, fn::String)
    println("# Reading metadata from ", fn)
    metadata_keep = read_metadata(fn)
    #if metadata_keep[:target] != "Solar" return df end
    calibration_target_substrings = ["Etalon","LFC","ThArS","LDLS","FlatBB"]

   if any(map(ss->contains( metadata_keep[:target],ss),calibration_target_substrings))
      return df
   end

    for k in propertynames(df)
        if typeof(metadata_keep[k]) == Missing continue end
        if typeof(metadata_keep[k]) == Nothing
            metadata_keep[k] = missing
            continue
        end
        if (eltype(df[!,k]) <: Union{Real,Union{<:Real,Missing}}) && !(typeof(metadata_keep[k]) <: Union{Real,Union{<:Real,Missing}})
            metadata_keep[k] = missing
        end
    end
    #println("metadata_keep = ", metadata_keep)
    if haskey(metadata_keep,:airmass) && typeof(metadata_keep[:airmass]) == String
       metadata_keep[:airmass] = parse(Float64,metadata_keep[:airmass])
    end
    push!(df, metadata_keep)
    return df
end

""" Read NEID data from FITS file, and return in a Spectra2DBasic object."""
function read_data   end

function read_data(f::FITS, metadata::Dict{Symbol,Any}; normalization::Symbol = :raw )
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: Union{FITSIO.ImageHDU,FITSIO.TableHDU},2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["SCIWAVE"]]), FITSIO.read(f[img_idx["SCIFLUX"]]), FITSIO.read(f[img_idx["SCIVAR"]])

    if normalization == :raw
        metadata[:normalization] = :raw
    elseif normalization == :blaze
        blaze = read_blaze(f)
        flux ./= blaze
        var ./= blaze.^2
        metadata[:normalization] = :blaze
    else
        @error "# Reading data directly with normalization " * string(normalization) * " is not implemented."
    end

    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR; normalization::Symbol = :raw ) where  { AR<:AbstractRange }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: Union{FITSIO.ImageHDU,FITSIO.TableHDU},2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["SCIWAVE"]],:,orders_to_read), FITSIO.read(f[img_idx["SCIFLUX"]],:,orders_to_read), FITSIO.read(f[img_idx["SCIVAR"]],:,orders_to_read)

    if normalization == :raw
        metadata[:normalization] = :raw
    elseif normalization == :blaze
        blaze = read_blaze(f,orders_to_read)
        flux ./= blaze
        var ./= blaze.^2
        metadata[:normalization] = :blaze
    else
        @error "# Reading data directly with normalization " * string(normalization) * " is not implemented."
    end
    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_data(fn::String, metadata::Dict{Symbol,Any}; normalization::Symbol = :raw  )
    f = FITS(fn)
    read_data(f, metadata, normalization=normalization)
end

function read_data(fn::String; normalization::Symbol = :raw )
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_data(f, metadata, normalization=normalization)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index}; normalization::Symbol = :raw )
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata, normalization=normalization)
end

function read_data(fn::String, orders_to_read::AR; normalization::Symbol = :raw  ) where  { AR<:AbstractRange }
    f = FITS(fn)
    metadata = read_metadata(f)
    read_data(f,metadata, orders_to_read, normalization=normalization)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index}, orders_to_read::AR; normalization::Symbol = :raw  ) where  { AR<:AbstractRange }
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata, orders_to_read, normalization=normalization)
end

function apply_barycentric_correction!(λ::AbstractArray{T,1}, z::Real) where { T<:Real }
    # TODO: WRITE ME
end

""" Read blaze for NEID science fiber from FITS file."""
function read_blaze end

function read_blaze(fn::String)
    f = FITS(fn)
    read_blaze(f)
end

function read_blaze(fn::String, orders_to_read::AR ) where  { AR<:AbstractRange }
    f = FITS(fn)
    read_blaze(f,orders_to_read)
end

function read_blaze(f::FITS )
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: Union{FITSIO.ImageHDU,FITSIO.TableHDU},2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    blaze  = FITSIO.read(f[img_idx["SCIBLAZE"]])
    return blaze
end

function read_blaze(f::FITS, orders_to_read::AR ) where  { AR<:AbstractRange }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: Union{FITSIO.ImageHDU,FITSIO.TableHDU},2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    blaze  = FITSIO.read(f[img_idx["SCIBLAZE"]],:,orders_to_read)
    return blaze
end

""" Read space delimited file with differential extinction corrections, interpolate to bjd's in df and insert into df[:,diff_ext_rv]. """
function read_differential_extinctions!(fn::String, df::DataFrame, df_time_col::Symbol = :bjd)
    df_diff_ext = CSV.read(fn, DataFrame, delim='\t')
    @assert any(isequal(:JD),propertynames(df_diff_ext))
    @assert any(isequal(:delta_vr),propertynames(df_diff_ext))
    diff_ext = LinearInterpolation(df_diff_ext.JD, df_diff_ext.delta_vr)
    df[!,:diff_ext_rv] = -diff_ext(df[!,df_time_col])/100 # cm/s -> m/s
    return df
end

""" Read normalization function for approximate solar continuum+blaze correction
TODO: Figure out mapping between these orders and what's stored in the FITS file before using. """
function read_solar_normalization(fn::String = joinpath(pkgdir(EchelleInstruments),"data","neid","NEID_normalization_function.csv") )
    CSV.read(fn, DataFrame)
end

""" Read calibrator data """
function read_cal_data(f::FITS, metadata::Dict{Symbol,Any}, orders_to_read::AR )  where { AR<:AbstractRange  }
    @assert length(f) >= 2
    @assert all(map(i->typeof(f[i]) <: Union{FITSIO.ImageHDU,FITSIO.TableHDU},2:length(f)))
    img_idx = Dict(map(i->first(read_key(f[i],"EXTNAME"))=>i,2:length(f)))
    λ, flux, var  = FITSIO.read(f[img_idx["CALWAVE"]],:,orders_to_read), FITSIO.read(f[img_idx["CALFLUX"]],:,orders_to_read), FITSIO.read(f[img_idx["CALVAR"]],:,orders_to_read)
    metadata[:normalization] = :raw
    spectrum = Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
    apply_doppler_boost!(spectrum,metadata)
    return spectrum
end

function read_cal_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    read_cal_data(f, metadata)
end

function read_cal_data(fn::String)
    f = FITS(fn)
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    read_cal_data(f, metadata)
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using SpectralKit: to_domain, from_domain

@testset "Chebyshev bounded" begin
    @test_throws DomainError BoundedLinear(-1.0, Inf)
    @test_throws DomainError BoundedLinear(-1.0, -2.0)

    f(x) = exp(-3*x)
    f′(x) = -3*f(x)
    A, B = 1, 5
    N = 11
    grid_dense = range(A, B, length = 100)
    trans = BoundedLinear(A, B)

    for grid_kind in (InteriorGrid(), EndpointGrid())
        basis0 = Chebyshev(grid_kind, N)
        basis = univariate_basis(Chebyshev, grid_kind, N, trans)
        @test dimension(basis) == N
        @test domain(basis) == (A, B)

        @testset "transformation" begin
            for i in 1:100
                x = rand_in_domain(i, A, B)
                @test from_domain(trans, basis0, to_domain(trans, basis0, x)) ≈ x
            end
        end

        g = @inferred grid(basis)
        @test length(g) == N
        C = collocation_matrix(basis, g)
        θ = C \ f.(g)

        if grid_kind ≡ InteriorGrid()
            @test all(A .< g .< B)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-5
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 5e-4
        else
            @test all(A .≤ g .≤ B)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-5
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 5e-4
        end
    end
end

@testset "Chebyshev semi-infinite" begin
    @test_throws DomainError SemiInfRational(-1.0, Inf)
    @test_throws DomainError SemiInfRational(-1.0, 0.0)
    @test_throws DomainError SemiInfRational(NaN, 2.0)

    f(x) = exp(-3*x)
    f′(x) = -3*f(x)
    grid_dense = range(3.1, 8, length = 100)

    A = 3.0
    L = 4.0
    N = 11
    trans = SemiInfRational(A, L)

    @testset "example mappings" begin
        basis0 = Chebyshev(InteriorGrid(), 10)
        @test from_domain(trans, basis0, 0) ≈ A + L
        @test from_domain(trans, basis0, -0.5) ≈ A + L / 3
        @test from_domain(trans, basis0, 0.5) ≈ A + 3 * L
    end

    for grid_kind in (InteriorGrid(), EndpointGrid())
        basis0 = Chebyshev(grid_kind, N)
        basis = univariate_basis(Chebyshev, grid_kind, N, trans)
        @test dimension(basis) == N
        @test domain(basis) == (A, Inf)

        @testset "transformation" begin
            for i in 1:100
                x = rand_in_domain(i, A, Inf)
                @test from_domain(trans, basis0, to_domain(trans, basis0, x)) ≈ x
            end
        end

        g = grid(basis)
        @test length(g) == N
        C = collocation_matrix(basis, g)
        θ = C \ f.(g)

        if grid_kind ≡ InteriorGrid()
            @test sum(g .< A + L) == 5
            @test all(A .< g .< Inf)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-8
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 1e-7
        else

            @test sum(g .< A + L) == 5
            @test all(A .≤ g .≤ Inf)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-8
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 1e-7
        end
    end
end

@testset "Chebyshev infinite" begin
    @test_throws DomainError InfRational(1.0, Inf)
    @test_throws DomainError InfRational(1.0, 0.0)
    @test_throws DomainError InfRational(1.0, -2.0)
    @test_throws DomainError InfRational(NaN, 1)

    f(x) = exp(-4*abs2(x))
    f′(x) = -8*x*f(x)
    grid_dense = range(-1, 2, length = 100)
    A = 0.0
    L = 1.0
    N = 20
    trans = InfRational(A, L)

    @testset "example mappings" begin
        basis0 = Chebyshev(InteriorGrid(), 10)
        @test from_domain(trans, basis0, 0) ≈ A
        @test from_domain(trans, basis0, -0.5) ≈ A - L / √3
        @test from_domain(trans, basis0, 0.5) ≈ A + L / √3
    end

    for grid_kind in (InteriorGrid(), EndpointGrid())
        basis0 = Chebyshev(grid_kind, N)
        basis = univariate_basis(Chebyshev, grid_kind, N, trans)
        @test dimension(basis) == N
        @test domain(basis) == (-Inf, Inf)

        @testset "transformation" begin
            for i in 1:100
                x = rand_in_domain(i, -Inf, Inf)
                @test from_domain(trans, basis0, to_domain(trans, basis0, x)) ≈ x
            end
        end

        g = grid(basis)
        @test length(g) == N
        C = collocation_matrix(basis, g)
        θ = C \ f.(g)

        if grid_kind ≡ InteriorGrid()
            @test sum(g .< A) == 10
            @test all(-Inf .< g .< Inf)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-4
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 2e-3
        else
            @test sum(g .< A) == 10
            @test all(-Inf .≤ g .≤ Inf)
            @test maximum(x -> abs(linear_combination(basis, θ, x) - f(x)),
                          grid_dense) ≤ 1e-4
            @test maximum(x -> abs(derivative(linear_combination(basis, θ), x) - f′(x)),
                          grid_dense) ≤ 2e-3
        end
    end
end


text:
using JuliaBlackBoard
using Documenter

makedocs(;
    modules=[JuliaBlackBoard],
    authors="jverzani <jverzani@gmail.com> and contributors",
    repo="https://github.com/jverzani/JuliaBlackBoard.jl/blob/{commit}{path}#L{line}",
    sitename="JuliaBlackBoard.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jverzani.github.io/JuliaBlackBoard.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jverzani/JuliaBlackBoard.jl",
)


text:
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

import Swifter: QueryResult, ResultInfo, App
import Swifter: path_of_image_from_param, save_image_url_to_file
import Swifter: look_like_view, image_scale

result = QueryResult(ResultInfo(:view,""),App(""),"/query",Dict("lhs"=>[]))
@test true == mimewritable(MIME"text/markdown", result)
@test contains(stringmime("text/markdown", result), "<img src=\"/image")

(path,simple) = path_of_image_from_param(Dict("lhs"=>Any[(:symbol,:vc),(:symbol,:view)]))
@test path == "%5B%5B%22symbol%22%2C%22vc%22%5D%2C%5B%22symbol%22%2C%22view%22%5D%5D"
@test simple == "vc.view"

ENV["SWIFTER_SAVE_IMAGE"] = false
@test "" == save_image_url_to_file("", "", "")

str = "<UILabel: 0x7fe29d447580; frame = (10 100; 500 200); "
(isview, address, xywh) = look_like_view(str)
@test (true, "0x7fe29d447580", [10,100,500,200]) == (isview, address, xywh)
@test "width:250.0px; height: 100.0px;" == image_scale(str)

str = "<UILabel: 0x7fe29d447580; frame = {{10, 100}, {500, 200}}; "
(isview, address, xywh) = look_like_view(str)
@test (true, "0x7fe29d447580", [10,100,500,200]) == (isview, address, xywh)
@test "width:250.0px; height: 100.0px;" == image_scale(str)

str = "string"
(isview, address, xywh) = look_like_view(str)
@test (false, nothing, nothing)== (isview, address, xywh)
@test "" == image_scale(str)


text:
function recurse(f, call, depth)
    f(call, depth)
    depth += 1
    map(call.children) do child
        recurse(f, child, depth)
    end
end

function recurse(f, trace::Callsite)
    recurse(f, trace, 1)
end

"""
    flatten(t::Callsite)

Helper function that flattens a trace generated by `concolic_execution` so that only the leaf nodes are present.
"""
function flatten(t::Callsite)
    stream = Tuple{Any, Any, Tuple}[]
    recurse(t) do call, _
        if isempty(call.children)
            push!(stream, (call.f, call.retval, call.args))
        end
        return
    end
    return stream
end

function filter(t::Callsite)
    stream = Tuple{Any, Any, Tuple}[]
    recurse(t) do call, _
        if isempty(call.children) && any(a->isa(a,Sym), call.args)
            if anything(call.retval) in call.args
                if call.f !== Base.convert
                    @warn """
                    Found method that assigns to an input argument
                    Skipping ($(call.f), $(call.retval), $(call.args))
                    """
                end
                return
            end
            push!(stream, (call.f, call.retval, call.args))
        end
        return
    end
    return stream
end

const FUNCTIONS_TO_IGNORE =(
    assert,
    prove,
)

"""
    verify(trace, merciless=false)

Use `verify(trace, true)` to see built-ins that cause concrete execution
"""
function verify(t::Callsite, merciless=false)
    topmost = t.f
    recurse(t) do call, depth
        if any(a->isa(a, Sym), call.args) && !isa(anything(call.retval), Sym)
            if typeof(call.f) <: Core.Builtin && !merciless
                return
            end
            if call.f == topmost || call.f ∈ FUNCTIONS_TO_IGNORE
                return
            end
            @warn "Function $(call.f) did not propagate taint, $((call.retval, call.args))"
        end
        return
    end
    return nothing
end

function Base.show(io::IO, t::Callsite)
    recurse(t) do call, depth
        base = "  "^depth
        write(io, base)
        print(io, call.retval)
        write(io, " ⇐ ")
        println(io, call.f, call.args)
        return
    end
    return nothing
end

function Base.show(io::IO, s::Sym)
    print(io, s.name)
end


text:
# Tests for distances
function t_affinityprop()

X = [randn(10,200)*2+rand(10,200) rand(10,200).+10 -5.0 .+ randn(10,400)]

Test.@test try
	W = [	j4pr.affinityprop(),
		j4pr.affinityprop(maxiter=200, tol=1e-5,damp=0.7)]
	Xd = j4pr.datacell(X)

	for w in W
		for dtr in [X,Xd], dts in [X,Xd]
			wt = dtr |> w
			dts |> wt
		end
	end
	true
catch
	false
end

end


text:
import CompEcon, BasisMatrices
using PyPlot
using PyCall

f1(x) = @. exp(-2x)
g1(x) = @. -2exp(-2x)

function univariate_interp(f=f1, g=g1)
    # fit approximant
    n = 10; a = -1; b = 1
    basis = CompEcon.fundefn(:cheb, n, a, b)
    c = CompEcon.funfitf(basis, f1)

    # graph approximation error for function and derivative
    x = BasisMatrices.nodeunif(1001, a, b)[1]
    yact = f(x)
    dact = g(x)
    yfit = CompEcon.funeval(c, basis, x)[1]
    dfit = CompEcon.funeval(c, basis, x, 1)[1]

    # Nice plot of function approximation error
    fig, ax = plt.subplots()
    ax[:plot](x, yfit-yact)
    ax[:plot](x, 0*x, "k--", linewidth=2)
    ax[:set_xlabel]("x")
    ax[:set_ylabel]("Error")
    ax[:set_title]("Error approximating function using Chebyshev")

    # Nice plot of derivative approximation error
    fig, ax = plt.subplots()
    ax[:plot](x, dfit-dact);
    ax[:plot](x, 0*x, "k--", linewidth=2)
    ax[:set_xlabel]("x")
    ax[:set_ylabel]("Error")
    ax[:set_title]("Error approximating derivative using Chebyshev")
end

f2(x) = cos.(x[:, 1]) ./ exp.(x[:, 2])
function bivaraite_interp(f=f2)
    n = [7 7]
    a = [0 0]
    b = [1 1]
    basis = CompEcon.fundefn(:cheb, n, a, b)
    c = CompEcon.funfitf(basis, f2)

    fig = plt.figure()
    nplot = [101 101]
    x, xcoord = BasisMatrices.nodeunif(nplot, a, b)
    yfit = CompEcon.funeval(c, basis, x)[1]
    error = reshape(yfit - f2(x), nplot...)
    poly3d = surf(xcoord[1], xcoord[2], error)

end

univariate_interp()
bivaraite_interp()


text:
import DataFrames
import Tables
import JuliaDB
import TypedTables
using Test

# push!(LOAD_PATH, "./src")
using Selections
const s = Selections

function _select_subset(tab, sub)
    t = NamedTuple{s.ensure_tuple(sub)}(Tables.columntable(tab))
    Tables.materializer(tab)(t)
end

function _select_rename(tab, renames)
    from = s.ensure_tuple(first.(renames))
    to = s.ensure_tuple(last.(renames))
    t = values(NamedTuple{from}(Tables.columntable(tab)))
    Tables.materializer(tab)(NamedTuple{to, typeof(t)}(t))
end

_pull(tab, col) = getproperty(tab, col)
_pull(tab::JuliaDB.IndexedTable, col) = JuliaDB.select(tab, col)

const df = DataFrames.DataFrame(a = 1:4, b = 'a':'d', c1 = [[float(i)] for i in 1:4])
const jdb = JuliaDB.table(df)
const tt = TypedTables.FlexTable(jdb)

@testset "Selections" begin
    include("test_select.jl")
end


text:
# Julia FizzBuzz #
function main()
    for i in 1:100
        if i % 3 == 0 && i % 5 == 0
            println("FizzBuzz")
        elseif i % 3 == 0
            println("Fizz")
        elseif i % 5 == 0
            println("Buzz")
        else
            println(i)
        end
    end
end
main()


text:
using PyPlotly
using Documenter

DocMeta.setdocmeta!(PyPlotly, :DocTestSetup, :(using PyPlotly); recursive=true)

makedocs(;
    modules=[PyPlotly],
    authors="SatoshiTerasaki <terasakisatoshi.math@gmail.com> and contributors",
    repo="https://github.com/AtelierArith/PyPlotly.jl/blob/{commit}{path}#{line}",
    sitename="PyPlotly.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://AtelierArith.github.io/PyPlotly.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/AtelierArith/PyPlotly.jl",
    devbranch="main",
)


text:
using TDMSReader
using Test

let a=TDMSReader.readtdms( TDMSReader._example_tdms )
    g=a.groups["Group"]
    c=g.channels["Channel1"]
    @testset "Basic Example File" begin
        @test isempty(a.props)
        @test g.props == Dict("prop" => "value", "num" => 10)
        @test isempty(c.props)
        @test isempty(c.data)
    end

    @testset "Indexing and Keys" begin
        @test Set(["Group"]) == keys(a)
        @test a.groups["Group"] == a["Group"] == a[1]
        @test keys(g.channels) == Set(["Channel1"])
        @test a["Group"] == a[1]
        @test a["Group"]["Channel1"]==a["Group","Channel1"]
        @test a["Group"][1]==a["Group",1]
        @test a[1]["Channel1"]==a[1,"Channel1"]
        @test a[1][1]==a[1,1]
        @test_throws BoundsError a[2]
        @test_throws BoundsError a[1,2]
        @test_throws KeyError a["BadKey"]
        @test_throws KeyError a[1,"BadKey"]
    end

    @testset "Equality Test" begin
        b=deepcopy(a)
        @test a[1,1]==b[1,1]
        @test a[1]==b[1]
        @test a==b

        a.props["NewVal"]="Val"
        b.props["NewVal"]="Val"
        @test a==a
        b.props["NewVal"]="Changed"
        @test a != b

        a[1].props["NewVal"]="Val"
        b[1].props["NewVal"]="Val"
        @test a[1]==a[1]
        b[1].props["NewVal"]="Changed"
        @test a[1] != b[1]

        a[1,1].props["NewVal"]="Val"
        b[1,1].props["NewVal"]="Val"
        @test a[1,1] == b[1,1]
        let b=deepcopy(b)
            b[1,1].props["NewVal"]="Changed"
            @test a[1,1] != b[1,1]
        end
        push!(a[1,1].data,4)
        @test a[1,1] != b[1,1]
        push!(b[1,1].data,4)
        @test a[1,1] == b[1,1]

        b.props["NewVal"]="Val"
        b[1].props["NewVal"]="Val"
        b[1,1].props["NewVal"]="Val"
        @test a[1,1]==b[1,1]
        @test a[1]==b[1]
        @test a==b
    end
end

let fn=TDMSReader._example_incremental
    a=TDMSReader.readtdms(first(fn))
    @testset "Incremental Files" begin
        @test a[1,1].data==[1:3;] && a[1,2].data==[4:6;]

        b=TDMSReader.readtdms(fn[2])
        @test a != b
        append!(a[1,1].data, [1:3;])
        append!(a[1,2].data, [4:6;])
        @test a==b

        b=TDMSReader.readtdms(fn[3])
        @test a != b
        append!(a[1,1].data,[1:3;])
        append!(a[1,2].data, [4:6;])
        @test a[1,1].data == b[1,1].data && a[1,2].data == b[1,2].data
        @test a[1,1].props["prop"] == "valid" && b[1,1].props["prop"] == "error"
        @test a[1,1] != a[1,2]
        a[1,1].props["prop"]="error"
        @test a == b

        #Add new voltage channel
        b=TDMSReader.readtdms(fn[4])
        @test collect(keys(b[1]))==["channel1","channel2","voltage"]
        append!(a[1,1].data,[1:3;])
        append!(a[1,2].data,[4:6;])
        a[1].channels["voltage"]=TDMSReader.Channel{Int}()
        append!(a[1,"voltage"].data, [7:11;])
        @test a == b

        b=TDMSReader.readtdms(fn[5])
        append!(a[1,1].data,[1:3;])
        append!(a[1,2].data,[1:27;])
        append!(a[1,3].data,[7:11;])
        @test a[1,1] == b[1,1]
        @test a[1,2] == b[1,2]
        @test a[1,3] == b[1,3]
        @test a == b

        # Stop appending channel #2
        b=TDMSReader.readtdms(fn[6])
        append!(a[1,1].data,[1:3;])
        append!(a[1,3].data,[7:11;])
        @test a == b

    end
end


text:
"""
$(TYPEDEF)

Newmark's integration scheme with given step-size and parameters `α` and `δ`.

### Fields

- `Δt` -- step-size
- `α`  -- parameter α of the method
- `δ`  -- parameter δ of the method

### References

See [[NEW59]](@ref).
"""
struct Newmark{N} <: AbstractSolver
    Δt::N
    α::N
    δ::N
end

function Newmark(; Δt::T, α::A, δ::D, check=true) where {T, A, D}
    N = promote_type(T, A, D)

    # consistency checks
    if check
        δ ≥ 0.5 || throw(ArgumentError("expected the parameter δ to be at least 1/2, but it is $δ"))
        μ = (1/2 + δ)^2 / 4
        α ≥ μ || throw(ArgumentError("expected the parameter α to be at least $μ, but it is $α"))
    end
    return Newmark(N(Δt), N(α), N(δ))
end

# notable special cases

"""
    Linear(Δt::N)

Linear integration scheme. Special case of Newmark with ``δ=1/2`` and ``α=1/6``.
"""
Linear(Δt::N) where N = Newmark(Δt=Δt, δ=1/2, α=1/6, check=false)
Linear(; Δt::N) where N = Linear(Δt)

"""
    Trapezoidal(Δt::N)

Trapezoidal integration scheme. Special case of Newmark with ``δ=1/2`` and ``α=1/4``.
"""
Trapezoidal(Δt::N) where N = Newmark(Δt=Δt, δ=1/2, α=1/4)
Trapezoidal(; Δt::N) where N = Trapezoidal(Δt)

function _build_solution(alg::Newmark{N}, U, U′, U′′, NSTEPS) where {N}
    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return Solution(alg, U, U′, U′′, t)
end

function _init(alg::Newmark)
    Δt = alg.Δt
    α = alg.α
    δ = alg.δ

    # compute integration constants
    ξ = 1 / (α * Δt)
    a₀ = ξ/Δt
    a₁ = δ*ξ
    a₂ = ξ
    a₃ = 1/(2α) - 1
    a₄ = δ/α - 1
    a₅ = (Δt/2) * (a₄-1)
    a₆ = Δt*(1-δ)
    a₇ = δ*Δt

    return a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇
end

function _solve(alg::Newmark{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int; kwargs...) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)

    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)

    a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇ = _init(alg)
    K̂ = K + a₀ * M + a₁ * C
    K̂⁻¹ = factorize(K̂)

    # initialize displacements, velocities and accelerations
    (U₀, U₀′) = initial_state(ivp)
    U₀′′ = M \ (R[1] - C * U₀′ - K * U₀)

    U = Vector{VT}(undef, IMAX)
    U′ = Vector{VT}(undef, IMAX)
    U′′ = Vector{VT}(undef, IMAX)
    U[1] = U₀
    U′[1] = U₀′
    U′′[1] = U₀′′

    @inbounds for i in 1:NSTEPS
        # calculate effective loads
        mᵢ = M * (a₀ * U[i] + a₂ * U′[i] + a₃ * U′′[i])
        cᵢ = C * (a₁ * U[i] + a₄ * U′[i] + a₅ * U′′[i])
        R̂ᵢ₊₁ = R[i+1] + mᵢ + cᵢ

        # solve for displacements
        U[i+1] = K̂⁻¹ \ R̂ᵢ₊₁

        # calculate accelerations and velocities
        U′′[i+1] = a₀ * (U[i+1] - U[i]) - a₂ * U′[i] - a₃ * U′′[i]
        U′[i+1] = U′[i] + a₆ * U′′[i] + a₇ * U′′[i+1]
    end

    return _build_solution(alg, U, U′, U′′, NSTEPS)
end

# case with possibly non-linear fi(x) term, uses N-R iteration scheme
function _solve(alg::Newmark{N},
                ivp::InitialValueProblem{ST, XT},
                NSTEPS::Int;
                reltol=1e-6,
                maxiter=10,
                t0=zero(N), kwargs...) where {N, VT, ST<:SecondOrderContinuousSystem, XT<:Tuple{VT, VT}}

    Δt = alg.Δt
    sys = system(ivp)
    (U₀, U₀′) = initial_state(ivp)

    IMAX = NSTEPS + 1
    @unpack M, C, fi, fe = sys

    # initial acceleration, obtained from the equilibrium
    # condition Mu'' + Cu' + fint(u) = fext(t), at time t = t0
    U₀′′ = M \ (fe(t0) - C * U₀′ - fi(U₀))

    # integration constants
    a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇ = _init(alg)

    # initialize displacements, velocities and accelerations
    U = Vector{VT}(undef, IMAX)
    U′ = Vector{VT}(undef, IMAX)
    U′′ = Vector{VT}(undef, IMAX)
    U[1] = U₀
    U′[1] = U₀′
    U′′[1] = U₀′′

    # tangent matrix function
    KT(x) = ForwardDiff.jacobian(fi, x)

    # times vector
    tvec = Vector{N}(undef, IMAX)
    tvec[1] = t0

    # time step index
    i = 1

    while i <= NSTEPS
        j = 1
        Ferr = +Inf
        utdt = U[i]

        # Newton-Raphson iteration for equilibrium at t + Δt
        while (j < maxiter) && (Ferr > reltol)

            # calculate effective loads
            mᵢ = M * (a₀ * U[i] + a₂ * U′[i] + a₃ * U′′[i])
            cᵢ = C * (a₁ * U[i] + a₄ * U′[i] + a₅ * U′′[i])
            q = a₀*M + a₁*C
            Feff = fe(tvec[i] + Δt) + mᵢ + cᵢ - fi(utdt) - q * utdt
            Keff = KT(utdt) + q

            # solve for displacements
            du = Keff \ Feff
            utdt += du

            # update stopping criterion
            Ferr = norm(Feff) / norm(fe(tvec[i] + Δt))

            j += 1
        end

        if j == maxiter
            @warn("maximum number of iterations reached")
            break
        end

        U[i+1]   = utdt
        U′′[i+1] = a₀ * (U[i+1] - U[i]) - a₂ * U′[i] - a₃ * U′′[i]
        U′[i+1]  = U′[i] + a₆ * U′′[i] + a₇ * U′′[i+1]

        tvec[i+1] = tvec[i] + Δt
        i += 1
    end

    return _build_solution(alg, U, U′, U′′, NSTEPS)
end

_next!(sampler, Unew) = nothing # ignored

function _solve_statistics(alg::Newmark{N},
                           ivp::InitialValueProblem{ST, XT},
                           NSTEPS::Int,
                           sampler, idx, vecsamples::Vector{Int}) where {N, VT, ST, XT<:Tuple{VT, VT}}

    sys = system(ivp)

    IMAX = NSTEPS + 1
    M, C, K, R = _unwrap(sys, IMAX)
    a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇ = _init(alg)
    K̂ = K + a₀ * M + a₁ * C
    K̂⁻¹ = factorize(K̂)

    # preallocate displacements, velocities and accelerations
    m = size(M, 1)
    U = [VT(undef, m) for _ in 1:IMAX]
    U′ = [VT(undef, m) for _ in 1:IMAX]
    U′′ = [VT(undef, m) for _ in 1:IMAX]

    # output vectors
    U_min = fill(Inf, IMAX)
    U_max = fill(-Inf, IMAX)

    U′_min = fill(Inf, IMAX)
    U′_max = fill(-Inf, IMAX)

    hist = Dict()
    MAXSAMPLES = vecsamples[end]

    runtimes = Dict()
    time0 = time_ns() * 1e-9

    U0new = Vector{N}(undef, 2*m)
    q = 1
    for k in 1:MAXSAMPLES

        # get next state
        _next!(sampler, U0new)
        copyto!(U[1], view(U0new, 1:m))
        copyto!(U′[1], view(U0new, m+1:2m))

        # initialization
        U₀′′ = M \ (R[1] - C * U′[1] - K * U[1])
        U′′[1] = U₀′′

        @inbounds for i in 1:NSTEPS
            # calculate effective loads
            mᵢ = M * (a₀ * U[i] + a₂ * U′[i] + a₃ * U′′[i])
            cᵢ = C * (a₁ * U[i] + a₄ * U′[i] + a₅ * U′′[i])
            R̂ᵢ₊₁ = R[i+1] + mᵢ + cᵢ

            # solve for displacements
            U[i+1] .= K̂⁻¹ \ R̂ᵢ₊₁

            # calculate accelerations and velocities
            U′′[i+1] .= a₀ * (U[i+1] - U[i]) - a₂ * U′[i] - a₃ * U′′[i]
            U′[i+1] .= U′[i] + a₆ * U′′[i] + a₇ * U′′[i+1]
        end

        # store observables
        @inbounds for i in 1:IMAX
            U_min[i] = min(U_min[i], U[i][idx])
            U_max[i] = max(U_max[i], U[i][idx])
            U′_min[i] = min(U′_min[i], U′[i][idx])
            U′_max[i] = max(U′_max[i], U′[i][idx])
        end

        if k == vecsamples[q]
            hist[:($k)] = Dict(:U_min=>copy(U_min), :U_max=>copy(U_max), :Udot_min=>copy(U′_min), :Udot_max=>copy(U′_max))
            runtimes[:($k)] = time_ns() * 1e-9 - time0
            q += 1
        end
    end

    t = range(zero(N), step=alg.Δt, length=(NSTEPS+1))
    return SolutionExtrema(alg, U_min, U_max, U′_min, U′_max, idx, hist, runtimes, t)
end


text:
#=
image_geom.jl
Methods related to an image geometry for image reconstruction
2017-10-02 Samuel Rohrer, University of Michigan
2019-03-05 Jeff Fessler, Julia 1.1 + tests
2019-06-23 Jeff Fessler, overhaul
=#

export MIRT_image_geom, image_geom, cbct
export image_geom_circle
export image_geom_ellipse

#using MIRT: jim, downsample2, downsample3
using Test: @test, @test_throws, @inferred
using ImageTransformations: imresize

"""
`MIRT_image_geom`

Image geometry "struct" with essential parameters
"""
struct MIRT_image_geom
	# options for 2D image geometry
	nx::Int				# image dimension 1
	ny::Int				# image dimension 2
	dx::Float32			# pixel size
	dy::Float32			# pixel size
	offset_x::Float32	# unitless
	offset_y::Float32	# unitless

	# options for 3D image geometry
	nz::Int				# image dimension 3
	dz::Float32			# voxel size
	offset_z::Float32	# unitless

	mask::Array{Bool}	# logical mask
end

"""
`image_geom_help( ; io)`
"""
function image_geom_help( ; io::IO = isinteractive() ? stdout : IOBuffer())
	print(io, "propertynames:\n")
	print(io, propertynames(image_geom(nx=1, dx=1)))

	print(io,
	"
	Derived values for 2D (and 3D) images:

	is3	is it 3d (nz > 0?)
	fovs [|dx|*nx |dy|*ny ...]
	np	sum(mask) = # pixels to be estimated

	dim	(nx, ny, [nz])
	x	1D x coordinates of each pixel
	y	1D y coordiantes of each pixel
	wx	(nx - 1/2) * dx + offset_x
	wy	(ny - 1/2) * dy + offset_y
	wz	(nz - 1/2) * dz + offset_z

	xg	x coordinates of each pixel as a grid (2D or 3D)
	yg	y coordinates of each pixel as a grid (2D or 3D)

	u	1D frequency domain coordinates of each pixel (kx)
	v	1D frequency domain coordinates of each pixel (ky)
	w	1D frequency domain coordinates of each pixel (kz)

	ug	2D or 3D grid of frequency domain coordinates
	vg	2D or 3D grid of frequency domain coordinates
	fg	NamedTuple of 2D or 3D frequency coordinates
	ones	ones(dim)
	zeros	zeros(dim)

	mask_outline	binary image showing 2D outline of mask

	Derived values for 3D images:

	z	z coordinates of each pixel
	zg	z coordinates of each pixel as a grid (3D only)
	wg	3D grid of frequency domain coordinates

	mask_or	[nx ny] logical 'or' of mask

	Methods:

	embed(x)	turn short column(s) into array(s)
	maskit(x)	opposite of embed
	shape(x)	reshape long column x to [nx ny [nz]] array
	unitv(...)	j=j | i=(ix,iy) | c=[cx cy]
				j: single index from 1 to length(z)
				i: [ix,iy[,iz]] index from 1 to nx,ny
				c: [cx,cy[,cz]] index from +/- n/2 center at floor(n/2)+1
	circ(rx=,ry=,cx=,cy=)	circle of given radius and center (cylinder in 3D)
	plot()		plot the image geometry

	Methods that return a new MIRT_image_geom:

	down(down::Int)		down-sample geometry by given factor
	over(over::Int)		over-sample geometry by given factor
	expand_nz(nz_pad)	expand image geometry in z by nz_pad on both ends
	\n")
end

"""
`MIRT_cbct_ig`

Structure suitable for passing to C routines `cbct_*`
based on the struct `cbct_ig` found in `cbct,def.h`
"""
struct MIRT_cbct_ig
	nx::Cint
	ny::Cint
	nz::Cint
	dx::Cfloat
	dy::Cfloat
	dz::Cfloat
	offset_x::Cfloat
	offset_y::Cfloat
	offset_z::Cfloat
	mask2::Ptr{Cuchar}	# [nx,ny] 2D support mask: 0 or 1 .. nthread
	iy_start::Ptr{Cint} # [nthread] for mask2
	iy_end::Ptr{Cint}	# [nthread] for mask2
end

"""
`cbct(ig::MIRT_image_geom; nthread=1)`
constructor for `MIRT_cbct_ig`
"""
function cbct(ig::MIRT_image_geom; nthread::Int=1)
	iy_start = [0]
	iy_end = [ig.ny]
	nthread != 1 && throw("only nthread=1 for now due to iy_start")
	return MIRT_cbct_ig(
		Cint(ig.nx), Cint(ig.ny), Cint(ig.nz), Cfloat(ig.dx),
		Cfloat(ig.dy), Cfloat(ig.dz),
		Cfloat(ig.offset_x), Cfloat(ig.offset_y), Cfloat(ig.offset_z),
			pointer(UInt8.(ig.mask_or)),
			pointer(Cint.(iy_start)), pointer(Cint.(iy_end)),
		)
end

# see https://docs.julialang.org/en/stable/manual/constructors/
"""
`ig = image_geom(...)`

Constructor for `MIRT_image_geom`

# Arguments
- `nx::Int = 128`
- `ny::Int = nx`
- `dx::Real = ?` (must specify one of `dx` or `fov`)
- `dy::Real = -dx`
- `offset_x::Real = 0` (unitless)
- `offset_y::Real = 0` (unitless)
- `fov::Real = ?` (if specified, then `nx*dx=ny*dy`)
- `nz::Int = 0`
- `dz::Real = ?` (need one of `dz` or `zfov` if `nz > 0`)
- `zfov::Real = ?` (if specified, then `nz*dz`)
- `offset_z::Real = 0` (unitless)
- `offsets::Symbol = :none` or :dsp
- `mask::Union{Symbol,AbstractArray{Bool}} = :all` | `:circ` | `:all_but_edge`
"""
function image_geom( ;
		nx::Int = 128,
		ny::Int = nx,
		dx::Real = NaN,
		dy::Real = NaN,
		offset_x::Real = 0,
		offset_y::Real = 0,
		fov::Real = NaN,
		nz::Int = 0,
		dz::Real = NaN,
		zfov::Real = NaN,
		offset_z::Real = 0,
		offsets::Symbol = :none,
		mask::Union{Symbol,AbstractArray{Bool}}	= :all,
	)

	# handle optional arguments

	is3 = !isnan(dz) || !isnan(zfov) || (nz > 0)

	# offsets
	if offsets === :dsp
		# check that no nonzero offset was specified
		(offset_x != 0 || offset_y != 0) && throw("offsets usage incorrect")
		offset_x = 0.5
		offset_y = 0.5
		if is3
			offset_z != 0 && throw("offset_z usage incorrect")
			offset_z = 0.5
		end
	elseif offsets != :none
		throw("offsets $offsets")
	end

	# transverse plane (2d) distances
	if true
		if isnan(fov)
			isnan(dx) && throw("dx or fov required")
		#	fov = nx * dx
		else
			!isnan(dx) && throw("dx and fov?")
			dx = fov / nx
			dy = fov / ny
		end

		dy = isnan(dy) ? -dx : dy # default dy
	end

	# 3D geometry
	if is3
		if isnan(zfov)
			isnan(dz) && throw("dz or zfov required")
		#	zfov = nz * dz
		else
			!isnan(dz) && throw("dz and zfov?")
			dz = zfov / nz
		end
	end

	# mask
	if mask === :all
		mask = is3 ? trues(nx, ny, nz) : trues(nx, ny)
	elseif mask === :circ
		mask = image_geom_circle(nx, ny, dx, dy)
	elseif mask === :all_but_edge
		mask = trues(nx,ny,nz)
		mask[  1,   :, :] .= false
		mask[end,   :, :] .= false
		mask[  :,   1, :] .= false
		mask[  :, end, :] .= false
	elseif isa(mask, Symbol)
		throw("mask symbol $mask")
	elseif size(mask,1) != nx || size(mask,2) != ny || (is3 && size(mask,3) != nz)
		throw("mask size $(size(mask)), nx=$nx ny=$ny nz=$nz")
	end

	# return the object
	return MIRT_image_geom(nx, ny, dx, dy, offset_x, offset_y,
		nz, dz, offset_z, mask)
end

"""
`image_geom_down_round()`
helper function needed to downsample `image_geom`
"""
function _down_round(val::Real=1, dd::Real=1, down::Real=1)
	out = val / down
	if out != round(out)
		out = 2 * round(out / 2) # keep it even if not an integer
	end
	dd = dd * down
	return Int(out), dd
end

"""
`ig_down = downsample(ig, down::Union{Int,Vector{Int}})`
down sample an image geometry by the factor `down`
cf `image_geom_downsample`
"""
function downsample(ig::MIRT_image_geom, down::Union{Int,Vector{Int}})

	# check if it is a scalar by seeing if the collection of size is empty
	if isempty(size(down))
		downv = ig.is3 ? [down,down,down] : downv = [down,down]
	else
		downv = down
	end

	# call the down round function
	down_nx, down_dx = _down_round(ig.nx, ig.dx, downv[1])
	down_ny, down_dy = _down_round(ig.ny, ig.dy, downv[2])
	# adjust to "pixel" units
	down_offset_x = ig.offset_x / downv[1]
	down_offset_y = ig.offset_y / downv[2]

	if ig.is3 # 3d case
		down_nz, down_dz = _down_round(ig.nz, ig.dz, downv[3])
		down_offset_z = ig.offset_z / downv[3]
	else
		down_nz = 0; down_dz = 0; down_offset_z = 0
	end

	# carefully down-sample the mask
	mdim_vec = [size(ig.mask)...] # tuple to vector
	if ig.is3
		if [down_nx, down_ny, down_nz] .* downv == mdim_vec
			down_mask = downsample3(ig.mask, downv) .> 0
		else
			throw("bug: bad mask size. need to address mask downsampling")
		end
	else
		if [down_nx,down_ny] .* downv == mdim_vec
			down_mask = downsample2(ig.mask, downv) .> 0
		else
			down_mask = imresize(ig.mask, (down_nx,down_ny)...) .> 0
		end
	end

	return MIRT_image_geom(
		down_nx, down_ny, down_dx, down_dy, down_offset_x, down_offset_y,
		down_nz, down_dz, down_offset_z, down_mask)
end

"""
`ig_new = image_geom_expand_nz(ig::MIRT_image_geom, nz_pad::Int)`
pad both ends
"""
function image_geom_expand_nz(ig::MIRT_image_geom, nz_pad::Int)
	!ig.is3 && throw("expand_nz only valid for 3D")
	out_nz = ig.nz + 2*nz_pad
	out_mask = cat(dims=3, repeat(ig.mask[:,:,1], 1, 1, nz_pad), ig.mask,
			repeat(ig.mask[:,:,end], 1, 1, nz_pad))
	return MIRT_image_geom(
		ig.nx, ig.ny, ig.dx, ig.dy, ig.offset_x, ig.offset_y,
		out_nz, ig.dz, ig.offset_z, out_mask)
end

"""
`ig_over = image_geom_over(ig::MIRT_image_geom, over::Int)`
over-sample an image geometry by the factor `over`
"""
function image_geom_over(ig::MIRT_image_geom, over::Int)
	if all(ig.mask .== true)
		mask_over = trues(ig.dim)
	else
		mask_over = imresize(ig.mask, ig.dim .* over) .> 0
	end
	return MIRT_image_geom(
		ig.nx*over, ig.ny*over, ig.dx/over, ig.dy/over,
		ig.offset_x*over, ig.offset_y*over,
		ig.nz*over, ig.dz/over, ig.offset_z*over,
		mask_over)
end

# ellipse that just inscribes the rectangle
# but keeping a 1 pixel border due to ASPIRE regularization restriction
function image_geom_ellipse(nx::Int, ny::Int, dx::Real, dy::Real ;
		rx::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
		ry::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
		cx::Real = 0, cy::Real = 0, over::Int=2)
	ig = image_geom(nx=nx, ny=ny, dx=dx, dy=dy)
	circ = ellipse_im(ig, [cx cy rx ry 0 1], oversample=over) .> 0
	return circ
end

# default is a circle that just inscribes the square
# but keeping a 1 pixel border due to ASPIRE regularization restriction
function image_geom_circle(nx::Int, ny::Int, dx::Real, dy::Real ;
		rx::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
		ry::Real = rx, cx::Real = 0, cy::Real = 0, nz::Int=0, over::Int=2)
	ig = image_geom(nx=nx, ny=ny, dx=dx, dy=dy)
	circ = ellipse_im(ig, [cx cy rx ry 0 1], oversample=over) .> 0
	if nz > 0
		circ = repeat(circ, 1, 1, nz)
	end
	return circ
end

"""
`out = image_geom_add_unitv(z::AbstractArray ; j=?, i=?, c=?)`

add a unit vector to an initial array `z` (typically of zeros)

# Arguments
- `j` single index from 1 to length(z)
- `i` [ix,iy[,iz]] index from 1 to nx,ny
- `c` [cx,cy[,cz]] index from +/- n/2 center at floor(n/2)+1

default with no arguments gives unit vector at center `c=[0,0]`
"""
function image_geom_add_unitv(
		z::AbstractArray{T} ; # starts with zeros()
		j::Int = 0,
		i::AbstractVector{Int} = zeros(Int, ndims(z)),
		c::AbstractVector{Int} = zeros(Int, ndims(z))
	) where {T <: Number}

	out = copy(z)

	if j > 0 && all(i .== 0)
		out[j] += one(T)
	elseif j == 0 && any(i .> 0)
		out[i...] += one(T)
	else
		tmp = c .+ Int.(floor.(size(out) ./ 2)) .+ 1
		out[tmp...] += one(T)
	end

	return out
end

"""
`image_geom_add_unitv_test()`
"""
function image_geom_add_unitv_test()
#	test != :test && throw(ArgumentError("test $test"))
	image_geom_add_unitv(zeros(3,4), j=2)
	image_geom_add_unitv(zeros(3,4), i=[2,3])
	image_geom_add_unitv(zeros(3,4))
	image_geom_add_unitv(zeros(3,4), c=[-1,-2])
	true
end

"""
`image_geom_plot(ig)`
"""
function image_geom_plot(ig::MIRT_image_geom; kwargs...)
	return ig.is3 ?
		jim(ig.x, ig.y, ig.mask_or,
			"(dx,dy,dz)=$(ig.dx),$(ig.dy),$(ig.dz)"; kwargs...) :
		jim(ig.x, ig.y, ig.mask, "(nx,ny)=$(ig.nx),$(ig.ny)"; kwargs...)
end

function Base.display(ig::MIRT_image_geom)
	ir_dump(ig)
end

# Extended properties

image_geom_fun0 = Dict([
	(:help, ig -> image_geom_help()),

	(:is3, ig -> ig.nz > 0),
	(:dim, ig -> ig.is3 ? (ig.nx, ig.ny, ig.nz) : (ig.nx, ig.ny)),
	(:fovs, ig -> ig.is3 ?
		[abs(ig.dx)*ig.nx, abs(ig.dy)*ig.ny, abs(ig.dz)*ig.nz] :
			[abs(ig.dx)*ig.nx, abs(ig.dy)*ig.ny]),

	(:zeros, ig -> zeros(Float32, ig.dim)), # (nx, ny, nz) : zeros(nx, ny)
	(:ones, ig -> ones(Float32, ig.dim)), # (nx, ny, nz) : zeros(nx, ny)

	(:x, ig -> ((0:ig.nx-1) .- ig.wx) * ig.dx),
	(:y, ig -> ((0:ig.ny-1) .- ig.wy) * ig.dy),
	(:z, ig -> ((0:ig.nz-1) .- ig.wz) * ig.dz),

	(:wx, ig -> (ig.nx - 1)/2 + ig.offset_x),
	(:wy, ig -> (ig.ny - 1)/2 + ig.offset_y),
	(:wz, ig -> (ig.nz - 1)/2 + ig.offset_z),

	(:xg, ig -> ig.is3 ? repeat(ig.x, 1, ig.ny, ig.nz) :
						repeat(ig.x, 1, ig.ny)),
	(:yg, ig -> ig.is3 ? repeat(ig.y', ig.nx, 1, ig.nz) :
						repeat(ig.y', ig.nx, 1)),
	(:zg, ig -> ig.is3 ? repeat(reshape(ig.z, 1, 1, ig.nz), ig.nx, ig.ny, 1) :
						zeros(ig.nx, ig.ny)),

	# DFT frequency sample grid
	(:u, ig -> (-ig.nx/2:(ig.nx/2-1)) / (ig.nx*ig.dx)),
	(:v, ig -> (-ig.ny/2:(ig.ny/2-1)) / (ig.ny*ig.dy)),
	(:w, ig -> (-ig.nz/2:(ig.nz/2-1)) / (ig.nz*ig.dz)),

	(:ug, ig -> ig.is3 ? repeat(ig.u, 1, ig.ny, ig.nz) :
						repeat(ig.u, 1, ig.ny)),
	(:vg, ig -> ig.is3 ? repeat(ig.v', ig.nx, 1, ig.nz) :
						repeat(ig.v', ig.nx, 1)),
	(:wg, ig -> ig.is3 ? repeat(reshape(ig.w, 1, 1, ig.nz), ig.nx, ig.ny, 1) :
						zeros(ig.nx, ig.ny)),
	(:fg, ig -> ig.is3 ? (ig.ug, ig.vg, ig.wg) : (ig.ug, ig.vg)),

	(:np, ig -> sum(ig.mask)),
	(:mask_or, ig -> mask_or(ig.mask)),
	(:mask_outline, ig -> mask_outline(ig.mask)),

	# simple functions

	(:plot, ig -> ((;kwargs...) -> image_geom_plot(ig; kwargs...))),
	(:embed, ig -> (x::AbstractArray{<:Number} -> embed(x, ig.mask))),
	(:maskit, ig -> (x::AbstractArray{<:Number} -> maskit(x, ig.mask))),
	(:shape, ig -> (x::AbstractArray{<:Number} -> reshape(x, ig.dim))),
	(:unitv, ig -> ((;kwargs...) -> image_geom_add_unitv(ig.zeros; kwargs...))),
	(:circ, ig -> ((;kwargs...) ->
		image_geom_circle(ig.nx, ig.ny, ig.dx, ig.dy, nz=ig.nz; kwargs...))),

	# functions that return new geometry:

	(:down, ig -> (down::Int -> downsample(ig, down))),
	(:over, ig -> (over::Int -> image_geom_over(ig, over))),
	(:expand_nz, ig -> (nz_pad::Int -> image_geom_expand_nz(ig, nz_pad))),

	])

# Tricky overloading here!

Base.getproperty(ig::MIRT_image_geom, s::Symbol) =
		haskey(image_geom_fun0, s) ? image_geom_fun0[s](ig) :
		getfield(ig, s)

Base.propertynames(ig::MIRT_image_geom) =
	(fieldnames(typeof(ig))..., keys(image_geom_fun0)...)

function image_geom_test2(ig::MIRT_image_geom)
	# test 2D functions provided by the constructor
	ig.dim
	ig.x
	ig.y
	ig.wx
	ig.wy
	ig.xg
	ig.yg
	ig.fovs
	ig.np
	ig.mask_outline
	ig.ones
	ig.zeros
	ig.u
	ig.v
	ig.ug
	ig.vg
	ig.fg
	@test ig.shape(ig.ones[:]) == ig.ones
	@test ig.embed(ig.ones[ig.mask]) == Float32.(ig.mask)
	@test ig.maskit(ig.ones) == Float32.(ones(ig.np))
	ig.unitv(j=4)
	ig.unitv(i=ones(Int, length(ig.dim)))
	ig.unitv(c=zeros(Int, length(ig.dim)))
#= todo-i: why do these fail?
	@inferred image_geom_ellipse(8, 10, 1, 2)
	@inferred ig.circ()
	@inferred ig.plot()
	@inferred ig.unitv()
=#
	ig.circ()
	ig.plot()
	ig.unitv()
	@inferred ig.down(2)
	@inferred ig.over(2)
	image_geom_ellipse(8, 10, 1, 2)
	true
end

function image_geom_test2()
	@inferred image_geom(nx=16, dx=2, offsets=:dsp, mask=:all_but_edge)
	@test_throws String image_geom(nx=16, dx=1, offsets=:bad)
	@test_throws String image_geom(nx=16, dx=1, mask=:bad) # mask type
	@test_throws String image_geom(nx=16, dx=1, mask=trues(2,2)) # mask size
	ig = image_geom(nx=16, dx=2)
	display(ig)
	image_geom_test2(ig)
	ig = image_geom(nx=16, dx=2, mask=:circ)
	ig.over(2)
	ig.down(3) # test both even and non-even factors
	ig.help
	true
end

function image_geom_test3(ig::MIRT_image_geom)
	@test image_geom_test2(ig)
	ig.wz
	ig.zg
	ig.mask_or
	@inferred ig.expand_nz(2)
	@inferred cbct(ig)
	true
end

function image_geom_test3()
	@inferred image_geom(nx=16, nz=4, dx=2, zfov=1)
	ig = image_geom(nx=16, nz=4, dx=2, dz=3)
	image_geom_test3(ig)
	true
end

"""
`image_geom(:test)`
self test
"""
function image_geom(test::Symbol)
	if test === :help
		image_geom_help()
		return true
	end
	test != :test && throw(ArgumentError("test $test"))
	@test image_geom_add_unitv_test()
	@test _down_round(4, 3, 2)[1] == 2
	@test image_geom_test2()
	@test image_geom_test3()
	@test image_geom(:help)
	true
end


text:
using ApproxFun, Base.Test

x=Fun(identity);
@test_approx_eq sqrt(cos(π/2*x))(.1) sqrt(cos(.1π/2))

x=Fun(identity,[-2.,2.])
u=sqrt(4-x.^2)/(2π)
@test_approx_eq u(.1) sqrt(4-.1^2)/(2π)
@test_approx_eq sum(u) 1

#this call threw an error, which we check
@test length(values(u))==1

f=Fun(x->x.*cot(π*x/2))
x=Fun(identity)
u=Fun((f./(1-x.^2)).coefficients,JacobiWeight(1.,1.,Interval()))
@test_approx_eq 1./(.1.*cot(π*.1/2)) (1./u)(.1)

@test_approx_eq (x./u)(.1) tan(π*.1/2)

f=Fun(x->exp(-x.^2),Line(0.,0.,-.5,-.5),400)
@test_approx_eq sum(f) sqrt(π)

f=Fun(x->exp(x)/sqrt(1-x.^2),JacobiWeight(-.5,-.5))
@test_approx_eq f(.1) (x->exp(x)/sqrt(1-x.^2))(.1)

S=JacobiWeight(-1.,-1.,Chebyshev([0.,1.]))
D=Derivative(S)

f=Fun(Fun(exp,[0.,1.]).coefficients,S)
x=.1
@test_approx_eq f(x) exp(x)*x^(-1).*(1-x)^(-1)/4
@test_approx_eq (D*f)(x) -exp(x)*(1+(x-3)*x)/(4*(x-1)^2*x^2)

S=JacobiWeight(-1.,0.,Chebyshev([0.,1.]))
D=Derivative(S)

f=Fun(Fun(exp,[0.,1.]).coefficients,S)
x=.1
@test_approx_eq f(x) exp(x)*x^(-1)/2
@test_approx_eq (D*f)(x) exp(x)*(x-1)/(2x^2)

## ODEs

println("    Bessel tests")

for ν in (1.,0.5,2.,3.5)
    println("        ν = $ν")
    S=JacobiWeight(-ν,0.,Chebyshev([0.,1.]))
    D=Derivative(S)
    x=Fun(identity,domain(S))
    L=(x^2)*D^2+x*D+(x^2-ν^2);
    u=linsolve([rdirichlet(S);rneumann(S);L],[bessely(ν,1.),.5*(bessely(ν-1.,1.)-bessely(ν+1.,1.))];
                tolerance=1E-13)
    @test_approx_eq_eps u(.1) bessely(ν,.1) eps(100000.)*max(abs(u(.1)),1)
    u=linsolve([rdirichlet(S),rneumann(S),L],[besselj(ν,1.),.5*(besselj(ν-1.,1.)-besselj(ν+1.,1.))];
                tolerance=1E-13)
    @test_approx_eq_eps u(.1) besselj(ν,.1) eps(100000.)*max(abs(u(.1)),1)

    u=Fun(x->bessely(ν,x),S)
    @test_approx_eq_eps u(.1) bessely(ν,.1) eps(10000.)*max(abs(u(.1)),1)
    u=Fun(x->besselj(ν,x),S)
    @test_approx_eq_eps u(.1) besselj(ν,.1) eps(10000.)*max(abs(u(.1)),1)
end

for ν in (1.,0.5,0.123,3.5)
    println("        ν = $ν")
    S=JacobiWeight(ν,0.,Chebyshev([0.,1.]))
    D=Derivative(S)
    x=Fun(identity,domain(S))
    L=(x^2)*D^2+x*D+(x^2-ν^2);

    u=linsolve([rdirichlet(S),rneumann(S),L],[besselj(ν,1.),.5*(besselj(ν-1.,1.)-besselj(ν+1.,1.))];
                tolerance=1E-10)
    @test_approx_eq_eps u(.1) besselj(ν,.1) eps(1000000.)*max(abs(u(.1)),1)

    u=Fun(x->besselj(ν,x),S)
    @test_approx_eq_eps u(.1) besselj(ν,.1) eps(10000.)*max(abs(u(.1)),1)
end

## f/g bugs

println("    Jacobi singularity tests")

x = Fun(identity)
f = exp(x)./(1-x.^2)

@test_approx_eq f(.1) exp(.1)./(1-.1^2)
f = exp(x)./(1-x.^2).^1
@test_approx_eq f(.1) exp(.1)./(1-.1^2)
f = exp(x)./(1-x.^2).^1.0
@test_approx_eq f(.1) exp(.1)./(1-.1^2)

## 1/f with poles

x=Fun(identity)
f=sin(10x)
g=1/f

@test_approx_eq g(.123) csc(10*.123)

println("    Ray tests")

## Ray

@test Inf in Ray()   # this was a bug

f=Fun(x->exp(-x),[0,Inf])
@test_approx_eq f'(.1) -f(.1)

x=Fun(identity,Ray())
f=exp(-x)
u=integrate(f)
@test_approx_eq (u(1.)-u(0)-1) -f(1)

x=Fun(identity,Ray())
f=x^(-0.123)*exp(-x)
@test_approx_eq integrate(f)'(1.) f(1.)

@test_approx_eq_eps sum(Fun(sech,[0,Inf])) sum(Fun(sech,[0,40.])) 1000000eps()

#Ei (Exp Integral)

y=Fun(Ray())
q=integrate(exp(-y)/y)
@test_approx_eq (q-last(q))(2.) (-0.04890051070806113)

## Line

f=Fun(x->exp(-x^2),Line())

@test_approx_eq f'(0.1) -2*0.1exp(-0.1^2)
@test_approx_eq (Derivative()*f)(0.1) -2*0.1exp(-0.1^2)

## PeriodicLine

d=PeriodicLine()
D=Derivative(d)
f=Fun(x->sech(x-.1),d)

@test_approx_eq_eps (D*f)(.2) -0.0991717226583897  100000eps()
@test_approx_eq_eps (D^2*f)(.2) -0.9752522555114987  1000000eps()

## LogWeight

x=Fun(identity,[-1.,1.])
f=exp(x+1)-1
@test_approx_eq log(f)(0.1) log(f(0.1))

x=Fun(identity,[0.,1.])
f=exp(x)-1
@test_approx_eq log(f)(0.1) log(f(0.1))

## Test divide sing

x=Fun(identity,[0,1])
@test_approx_eq Fun(exp(x)/x-1/x,Chebyshev)(0.1) (exp(0.1)-1)/0.1

x=Fun(identity,[0,1])
f=1/x
p=integrate(f)
@test_approx_eq (p-p(1.))(0.5) log(0.5)

f=1/(1-x)
p=integrate(f)
@test_approx_eq (p-p(0.))(0.5) -log(1-0.5)

y=Fun(Ray())
@test_approx_eq (y^2)(10.) y(10.)^2
@test_approx_eq 1/y(10.) (1/y)(10.)
@test_approx_eq (1/y^2)(10.) 1/y(10.)^2
@test_approx_eq (-1/y^2)'(10.) 2/(y(10.)^3)
@test_approx_eq exp(-1/y^2)(5.) exp(-1/y(5.)^2)

# catch bug from SIE

a=1+10*im;b=2-6*im
d=Curve(Fun(x->1+a*x+b*x^2))

x=Fun(d)
w=sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x))

@test_approx_eq sum(w/(x-2.))/(2π*im) (-4.722196879007759+2.347910413861846im)
@test_approx_eq linesum(w*log(abs(x-2.)))/π (88.5579588360686)

## Dirac Space

a,b=DiracDelta(0.),DiracDelta(1.)
f=Fun(exp)
g=a+0.2b+f
@test_approx_eq pieces(g)[2](0.) 1.
@test_approx_eq g(.1) exp(.1)
@test_approx_eq sum(g) (sum(f)+1.2)

#Checks prevoius bug
δ=DiracDelta()
x=Fun()
w=sqrt(1-x^2)
w+δ

## multiplicities
x=Fun(identity,[-1,1.])
@test_approx_eq (1/x^2)(0.1) 100.
@test_approx_eq (1/x^2)(-0.1) 100.

fc=x*(1+x)^2
@test_approx_eq (1/fc)(0.1) 1/fc(0.1)

fc=x*(1-x)^2
@test_approx_eq (1/fc)(0.1) 1/fc(0.1)

## erf(sqrt(x))

x=Fun([0.,1.])
@test_approx_eq erf(sqrt(x))(0.1) erf(sqrt(0.1))
@test_approx_eq erfc(sqrt(x))(0.1) erfc(sqrt(0.1))

## norm(u-x)

@test_approx_eq norm(Fun(exp,Legendre([0,1]))+sqrt(x)) 2.491141949903508

## Test Jacobi special conversions

S1,S2=JacobiWeight(3.,1.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(1.,0.))
f=Fun([1,2,3.],S1)
C=Conversion(S1,S2)
Cf=C*f
@test_approx_eq Cf(0.1) f(0.1)

S1,S2=JacobiWeight(3.,2.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(0.,0.))
f=Fun([1,2,3.],S1)
C=Conversion(S1,S2)
Cf=C*f
@test_approx_eq Cf(0.1) f(0.1)


text:
YoungsModulus(p::RectilinearPointLoad) = p.E
YoungsModulus(inp::InpLinearElasticity) = inp.inp_content.E
PoissonRatio(p::RectilinearPointLoad) = p.ν
PoissonRatio(inp::InpLinearElasticity) = inp.inp_content.ν


text:
using DashBootstrapComponents, DashHtmlComponents

top_card = dbc_card(
    [
        dbc_cardimg(src="/static/images/placeholder286x180.png", top=true),
        dbc_cardbody(html_p("This card has an image at the top", className="card-text")),
    ],
    style=Dict("width" => "18rem"),
);

bottom_card = dbc_card(
    [
        dbc_cardbody(html_p("This has a bottom image", className="card-text")),
        dbc_cardimg(src="/static/images/placeholder286x180.png", bottom=true),
    ],
    style=Dict("width" => "18rem"),
);

cards = dbc_row([dbc_col(top_card, width="auto"), dbc_col(bottom_card, width="auto")]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Test, LinearAlgebra, Statistics
using StableRNGs
using Lighthouse
using Lighthouse: plot_reliability_calibration_curves, plot_pr_curves,
                  plot_roc_curves, plot_kappas, plot_confusion_matrix,
                  evaluation_metrics_plot
using Base.Threads
using CairoMakie

# Needs to be set for figures
# returning true for showable("image/png", obj)
# which TensorBoardLogger.jl uses to determine output
CairoMakie.activate!(type="png")
plot_results = joinpath(@__DIR__, "plot_results")
# Remove any old plots
isdir(plot_results) && rm(plot_results; force=true, recursive=true)
mkdir(plot_results)

macro testplot(fig_name)
    path = joinpath(plot_results, string(fig_name, ".png"))
    return quote
        fig = $(esc(fig_name))
        @test fig isa Makie.FigureLike
        save($(path), fig)
    end
end

include("plotting.jl")
include("metrics.jl")
include("learn.jl")
include("utilities.jl")
include("logger.jl")


text:
export catalogquery

"""
    read_catalog(catalog)

Read a SCEDC-format listing of catalog data from the Southern California Earthquake Data
Center. For a format description, see http://www.data.scec.org/ftp/catalogs/SCEC_DC.

Catalog maybe a filepath or an `IOStream`.
"""
function read_catalog(catalog)
    cols = [
        "YYYY/MM/DD",
        "HH:mm:SS.ss",
        "ET",
        "GT",
        "MAG",
        "M",
        "LAT",
        "LON",
        "DEPTH",
        "Q",
        "EVID",
        "NPH",
        "NGRM",
    ]
    types = [
        String,
        String,
        String,
        String,
        Float64,
        String,
        Float64,
        Float64,
        Float64,
        String,
        Int,
        Int,
        Int,
    ]

    # read into dataframe
    df = CSV.File(
        catalog,
        header=cols,
        datarow=11,
        types=types,
        comment="#",
        delim=' ',
        ignorerepeated=true,
        silencewarnings=true,
    ) |> DataFrame
    dropmissing!(df)

    # format dates
    dateformater = Dates.DateFormat("y/m/dTH:M:S.s")
    df[!,"EVENTTIME"] = DateTime.(parse_hms.(df[!,"YYYY/MM/DD"],df[!,"HH:mm:SS.ss"]),dateformater)
    select!(df,Not([Symbol("YYYY/MM/DD"),Symbol("HH:mm:SS.ss")]))
    return df
end

"""
    catalogquery()

Use S3 to query SCEDC-pds event catalog. Returns a DataFrame of events.

# Keywords
- `starttime::TimeType`: The start time/day of the event query.
- `endtime::TimeType`: The end time/day of the event query.
- `minlatitude::Real`: Minimum latitude in event query.
- `maxlatitude::Real`: Maximum latitude in event query.
- `minlongitude::Real`: Minimum longitude in event query.
- `maxlongitude::Real`: Maximum longitude in event query.
- `minmagnitude::Real`: Minimum magnitude in event query.
- `maxmagnitude::Real`: Maximum magnitude in event query.
- `mindepth::Real`: Minimum event depth in km.
- `maxdepth::Real`: Maximum event depth in km.
- `eventtype::String`: Event type: {"eq"=>"earthquake", "qb"=> "quarry blast",
    "sn"=>"sonic boom', "nt"=>"nuclear blast", "uk"=>"not reported","*"=>"all events"}
"""
function catalogquery(aws::AWSConfig;
    starttime::TimeType=today(),
    endtime::TimeType=today(),
    minlatitude::Real=-90,
    maxlatitude::Real=90,
    minlongitude::Real=-180,
    maxlongitude::Real=180,
    minmagnitude::Real=-10,
    maxmagnitude::Real=10,
    mindepth::Real=-10,
    maxdepth::Real=100,
    eventtype::String="*",
)
    all_types = ["eq","qb","sn","nt","uk","*"]
    # do filein
    @assert starttime <= endtime "Starttime ($starttime) must be less than endtime ($endtime)."
    @assert minlatitude >= -90 "Minimum latitude must be greater than -90"
    @assert maxlatitude <= 90 "Maximum latitude must be less than 180"
    @assert minlongitude >= -180 "Minimum longitude must be greater than -180"
    @assert maxlongitude <= 180 "Maximum longitude must be less than 180"
    @assert minmagnitude >= -10 "Minimum magnitude must be greater than -10"
    @assert maxmagnitude <= 10 "Maximum magnitude must be less than 10"
    @assert in(lowercase(eventtype), all_types) "eventtype '$eventtype' must be one of $all_types"

    # get each year
    years = year(starttime):year(endtime)
    alldf = DataFrame()
    for yr in years
        yrstream = s3_get(aws,"scedc-pds","earthquake_catalogs/SCEC_DC/$yr.catalog")
        df = read_catalog(yrstream)
        filter!(row -> starttime <= row.EVENTTIME <= endtime, df)
        filter!(row -> minlatitude <= row.LAT <= maxlatitude,df)
        filter!(row -> minlongitude <= row.LON <= maxlongitude,df)
        filter!(row -> minmagnitude <= row.MAG <= maxmagnitude,df)
        filter!(row -> mindepth <= row.DEPTH <= maxdepth,df)
        if eventtype != "*"
            filter!(row -> row.ET == eventtype, df)
        end
        if !isempty(df)
            append!(alldf,df)
        end
    end

    if isempty(alldf)
        throw(DomainError("No events for query."))
    end
    return alldf
end
catalogquery(a...;b...) = catalogquery(global_aws_config(region="us-west-2"), a...; b...)
catalogquery(d::Dict, a...;b...) = catalogquery(global_aws_config(region=d[:region]), a...; b...)

function parse_hms(d::String,hms::String)
    # check that date matches HH:MM:SS.s syntax
    r = r"^[0-2][0-3]:[0-5][0-9]:[0:5][0-9].[0-9][0-9]"
    if occursin(r,hms)
        return d * "T" * hms
    end

    # string does not match
    strm=match(r"(?<hour>\d+):(?<minute>\d+):(?<second>\d+).(?<dec>\d+)",hms)
    dec = parse(Int,strm[:dec])
    h = parse(Int,strm[:hour])
    m = parse(Int,strm[:minute])
    s = parse(Int,strm[:second])

    if s > 59
        m +=  s ÷ 60
        s = s % 60
    end

    if m > 59
        h += m ÷ 60
        m = m % 60
    end

    if h > 23
        newday = h ÷ 24
        h = h % 24
        dateformater = DateFormat("yyyy/mm/dd")
        d = Dates.format(Date(d,dateformater)+ Day(newday),dateformater)
    end

    h = lpad(string(h),2,'0')
    m = lpad(string(m),2,'0')
    s = lpad(string(s),2,'0')
    return d * "T$h:$m:$s.$dec"
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export checkbox

"""
    checkbox(;kwargs...)

A Checkbox component.
Capture boolean input from user. For more information, see: https://mantine.dev/core/checkbox/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `checked` (Bool; optional): State of check box
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Checkbox color
- `disabled` (Bool; optional): A checkbox can show it is currently unable to be interacted with
- `label` (String; optional): Checkbox label
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined label font-size and checkbox width and height in px
- `style` (Dict; optional): Inline style override
- `transitionDuration` (Real; optional): Check/uncheck transition duration, set to 0 to disable all transitions
"""
function checkbox(; kwargs...)
        available_props = Symbol[:id, :checked, :className, :color, :disabled, :label, :size, :style, :transitionDuration]
        wild_props = Symbol[]
        return Component("checkbox", "Checkbox", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
# Completed: 69/69

import Base: acos, acosd, acosh, acot, acotd, acoth, acsc, acscd, acsch, asec, asecd, asech, asin, asind, asinh, atan, atand, atanh, cbrt, clamp, cos, cosc, cosd, cosh, cospi, cot, cotd, coth, csc, cscd, csch, deg2rad, exp, exp10, exp2, expm1, exponent, hypot, ldexp, log, log10, log1p, log2, max, min, mod2pi, rad2deg, rem2pi, sec, secd, sech, significand, sin, sinc, sind, sinh, sinpi, sqrt, tan, tand, tanh

# julia> for n in names(Base.Math); println(n); end
# @evalpoly(z,c...): Evaluate the polynomial \sum_k c[k] z^{k-1}.
# Math: package name
# ^: Handled in base.jl
@primitive acos(x),dy (dy.*(-1 ./ sqrt.(1 .- abs2.(x))))
@primitive acosd(x),dy (dy.*(-(180 / pi) ./ sqrt.(1 .- abs2.(x))))
@primitive acosh(x),dy (dy.*(1 ./ sqrt.(abs2.(x) .- 1)))
@primitive acot(x),dy (dy.*(-1 ./ (1 .+ abs2.(x))))
@primitive acotd(x),dy (dy.*(-(180 / pi) ./ (1 .+ abs2.(x))))
@primitive acoth(x),dy (dy.*(1 ./ (1 .- abs2.(x))))
@primitive acsc(x),dy (dy.*(-1 ./ sqrt.((abs2.(x) .* (x .- 1)) .* (x .+ 1))))
@primitive acscd(x),dy (dy.*(-(180 / pi) ./ sqrt.((abs2.(x) .* (x .- 1)) .* (x .+ 1))))
@primitive acsch(x),dy (dy.*(-1 ./ sqrt.(abs2.(abs2.(x)) + abs2.(x))))
@primitive asec(x),dy (dy.*(1 ./ sqrt.(abs2.(abs2.(x)) - abs2.(x))))
@primitive asecd(x),dy (dy.*((180 / pi) ./ sqrt.(abs2.(abs2.(x)) - abs2.(x))))
@primitive asech(x),dy (dy.*(-1 ./ sqrt.(abs2.(x) - abs2.(abs2.(x)))))
@primitive asin(x),dy (dy.*(1 ./ sqrt.(1 .- abs2.(x))))
@primitive asind(x),dy (dy.*((180 / pi) ./ sqrt.(1 .- abs2.(x))))
@primitive asinh(x),dy (dy.*(1 ./ sqrt.(1 .+ abs2.(x))))
@primitive atan(x),dy (dy.*(1 ./ (1 .+ abs2.(x)))); @primitive atan(x1,x2),dy,y  unbroadcast(x1,dy.*(x2./(abs2.(x1)+abs2.(x2))))  unbroadcast(x2,dy.*(-x1./(abs2.(x1)+abs2.(x2))))
@primitive atand(x),dy (dy.*((180 / pi) ./ (1 .+ abs2.(x))))
@primitive atanh(x),dy (dy.*(1 ./ (1 .- abs2.(x))))
@primitive cbrt(x),dy,y (dy.*(1 ./ (3 .* abs2.(y))))
@primitive clamp(x,lo,hi),dy,y  unbroadcast(x,dy.*(lo .<= x .<= hi))
# clamp! overwrites.
@primitive cos(x),dy (dy.*(-(sin.(x))))
@primitive cosc(x),dy,y (dy.*((-2 .* y) ./ x - sinc.(x) .* (pi ^ 2)))
@primitive cosd(x),dy (dy.*((-(sind.(x)) .* pi) / 180))
@primitive cosh(x),dy (dy.*(sinh.(x)))
@primitive cospi(x),dy (dy.*(-(sinpi.(x)) .* pi))
@primitive cot(x),dy (dy.*(-(abs2.(csc.(x)))))
@primitive cotd(x),dy (dy.*((-(abs2.(cscd.(x))) .* pi) / 180))
@primitive coth(x),dy (dy.*(-(abs2.(csch.(x)))))
@primitive csc(x),dy,y (dy.*(-1 .* y .* cot.(x)))
@primitive cscd(x),dy,y (dy.*(((-1 .* y .* cotd.(x)) .* pi) / 180))
@primitive csch(x),dy,y (dy.*(-1 .* y .* coth.(x)))
@primitive deg2rad(x),dy (dy.*(pi / 180))
@primitive exp(x),dy,y (dy.*(y))
@primitive exp10(x),dy,y (dy.*(y .* log(10)))
@primitive exp2(x),dy,y (dy.*(y .* log(2)))
@primitive expm1(x),dy,y (dy.*(1 .+ y))
@zerograd exponent(x)
# frexp # returns 2 values
@primitive hypot(x1,x2),dy,y  unbroadcast(x1,dy.*(x1./y))  unbroadcast(x2,dy.*(x2./y))
@primitive ldexp(x,n),dy  (dy.*exp2.(n))
@primitive log(x),dy (dy.*(1 ./ x)); @primitive log(x1,x2),dy,y  unbroadcast(x1,dy.*(-log.(x2)./(x1.*abs2.(log.(x1)))))  unbroadcast(x2,dy.*(1 ./ (x2.*log.(x1))))
@primitive log10(x),dy (dy.*(1 ./ (log(10) .* x)))
@primitive log1p(x),dy (dy.*(1 ./ (1 .+ x)))
@primitive log2(x),dy (dy.*(1 ./ (log(2) .* x)))
@primitive max(x1,x2),dy,y  unbroadcast(x1,dy.*(y.==x1))  unbroadcast(x2,dy.*(y.==x2))
@primitive min(x1,x2),dy,y  unbroadcast(x1,dy.*(y.==x1))  unbroadcast(x2,dy.*(y.==x2))
# minmax # returns 2 values
@primitive mod2pi(x),dy dy
# modf # returns 2 values
@primitive rad2deg(x),dy (dy.*(180 / pi))
@primitive rem2pi(x,r),dy  dy
@primitive sec(x),dy,y (dy.*(y .* tan.(x)))
@primitive secd(x),dy,y (dy.*(((y .* tand.(x)) .* pi) / 180))
@primitive sech(x),dy,y (dy.*(-1 .* y .* tanh.(x)))
@primitive significand(x),dy (dy.*exp2.(-exponent.(x)))
@primitive sin(x),dy (dy.*(cos.(x)))
@primitive sinc(x),dy (dy.*(cosc.(x)))
# sincos # returns 2 values
@primitive sind(x),dy (dy.*((cosd.(x) .* pi) / 180))
@primitive sinh(x),dy (dy.*(cosh.(x)))
@primitive sinpi(x),dy (dy.*(cospi.(x) .* pi))
@primitive sqrt(x),dy,y (dy.*(1 ./ (2 .* y)))
@primitive tan(x),dy,y (dy.*(1 .+ abs2.(y)))
@primitive tand(x),dy,y (dy.*(((1 .+ abs2.(y)) .* pi) / 180))
@primitive tanh(x),dy,y (dy.*(1 .- abs2.(y)))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SassC_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SassC")
JLLWrappers.@generate_main_file("SassC", UUID("e813748c-f071-5ac4-9ddd-9e9e72545407"))
end  # module SassC_jll


text:
###############################################################################
#
# Set X and Y ticks for the axis
#
###############################################################################
"""
    set_xticks!(axs::Array, xticks::Array)
    set_xticks!(axs::Array,
                xticks::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
    )

Set X-axis ticks for a given
- `axs` Array of axis
- `xticks` Array of X-axis ticks
"""
function set_xticks!(axs::Array, xticks::Array)
    for indx in eachindex(axs)
        axs[indx].set_xticks(xticks[indx]);
    end

    return nothing
end

function set_xticks!(
            axs::Array,
            xticks::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
)
    for indx in eachindex(axs)
        axs[indx].set_xticks(xticks);
    end

    return nothing
end

"""
    set_yticks!(axs::Array, yticks::Array)
    set_yticks!(axs::Array,
                yticks::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
    )

Set Y-axis ticks for a given
- `axs` Array of axis
- `yticks` Array of Y-axis ticks
"""
function set_yticks!(axs::Array, yticks::Array)
    for indx in eachindex(axs)
        axs[indx].set_yticks(yticks[indx]);
    end

    return nothing
end

function set_yticks!(
            axs::Array,
            yticks::Union{Array{Int,1},Array{Float32,1},Array{Float64,1}}
)
    for indx in eachindex(axs)
        axs[indx].set_yticks(yticks);
    end

    return nothing
end

"""
    set_xyticks!(axs::Array, xticks::Array, yticks::Array)

Set Y-axis ticks for a given
- `axs` Array of axis
- `xticks` Array of X-axis ticks
- `yticks` Array of Y-axis ticks
"""
function set_xyticks!(axs::Array, xticks::Array, yticks::Array)
    set_xticks!(axs, xticks);
    set_yticks!(axs, yticks);

    return nothing
end

###############################################################################
#
# Set ticklabels for the axis
#
###############################################################################
"""
    set_xticklabels!(axs::Array, xticklabs::Array)
    set_xticklabels!(axs::Array, xticklabs::Array{String,1})

Set X-axis ticks for a given
- `axs` Array of axis
- `xticklabs` Array of X-axis tick labels
"""
function set_xticklabels!(axs::Array, xticklabs::Array)
    for indx in eachindex(axs)
        axs[indx].set_xticklabels(xticklabs[indx]);
    end

    return nothing
end

function set_xticklabels!(axs::Array, xticklabs::Array{String,1})
    for indx in eachindex(axs)
        axs[indx].set_xticklabels(xticklabs);
    end

    return nothing
end

"""
    set_yticklabels!(axs::Array, yticklabs::Array)
    set_yticklabels!(axs::Array, yticklabs::Array{String,1})

Set Y-axis ticks for a given
- `axs` Array of axis
- `yticklabs` Array of Y-axis tick labels
"""
function set_yticklabels!(axs::Array, yticklabs::Array)
    for indx in eachindex(axs)
        axs[indx].set_yticklabels(yticklabs[indx]);
    end

    return nothing
end

function set_yticklabels!(axs::Array, yticklabs::Array{String,1})
    for indx in eachindex(axs)
        axs[indx].set_yticklabels(yticklabs);
    end

    return nothing
end

"""
    set_xyticklabels!(axs::Array, xticklabs::Array, yticklabs::Array)

Set X-axis and Y-axis ticks for a given
- `axs` Array of axis
- `xticklabs` Array of X-axis tick labels
- `yticklabs` Array of Y-axis tick labels
"""
function set_xyticklabels!(axs::Array, xticklabs::Array, yticklabs::Array)
    set_xticklabels!(axs, xticklabs);
    set_yticklabels!(axs, yticklabs);

    return nothing
end


text:
#=##############################################################################
# DESCRIPTION
    Linear algebra tools

# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2021
  * License   : MIT License
=###############################################################################

dot(A, B) = sum(a*b for (a,b) in zip(A, B))
norm(A) = sqrt(mapreduce(x->x^2, +, A))

cross1(A, B) = A[2]*B[3] - A[3]*B[2]
cross2(A, B) = A[3]*B[1] - A[1]*B[3]
cross3(A, B) = A[1]*B[2] - A[2]*B[1]

function cross!(out, A, B)
    out[1] = cross1(A, B)
    out[2] = cross2(A, B)
    out[3] = cross3(A, B)
end

function cross(A, B)
    out = zero(A)
    cross!(out, A, B)
    return out
end


text:
using GalacticNonconvex, GalacticOptim, Zygote, Pkg
using Test

@testset "GalacticNonconvex.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())

    GalacticNonconvex.Nonconvex.@load MMA
    GalacticNonconvex.Nonconvex.@load Ipopt
    GalacticNonconvex.Nonconvex.@load NLopt
    GalacticNonconvex.Nonconvex.@load BayesOpt
    GalacticNonconvex.Nonconvex.@load Juniper
    GalacticNonconvex.Nonconvex.@load Pavito
    GalacticNonconvex.Nonconvex.@load Hyperopt
    GalacticNonconvex.Nonconvex.@load MTS
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[1.5, 1.5])

    sol = solve(prob, MMA02())
    @test 10 * sol.minimum < l1

    sol = solve(prob, GCMMA())
    @test 10 * sol.minimum < l1

    sol = solve(prob, IpoptAlg())
    @test 10 * sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LN_NELDERMEAD))
    @test 10 * sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LD_LBFGS))
    @test 10 * sol.minimum < l1

    sol = solve(prob, MTSAlg())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler())
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    sol = solve(prob, MMA02(), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, GCMMA(), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, IpoptAlg(), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LN_NELDERMEAD), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LD_LBFGS), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, MTSAlg(), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10 * sol.minimum < l1

    ### suboptions
    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(; max_iter=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(; maxeval=100))
    @test 10 * sol.minimum < l1
end


text:
module Day25
__revise_mode__ = :eval

using Advent: save_input

const m = 20201227

function loopsize(target)
    x = 1
    for i in Iterators.countfrom()
        x = rem(x * 7, m)
        x == target && return i
    end
end

function solve_a(filename)
    cpub, dpub = parse.(Int, readlines(filename))
    cloops = loopsize(cpub)
    dloops = loopsize(dpub)
    resd = powermod(dpub, cloops, m)
    resc = powermod(cpub, dloops, m)
    @assert resd == resc
    resd
end

solve_b(filename) = nothing

solve(filename = save_input(2020, 25)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
#=
Created on Thursday 2 April 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Simple simulation to illustrate OTSIN.
=#

using OTSIN
using Plots
using StatsPlots
using Distributions
using DataFrames

# GENERATE DATA
# -------------

top_traits = [Normal(1.6, 0.4), Normal(2.2, 0.6), Normal(2.7, 0.3), Normal(5, 0.9)]
top_names = ["A" "B" "C" "D"]

bottom_traits = [1.5, 2, 2.5, 4, 5.4]
bottom_names = ["a", "b", "c", "d", "e"]

gen_inv_trait = Normal(3.7, 1.4)
spec_inv_trait = Normal(2.85, 0.25)

a, b = [0.20, 0.20, 0.1, 0.2, 0.3], [0.30, 0.4, 0.1, 0.2]

@assert sum(a) ≈ sum(b)

# PLOTTING
# --------

xvals = 0.5:0.01:7

top_prefs = [pdf.(top_traits[1], xvals) pdf.(top_traits[2], xvals) pdf.(top_traits[3], xvals) pdf.(top_traits[4], xvals)]

p_traits = plot(xvals, top_prefs, color=:orange, linewidth=2, label="", linestyle=:dash)
vline!(bottom_traits, color=:green, linewidth=2, label="bottom sp.")
vline!(mean.(top_traits), color=:orange, linewidth=2, alpha=0.5, linestyle=:dash, label="top sp.")

# gen invader
plot!(xvals, pdf.(gen_inv_trait, xvals), color=:purple, linewidth=2, label="", linestyle=:dot)
vline!([mean(gen_inv_trait)], color=:purple, linewidth=2, alpha=0.5, linestyle=:dot, label="top gen. invader sp.")

# spec invader
plot!(xvals, pdf.(spec_inv_trait, xvals), color=:red, linewidth=2, label="", linestyle=:dashdot)
vline!([mean(spec_inv_trait)], color=:red, linewidth=2, alpha=0.5, linestyle=:dashdot, label="top spec. invader sp.")

xticks!([bottom_traits; mean.(top_traits); mean.([gen_inv_trait, spec_inv_trait])],
                        [bottom_names; vec(top_names); "E"; "F"])
ylabel!("utility")

p_bar_bottom = bar(bottom_names, a, color=:green, label="", ylabel="frequency", title="Density bottom sp.")
p_bar_top = bar(vec(top_names), b, color=:orange, label="", ylabel="frequency", title="Density top sp.")

# OT SIMULATIONS
# --------------

function make_M(bottom_traits, top_traits; α=10)
    n, m = size(bottom_traits, 1), size(top_traits, 1)
    M = zeros(n, m)

    for i in 1:n
        for j in 1:m
            M[i,j] = pdf(top_traits[j], bottom_traits[i])
        end
    end
    return M
end

M = make_M(bottom_traits, top_traits, α=5)

p_heatmap_M = heatmap(M, xticks=(1:4, top_names), yticks=(1:5, bottom_names), title="Utility matrix M")

l = @layout [a{0.7w}  [b; c; d]]
plot(p_traits, p_heatmap_M, p_bar_top, p_bar_bottom, layout=l, size=(800, 500))
savefig("figures/simulation_outline.svg")

av_utilities = Matrix{Union{Nothing, Float32}}(nothing, 6, 6)

Ql = optimaltransport(M, a, b, λ=1)
Qh = optimaltransport(M, a, b, λ=20)

p_ot_neutral = groupedbar(bottom_names, Ql, labels=top_names, bar_position = :stack,
        title="Neutrally-driven", legend=:none)
p_ot_utility = groupedbar(bottom_names, Qh, labels=top_names, bar_position = :stack,
        title="Utility-driven", legend=:none)

av_utilities[1,1:4] = sum(Qh .* M, dims=1)[:]
av_utilities[2,1:4] = sum(Ql .* M, dims=1)[:]

plot(p_ot_neutral, p_ot_utility)

report_utility(Ql, M, bottom_names, top_names)

report_utility(Qh, M, bottom_names, top_names)

a_removal = copy(a)
a_removal[2] = 0.0
a_removal ./ sum(a_removal)

Q_extinction = optimaltransport(M, a_removal, b, λ=20)

av_utilities[4,1:4] = sum(Q_extinction .* M, dims=1)[:]

p_extincition = groupedbar(bottom_names, Q_extinction, labels=top_names, bar_position = :stack,
        title="Extinction of sp. b", legend=:none)

Q_extinction

M_gen_inv = make_M(bottom_traits, [top_traits; gen_inv_trait])

top_names_inv = [top_names..., "E"]

heatmap(M_gen_inv, xticks=(1:5, top_names_inv), yticks=(1:5, bottom_names))

b_inv = OTSIN.normalize([b; 0.1])

Q_gen_inv = optimaltransport(M_gen_inv, a, b_inv, λ=20.0)

av_utilities[5,1:5] = sum(Q_gen_inv .* M_gen_inv, dims=1)[:]

p_gen_inv = groupedbar(bottom_names, Q_gen_inv, labels=[top_names "E"], bar_position = :stack,
        title="Effect of generalist invader", legend=:top)

report_utility(Q_gen_inv, M_gen_inv, bottom_names, top_names_inv)

M_spec_inv = make_M(bottom_traits, [top_traits; spec_inv_trait])

top_names_inv = [top_names "F"]

heatmap(M_spec_inv, xticks=(1:5, top_names_inv), yticks=(1:5, bottom_names))

Q_spec_inv = optimaltransport(M_spec_inv, a, b_inv, λ=20.0)

av_utilities[6,[1, 2, 3, 4, 6]] = sum(Q_spec_inv .* M_spec_inv, dims=1)[:]

p_spec_inv =  groupedbar(bottom_names, Q_spec_inv, labels=[top_names "F"], bar_position = :stack,
        title="Effect of specialist invader", legend=:top)

plot(p_ot_utility, p_ot_neutral, p_optimal_free, p_extincition, p_gen_inv, p_spec_inv,
        layout=(2, 3), size=(1500, 1.5*500))

savefig("figures/simulation_settings.svg")

p_ut_a, p_ut_b = pert_utility(M, a, b, λ=20)
p_entr_a, p_entr_b = pert_entropy(M, a, b, λ=20)

plot_pert_bottom = scatter(bottom_names, [p_ut_a p_entr_a],
    label=["relative effect utility", "relative effect entropy"],
    title="Perturbations bottom sp.", markershape=[:circle :utriangle])
plot_pert_top = scatter(top_names, [p_ut_b p_entr_b],
    label=["relative effect utility", "relative effect entropy"],
    legend=:none, title="Perturbations top sp.", markershape=[:circle :utriangle])

plot(plot_pert_bottom, plot_pert_top, layout=grid(2, 1), size=(600, 400))
savefig("figures/perturbation.svg")

results = DataFrame(av_utilities, [:A, :B, :C, :D, :E, :F])

println(results)

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

"""
    pagerank(g::ADiGraph, α=0.85, n=100, ϵ = 1.0e-6)

Calculates the [PageRank](https://en.wikipedia.org/wiki/PageRank) of the graph
`g`. Can optionally specify a different damping factor (`α`), number of
iterations (`n`), and convergence threshold (`ϵ`). If convergence is not
reached within `n` iterations, an error will be returned.
"""
function pagerank(g::ADiGraph, α=0.85, n=100, ϵ = 1.0e-6)
    A = adjacency_matrix(g,:in,Float64)
    S = 1 ./ vec(sum(A, dims=1))
    S[S .== Inf] .= 0.0
    M = A' # need a separate line due to bug #17456 in julia
    M = (Diagonal(S) * M)'
    N = nv(g)
    x = fill(1/N, N)
    p = fill(1/N, N)
    dangling_weights = p
    is_dangling = findall(S .== 0)

    for _ in 1:n
        xlast = x
        x = α * (M * x + sum(x[is_dangling]) * dangling_weights) + (1 - α) * p
        err = sum(abs, x - xlast)
        err < N * ϵ && return x
    end
    error("Pagerank did not converge after $n iterations.")
end