text:
import Base.convert

# How to convert a Dual to a Float
convert{N,T}(Float64, x::Dual{N,T}) = value(x)

"""
    Assemble Jacobian and Residual
"""
function assembleResidualAndJacobian{T}(solver::Solver, sys::System{T})
    startLog(main_perf_log, "assembleResidualAndJacobian()")

    mesh = sys.mesh

    updateGhostedSolution!(solver)

    solution = solver.ghosted_solution

    vars = sys.variables
    n_vars = length(vars)

    # Reset the Residual and Jacobian
    zero!(solver.rhs)
    zero!(solver.mat)

    # Each inner array is of length n_dofs for each var
    var_residuals = Array{Array{T}}(n_vars)
    for i in 1:n_vars
        var_residuals[i] = Array{T}(0)
    end

    # Each inner matrix is i_n_dofs * j_n_dofs
    var_jacobians = Matrix{Matrix{Float64}}((n_vars, n_vars))
    for i in 1:n_vars
        for j in 1:n_vars
            var_jacobians[i,j] = Matrix{Float64}((0,0))
        end
    end

    # Execute the element loop and accumulate Kernel contributions
    for elem in mesh.local_elements
        reinit!(sys, elem, solution)

        # Resize and zero residual vectors and jacobian matrices
        for i_var in sys.variables
            resize!(var_residuals[i_var.id], i_var.n_dofs)
            fill!(var_residuals[i_var.id], 0.)

            for j_var in sys.variables
                # Because for some inexplicable reason Julia doesn't provide a resize!() for Matrices...
                if size(var_jacobians[i_var.id, j_var.id]) != (i_var.n_dofs, j_var.n_dofs)
                    var_jacobians[i_var.id, j_var.id] = Matrix{Float64}((i_var.n_dofs, j_var.n_dofs))
                end

                fill!(var_jacobians[i_var.id, j_var.id], 0.)
            end
        end

        # Get the Residual/Jacobian contributions from all Kernels
        computeResidualAndJacobian!(var_residuals, var_jacobians, vars, sys.kernels)

        # Scatter those entries back out into the Residual and Jacobian
        for i_var in sys.variables
            plusEquals!(solver.rhs, var_residuals[i_var.id], i_var.dofs)

            for j_var in sys.variables
                plusEquals!(solver.mat, var_jacobians[i_var.id, j_var.id], i_var.dofs, j_var.dofs)
            end
        end
    end

    assemble!(solver.mat)
    assemble!(solver.rhs)

    # Now apply BCs

    boundary_info = mesh.boundary_info
    bcs = sys.bcs

    # First: get the set of boundary IDs we need to operate on:
    bids = Set{Int64}()
    for bc in bcs
        union!(bids, bc.bids)
    end

    # Reusable storage for calling residual and jacobian calculations on NodalBCs
    temp_residual = Array{T}(1)
    temp_jacobian = Array{Float64}(n_vars)

    rows_to_zero = []

    entries_in_bc_rows= []

    proc_id = MPI.Comm_rank(MPI.COMM_WORLD)

    # Now go over each nodeset and apply the BCs
    for bid in bids

        # Grab the nodeset for this bid
        node_list = boundary_info.node_list[bid]

        # Iterate over each node and apply the boundary conditions
        for node in node_list
            if node.processor_id == proc_id
                reinit!(sys, node, solution)

                # Apply all of the BCs that should be applied here
                for bc in bcs
                    if bid in bc.bids
                        computeResidualAndJacobian!(temp_residual, temp_jacobian, vars, bc)

                        # First set the residual
                        solver.rhs[bc.u.nodal_dof] = temp_residual[1]

                        # Now - we need to zero out the row in the matrix corresponding to this dof
                        push!(rows_to_zero, bc.u.nodal_dof)

                        # And put this piece in place
                        for v in vars
                            #solver.mat[bc.u.nodal_dof,v.nodal_dof] = temp_jacobian[v.id]
                            push!(entries_in_bc_rows, (bc.u.nodal_dof,v.nodal_dof,temp_jacobian[v.id]))
                        end
                    end
                end
            end
        end
    end

    zeroRows!(solver.mat, rows_to_zero)
    assemble!(solver.mat)

    for entry in entries_in_bc_rows
        solver.mat[entry[1], entry[2]] = entry[3]
    end

    assemble!(solver.mat)
    assemble!(solver.rhs)

    stopLog(main_perf_log, "assembleResidualAndJacobian()")
end


text:
using StatFiles, StatsModels, DataFrames, DataFramesMeta, GLM
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))[:,[:wage; :tenure; :ttl_exp; :collgrad; :industry]]
dropmissing!(df)
desc = describe(df, :eltype)
for i in axes(desc, 1)  # needed only for lm()
  if desc[i,:eltype] == Float32
    sym = desc[i,:variable]
    @transform!(df, @byrow $sym=Float64($sym))
  end
end
f = @formula(wage ~ 1 + tenure + ttl_exp * collgrad)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

fit = lm(f, df)

isa(fit.model, LinearModel)
_terms = fit.mf.f.rhs.terms
k = length(_terms)
_I = zeros(k,k); _I[diagind(_I)] .= 1.

termcols = Dict(t => _I[:,i] for (i,t) in enumerate(_terms))
newtermcols = Dict{AbstractTerm, Vector{Float64}}()
for pair in termcols
  if isa(pair.first, InteractionTerm)
    newtermcols[InteractionTerm((pair.first.terms[2], pair.first.terms[1]))] = pair.second
  end
end
merge!(termcols, newtermcols)

Hstr = "wage + 2 * tenure ==3"
H = Meta.parse(Hstr)

dump(Meta.parse("1 + ttl_exp & collgrad"))

function processH(H::Expr)
  if H.head==:call
    if H.args[1] == :(==)
      return Expr(:(-), H.args[2], H.args[3])
    end
  end
end

ext: jl


text:
module Burgers

using DataDeps, MAT, MLUtils
using NeuralOperators, Flux
using CUDA, FluxTraining, BSON

include("Burgers_deeponet.jl")

function register_burgers()
    register(DataDep(
        "Burgers",
        """
        Burgers' equation dataset from
        [fourier_neural_operator](https://github.com/zongyi-li/fourier_neural_operator)
        """,
        "http://www.med.cgu.edu.tw/NeuralOperators/Burgers_R10.zip",
        "9cbbe5070556c777b1ba3bacd49da5c36ea8ed138ba51b6ee76a24b971066ecd",
        post_fetch_method=unpack
    ))
end

function get_data(; n=2048, Δsamples=2^3, grid_size=div(2^13, Δsamples), T=Float32)
    file = matopen(joinpath(datadep"Burgers", "burgers_data_R10.mat"))
    x_data = T.(collect(read(file, "a")[1:n, 1:Δsamples:end]'))
    y_data = T.(collect(read(file, "u")[1:n, 1:Δsamples:end]'))
    close(file)

    x_loc_data = Array{T, 3}(undef, 2, grid_size, n)
    x_loc_data[1, :, :] .= reshape(repeat(LinRange(0, 1, grid_size), n), (grid_size, n))
    x_loc_data[2, :, :] .= x_data

    return x_loc_data, reshape(y_data, 1, :, n)
end

function get_dataloader(; ratio::Float64=0.9, batchsize=100)
    𝐱, 𝐲 = get_data(n=2048)
    data_train, data_test = splitobs((𝐱, 𝐲), at=ratio)

    loader_train = DataLoader(data_train, batchsize=batchsize, shuffle=true)
    loader_test = DataLoader(data_test, batchsize=batchsize, shuffle=false)

    return loader_train, loader_test
end

__init__() = register_burgers()

function train(; cuda=true, η₀=1f-3, λ=1f-4, epochs=500)
    if cuda && CUDA.has_cuda()
        device = gpu
        CUDA.allowscalar(false)
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    model = FourierNeuralOperator(ch=(2, 64, 64, 64, 64, 64, 128, 1), modes=(16, ), σ=gelu)
    data = get_dataloader()
    optimiser = Flux.Optimiser(WeightDecay(λ), Flux.ADAM(η₀))
    loss_func = l₂loss

    learner = Learner(
        model, data, optimiser, loss_func,
        ToDevice(device, device),
    )

    fit!(learner, epochs)

    return learner
end

end


text:
function show_op(o::Operator)
    o == crossover && return 'c'
    o == mutation && return 'm'
    return '-'
end

function show_ref(e::Entry)
    e.index in keys(ref_counter) && return ref_counter[e.index]
    return -1
end

function Base.show(io::IO, e::Entry)
    println(io, "  ID Index  Op    Entry(id)         Entry(index)    ref")
    @printf(io, "%4d  %4d   %c ", e.id, e.index, show_op(e.operator))
    @printf(io, "(%4d %4d %4d)   ",  e.id_entry[1], e.id_entry[2], e.id_entry[3])
    @printf(io, "(%4d %4d %4d)", e.index_entry[1], e.index_entry[2], e.index_entry[3])
    @printf(io, "  %3d", show_ref(e))
end

function Base.show(io::IO, x::MIME"text/plain", ve::Vector{Entry})
    println(io, "  ID Index  Op    Entry(id)         Entry(index)    ref")
    for (i, e) in enumerate(ve)
        @printf(io, "%4d  %4d   %c ", e.id, e.index, show_op(e.operator))
        @printf(io, "(%4d %4d %4d)   ", e.id_entry[1], e.id_entry[2], e.id_entry[3])
        @printf(io, "(%4d %4d %4d)", e.index_entry[1], e.index_entry[2], e.index_entry[3])
        @printf(io, "   %3d", show_ref(e))
        i  == length(ve) || println(io)
    end
end

function print_latest_fitness(rec, g)
    if g == 1
        println("     train_fitness  test_fitness")
    end
    @printf("%4d %10.4f     %10.4f\n", g, rec.data["train_fitness"][end], rec.data["test_fitness"][end])
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
export
    COLOR_ASPHALT,
    COLOR_LANE_MARKINGS_WHITE,
    COLOR_LANE_MARKINGS_YELLOW,
    COLOR_CAR_EGO,
    COLOR_CAR_OTHER,
    MONOKAI

const COLOR_ASPHALT       = colorant"0x708090FF"
const COLOR_LANE_MARKINGS_WHITE = colorant"0xDFDFDFFF"
const COLOR_LANE_MARKINGS_YELLOW = colorant"0xFFEF00FF"
const COLOR_CAR_EGO       = colorant"0x79ABFFFF" # bluish
const COLOR_CAR_OTHER     = colorant"0xFF007FFF" # reddish

const MONOKAI = Dict(
            "foreground" => colorant"0xCFBFADFF",
            "background" => colorant"0x272822FF",
            "color1"     => colorant"0x52E3F6FF", # light blue
            "color2"     => colorant"0xA7EC21FF", # light green
            "color3"     => colorant"0xFF007FFF", # red
            "color4"     => colorant"0xF9971FFF", # orange
            "color5"     => colorant"0x79ABFFFF", # cobalt
        )

function Vec.lerp(a::Colorant, b::Colorant, t::Real)

    ra = red(a)
    rb = red(b)
    ga = green(a)
    gb = green(b)
    ba = blue(a)
    bb = blue(b)
    aa = alpha(a)
    ab = alpha(b)

    r = ra + (rb - ra)*t
    g = ga + (gb - ga)*t
    b = ba + (bb - ba)*t
    a = aa + (ab - aa)*t

    RGBA(r,g,b,a)
end


text:
using Test, SafeTestsets

const is_WINDOWS = Sys.iswindows()

println("importing BellScenario.jl")
@time using BellScenario

function _test_runner()
    @time @testset "BellScenario.jl" begin
        @testset "unit tests:" begin
            println("running unit tests.")
            for test in readdir("./test/unit/")
                # run only julia files in test directory
                if occursin(r"^.*\.jl$", test)
                    println("./unit/$test")
                    @time include("./unit/$test")
                end
            end
        end
        @time @testset "integration tests:" begin
            println("running integration tests.")

            println("./test/integration/file_io.jl")
            @time @safetestset "./test/integration/file_io.jl" begin include("integration/file_io.jl") end

            if !is_WINDOWS
                println("./test/integration/LocalPolytope/adjacency_decomposition.jl")
                @time @safetestset "./test/integration/LocalPolytope/adjacency_decomposition.jl" begin include("integration/LocalPolytope/adjacency_decomposition.jl") end
            end

            # for test in readdir("./test/integration/")
            #     # run only julia files in test directory
            #     if occursin(r"^.*\.jl$", test)
            #         println("./integration/$test")
            #         @time include("./integration/$test")
            #     end
            # end
        end
        @testset "regression tests:" begin
            println("running regression tests.")
            for test in readdir("./test/regression/")
                # run only julia files in test directory
                if occursin(r"^.*\.jl$", test)
                    println("./regression/$test")
                    @time include("./regression/$test")
                end
            end
        end
    end
end

# Pkg.test("BellScenario") runs from ./test directory. Development tests from root.
dir = pwd()
if occursin(r".*test$", dir)
    cd(_test_runner, "../")
elseif occursin(r".*BellScenario", dir)
    _test_runner()
else
    error("runtests.jl is currently running from the $(pwd()) directory with contents $(readdir()). runtests.jl must be run from the ./BellScenario.jl or ./BellScenario.jl/test directories.")
end


text:
mutable struct UnfoldProposeState{T}
    choices::DynamicChoiceMap # NOTE: could be read-only vector assignment
    weight::Float64
    retvals::Vector{T}
    state::T
end

function process_new!(gen_fn::Unfold{T,U}, params::Tuple, key::Int,
                      state::UnfoldProposeState{T}) where {T,U}
    local new_state::T
    kernel_args = (key, state.state, params...)
    (submap, weight, new_state) = propose(gen_fn.kernel, kernel_args)
    set_submap!(state.choices, key, submap)
    state.weight += weight
    state.retvals[key] = new_state
    state.state = new_state
end

function propose(gen_fn::Unfold{T,U}, args::Tuple) where {T,U}
    len = args[1]
    init_state = args[2]
    params = args[3:end]
    choices = choicemap()
    state = UnfoldProposeState{T}(choices, 0., Vector{T}(undef,len), init_state)
    for key=1:len
        process_new!(gen_fn, params, key, state)
    end
    (state.choices, state.weight, PersistentVector{T}(state.retvals))
end


text:
"""
    smoke_anim(sname::String, imgs::AbstractArray, bgimg::Union{AbstractArray,Nothing}; figtitle=" ", cb_title=" ", fps=10, ksize= 3, inverse=false, col=:rainbow, alpha=1.0, clim=(-Inf,Inf))

## Description
Function to create an animation heatmap
## Arguments
- sname: name of the file to be saved (.gif)
- imgs: array of original images
- bgimg: background image
- figtitle: figure title
- fps: frames per second
- ksize: size of the kernel (Kernel.gaussian)
- inverse: order of subtraction (background subtraction)
- col: color palette
- cb_title: colorbar title
- alpha: opacity
- clim: fixed limits for the heatmap
saves a gif animation
## Example
```jldoctest
julia> using VisualEFM

julia> fnames = "frame" .* string.(collect(1:100)) .* ".jpg";

# (Not run)
julia> im = read_image.(fnames);

julia> smoke_anim("anim_smoke.gif", im, nothing, figtitle="Animation", cb_title="Grayscale", fps=20)
```
"""
function smoke_anim(sname::String, imgs::AbstractArray, bgimg::Union{AbstractArray,Nothing};
                    figtitle=" ", cb_title=" ", fps=10, ksize= 3, inverse=false,
                    col=:rainbow, alpha=1.0, clim=(-Inf,Inf))

    n = length(imgs)

    if isnothing(bgimg)
        anim = @animate for i in 1:n
            img_gs = gauss_grayscale(imgs[i], ksize=ksize)
            img_gs = Float64.(img_gs)
            plot(imgs[i], axis=nothing, bordercolor="white")
            heatmap!(img_gs, color=col, colorbar_title=string(cb_title), clim=clim, alpha=alpha)
            title!(string(figtitle))
        end

    else
        bgimg_gs = gauss_grayscale(bgimg, ksize=ksize)
        anim = @animate for i in 1:n
            img_gs = gauss_grayscale(imgs[i], ksize=ksize)
            img_diff = backsubtraction(img_gs,bgimg_gs, inverse=inverse)
            img_diff = Float64.(img_diff)
            plot(imgs[i], axis=nothing, bordercolor="white")
            heatmap!(img_diff, color=col, colorbar_title=string(cb_title), clim=clim, alpha=alpha)
            title!(string(figtitle))
        end
    end

    gif(anim, sname, fps=fps)

end

"""
    smoke_statsmap(statfun, imgs::AbstractArray, bgimg::Union{AbstractArray,Nothing}; figtitle=" ", cb_title=" ", ksize = 3, inverse = false, col=:rainbow, alpha=1.0, clim=(-Inf,Inf), showPlot=true)

## Description
Plot statistics of an array of images.
## Arguments
- statfun: function to be applied to an array of images (works with mean and std)
- imgs: array of original images
- bgimg: background image
- figtitle: figure title
- ksize: size of the kernel (Kernel.gaussian)
- inverse: order of subtraction (background subtraction)
- col: color palette
- cb_title: colorbar title
- alpha: opacity
- clim: fixed limits for the heatmap
- showPlot: true to return a plot
returns a plot
## Example
```jldoctest
julia> using VisualEFM

julia> fnames = "frame" .* string.(collect(1:100)) .* ".jpg";

# (Not run)
julia> im = read_image.(fnames);

julia> smoke_statsmap(Statistics.mean, im, nothing, cb_title="Mean Grayscale Values")
```
"""
function smoke_statsmap(statfun, imgs::AbstractArray, bgimg::Union{AbstractArray,Nothing};
                        figtitle=" ", cb_title=" ", ksize = 3, inverse = false,
                        col=:rainbow, alpha=1.0, clim=(-Inf,Inf), showPlot=true)

    imgs_gray = gauss_grayscale.(imgs, ksize=ksize)

    if !isnothing(bgimg)
        bgim_gray = gauss_grayscale(bgimg, ksize=ksize)
        imgs_gray = [backsubtraction(i,bgim_gray,inverse=inverse) for i in  imgs_gray]
    end

    imgs_array = [Float64.(j) for j in imgs_gray]

    res = statfun(imgs_array)

    if showPlot
        fig = heatmap(res, color=col, alpha=alpha, clim=clim, axis=nothing,
               	     title=string(figtitle), bordercolor="white", colorbar_title=string(cb_title))
        display(fig)
    end

    return [imgs_array, res]
end


text:
module MagneticResonanceSignals

using StaticArrays
using DataStructures
using AxisArrays

using Statistics
using Unitful
using FFTW

using LinearAlgebra

# Plotting
using RecipesBase
using UnitfulRecipes
using Colors

include("core.jl")
include("LCOSY.jl")
include("PRESS.jl")

include("fixedstring.jl")  # IO Util

include("io_twix.jl")
include("io_felix.jl")
include("io_rda.jl")
include("io_nmrpipe.jl")

# TODO: Kill off MRAxis
include("MRAxis.jl")

include("coil_combination.jl")
include("processing.jl")
include("hsvd.jl")

include("windows.jl")
include("plotting.jl")

# Data structures and metadata access for MR experiments
export
    MRExperiment,
    meta_search,
    standard_metadata,
    scanner_time,
    mr_load,
    sampledata,
    count_cycles

# IO
export load_rda,
    save_rda,
    load_twix,
    save_felix,
    save_nmrpipe

# Conversion utils
export twix_to_nmrpipe

# Axes
# TODO: Deprecate / remove all this in favour of AxisArray
export
    MRAxis,
    hertz_to_ppm,
    ppm_to_hertz,
    time_axis,
    frequency_axis,
    frequency_axis_ppm

# Low level signal processing
export
    # Channel combination
    pca_channel_combiner,
    combine_channels,
    # Spectral windows
    zeropad,
    apply_window, apply_window!,
    sinebell

# High level signal processing
export
    # Raw FIDs
    extract_fids,
    # Averaging of repeats
    simple_averaging,
    # Fourier transform
    spectrum

# Plotting
export
    felix_colors

const _local_basefactors = Unitful.basefactors

function __init__()
    Unitful.register(MagneticResonanceSignals)
    merge!(Unitful.basefactors, _local_basefactors)
end

end # module

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using MCMC

println("    Testing basic EmpMCTuner constructors...")

mctuners = [EmpMCTuner(0.85),
  EmpMCTuner(0.85, adaptStep=50),
  EmpMCTuner(0.85, maxStep=100),
  EmpMCTuner(0.85, targetPath=0.75),
  EmpMCTuner(0.85, verbose=true)]

println("    Testing that EmpMCTuner tuners works with all samplers...")

npars = 3
mcmodel = model(v->-dot(v,v), grad=v->-2v, tensor=v->-2*ones(npars), dtensor=v->zeros(npars, npars, npars),
  init=ones(npars))

for mctuner in mctuners
  mcsamplers = [HMC(0.75, mctuner),
    RMHMC(0.75, mctuner),
    MALA(0.75, mctuner),
    SMMALA(0.75, mctuner),
    PMALA(0.75, mctuner)]

  for mcsampler in mcsamplers
    mcchain = run(mcmodel, HMC(0.75, mctuner), SerialMC(steps=5000, burnin=1000))
  end
end


text:
@testset "Testing MNIST Conv" begin
    model = Chain(
        # First convolution, operating upon a 28x28 image
        CrossCor((3, 3), 1=>16, pad=(1,1,1,1), relu),
        MaxPool((2,2)),

        # Second convolution, operating upon a 14x14 image
        CrossCor((3, 3), 16=>32, pad=(1,1,1,1), relu),
        MaxPool((2,2)),

        # Third convolution, operating upon a 7x7 image
        CrossCor((3, 3), 32=>32, pad=(1,1,1,1), relu),
        MaxPool((2,2)),

        # Reshape 3d tensor into a 2d one, at this point it should be (3, 3, 32, N)
        # which is where we get the 288 in the `Dense` layer below:
        x -> reshape(x, :, size(x, 4)),
        Dense(288, 10),

        # Finally, softmax to get nice probabilities
        softmax,
    )

    backend = nGraph.Backend()
    x = rand(Float32, 28, 28, 1, 100)
    expected = model(x)
    f = nGraph.compile(backend, model, x)

    # Test that the forward passes match
    @test isapprox(expected, parent(f()))

    @time model(x)
    @time model(x)
    @time f()
    @time f()
end


text:
include("signature.jl")
include("domain.jl")
include("problem.jl")
include("state.jl")
include("action.jl")
include("diff.jl")
include("utils.jl")

"""
    satisfy(domain::Domain, state::State, term::Term)
    satisfy(domain::Domain, state::State, terms::AbstractVector{<:Term})

Returns whether the queried `term` or `terms` can be satisfied in the given
`domain` and `state`.
"""
satisfy(domain::Domain, state::State, term::Term) =
    satisfy(domain, state, [term])
satisfy(domain::Domain, state::State, terms::AbstractVector{<:Term}) =
    error("Not implemented.")

"""
    satisfiers(domain::Domain, state::State, term::Term)
    satisfiers(domain::Domain, state::State, terms::AbstractVector{<:Term})

Returns a list of satisfying substitutions of the queried `term` or `terms`
within the given `domain` and `state`.
"""
satisfiers(domain::Domain, state::State, term::Term) =
    satisfiers(domain, state, [term])
satisfiers(domain::Domain, state::State, terms::AbstractVector{<:Term}) =
    error("Not implemented.")

"""
    evaluate(domain::Domain, state::State, term::Term)

Evaluates a grounded `term` in the given `domain` and `state`. If `term`
refers to a numeric fluent, the value of the fluent is returned.
For logical predicates, `evaluate` is equivalent to `satisfy`.
"""
evaluate(domain::Domain, state::State, term::Term) =
    error("Not implemented.")

"""
    initstate(domain::Domain, problem::Problem)
    initstate(domain::Domain, objtypes[, fluents])

Construct the initial state for a given planning `domain` and `problem`, or
from a `domain`, a map of objects to their types (`objtypes`), and an optional
list of `fluents`.

Fluents can either be provided as a list of `Term`s representing the initial
fluents in a PDDL problem, or as a map from fluent names to fluent values.
"""
initstate(domain::Domain, problem::Problem) =
    error("Not implemented.")
initstate(domain::Domain, objtypes::AbstractDict, fluents=Term[]) =
    error("Not implemented.")

"""
    goalstate(domain::Domain, problem::Problem)
    goalstate(domain::Domain, objtypes, terms)

Construct a (partial) goal state from a `domain` and `problem`, or from
a `domain`, a map of objects to their types (`objtypes`), and goal `terms`.
"""
goalstate(domain::Domain, problem::Problem) =
    error("Not implemented.")
goalstate(domain::Domain, objtypes::AbstractDict, terms) =
    error("Not implemented.")

"""
    transition(domain::Domain, state::State, action::Term)
    transition(domain::Domain, state::State, actions)

Returns the successor to `state` in the given `domain` after applying a single
`action` or a set of `actions` in parallel.
"""
transition(domain::Domain, state::State, action::Term) =
    error("Not implemented.")
transition(domain::Domain, state::State, actions) =
    error("Not implemented.")

"""
    transition!(domain::Domain, state::State, action::Term)
    transition!(domain::Domain, state::State, actions)

Variant of [`transition`](@ref) that modifies `state` in place.
"""
transition!(domain::Domain, state::State, action::Term; options...) =
    error("Not implemented.")
transition!(domain::Domain, state::State, actions; options...) =
    error("Not implemented.")

"""
    available(domain::Domain, state::State, action::Action, args)
    available(domain::Domain, state::State, action::Term)

Check if an `action` parameterized by `args` can be executed in the given
`state` and `domain`. Action parameters can also be specified as the arguments
of a compound `Term`.
"""
available(domain::Domain, state::State, action::Action, args) =
    error("Not implemented.")
available(domain::Domain, state::State, action::Term) =
    available(domain, state, get_actions(domain)[action.name], action.args)

"""
    available(domain::Domain, state::State)

Return an iterator over available actions in a given `state` and `domain`.
"""
available(domain::Domain, state::State) =
    error("Not implemented.")

"""
    execute(domain::Domain, state::State, action::Action, args)
    execute(domain::Domain, state::State, action::Term)

Execute an `action` parameterized by `args` in the given `state`, returning
the resulting state.  Action parameters can also be specified as the arguments
of a compound `Term`.
"""
execute(domain::Domain, state::State, action::Action, args; options...) =
    error("Not implemented.")
execute(domain::Domain, state::State, action::Term; options...) =
    if action.name == get_name(no_op)
        execute(domain, state, no_op, (); options...)
    else
        execute(domain, state, get_actions(domain)[action.name],
                action.args; options...)
    end

"""
    execute!(domain::Domain, state::State, action::Action, args)
    execute!(domain::Domain, state::State, action::Term)

Variant of [`execute`](@ref) that modifies `state` in-place.
"""
execute!(domain::Domain, state::State, action::Action, args; options...) =
    error("Not implemented.")
execute!(domain::Domain, state::State, action::Term; options...) =
    if action.name == get_name(no_op)
        execute!(domain, state, no_op, (); options...)
    else
        execute!(domain, state, get_actions(domain)[action.name],
                 action.args; options...)
    end

"""
    relevant(domain::Domain, state::State, action::Action, args)
    relevant(domain::Domain, state::State, action::Term)

Check if an `action` parameterized by `args` is relevant (can lead to) a `state`
in the given `domain`. Action parameters can also be specified as the arguments
of a compound `Term`.
"""
relevant(domain::Domain, state::State, action::Action, args) =
    error("Not implemented.")
relevant(domain::Domain, state::State, action::Term) =
    relevant(domain, state, get_actions(domain)[action.name], action.args)

"""
    relevant(domain::Domain, state::State)

Return an iterator over relevant actions in a given `state` and `domain`.
"""
relevant(domain::Domain, state::State) =
    error("Not implemented.")

"""
    regress(domain::Domain, state::State, action::Action, args)
    regress(domain::Domain, state::State, action::Term)

Compute the pre-image of an `action` parameterized by `args` with respect to
a `state`. Action parameters can also be specified as the arguments of a
compound `Term`.
"""
regress(domain::Domain, state::State, action::Action, args; options...)=
    error("Not implemented.")
regress(domain::Domain, state::State, action::Term; options...) =
    if action.name == get_name(no_op)
        regress(domain, state, no_op, (); options...)
    else
        regress(domain, state, get_actions(domain)[action.name],
                action.args; options...)
    end

"""
    regress!(domain::Domain, state::State, action::Action, args)
    regress!(domain::Domain, state::State, action::Term)

Variant of [`regress`](@ref) that modifies `state` in-place.
"""
regress!(domain::Domain, state::State, action::Action, args; options...)=
    error("Not implemented.")
regress!(domain::Domain, state::State, action::Term; options...) =
    if action.name == get_name(no_op)
        regress!(domain, state, no_op, (); options...)
    else
        regress!(domain, state, get_actions(domain)[action.name],
                 action.args; options...)
    end

"""
    update(domain::Domain, state::State, diff::Diff)

Updates a PDDL `state` with a state difference.
"""
update(domain::Domain, state::State, diff::Diff) =
    update(domain, copy(state), diff)

"""
    update!(domain::Domain, state::State, diff::Diff)

Updates a PDDL `state` in-place with a state difference.
"""
update!(domain::Domain, state::State, diff::Diff) =
    error("Not implemented.")


text:
using Random;

export contour_block_SS

"""
    contour_block_SS([eltype,] nep [,mintegrator];[tol,][logger,][σ,][radius,][linsolvercreator,][N,][neigs,][k,][L])

This is an implementation of the block_SS contour integral method which is
based on the computation of higher order moments.
The contour is an ellipse centered at `σ` with radii given in `radius`, or if only one `radius` is given, the contour is a circle. The numerical quadrature method is specified in `mintegrator`,
which is a type inheriting from `MatrixIntegrator`, by default
`MatrixTrapezoidal`. For a parallell implementation of the
integrator use `MatrixTrapezoidalParallel`.
 The integer `k`
specifies size of the probe subspace. `N` corresponds to the
number of quadrature points. The integer L specifies the number of moments.
Ellipses are the only supported contours. The
`linsolvercreator` must create a linsolver that can handle (rectangular) matrices
as right-hand sides, not only vectors. We integrate in complex arithmetic so
`eltype` must be complex type.

# Example

```julia-repl
julia> nep=SPMF_NEP([[0 1 ; 1 1.0], [1 0 ; 0 0]], [s->one(s),s->exp(1im*s^2)]);
julia> λ,V=contour_assu(nep,radius=3,neigs=6)
julia> @show λ
6-element Array{Complex{Float64},1}:
  4.496403249731884e-15 + 2.506628274630998im
        -2.506628274631 - 2.8727020762175925e-15im
  3.219972424519104e-16 - 2.5066282746310034im
     2.5066282746310096 - 1.1438072192922029e-15im
 -2.3814273710772784e-7 - 7.748469160458366e-8im
   2.381427350935646e-7 + 7.748467479992284e-8im
```

# References

* Asakura, Sakurai, Tadano, Ikegami, Kimura, A numerical method for nonlinear eigenvalue problems using contour integrals, JSIAM Letters, 2009 Volume 1 Pages 52-55
* Van Beeumen,  Meerbergen, Michiels. Connections between contour integration and rational Krylov methods for eigenvalue problems, 2016, TW673, https://lirias.kuleuven.be/retrieve/415487/
"""
contour_block_SS(nep::NEP;params...)=contour_block_SS(ComplexF64,nep;params...);
contour_block_SS(nep::NEP,MIntegrator;params...)=contour_block_SS(ComplexF64,nep,MIntegrator;params...);
function contour_block_SS(
    ::Type{T},
    nep::NEP,
    ::Type{MIntegrator}=MatrixTrapezoidal;
    tol::Real=sqrt(eps(real(T))), # Note tol is quite high for this method
    σ::Number=zero(complex(T)),
    logger=0,
    linsolvercreator=BackslashLinSolverCreator(),
    neigs=Inf, # Number of wanted eigvals (currently unused)
    k::Integer=3, # Columns in matrix to integrate
    radius::Union{Real,Tuple,Array}=1, # integration radius
    N::Integer=1000,  # Nof quadrature nodes
    K::Integer=3, # Nof moments
    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
    sanity_check=true,
    Shat_mode=:native, # native or JSIAM-mode
    rank_drop_tol=tol # Used in sanity checking
)where{T<:Number, MIntegrator<:MatrixIntegrator}

    @parse_logger_param!(logger)

    n = size(nep,1);

    # Notation: L in JSIAM-paper corresponds to k in Beyn's paper.
    # Input params the same as contourbeyn, but
    # the code is like JSIAM-paper
    L=k

    Random.seed!(10); # Reproducability (not really)
    U = rand(T,n,L);
    V = rand(T,n,L);

    function local_linsolve(λ::TT,V::Matrix{TT}) where {TT<:Number}
        local M0inv::LinSolver = create_linsolver(linsolvercreator, nep, λ+σ);
        # This requires that lin_solve can handle rectangular
        # matrices as the RHS
        return lin_solve(M0inv,V);
    end

    # The step-references refer to the JSIAM-paper

    push_info!(logger,"Computing integrals")

    local Shat
    push_info!(logger,"Forming Mhat and Shat")
    Shat = zeros(T,n,L,2*K)
    Mhat = zeros(T,L,L,2*K)

    length(radius)==1 ? radius1=(radius,radius) : radius1=radius

    if (Shat_mode==:JSIAM)
        # This is the way the JSIAM-paper proposes to compute Shat
        if (length(radius)>1)
            error("JSIAM Shat_mode does not support ellipses");
        end

        # Quadrature points: Only circle supported
        w = exp.(2im*pi*(0.5 .+ (0:(N-1)))/N);
        omega = radius*w;

        push_info!(logger,"Forming all linear systems F(s)^{-1}V:",
                   continues=true)
        # Step 2: Precompute all the linear systems
        FinvV =zeros(T,n,L,N);
        for k = 1:N
            FinvV[:,:,k]=local_linsolve(omega[k],V);
        end
        push_info!(logger,"");

        # Step 3-4: Compute all the integrals and store in Shat (

        for k=0:(2*K-1)
            for j=0:N-1
                d=((omega[j+1])/radius)^(k+1)
                Shat[:,:,k+1] += d*FinvV[:,:,j+1]/N;
            end
        end

    elseif (Shat_mode==:native)

        # This deviates from the JSIAM-paper description, since
        # we do not precompute linear systems, but instead
        # compute linear system in combination with the quadrature.
        # It handles the scaling differently.
        # This version is also more extendable.

        g(t) = complex(radius1[1]*cos(t),radius1[2]*sin(t)) # ellipse
        gp(t) = complex(-radius1[1]*sin(t),radius1[2]*cos(t)) # derivative
        Tv(λ) = local_linsolve(T(λ),V)
        f(t) = Tv(g(t))*gp(t)/(2im*pi)

        gv=Vector{Function}(undef,2*K)
        for k=0:(2*K-1)
            gv[k+1]= s -> g(s)^k;
        end

        # Call the integrator
        Shat=integrate_interval(MIntegrator, ComplexF64,
                                f,gv,0,2*pi,N,logger )

    else
        error("Unknown Shat_mode: $Shat_mode")
    end

    for k=0:(2*K-1)
        Mhat[:,:,k+1]=U'*Shat[:,:,k+1] # Step 4: Mhat=U'*Shat
    end

    # Construct H-matrices:
    push_info!(logger,"Computing Hhat and Hhat^{<}")
    m=K*L;
    Hhat=zeros(T,m,m)   # Hhat
    Hhat2=zeros(T,m,m)  # Hhat^{<}
    for i=1:K
        for j=1:K
            Hhat[(i-1)*L .+ (1:L), (j-1)*L .+ (1:L)] = Mhat[:,:,i+j-2+1];
            Hhat2[(i-1)*L .+ (1:L), (j-1)*L .+ (1:L)] = Mhat[:,:,i+j-1+1];
        end
    end

    # Extraction more similar to Algorithm 1
    # in https://arxiv.org/pdf/1510.02572.pdf
    push_info!(logger,"Computing SVD prepare for eigenvalue extraction ",continues=true)
    F=svd(Hhat)
    UU=F.U;
    SS=F.S;
    VV=F.V

    # rank_drop_tol = δ in reference
    pp =   count( SS/SS[1] .> rank_drop_tol);
    mprime=pp; # To make closer to notation in paper

    push_info!(logger," mprime=$mprime");

    # Pick relevant eigvecs
    UU_H1 = UU[:,1:mprime]
    VV_H1 = VV[:,1:mprime]

    # Step 7: Project the moment matrices:
    Hhat_mprime = UU_H1'*Hhat*VV_H1;
    Hhat2_mprime = UU_H1'*Hhat2*VV_H1;

    # Step 8:
    (xi,X)=eigen(Hhat2_mprime,Hhat_mprime)

    # Step 10: Extract eigpair

    # Compute S-matrix (by reshaping parts of Shat-tensor)
    S=zeros(T,n,L*K)
    for j=0:(K-1)
        S[:,j*L .+ (1:L)]=Shat[:,:,j+1]
    end
    V=S*VV_H1*X;

    # Reverse the shift

    Shat_mode == :JSIAM  ? factor=radius : factor=1
    λ=σ .+ factor*xi

    return λ,V
end


text:
export getlb, getub, sosdata, getsmp, hassmp, unstable_periodic_switchings

mutable struct Lyapunov{PT <: HybridSystems.StateProperty,
                        DT <: HybridSystems.TransitionProperty}
    d::Int
    soslb::Float64
    dual::DT
    sosub::Float64
    primal::PT
end

mutable struct SOSData{S, TT,
                       XT <: HybridSystems.StateProperty,
                       PT <: HybridSystems.StateProperty,
                       DT <: HybridSystems.TransitionProperty}
    x::XT
    lb::Float64
    ub::Float64
    # There will typically only be lyapunov for small d so a dictionary would be overkill
    lyaps::Vector{Union{Nothing, Lyapunov{PT, DT}}}
    smp::Union{Nothing, DiscretePeriodicSwitching{S, TT}}
    # Unstable periodic switchings
    upsw::Union{Nothing, Set{DiscretePeriodicSwitching{S, TT}}}
end
function SOSData{S, TT, XT, PT, DT}(s::S) where {S, TT,
                                                 XT <: HybridSystems.StateProperty,
                                                 PT <: HybridSystems.StateProperty,
                                                 DT <: HybridSystems.TransitionProperty}
    y = HybridSystems.state_property(s, Vector{PolyVar{true}})::XT
    for st in states(s)
        @polyvar x[1:statedim(s, st)]
        y[st] = x
    end
    lyaps = Union{Nothing, Lyapunov{PT, DT}}[]
    SOSData{S, TT, XT, PT, DT}(y, 0, Inf, lyaps, nothing, nothing)
end

const PolynomialLyapunov{T} = SetProg.Sets.PolySet{T, SetProg.Sets.MonoBasis, T}
const MeasureLyapunov{T} = MultivariateMoments.MomentMatrix{T, SetProg.Sets.MonoBasis}

function SOSData(s::AbstractDiscreteSwitchedSystem)
    S = typeof(s)
    TT = transitiontype(s)
    XT = HybridSystems.state_property_type(S, Vector{PolyVar{true}})
    PT = HybridSystems.state_property_type(S, PolynomialLyapunov{Float64})
    DT = HybridSystems.transition_property_type(S, MeasureLyapunov{Float64})
    SOSData{S, TT, XT, PT, DT}(s)
end

const sosdatakey = :SwitchOnSafetyData

function sosdata(s::AbstractSwitchedSystem)
    if !haskey(s.ext, sosdatakey)
        s.ext[sosdatakey] = SOSData(s)
    end
    s.ext[sosdatakey]
end

#states(s::AbstractSwitchedSystem) = 1:1
#nstates(s::AbstractSwitchedSystem) = 1
#statedim(s::AbstractSwitchedSystem, i::Int) = s.n[i]
# using MP.variables triggers the invalid range update Julia v0.6 bug
MultivariatePolynomials.variables(s::SOSData, state::Int) = s.x[state]
MultivariatePolynomials.variables(s::AbstractSwitchedSystem, state::Int) = variables(sosdata(s), state)

getlyaps(s::SOSData) = s.lyaps

getlb(s::SOSData) = s.lb
function updatelb!(s::SOSData, lb, smp=nothing)
    s.lb = max(s.lb, lb)
    lb
end

getub(s::SOSData) = s.ub
function updateub!(s::SOSData, ub)
    s.ub = min(s.ub, ub)
    ub
end

function updateb!(s::AbstractSwitchedSystem, lb, ub)
    updatelb!(s, lb), updateub!(s, ub)
end

function updatesmp!(s::SOSData, smp::AbstractPeriodicSwitching)
    updatelb!(s, smp.growthrate)
    if s.smp === nothing || isbetter(smp, s.smp)
        s.smp = smp
    end
    smp
end

function notifyperiodic!(s::SOSData, psw::AbstractPeriodicSwitching, tol=sqrt(eps(Float64)))
    if s.upsw !== nothing && psw.growthrate+tol >= 1
        if !(psw in s.upsw)
            push!(s.upsw, psw)
        end
    end
end

function compute_unstable_periodic_switchings!(s::SOSData{S, TT}) where {S, TT}
    s.upsw = Set{DiscretePeriodicSwitching{S, TT}}()
end
unstable_periodic_switchings(s::SOSData) = s.upsw

hassmp(s::SOSData) = s.smp !== nothing
function getsmp(s::SOSData)
    if !hassmp(s)
        error("No smp found")
    end
    return s.smp
end

for f in (:getsmp, :hassmp, :unstable_periodic_switchings, :notifyperiodic!, :updatesmp!, :updateub!, :getub, :updatelb!, :getlyaps, :getlb)
    @eval begin
        $f(s::AbstractSwitchedSystem, args...) = $f(sosdata(s), args...)
    end
end


text:
# This file contains code that helps to query from the registry to determine the format

## Code dealing with Blocks
@with_kw struct Block
    IMAX::Int64
    JMAX::Int64
    KMAX::Int64
    X::Array{Float64,3}
    Y::Array{Float64,3}
    Z::Array{Float64,3}
end

""" Reads a chunk X, Y, or Z variable from a plot3d file
    Args
        f: IO
        IMAX:
        JMAX:
        KMAX
    Returns
        Array 3 dimensions for either X,Y, or Z
    """
function read_plot3D_chunk(f::IO,IMAX::Int64,JMAX::Int64,KMAX::Int64)::Array{Float64,3}
    A = Array{Float64,3}(undef,IMAX,JMAX,KMAX)
    for k in 1:KMAX
        for j in 1:JMAX
            for i in 1:IMAX
                A[i,j,k] = Float64(read(f,Float32))
            end
        end
    end
    return A
end

""" Reads a multi-block Plot3D File
    Args
        filename (string)
"""
function read_multi_block(filename::String)

    blocks = Block[]
    if isfile(filename)
        io = open(filename,"r")
        IMAX = Int64[]; JMAX = Int64[]; KMAX = Int64[]
        NBLOCKS = Int64(read(io,Int32))
        for i in 1:NBLOCKS
            append!(IMAX,Int64(read(io,Int32)))
            append!(JMAX,Int64(read(io,Int32)))
            append!(KMAX,Int64(read(io,Int32)))
        end
        for i in 1:NBLOCKS
            X = read_plot3D_chunk(io,IMAX[i],JMAX[i],KMAX[i])
            Y = read_plot3D_chunk(io,IMAX[i],JMAX[i],KMAX[i])
            Z = read_plot3D_chunk(io,IMAX[i],JMAX[i],KMAX[i])
            b = Block(IMAX=IMAX[i],JMAX=JMAX[i],KMAX=KMAX[i],X=X,Y=Y,Z=Z)
            push!(blocks,b)
        end
        close(io)
    end
    return blocks
end
## End Blocks


text:
using Fronts
using Documenter

makedocs(;
    modules=[Fronts],
    authors="Gabriel S. Gerlero",
    repo="https://github.com/gerlero/Fronts.jl/blob/{commit}{path}#L{line}",
    sitename="Fronts.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gerlero.github.io/Fronts.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Equations" => "equations.md",
        "Problems" => "problems.md",
        "Solving" => "solvers.md",
        "Solutions" => "solution.md",
        "Boltzmann transformation" => "boltzmann.md",
        "Inverse problems" => "inverse.md",
        "Unsaturated flow models" => "PorousModels.md",
    ],
)

deploydocs(;
    repo="github.com/gerlero/Fronts.jl",
)


text:
export inpolygon

"""
Return whether the given point (px, py) is inside or on the edge of the
polygon defined by vertices (polyx, polyy). Uses Hormann-Agathos
"Point in Polygon" algorithm to support arbitrary polygons.
Implementation inspired by 'isinside' in Luxor
(see https://github.com/JuliaGraphics/Luxor.jl/blob/master/src/polygons.jl)
"""
function inpolygon(px, py, polyx, polyy)
    if !isfinite(px) || !isfinite(py)
        return false
    end

    c = false
    @inbounds for idx in 1:length(polyx)
        q1x = polyx[idx]
        q1y = polyy[idx]

        nextidx = 1 + mod(idx, length(polyx))
        q2x = polyx[nextidx]
        q2y = polyy[nextidx]

        if q1x == px && q1y == py
            return true # on vertex
        end
        if q2y == py
            if q2x == px
                return true # on vertex
            elseif q1y == py && (q2x > px) == (q1x < px)
                return true # on edge
            end
        end
        if (q1y < py) != (q2y < py) # crossing
            if q1x >= px
                if q2x > px
                    c = !c
                else
                    det = det3p(q1x, q1y, q2x, q2y, px, py)
                    if isapprox(0, det)
                        return true # on edge
                    elseif (det > 0) == (q2y > q1y)
                        c = !c
                    end
                end
            elseif q2x > px
                det = det3p(q1x, q1y, q2x, q2y, px, py)
                if isapprox(0, det)
                    return true # on edge
                elseif (det > 0) == (q2y > q1y)
                    c = !c
                end
            end
        end
    end
    return c
end

function det3p(q1x, q1y, q2x, q2y, px, py)
    (q1x - px) * (q2y - py) - (q2x - px) * (q1y - py)
end


text:
using Twitter, App

twitterauth(App.TWITTER_TOKENS...)

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
# vector kernels for mapping

function vecmap!(n::Int, dst::Array, idst::Int, fun::Functor{1}, a::Array, ia::Int)
    idst_end = idst + n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, a[ia])
        idst += 1
        ia += 1
    end
end

function vecmap!(n::Int, dst::Array, idst::Int, sdst::Int, fun::Functor{1}, a::Array, ia::Int, sa::Int)
    idst_end = idst + sdst * n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, a[ia])
        idst += sdst
        ia += sa
    end
end

function vecmap!(n::Int, dst::Array, idst::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    idst_end = idst + n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        idst += 1
        ia += 1
        ib += 1
    end
end

function vecmap!(n::Int, dst::Array, idst::Int, sdst::Int, fun::Functor{2}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                            b::ArrayOrNum, ib::Int, sb::Int)
    idst_end = idst + sdst * n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia), getvalue(b, ib))
        idst += sdst
        ia += sa
        ib += sb
    end
end

function vecmap!(n::Int, dst::Array, idst::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int,
                                                                 b::ArrayOrNum, ib::Int,
                                                                 c::ArrayOrNum, ic::Int)
    idst_end = idst + n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        idst += 1
        ia += 1
        ib += 1
        ic += 1
    end
end

function vecmap!(n::Int, dst::Array, idst::Int, sdst::Int, fun::Functor{3}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                            b::ArrayOrNum, ib::Int, sb::Int,
                                                                            c::ArrayOrNum, ic::Int, sc::Int)
    idst_end = idst + sdst * n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia), getvalue(b, ib), getvalue(c, ic))
        idst += sdst
        ia += sa
        ib += sb
        ic += sc
    end
end

function vecmapdiff!(n::Int, dst::Array, idst::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, b::ArrayOrNum, ib::Int)
    idst_end = idst + n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        idst += 1
        ia += 1
        ib += 1
    end
end

function vecmapdiff!(n::Int, dst::Array, idst::Int, sdst::Int, fun::Functor{1}, a::ArrayOrNum, ia::Int, sa::Int,
                                                                                b::ArrayOrNum, ib::Int, sb::Int)
    idst_end = idst + sdst * n
    @inbounds while idst < idst_end
        dst[idst] = evaluate(fun, getvalue(a, ia) - getvalue(b, ib))
        idst += sdst
        ia += sa
        ib += sb
    end
end


text:
@assert [[1;2];[3;4]] == [1;2;3;4]
s = """
r = [[1;2];[3;4]]
"""
msgs = lintstr(s)
@test msgs[1].code == :W444
@test contains(msgs[1].message, "nested vcat is treated as a 1-dimensional array")

s = """
r = [[1,2],[3,4]]
"""
msgs = lintstr(s)
@test isempty(msgs)

@assert [[1 2] [3 4]] == [1 2 3 4]
s = """
r = [[1 2]  [3 4]]
"""
msgs = lintstr(s)
@test msgs[1].code == :W445
@test contains(msgs[1].message, "nested hcat is treated as a 1-row horizontal array of " *
    "dim=2")

s = """
x = Any[[1,2],[7,8]]
y = Array[[1,2],[3,4]]
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Array{Float64,2})
    y = x[1, 2, 3]
    y
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E436
@test contains(msgs[1].message, "more indices than dimensions")

s = """
function f(x::Array{Float64,2})
    x[1, 2, 3]
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E436
@test contains(msgs[1].message, "more indices than dimensions")

s = """
function f(x::Array{Float64,2})
    y = x[Colon(), 1]
    for i in y
        println(i)
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x::Array{Float64,2})
    y = x[1,1]
    @lintpragma("Info type y")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(y) == Float64")

s = """
function f(t)
    x1 = zeros(1, 2)
    x2 = zeros(Int64, 2, 2)
    x3 = zeros(t, 2, 2)
    x4 = zeros(x1)
    @lintpragma("Info type x1")
    @lintpragma("Info type x2")
    @lintpragma("Info type x3")
    @lintpragma("Info type x4")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x1) == Array{Float64,2}")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(x2) == Array{Int64,2}")
if VERSION ≥ v"0.6-"
    @test msgs[3].code == :I271
    @test contains(msgs[3].message, "typeof(x3) == $Array")
end
@test msgs[4].code == :I271
@test contains(msgs[4].message, "typeof(x4) == Array{Float64,2}")

# more array function
s = """
function f(t::Array{Int64,2}, m, n)
    x2 = reshape(t, 1)
    x3 = reshape(t, (1,2))
    x4 = reshape(m, (1,2))
    x6 = reshape(t, 1, 2)
    x7 = t'
    x8 = (1, 2)
    @lintpragma("Info type x2")
    @lintpragma("Info type x3")
    @lintpragma("Info type x4")
    @lintpragma("Info type x6")
    @lintpragma("Info type x7")
    @lintpragma("Info type x8")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x2) == Array{Int64,1}")
@test msgs[2].code == :I271
@test contains(msgs[2].message, "typeof(x3) == Array{Int64,2}")
@test msgs[3].code == :I271
@test contains(msgs[3].message, "typeof(x4) == Any")
@test msgs[4].code == :I271
@test contains(msgs[4].message, "typeof(x6) == Array{Int64,2}")
@test msgs[5].code == :I271
@test contains(msgs[5].message, "typeof(x7) == Array{Int64,2}")

s = """
function f(a::Array{Float64})
    x = a[1,2]
    @lintpragma("Info type x")
    return x
end
"""
msgs = lintstr(s)
# it could be Float64, or it could be an array still!
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(x) == Any")

s = """
s = "abcdef"
s = s[chr2ind(s,2) :end]
"""
msgs = lintstr(s)
@test msgs[1].code == :I681
@test contains(msgs[1].message, "ambiguity of :end as a symbol vs as part of a range")

s = """
s = "abcdef"
sndlast = s[end -1]
"""
msgs = lintstr(s)
@test msgs[1].code == :I682
@test contains(msgs[1].message, "ambiguity of `[end -n]` as a matrix row vs index [end-n]")

s = """
function f()
    x1 = zeros(100, 100)
    x2 = Array(Float64, (100, 100))
    x1[1, 1]
    x2[1, 1]
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(y::Array{Float64, 3}, x1::Int64)
    reshape(y[Colon(), x1, Colon()]', size(y, 1), size(y, 3)')
end
"""
msgs = lintstr(s)
@test isempty(msgs)


text:
module JuliaChem

Base.include(@__MODULE__,"../deps/src/jeri.jl")
Base.include(@__MODULE__, "modules/JCModules.jl")

Base.include(@__MODULE__, "basis/JCBasis.jl")
Base.include(@__MODULE__, "io/JCInput.jl")
Base.include(@__MODULE__, "molecule/JCMolecule.jl")
Base.include(@__MODULE__, "rhf/JCRHF.jl")
Base.include(@__MODULE__, "runtime/JCRuntime.jl")

end

#================================================#
#== we want to precompile all involved modules ==#
#================================================#
#=
if (isfile("../snoop/precompile_Base.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Base.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Blosc.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Blosc.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Blosc_jll.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Blosc_jll.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Core.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Core.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_DocStringExtensions.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_DocStringExtensions.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_HDF5.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_HDF5.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_HDF5_jll.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_HDF5_jll.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Libdl.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Libdl.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Lz4_jll.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Lz4_jll.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_MPI.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_MPI.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Parsers.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Parsers.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Pkg.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Pkg.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Requires.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Requires.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_unknown.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_unknown.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Zlib_jll.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Zlib_jll.jl")
    _precompile_()
end
if (isfile("../snoop/precompile_Zstd_jll.jl"))
    Base.include(@__MODULE__, "../snoop/precompile_Zstd_jll.jl")
    _precompile_()
end
=#


text:
transition(domain::GroundDomain, state::State, action::Term; options...) =
    execute(domain, state, action; options...)
transition!(domain::GroundDomain, state::State, action::Term; options...) =
    execute!(domain, state, action; options...)


text:
import YAML

include(joinpath(@__DIR__, "common.jl"))

function generate(platform::Platform)
    commands = """
    echo "--- Print kernel information"
    uname -a

    echo "--- Print CPU information"
    # These machines have multiple cores. However, it should be sufficient to
    # just print the information for one of the cores.
    sed -n '1,/^\$\$/p' /proc/cpuinfo

    if [[ "$(platform.arch)" == "aarch64" ]]; then
      echo "--- Patch glibc host environment"
      curl -LO https://github.com/JuliaBinaryWrappers/DebianGlibc_jll.jl/releases/download/DebianGlibc-v2.33.0%2B1/DebianGlibc.v2.33.0.aarch64-linux-gnu.tar.gz
      tar -C / -xf DebianGlibc.v2.33.0.aarch64-linux-gnu.tar.gz
    fi

    echo "--- Generate build environment"
    cmake --version
    rm -rf obj
    mkdir obj
    cd obj
    cmake ..

    echo "--- Build"
    make --output-sync -j\$\${JULIA_CPU_THREADS:?}

    echo "--- Test"
    mkdir -p Testing/Temporary
    mv ../.buildkite/CTestCostData.txt Testing/Temporary
    julia ../.buildkite/capture_tmpdir.jl ctest --output-on-failure -j\$\$(expr \$\${JULIA_CPU_THREADS:?} - 2)
    """
    job_label = "Test $(platform.arch)"
    job_key = "test-$(platform.arch)"
    yaml = Dict(
        "steps" => [
            Dict(
                "label" => job_label,
                "key" => job_key,
                "timeout_in_minutes" => 45,
                "agents" => Dict(
                    "sandbox.jl" => "true",
                    "queue" => "juliaecosystem",
                    "arch" => "$(platform.arch)",
                    "os" => "linux",
                ),
                "commands" => commands,
                "plugins" => [
                    Dict(
                        "JuliaCI/julia#v1" => Dict(
                            "persist_depot_dirs" => "packages,artifacts,compiled",
                            "version" => "1.7",
                        ),
                    ),
                    Dict(
                        "staticfloat/sandbox#v1" => Dict(
                            "rootfs_treehash" => "$(platform.rootfs_treehash)",
                            "verbose" => true,
                            "rootfs_url" => "https://github.com/JuliaCI/rootfs-images/releases/download/$(platform.rootfs_tag)/rr.$(platform.arch).tar.gz",
                            "workspaces" => ["/cache:/cache"],
                        ),
                    ),
                ],
                "soft_fail" => "$(platform.allow_fail)",
            ),
        ],
    )
    if platform.commit_status
        let
            notify = [
                Dict(
                    "github_commit_status" => Dict(
                        "context" => job_key,
                    ),
                ),
            ]
            only(yaml["steps"])["notify"] = notify
        end
    end
    yaml_path = get_yaml_path(platform)
    rm(yaml_path; force = true)
    YAML.write_file(yaml_path, yaml)
end

generate.(platforms)

max_stars_repo_name: mozilla/rr

max_issues_repo_name: rocallahan/rr

max_forks_repo_name: rocallahan/rr


text:
include("../src/components/welfare_component.jl")

@testset "welfare" begin

    Precision = 1.0e-11
    T = 60
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, welfare, :welfare)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :welfare, :CPC, read_params(f, "B126:BI126", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :welfare, :l, p[:l])
    set_param!(m, :welfare, :elasmu, p[:elasmu])
    set_param!(m, :welfare, :rr, p[:rr])
    set_param!(m, :welfare, :scale1, p[:scale1])
    set_param!(m, :welfare, :scale2, p[:scale2])

# Run the one-component model
    run(m)

# Extract the generated variables
    CEMUTOTPER  = m[:welfare, :CEMUTOTPER]
    PERIODU     = m[:welfare, :PERIODU]
    UTILITY     = m[:welfare, :UTILITY]

# Extract the true values
    True_CEMUTOTPER = read_params(f, "B129:BI129", T)
    True_PERIODU    = read_params(f, "B128:BI128", T)
    True_UTILITY    = read_params(f, "B130")

# Test that the values are the same
    @test maximum(abs, CEMUTOTPER .- True_CEMUTOTPER) ≈ 0. atol = Precision
    @test maximum(abs, PERIODU .- True_PERIODU) ≈ 0. atol = Precision
    @test abs(UTILITY - True_UTILITY) ≈ 0. atol = Precision

end


text:
using LightGraphs
using GraphPlot
using Plots
using GraphRecipes

g = LightGraphs.grid([20, 1], periodic=true)

n_start = copy(nv(g))
for v in 1:n_start
    add_vertex!(g)
    add_edge!(g, n_start + v, v)
end
n_new = copy(nv(g))
for v in (n_start + 1):n_new
    add_edge!(g, v, ((v + 1) % n_start) + n_start)
end

gplot(g, layout = spring_layout)
graphplot(g, dim=3, markersize=2, curvature_scalar=0.01)


text:
using RPOMDPModels, RPOMDPs, RPOMDPToolbox
using RobustValueIteration
using SimpleProbabilitySets
using DataFrames, ProgressMeter, CSV, BenchmarkTools, StatsBase
const RPBVI = RobustValueIteration
TOL = 1e-6

# set experiment parameters
sname = "assessment"
sversion = "7.1"
disc = 0.9
maxiter = 45
nbrand = 40
nreps = 10
nsteps = 45
nrows = 3

# intialize problems
ip = CyberIPOMDP(disc)
rip = CyberRIPOMDP(disc)
tip = CyberTestIPOMDP(disc)

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    (m - hw, m + hw)
end

# # select belief
# srand(7971023)
# s1 = zeros(27); s1[14] = 1;
# s2 = zeros(27); s2[1] = 1;
# s32 = zeros(27); s32[1] = 0.5; s32[2] = 0.5;
# s42 = zeros(27); s42[1] = 0.5; s42[4] = 0.5;
# s52 = zeros(27); s52[1] = 0.5; s52[10] = 0.5;
# s33 = zeros(27); s33[1] = 0.5; s33[3] = 0.5;
# s43 = zeros(27); s43[1] = 0.5; s43[7] = 0.5;
# s53 = zeros(27); s53[1] = 0.5; s53[19] = 0.5;
# s6 = zeros(27); s6[1] = 0.5; s6[6] = 0.5;
# s7 = zeros(27); s7[1] = 0.5; s7[22] = 0.5;
# ss = [s1, s2, s32, s42, s52, s33, s43, s53, s6, s7]
# nS = length(states(ip))
# bs = Vector{Vector{Float64}}(length(ss) + nbrand)
# for i = 1:length(ss)
#     bs[i] = ss[i]
# end
# for i = (length(ss)+1):(length(ss) + nbrand)
#     @show i
#     bs[i] = psample(zeros(nS), ones(nS))
# end
# push!(bs, vcat(fill(0.0, nS - 1), 1.0))
# push!(bs, fill(1/nS, nS))
# bsmatrix = Array{Float64}(length(bs), length(bs[1]))
# for i = 1:length(bs)
#     for j = 1:length(bs[i])
#         bsmatrix[i,j] = bs[i][j]
#     end
# end
# pathbp = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
# fnbp = string("belief_points_", sname, "_", sversion, ".csv")
# CSV.write(joinpath(pathbp, fnbp), DataFrame(bsmatrix))

# load bs (if previous bs was large)
pathbp = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
fnbp = string("belief_points_", sname, "_7.0.csv")
bsdf = CSV.read(joinpath(pathbp, fnbp))
bsmatrix = convert(Array, bsdf[:,:])
numrowbs, numcolbs = size(bsmatrix)
bs = Vector{Vector{Float64}}(numrowbs)
for i = 1:numrowbs
    bs[i] = bsmatrix[i,:]
end
length(bs) == nbrand + 12

# intialize solver
solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)

# solve
srand(5917293)
@time solip = RPBVI.solve(solver, ip, verbose = true);
@time solrip = RPBVI.solve(solver, rip, verbose = true);
@time soltip = RPBVI.solve(solver, tip, verbose = true);

# calculate values
e1 = zeros(27); e1[1] = 1
println("Standard Value: ", policyvalue(solip, e1))
println("Robust Value: ", policyvalue(solrip, e1))
println("Off Nominal Precise Value: ", policyvalue(soltip, e1))

# ipomdp and ripomdp actions for some interesting states
actionind = ["2,2,2", "1,1,1", "1,1,1 - 1,1,2", "1,1,1 - 1,2,1", "1,1,1 - 2,1,1",
        "1,1,1 - 1,1,3" ,"1,1,1 - 1,3,1", "1,1,1 - 3,1,1",
        "1,1,1 - 1,2,3", "1,1,1 - 3,2,1"]
dbsip = [DiscreteBelief(ip, states(ip), s) for s in ss]
asip = [action(solip, db) for db in dbsip]
dbsrip = [DiscreteBelief(rip, states(rip), s) for s in ss]
asrip = [action(solrip, db) for db in dbsrip]
dbstip = [DiscreteBelief(tip, states(tip), s) for s in ss]
astip = [action(soltip, db) for db in dbstip]
actiondata = DataFrame(Belief = actionind, StdAction = asip,
    RobustAction = asrip, OffNominalAction = astip)
# actiondata = DataFrame(Belief = actionind, StdAction = asip, OffNominalAction = astip)
@show actiondata

# # save policies
# policyname = "off_nominal"
# bpversion = "7.0"
# policyversion = sversion
# sol = soltip
# prob_policy = tip
# alphamatrix = Array{Float64}(length(bs), length(bs[1]))
# for i = 1:length(bs)
#     for j = 1:length(bs[i])
#         alphamatrix[i,j] = sol.alphas[i][j]
#     end
# end
# actionmatrix = action_index.(prob_policy, sol.action_map)
# pathbp = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
# fnalphas = string("policy_alphas_", policyname, "_bp", bpversion, "_results", sversion, ".csv")
# fnactions = string("policy_actions_", policyname, "_bp", bpversion, "_results", sversion, ".csv")
# CSV.write(joinpath(pathbp, fnalphas), DataFrame(alphamatrix))
# CSV.write(joinpath(pathbp, fnactions), DataFrame(action_ind = actionmatrix))

# set sim structures
nreps = 25
nsteps = 100
bus = [updater(solip), updater(solrip), updater(soltip)]
binits = [SparseCat(states(ip), initial_belief(ip)),
    SparseCat(states(rip), initial_belief(rip)),
    SparseCat(states(tip), initial_belief(tip))]
sols = [solip, solrip, soltip]
psim = RolloutSimulator(max_steps = nsteps)
simvals = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated values
simps = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated percent correct
ves = Array{Float64}(nrows) # expected values
vms = Array{Float64}(nrows) # mean of sim values
vss = Array{Float64}(nrows) # std dev of sim values
vci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim values
pms = Array{Float64}(nrows) # mean percent correct
pss = Array{Float64}(nrows) # std percent correct
pci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim percent correct

# run simluation
simprob = tip
if simprob == tip
    simprobnames = fill("Off Nominal", nrows)
elseif simprob == rip
    simprobnames = fill("Robust", nrows)
elseif simprob == ip
    simprobnames = fill("Nominal", nrows)
else
    println("Error: Unexpected simprob")
    simprobnanmes = fill("Undefined", nrows)
end
simdynamics = (simprob == rip) ? :worst : :precise
rseed = 92378432
sinit = rand(psim.rng, initial_state_distribution(simprob))
for i in 1:nrows
    srand(rseed)
    println("Run $i...")
    @showprogress 1 "Simulating value..." for j = 1:nreps
        if simdynamics == :worst
            sv, sp = simulate_worst(psim, simprob, sols[i], bus[i], binits[i], sinit, solrip.alphas)
        else
            sv, sp = simulate(psim, simprob, sols[i], bus[i], binits[i], sinit)
        end
        simvals[i][j] = sv
        simps[i][j] = sp
    end
    vms[i] = mean(simvals[i])
    vss[i] = std(simvals[i])
    vci[i] = meanci(simvals[i])
    pms[i] = mean(simps[i])
    pss[i] = std(simps[i])
    pci[i] = meanci(simps[i])
end

ves = [policyvalue(solip, e1), policyvalue(solrip, e1),
        policyvalue(soltip, e1)]
rdata = DataFrame(ID = ["Nominal", "Robust","OffNominal"], ExpValue = ves,
            SimMean = vms, SimStd = vss, SimCI = vci,
            SimPercentMean = pms, SimPercentStd = pss, SimPercentCI = pci,
            SimProb = simprobnames)
simdata = DataFrame(NominalSim = simvals[1], RobustSim = simvals[2],
            OffNominal = simvals[3])
@show rdata

sversion = "8.1"
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
fnactions = string("exp_actions_", sname, "_", sversion, ".csv")
fnresults = string("exp_results_", sname, "_", sversion, ".csv")
fnsim = string("exp_sim_values_", sname, "_", sversion, ".csv")
CSV.write(joinpath(path, fnactions), actiondata)
CSV.write(joinpath(path, fnresults), rdata)
CSV.write(joinpath(path, fnsim), simdata)

println("Data Collection Complete.")

max_stars_repo_name: ajkeith/RobustInfoPOMDP

max_issues_repo_name: ajkeith/RobustInfoPOMDP

max_forks_repo_name: ajkeith/RobustInfoPOMDP


text:
module UInt24s

using BitIntegers
#using ..UInt12s

export UInt24

BitIntegers.@define_integers 24

const DEFAULT_TYPE = UInt16

# upper_and_lower(x::UInt24, ::Type{T} ) where T = ( T(x & 0xfff), T(x >> 12) )
first(::Type{T}, x::UInt24) where T = T(x & 0xfff)
last(::Type{T}, x::UInt24) where T = T(x >> 12)
Base.first(x::UInt24) = first(DEFAULT_TYPE, x)
Base.last(x::UInt24) = last(DEFAULT_TYPE, x)

end


text:
function docs(args::Dict)
  docs()
end

function docs()

  package_name = rsplit(pwd(), "/"; limit=2)[2]

  package_dir = rsplit(dirname(@__FILE__), "/"; limit=9)[1]

  target_name = replace(package_name, ".jl", "")

  target_path = "$package_dir/$target_name"

  if target_path != pwd()
    initial_dir = pwd()

    cd("..")
    cp(package_name, target_path, remove_destination=true)
    cd(initial_dir)
  end

  include("docs/make.jl")

end


text:
export LinearBackwardEulerSolver, AbstractBackwardEulerSolver

"""
    op! = EulerOperator(f!, ϵ)

Construct a linear operator which performs an explicit Euler step ``Q + α
f(Q)``, where `f!` and `op!` both operate inplace, with extra arguments passed
through, i.e.
```
op!(LQ, Q, args...)
```
is equivalent to
```
f!(dQ, Q, args...)
LQ .= Q .+ ϵ .* dQ
```
"""
mutable struct EulerOperator{F, FT}
    f!::F
    ϵ::FT
end

function (op::EulerOperator)(LQ, Q, args...)
    op.f!(LQ, Q, args..., increment = false)
    @. LQ = Q + op.ϵ * LQ
end

abstract type AbstractBackwardEulerSolver end

"""
    (be::AbstractBackwardEulerSolver)(Q, Qhat, α, param, time)

Each concrete implementations of `AbstractBackwardEulerSolver` should provide a
callable version which solves the following system for `Q`
```
    Q = Qhat + α f(Q, param, time)
```
where `f` is the ODE tendency function, `param` are the ODE parameters, and
`time` is the current ODE time. The arguments `Q` should be modified in place
and should not be assumed to be initialized to any value.
"""
(be::AbstractBackwardEulerSolver)(Q, Qhat, α, p, t) =
    throw(MethodError(be, (Q, Qhat, α, p, t)))

"""
    Δt_is_adjustable(::AbstractBackwardEulerSolver)

Return `Bool` for whether this backward Euler solver can be updated. default is
`false`.
"""
Δt_is_adjustable(::AbstractBackwardEulerSolver) = false

"""
    update_backward_Euler_solver!(::AbstractBackwardEulerSolver, α)

Update the given backward Euler solver for the parameter `α`; see
['AbstractBackwardEulerSolver'](@ref). Default behavior is no change to the
solver.
"""
update_backward_Euler_solver!(::AbstractBackwardEulerSolver, Q, α) = nothing

"""
    setup_backward_Euler_solver(solver, Q, α, tendency!)

Returns a concrete implementation of an `AbstractBackwardEulerSolver` that will
solve for `Q` in systems of the form of
```
    Q = Qhat + α f(Q, param, time)
```
where `tendency!` is the in-place tendency function. Not the array `Q` is just
passed in for type information, e.g., `Q` the same `Q` will not be used for all
calls to the solver.
"""
setup_backward_Euler_solver(solver::AbstractBackwardEulerSolver, _...) = solver

"""
    LinearBackwardEulerSolver(::AbstractLinearSolver; isadjustable = false)

Helper type for specifying building a backward Euler solver with a linear
solver.  If `isadjustable == true` then the solver can be updated with a new
time step size.
"""
struct LinearBackwardEulerSolver{LS}
    solver::LS
    isadjustable::Bool
    LinearBackwardEulerSolver(solver; isadjustable = false) =
        new{typeof(solver)}(solver, isadjustable)
end

"""
    LinBESolver

Concrete implementation of an `AbstractBackwardEulerSolver` to use linear
solvers of type `AbstractLinearSolver`. See helper type
[`LinearBackwardEulerSolver`](@ref)
"""
mutable struct LinBESolver{FT, FAC, LS, F} <: AbstractBackwardEulerSolver
    α::FT
    factors::FAC
    solver::LS
    isadjustable::Bool
    rhs!::F
end
Δt_is_adjustable(lin::LinBESolver) = lin.isadjustable

function setup_backward_Euler_solver(lin::LinearBackwardEulerSolver, Q, α, rhs!)
    FT = eltype(α)
    factors =
        prefactorize(EulerOperator(rhs!, -α), lin.solver, Q, nothing, FT(NaN))
    LinBESolver(α, factors, lin.solver, lin.isadjustable, rhs!)
end

function update_backward_Euler_solver!(lin::LinBESolver, Q, α)
    lin.α = α
    FT = eltype(Q)
    lin.factors = prefactorize(
        EulerOperator(lin.rhs!, -α),
        lin.solver,
        Q,
        nothing,
        FT(NaN),
    )
end

function (lin::LinBESolver)(Q, Qhat, α, p, t)
    if lin.α != α
        @assert lin.isadjustable
        update_backward_Euler_solver!(lin, Q, α)
    end
    linearsolve!(lin.factors, lin.solver, Q, Qhat, p, t)
end


text:
using MathepiaData
using Test

@testset "MathepiaData.jl" begin
    # Write your tests here.
end


text:
export phase, unphase, findploidy

"""
    isbinary(filepath::String)
Returns `true` if the `filepath` is a binary file.
"""
function isbinary(filepath::String)
    !isvalid(String(Base.read(filepath, 1000)))
end

## Utilities for file reading/writing ##

#TODO change in docs
"""
    findploidy(genotypes::T where T<:AbstractVector)
Used internally in the `genepop` and `delimited` file parsers to scan the genotypes
of a sample and return the ploidy of the first non-missing locus.
"""
@inline function findploidy(genotypes::T) where T<:AbstractVector
    @inbounds Int8(length(first(skipmissing(genotypes))))
end

"""
    phase(loc::Union{String, Int}, type::DataType, digit::Int)
Takes a String of numbers or Integers and returns a typed locus appropriate for PopGen.jl as used in the
`genepop` and `csv` file parsers. Use `type` to specify output type (`Int8` or `Int16`),
and `digit` to specify the number of digits/characters used per allele in a locus.

## Examples
```
ph_locus = phase("128114", Int16, 3)
map(i -> phase(i, Int16, 3), ["112131", "211112", "001003", "516500"])
# or #
[phase(i, Int8, 2) for i in ["0101", "0103", "0202", "0103"]]
```
"""
@inline function phase(loc::T, type::DataType, digit::Int) where T<:AbstractString
    loc == "-9" || iszero(parse(Int, loc)) && return missing
    phased = map(i -> parse(type, join(i)), Iterators.partition(loc, digit))
    NTuple{length(phased), type}(sort(phased))
end

phase(loc::Missing, type::DataType, digit::Int) = missing

@inline function phase(loc::T, type::DataType, digits::T) where T<:Integer
    loc == -9 || iszero(loc) && return missing
    out = type[]
    units = 10^digits
    d,r = divrem(loc, units)
    out = type[r]
    #@inbounds push!(out, r)
    @inbounds while d != 0
        d,r = divrem(d, units)
        @inbounds push!(out, r)
    end
    return NTuple{length(out), type}(sort(out))
end

"""
    unphase(geno::T; digits::Int = 3, ploidy::Int = 2, miss::Int = 0) where T <: Genotype
Takes a `Genotype` (e.g. `(131, 94)`) and returns a string of concatenated
alleles padded with *n* number of zeroes, where *n* is given by `digits = `.
`missing` values are returned as either a string of 'digits × ploidy' zeroes (`miss = 0`)
or `"-9"` (`miss = -9`). The `ploidy` flag is only relevant for unphasing `missing` genotypes
and not used otherwise.

### Example
```
unphase((1,2,3,4), digits = 3)
"001002003004"

unphase(missing, digits = 2, ploidy = 2, miss = -9)
"-9"

unphase(missing, digits = 2, ploidy = 2, miss = 0)
"0000"
```
"""
function unphase(geno::T; digits::Int = 3, ploidy::Int = 2, miss::Int = 0) where T <: Genotype
    join(map(i -> lpad(i, digits, "0"), geno))
end

function unphase(geno::Missing; digits::Int = 3, ploidy::Int, miss::Int = 0)
    if miss == 0
        return "0"^(digits * ploidy)
    else
        return "$miss"
    end
end