text:
file = matopen("temp.mat","w")
#title:"Lock Prediction"
#write(file,"title",title)

#length:6 -> {"Actual", "Our contention model", "LR+class", "quad+class", "Dec. tree regression", "Orig. Thomasian"};
#length:5 -> {"Our contention model", "LR+class", "quad+class", "Dec. tree regression", "Orig. Thomasian"};
lenOfLegends=length(legends)
write(file,"lenOfLegends",lenOfLegends)

#=
for i=1:lenOfLegends
	write(file,string("legends",i),legends[i])
end
=#

lenOfXdata=length(Xdata)
write(file,"lenOfXdata",lenOfXdata)
for i=1:lenOfXdata
	write(file,string("Xdata",i),float(Xdata[i]))
end

lenOfYdata=length(Ydata)
write(file,"lenOfYdata",lenOfYdata)
for i=1:lenOfYdata
	write(file,string("Ydata",i),float(Ydata[i]))
end

#Xlabel = "TPS";
#Ylabel = "Total time spent acquiring row locks (seconds)";
#=
write(file,"Xlabel",Xlabel)
write(file,"Ylabel",Ylabel)
=#

lenOfMeanAbsError=length(meanAbsError)
write(file,"lenOfMeanAbsError",lenOfMeanAbsError)
for i=1:lenOfMeanAbsError
	write(file,string("meanAbsError",i),float(meanAbsError[i]))
end

lenOfMeanRelError=length(meanRelError)
write(file,"lenOfMeanRelError",lenOfMeanRelError)
for i=1:lenOfMeanRelError
	write(file,string("meanRelError",i),float(meanRelError[i]))
end

#length>0 -> {"Our contention model", "LR+class", "quad+class", "Dec. tree regression", "Orig. Thomasian"};
lenOfErrorHeader=length(errorHeader)
write(file,"lenOfErrorHeader",lenOfErrorHeader)
#=
for i=1:lenOfErrorHeader
	write(file,string("errorHeader",i),errorHeader[i])
end
=#

lenOfExtra=length(extra)
write(file,"lenOfExtra",lenOfExtra)
for i=1:lenOfExtra
	write(file,string("extra",i),extra[i])
end
close(file)

#=
file = matopen("lockP2.mat","w")
write(file,"title",title)
write(file,"legends",legends)
write(file,"Xdata",Xdata)
write(file,"Ydata",Ydata)
write(file,"Xlabel",Xlabel)
write(file,"Ylabel",Ylabel)
write(file,"meanAbsError",meanAbsError)
write(file,"meanRelError",meanRelError)
write(file,"errorHeader",errorHeader)
write(file,"extra",extra)
close(file)
=#

#=
matwrite("lockP.mat",{
	"title" => title,
	"legends" => legends,
	"Xdata" => Xdata,
	"Xlabel" => Xlabel,
	"Ylabel" => Ylabel,
	"meanAbsError" => meanAbsError,
	"meanRelError" => meanRelError,
	"errorHeader" => errorHeader,
	"extra" => extra
})
=#

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
module MultiWriters
export MultiWriter, mwprint, mwwrite, mwprintln, mwshow

struct MultiWriter <: IO
    writers::Tuple{Vararg{IO}}

    function MultiWriter(writers::IO...)
        new(writers)
    end
end

function _write(f::Function, ios, x...)
    for io in ios.writers
        f(io, x...)
    end
end

for f in [:write, :print, :println, :show]
    fnname = Symbol("mw", f)
    eval(quote
        function $(fnname)(io::MultiWriter, x...)
            _write($(f), io, x...)
        end
    end)
end

# argtypes(m::Method) = fieldtypes(m.sig)
#
# for f in [:write, :print, :show]
#     ms = methods(eval(f)).ms
#     args_list = sizehint!(Set{Symbol}(), length(ms))
#     for m in ms
#         nargs = m.nargs
#         println(m)
#         args = argtypes(m)
#         fn_args = Symbol(join(["a_$(i)::$type" for (i, type) in enumerate(args[3:end])], ","))
#         push!(args_list, fn_args)
#     end

#     for arg in args_list
#         eval(quote
#             Base.$(f)(io::MultiWriter, $(arg)) = _write($(f), io, $(arg))
#         end)
#     end
# end

end # module

ext: jl


text:
# This file defines the State type. This type defines a state that a personnel
#   member can be in.

# The State type requires the Attrition type.
requiredTypes = [ "attrition" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export State
"""
This type defines a state that a personnel member in the simulation can have. A
state in this instance is a collection of attributes and values that must be
satisfied before a personnel member is considered to be in that state.

The type contains the following fields:
* `name::String`: the name of the state.
* `requirements::Dict{String, Vector{String}}`: the attributes (key) and the
  value(s) these attributes must have (value) for the personnel member to have
  this particular state.
* `isInitial::Bool`: this flag states whether the state is an initial state or
  not.
* `stateTarget::Int`: the target number of personnel members in this state.
* `stateRetAge::Float64`: the retirement age of this state. If this is zero, the
  default retirement age is used.
* `attrScheme::Attrition`: The attrition scheme attached to this state.
* `inStateSince::Dict{String, Float64}`: a record of all the persons in the
  state, including the time they last entered the state.
* `isLockedForTransition::Dict{String, Bool}`: a record of all the persons ready
  to undergo a transition.
"""
type State

    name::String
    requirements::Dict{String, Vector{String}}
    isInitial::Bool
    stateTarget::Int
    stateRetAge::Float64
    attrScheme::Attrition
    inStateSince::Dict{String, Float64}
    isLockedForTransition::Dict{String, Bool}

    # Basic constructor.
    function State( name::String, isInitial::Bool = false )

        newState = new()
        newState.name = name
        newState.requirements = Dict{String, Vector{String}}()
        newState.isInitial = isInitial
        newState.stateTarget = 0
        newState.stateRetAge = 0.0
        newState.attrScheme = Attrition()
        newState.inStateSince = Dict{String, Float64}()
        newState.isLockedForTransition = Dict{String, Bool}()
        return newState

    end  # State( name, isInitial )

end  # type State


text:
using ProgressMeter

function main()
    out_dir  = "result"
    PARAMDAT = "PARAMDAT.json"
    fwrite   = "write"

    nch  = 2    # rho,u,v,p
    nval = 6    # rho,u,v,p + N2, N
    nre  = 2    # N2+N2 -> 2N+N2
                # N2+N  -> 2N+N
    Rd   = 287.0
    R    = 8.314

    xmax, ymax, nodes, vecAx, vecAy = read_allgrid()
    out_file_front, out_ext, restartnum, restart_file, init_small, norm_ok,
    time_integ, nt, dt, every_outnum, in_nt, dtau,
    init_rho, init_u, init_v, init_p, init_N2, init_N, specific_heat_ratio, Rd, bdcon = input_para(PARAMDAT)

    Qbase, cellxmax, cellymax, restartnum = set_initQbase(xmax, ymax, restart_file, init_rho, init_u, init_v, init_p, init_N2, init_N,
                                                        specific_heat_ratio, out_file_front, out_ext, out_dir, restartnum, Rd, nval)

    # init Delta_Qcon_hat
    volume = set_volume(nodes, cellxmax, cellymax)
    Minf   = set_Minf(bdcon, specific_heat_ratio)
    reset_write(fwrite)

    # main loop
    print("threads num : ")
    println(Threads.nthreads())

    # check bd
    check_bd(bdcon)

    #throw(UndefVarError(:x))
    # main loop
    prog = Progress(nt,1)
    for t in 1:nt
        next!(prog)

        evalnum = t+restartnum
        if time_integ == "1"
            # exlicit scheme

            #println(Qbase[1,3,:])

            Qbase    = set_boundary(Qbase, cellxmax, cellymax, vecAx, vecAy, bdcon, nval)
            Qcon     = base_to_conservative(Qbase, cellxmax, cellymax, specific_heat_ratio, nval, nch)
            Qcon_hat = setup_Qcon_hat(Qcon, cellxmax, cellymax, volume, nval, nch)

            # initial_setup
            chmu, chlambda_tr, chD, chi = chemical_value(Qbase, cellxmax, cellymax, Rd, nval, nch)

            # mu     = set_mu(Qbase,cellxmax,cellymax,specific_heat_ratio,Rd)
            # lambda = set_lambda(Qbase,cellxmax,cellymax,mu,specific_heat_ratio,Rd)

            # RHS
            # sorce term
            wdot = set_wdot(Qbase, cellxmax, cellymax, Rd, nch, nre, nval)

            #println(Qbase[1,3,:])
            #println(Qcon[1,3,:])
            # advection_term
            #E_adv_hat, F_adv_hat = AUSM(Qbase,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio)
            E_adv_hat, F_adv_hat = AUSM_plus(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, Minf, volume, nval, nch)

            #println("AUSM")

            #println(Qcon[40,101,:])
            #println(Qcon[40,102,:])
            #println(E_adv_hat[2,10,:])
            #println(E_adv_hat[35,3,:])
            #println(F_adv_hat[2,100,:])
            #println(wdot[2,2,:])
            #println(F_adv_hat[35,3,:])

            #=
            cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus = cal_jacobi(Qbase,Qcon,cellxmax,cellymax,specific_heat_ratio,vecAx,vecAy,volume)
            cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus = setup_cell_flux_hat(Qcon,cellxmax,cellymax,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus)
            E_adv_hat,F_adv_hat = FVS(cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus,cellxmax,cellymax)
            =#

            #println("FVS")
            #println(E_adv_hat[52,50,:])
            #println(E_adv_hat[35,3,:])
            #println(F_adv_hat[52,50,:])
            #println(F_adv_hat[35,3,:])

            # viscos_term
            E_vis_hat, F_vis_hat = central_diff(Qbase, Qcon, cellxmax, cellymax, chmu, chlambda_tr, chD, chi,
                                                vecAx, vecAy, specific_heat_ratio, volume, Rd, nval, nch)
            #E_vis_hat = zeros(cellxmax+1, cellymax+1, nval)
            #F_vis_hat = zeros(cellxmax+1, cellymax+1, nval)

            RHS = setup_RHS(cellxmax, cellymax, E_adv_hat, F_adv_hat, E_vis_hat, F_vis_hat, wdot, nval, volume)
            #=
            println(E_adv_hat[2,50,:])
            println(F_adv_hat[2,50,:])
            println(E_vis_hat[2,50,:])
            println(F_vis_hat[2,50,:])
            println(wdot[2,50,:])
            println(RHS[2,50,:])
            =#

            #=
            println("          ")
            println(E_adv_hat[100,2,:])
            println(F_adv_hat[100,2,:])
            println(E_adv_hat[100,3,:])
            println(F_adv_hat[100,3,:])
            println(Qcon_hat[100,2,:])
            println(RHS[100,2,:])
            =#

            #println(RHS[2,4,:])
            # time integral
            Qcon_hat = time_integration_explicit(dt, Qcon_hat, RHS, cellxmax, cellymax, nval)

            Qcon  = Qhat_to_Q(Qcon_hat, cellxmax, cellymax, volume, nval, nch)
            Qbase = conservative_to_base(Qcon, cellxmax, cellymax, specific_heat_ratio, nval, nch)

        elseif time_integ == "2"
            # Qcn
            Qbasen = copy(Qbase)
            Qconn = base_to_conservative(Qbasen,cellxmax,cellymax,specific_heat_ratio)
            Qconn_hat = setup_Qcon_hat(Qconn,cellxmax,cellymax,volume)

            Qbasem = copy(Qbase)

            for tau in 1:in_nt
                # LHS (A_adv_hat=jacobian)
                Qbasem = set_boundary(Qbasem,cellxmax,cellymax,vecAx,vecAy,bdcon)

                Qcon = base_to_conservative(Qbasem,cellxmax,cellymax,specific_heat_ratio)
                Qcon_hat = setup_Qcon_hat(Qcon,cellxmax,cellymax,volume)

                # initial_setup
                mu = set_mu(Qbasem,cellxmax,cellymax,specific_heat_ratio,Rd)
                lambda = set_lambda(Qbasem,cellxmax,cellymax,mu,specific_heat_ratio,Rd)

                # RHS
                #advection_term
                #E_adv_hat, F_adv_hat = AUSM(Qbasem,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio)
                cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus = cal_jacobi(Qbasem,Qcon,cellxmax,cellymax,specific_heat_ratio,vecAx,vecAy,volume)
                cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus = setup_cell_flux_hat(Qcon,cellxmax,cellymax,cell_Ahat_plas,cell_Ahat_minus,cell_Bhat_plas,cell_Bhat_minus)
                E_adv_hat,F_adv_hat = FVS(cell_E_hat_plas,cell_E_hat_minus,cell_F_hat_plas,cell_F_hat_minus,cellxmax,cellymax)

                # viscos_term
                E_vis_hat, F_vis_hat = central_diff(Qbasem,Qcon,cellxmax,cellymax,mu,lambda,vecAx,vecAy,specific_heat_ratio,volume,Rd)

                RHS = setup_RHS(cellxmax,cellymax,E_adv_hat, F_adv_hat,E_vis_hat, F_vis_hat)

                #println(E_adv_hat[:,50,:])
                #println(E_adv_hat[:,50,:])

                #println(RHS[:,50,:])

                #throw(UndefVarError(:x))

                # lusgs_advection_term
                A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_shig, B_beta_shig = one_wave(Qbasem,Qcon,cellxmax,cellymax,vecAx,vecAy,specific_heat_ratio,volume)
                # lusgs_viscos_term
                jalphaP, jbetaP = central_diff_jacobian(Qbasem,Qcon,cellxmax,cellymax,mu,lambda,vecAx,vecAy,specific_heat_ratio,volume)
                #jalphaP = zeros(cellxmax, cellymax)
                #jbetaP = zeros(cellxmax, cellymax)

                # LUSGS
                delta_Q = zeros(cellxmax,cellymax,4)

                ite=0
                while true
                    delta_Q_temp = copy(delta_Q)
                    delta_Q = lusgs(dt,Qcon_hat,Qconn_hat,delta_Q,A_adv_hat_p, A_adv_hat_m, B_adv_hat_p, B_adv_hat_m, A_beta_shig, B_beta_shig,jalphaP, jbetaP,RHS,cellxmax,cellymax,volume)

                    # println(RHS[:,50,:])
                    #println(delta_Q[:,50,:])

                    res = set_res(delta_Q,delta_Q_temp,cellxmax,cellymax)
                    norm2 = check_converge(res,RHS,cellxmax,cellymax,init_small)

                    ite += 1
                    if norm2[1] < norm_ok && norm2[2] < norm_ok && norm2[3] < norm_ok && norm2[4] < norm_ok
                        break
                    end
                    if ite % 100 ==0
                        println(norm2)
                    end
                end
                #println(ite)
                #println(delta_Q[:,50,:])

                for i in 2:cellxmax-1
                    for j in 2:cellymax-1
                        for k in 1:4
                            Qcon_hat[i,j,k] = Qcon_hat[i,j,k] + delta_Q[i,j,k]
                        end
                    end
                end

                Qcon = Qhat_to_Q(Qcon_hat,cellxmax,cellymax,volume)
                Qbasem = conservative_to_base(Qcon,cellxmax,cellymax,specific_heat_ratio)

                #println("xxx")
            end
            Qbase = copy(Qbasem)
        end

        if round(evalnum) % every_outnum == 0
            println("\n")
            println("nt_______________________________"*string(round(evalnum)))
            output_result(evalnum, Qbase, cellxmax, cellymax, specific_heat_ratio, out_file_front, out_ext, out_dir, Rd, nval)
        end

        check_divrege(Qbase, cellxmax, cellymax, Rd, fwrite)

    end
end

# -- main --
main()


text:
#!/usr/bin/env julia

using FromFile

@from "Helper.jl" import Helper

function take_apart(right::AbstractString)::Vector{Int}
    result::Vector{Int} = []

    parts = split(right, ",")
    for part in parts
        _, r = split(part, "=")
        a, b = split(r, "..")
        append!(result, parse(Int, a))
        append!(result, parse(Int, b))
    end
    #
    return result
end

function main()
    # fname = "example1.txt"
    # fname = "example2.txt"
    fname = "input.txt"

    bag::Set{Tuple{Int, Int, Int}} = Set()
    lines = Helper.read_lines(fname)
    for line in lines
        left, right = split(line)
        x1, x2, y1, y2, z1, z2 = take_apart(right)
        for x in x1:x2
            if !(-50 <= x <= 50)
                continue
            end
            for y in y1:y2
                if !(-50 <= y <= 50)
                    continue
                end
                for z in z1:z2
                    if !(-50 <= z <= 50)
                        continue
                    end
                    t = (x, y, z)
                    if left == "on"
                        push!(bag, t)
                    else
                        delete!(bag, t)
                    end
                end
            end
        end
    end
    println(length(bag))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
facts("Non-negative CANDECOMP") do

r = 2
T = rand_kruskal3(r, (10, 20, 30), true)

@time factors = nncp(T, r, compute_error=true)
@fact length(factors.factors) --> ndims(T)
@fact map(size, factors.factors) --> (collect(zip(size(T), (r, r, r)))...)
@fact rank(factors) --> r
@fact rel_residue(factors) --> less_than(0.05)

end


text:
module MetalheadUtils

using Images, Metalhead, Flux

export imgFromURL, normalize!
export postprocess, impreprocess, vggBlockLen
export Features

function imgFromURL(url)
	mktemp() do fn,f
		download(url, fn)
    load(fn)
	end
end

## the reverse of Metalhead.preprocess
function postprocess(x)
	μ = [0.485, 0.456, 0.406]
	σ = [0.229, 0.224, 0.225]
	xx = permutedims(x[:,:,:,1], (3, 2, 1))./255
	xxx = N0f8.(xx.*σ .+ μ)
	[RGB(xxx[1,i,j], xxx[2,i,j], xxx[3,i,j]) for i = 1:size(xxx,2), j = 1:size(xxx,3)]
end

function normalize!(x, lower, upper)
	x .= min.(x, upper)
	x .= max.(x, lower)
	x
end

## from Metalhead.preprocess
function impreprocess(im)
    # Resize such that smallest edge is 256 pixels long
    im = Metalhead.resize_smallest_dimension(im, 256)
    # Center-crop to 224x224
    im = Metalhead.center_crop(im, 224)
    im
end

const vggBlockLen = [3, 3, 5, 5, 5]

function Features(chain::Chain, X, layernums)
	[chain[1:k](X) for k in layernums]
end

function Features(nn, X, layernums)
	Features(nn.layers, X, layernums)
end

end


text:
"""
Current elements:
    :times: (quantity array 1D, numpy array 1D, or list) The times of each spike.
    :units: (quantity units) Required if :attr:`times` is a list or :class:`~numpy.ndarray`, not if it is a :class:`~quantities.Quantity`.
    :t_stop: (quantity scalar, numpy scalar, or float) Time at which :class:`SpikeTrain` ended. This will be converted to the same units as :attr:`times`. This argument is required because it specifies the period of time over which spikes could have occurred.  Note that :attr:`t_start` is highly recommended for the same reason.
    :t_start: (quantity scalar, numpy scalar, or float) Time at which :class:`SpikeTrain` began. This will be converted to the same units as :attr:`times`.  Default: 0.0 seconds.
    :name: (str) A label for the dataset.
    :description: (str) Text description.
    :file_origin: (str) Filesystem path or URL of the original data file.
    :waveforms: (quantity array 3D (spike, channel, time)) The waveforms of each spike.
    :sampling_rate: (quantity scalar) Number of samples per unit time for the waveforms.
    :left_sweep: (quantity array 1D) Time from the beginning of the waveform to the trigger time of the spike.
    :sort: (bool) If True, the spike train will be sorted by time.

"""
struct SpikeTrain
    times::Array
    t_stop::Number
    t_sart::Number
		name::String
		description::String
    waveforms::Array
    sampling_rate::Number
    left_sweep::Number
    right_sweep::Number
    sort::Bool
end


text:
##################################################
#   Override default handling (default = inline) #
##################################################

get_latex_function(r::DiffEqBase.AbstractReactionNetwork) = latexalign

get_md_function(args::DiffEqBase.AbstractReactionNetwork) = mdalign

###############################################
#         Overload environment functions      #
###############################################

"""
latexalign(r::AbstractReactionNetwork; noise=false, symbolic=true)

Generate an align environment from a reaction network.

### kwargs
- noise::Bool - output the noise function?
- symbolic::Bool - use symbolic calculation to reduce the expression?
- bracket::Bool - Surround the variables with square brackets to denote concentrations.
- clean::Bool - Clean out redundant "1*". Only useful for DiffEqBiological@v3.4.2 or earlier.
"""
function latexalign(r::DiffEqBase.AbstractReactionNetwork; bracket=false, noise=false, symbolic=false, clean=false, kwargs...)
    lhs = [Meta.parse("d$x/dt") for x in r.syms]
    if !noise
        if symbolic
            rhs = r.f_symfuncs
        else
            if clean
                rhs = clean_subtractions.(r.f_func)
            else
                rhs = r.f_func
            end
        end
    else
        vec = r.g_func
        M = reshape(vec, :, length(r.syms))
        M = permutedims(M, [2,1])
        expr_arr = Meta.parse.([join(M[i,:], " + ") for i in 1:size(M,1)])

        if symbolic
            rhs = [SymEngine.Basic(ex) for ex in expr_arr]
        else
            for i in 1:length(expr_arr)
                filter!(x -> x != 0, expr_arr[i].args)
            end
            rhs = expr_arr
        end
    end
    if bracket
        rhs = add_brackets(rhs, r.syms)
        lhs = [:(d[$x]/dt) for x in r.syms]
    end
    return latexalign(lhs, rhs; kwargs...)
end

function chemical_arrows(rn::DiffEqBase.AbstractReactionNetwork;
    expand = true, double_linebreak=false, mathjax=true, starred=false, kwargs...)
    str = starred ? "\\begin{align*}\n" : "\\begin{align}\n"
    eol = double_linebreak ? "\\\\\\\\\n" : "\\\\\n"

    mathjax && (str *= "\\require{mhchem}\n")

    backwards_reaction = false
    for (i, r) in enumerate(rn.reactions)
        if backwards_reaction
            backwards_reaction = false
            continue
        end
        str *= "\\ce{ "

        ### Expand functions to maths expressions
        rate = deepcopy(r.rate_org)
        expand && (rate = DiffEqBiological.recursive_clean!(rate))
        expand && (rate = DiffEqBiological.recursive_clean!(rate))

        ### Generate formatted string of substrates
        substrates = [latexraw("$(substrate.stoichiometry== 1 ? "" : "$(substrate.stoichiometry) * ") $(substrate.reactant)"; kwargs...) for substrate in r.substrates ]
        isempty(substrates) && (substrates = ["\\varnothing"])

        str *= join(substrates, " + ")

        ### Generate reaction arrows
        if i + 1 <= length(rn.reactions) && r.products == rn.reactions[i+1].substrates && r.substrates == rn.reactions[i+1].products
            ### Bi-directional arrows
            rate_backwards = deepcopy(rn.reactions[i+1].rate_org)
            expand && (rate_backwards = DiffEqBiological.recursive_clean!(rate_backwards))
            expand && (rate_backwards = DiffEqBiological.recursive_clean!(rate_backwards))
            str *= " &<=>"
            str *= "[" * latexraw(rate; kwargs...) * "]"
            str *= "[" * latexraw(rate_backwards; kwargs...) * "] "
            backwards_reaction = true
        else
            ### Uni-directional arrows
            str *= " &->"
            str *= "[" * latexraw(rate; kwargs...) * "] "
        end

        ### Generate formatted string of products
        products = [latexraw("$(product.stoichiometry== 1 ? "" : "$(product.stoichiometry) * ") $(product.reactant)"; kwargs...) for product in r.products ]
        isempty(products) && (products = ["\\varnothing"])
        str *= join(products, " + ")
        str *= "}$eol"
    end
    str = str[1:end-length(eol)] * "\n"

    str *= starred ? "\\end{align*}\n" : "\\end{align}\n"

    latexstr = LaTeXString(str)
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

"""
clean_subtractions(ex::Expr)

Replace additions of negative terms with subtractions.

This is a fairly stupid function which is designed for a specific problem
with reaction networks. It is neither recursive nor very general.

Note: this function was moved to DiffEqBiological and is only retained here for
compatibility with older versions.

Return :: cleaned out expression
"""
function clean_subtractions(ex::Expr)
    ex.args[1] != :+ && return ex

    term = ex.args[2]

    ### Sort out the first term
    if term isa Expr && length(term.args) >= 3 && term.args[1:2] == [:*, -1]
        result = :(- *($(term.args[3:end]...)))
    else
        result = :($term)
    end

    ### Sort out the other terms
    for term in ex.args[3:end]
        if term isa Expr && length(term.args) >= 3 && term.args[1:2] == [:*, -1]
            result = :($result - *($(term.args[3:end]...)))
        else
            result = :($result + $term)
        end
    end
    return result
end

function clean_subtractions(arg)
    @warn "It appears that you are using a version of DiffEqBiological which does not require the latexify `clean` keyword argument. The `clean=true` specification will be ignored."
    return arg
end


text:
using LightGraphs

abstract type Constraint end
abstract type OnePropagationConstraint <: Constraint end

include("variables/variables.jl")
include("core/core.jl")
include("constraints/constraints.jl")


text:
"""
    RegionGrid(
        geo :: RectRegion{ST,FT},
        lon :: Vector{<:Real},
        lat :: Vector{<:Real}
    ) -> RectGrid{FT}

Creates a rectilinear GeoRegion `RegID`.

Arguments
=========

- `geo` : A RectRegion struct type
- `lon` : A vector containing the longitude points
- `lat` : A vector containing the latitude points
"""
function RegionGrid(
    geo :: RectRegion{ST,FT},
    lon :: Vector{<:Real},
    lat :: Vector{<:Real}
) where {FT <: Real, ST <: AbstractString}

    @info "$(now()) - GeoRegions.jl - Creating a RegionGrid for the $(geo.name) GeoRegion"

    @debug "$(now()) - GeoRegions.jl - Determining indices of longitude and latitude boundaries in the given dataset ..."
    N = geo.N
    S = geo.S
    E = geo.E
    W = geo.W

    igrid = regiongrid([N,S,E,W],lon,lat);
    iN = igrid[1]; iS = igrid[2]; iE = igrid[3]; iW = igrid[4]
    nlon = deepcopy(lon)

    @info "$(now()) - GeoRegions.jl - Creating vector of latitude indices to extract ..."
    if     iN < iS; iNS = iN : iS
    elseif iS < iN; iNS = iS : iN
    else;           iNS = iN;
    end

    @info "$(now()) - GeoRegions.jl - Creating vector of longitude indices to extract ..."
    if     iW < iE; iWE = vcat(iW:iE)
    elseif iW > iE || (iW == iE && bounds[3] != bounds[4])
          iWE = vcat(iW:length(lon),1:iE); nlon[1:(iW-1)] .+= 360
    else; iWE = [iW];
    end

    while maximum(nlon) > 360; nlon .-= 360 end
    nlon = nlon[iWE]
    nlat =  lat[iNS]

    return RectGrid{FT}(igrid,iWE,iNS,nlon,nlat)

end

"""
    RegionGrid(
        geo :: PolyRegion{ST,FT},
        lon :: Vector{<:Real},
        lat :: Vector{<:Real}
    ) -> PolyGrid{FT}

Creates a rectilinear GeoRegion `RegID`.

Arguments
=========

- `geo` : A PolyRegion struct type
- `lon` : A vector containing the longitude points
- `lat` : A vector containing the latitude points
"""
function RegionGrid(
    geo :: PolyRegion{ST,FT},
    lon :: Vector{<:Real},
    lat :: Vector{<:Real}
) where {FT <: Real, ST <: AbstractString}

    @info "$(now()) - GeoRegions.jl - Creating a RegionGrid for the $(geo.name) GeoRegion"

    @debug "$(now()) - GeoRegions.jl - Determining indices of longitude and latitude boundaries in the given dataset ..."
    N = geo.N
    S = geo.S
    E = geo.E
    W = geo.W

    igrid = regiongrid([N,S,E,W],lon,lat);
    iN = igrid[1]; iS = igrid[2]; iE = igrid[3]; iW = igrid[4]
    nlon = deepcopy(lon)

    @info "$(now()) - GeoRegions.jl - Creating vector of latitude indices to extract ..."
    if     iN < iS; iNS = vcat(iN:iS)
    elseif iS < iN; iNS = vcat(iS:iN)
    else;           iNS = [iN];
    end

    @info "$(now()) - GeoRegions.jl - Creating vector of longitude indices to extract ..."
    if     iW < iE; iWE = vcat(iW:iE)
    elseif iW > iE || (iW == iE && bounds[3] != bounds[4])
          iWE = vcat(iW:length(lon),1:iE); nlon[1:(iW-1)] .+= 360
    else; iWE = [iW];
    end

    while maximum(nlon) > 360; nlon .-= 360 end

    @info "$(now()) - GeoRegions.jl - Since the $(geo.name) GeoRegion is a PolyRegion, we need to defined a mask as well ..."
    nlon = nlon[iWE]
    nlat =  lat[iNS]
    mask = Array{FT,2}(undef,length(nlon),length(nlat))
    for ilat = 1 : length(nlat), ilon = 1 : length(nlon)
        ipnt = Point2(nlon[ilon],nlat[ilat])
        if isinGeoRegion(ipnt,geo,throw=false)
              mask[ilon,ilat] = 1
        else; mask[ilon,ilat] = NaN
        end
    end

    return PolyGrid{FT}(igrid,iWE,iNS,nlon,nlat,mask)

end

function regionpoint(
    plon::Real, plat::Real,
    rlon::Vector{<:Real}, rlat::Vector{<:Real}
)

    plon = mod(plon,360); ispointinregion(plon,plat,rlon,rlat); rlon = mod.(rlon,360);
    ilon = argmin(abs.(rlon.-plon)); ilat = argmin(abs.(rlat.-plat));

    return [ilon,ilat]

end

function regiongrid(gridbounds::Vector{<:Real},rlon::Vector{<:Real},rlat::Vector{<:Real})

    N,S,E,W = gridbounds; isgridinregion(gridbounds,rlon,rlat)

    E = mod(E,360); W = mod(W,360); nlon = mod.(rlon,360);
    iN = argmin(abs.(rlat.-N)); iS = argmin(abs.(rlat.-S)); iW = argmin(abs.(nlon.-W));
    if E == W;
        if gridbounds[3] != gridbounds[4]
              if iW != 1; iE = iW - 1; else; iE = length(nlon); end
        else; iE = iW
        end
    else; iE = argmin(abs.(nlon.-E));
    end

    return [iN,iS,iE,iW]

end

function regioninfo(gridbounds::Vector{<:Real},rlon::Vector{<:Real},rlat::Vector{<:Real})

    iN,iS,iE,iW = regiongrid(gridbounds,rlon,rlat);
    nlon = deepcopy(rlon)

    if     iN < iS; iNS = vcat(iN:iS)
    elseif iS < iN; iNS = vcat(iS:iN)
    else;           iNS = [iN];
    end

    if     iW < iE; iWE = vcat(iW:iE)
    elseif iW > iE || (iW == iE && bounds[3] != bounds[4])
          iWE = vcat(iW:length(lon),1:iE); nlon[1:(iW-1)] .+= 360
    else; iWE = [iW];
    end

    while maximum(nlon); nlon .-= 360 end
    nlon = nlon[iWE]
    nlat =  lat[iNS]

    return RectGrid{FT}(igrid,iWE,iNS,nlon,nlat)

end


text:
mutable struct Snake
    head::Point
    body::Vector{Point}
    last_move::Direction

    function Snake(head::Point=Point(WIDTH ÷ 2, HEIGHT ÷ 2))

        # Create the initial body
        body = [
            head,
            Point(head.x - BLOCK_SIZE, head.y),
            Point(head.x - 2 * BLOCK_SIZE, head.y)
        ]

        return new(head, body, RIGHT)
    end
end

function Base.show(io::IO, s::Snake)
    print("$(s.head)[HEAD]")
    for p in s.body[2:end]
        print(" <- ", p)
    end
end

function move!(s::Snake, d::Direction)
    x = s.head.x
    y = s.head.y

    # Move the head
    if d == UP
        y -= BLOCK_SIZE
    elseif d == DOWN
        y += BLOCK_SIZE
    elseif d == LEFT
        x -= BLOCK_SIZE
    elseif d == RIGHT
        x += BLOCK_SIZE
    end

    # Adding the new head to the snake
    head = Point(x, y)
    s.head = head
    pushfirst!(s.body, head)

    # Updating the last move of the snake
    s.last_move = d

    # Return the old tail
    return pop!(s.body)
end

function is_collision(s::Snake, p::Union{Point, Nothing}=nothing)
    # Defaults to the snake's head if no point is given
    if p === nothing
        p = s.head
    end

    # Hits screen edges
    if p.x > WIDTH - BLOCK_SIZE ||
       p.x < 0 ||
       p.y > HEIGHT - BLOCK_SIZE ||
       p.y < 0

        return true
    end

    # Check if snake hits itself
    return p in s.body[2:end]
end


text:
using StochasticDiffEq, Test, Random, DiffEqDevTools
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich
Random.seed!(100)
dts = 1 ./2 .^(10:-1:2) #14->7 good plot

prob = prob_sde_linear_stratonovich
sim  = test_convergence(dts,prob,EulerHeun(),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,LambaEulerHeun(),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ISSEulerHeun(),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim = test_convergence(dts,prob,RKMil_General(interpretation=:Stratonovich),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim  = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich),trajectories=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,SROCK1(interpretation=:Stratonovich),trajectories=Int(2e2))
@test abs(sim.𝒪est[:l2]-1) < 0.15

sim  = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2))
@test abs(sim.𝒪est[:final]-2) < 0.20

println("Now 2D")

prob = prob_sde_2Dlinear_stratonovich

sim  = test_convergence(dts,prob,EulerHeun(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,LambaEulerHeun(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ISSEulerHeun(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

println("RKMils")

sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim = test_convergence(dts,prob,RKMil_General(interpretation=:Stratonovich),trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim  = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich),
                        trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=:Stratonovich),
                        trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=:Stratonovich),
                        trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,SROCK1(interpretation=:Stratonovich),trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,KomBurSROCK2(),trajectories=Int(2e2))
@test abs(sim.𝒪est[:final]-2) < 0.20

Random.seed!(200)
sol = solve(prob,EulerHeun(),dt=1/4)


text:
dy = iy - ix
dz = iz - ix
@inbounds for i = ix:ix-1+n
    x[i] = y[dy+i] + z[dz+i]*a
end


text:
using ChainRules: level2partition, level3partition, chol_blocked_rev, chol_unblocked_rev

@testset "Factorizations" begin
    @testset "svd" begin
        for n in [4, 6, 10], m in [3, 5, 10]
            X = randn(n, m)
            F, dX_pullback = rrule(svd, X)
            for p in [:U, :S, :V]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = randn(size(Y)...)

                dself1, dF, dp = dF_pullback(Ȳ)
                @test dself1 === NO_FIELDS
                @test dp === DoesNotExist()

                ΔF = unthunk(dF)
                dself2, dX = dX_pullback(ΔF)
                @test dself2 === NO_FIELDS
                X̄_ad = unthunk(dX)
                X̄_fd = only(j′vp(central_fdm(5, 1), X->getproperty(svd(X), p), Ȳ, X))
                @test all(isapprox.(X̄_ad, X̄_fd; rtol=1e-6, atol=1e-6))
            end
            @testset "Vt" begin
                Y, dF_pullback = rrule(getproperty, F, :Vt)
                Ȳ = randn(size(Y)...)
                @test_throws ArgumentError dF_pullback(Ȳ)
            end
        end

        @testset "+" begin
            X = [1.0 2.0; 3.0 4.0; 5.0 6.0]
            F, dX_pullback = rrule(svd, X)
            X̄ = Composite{typeof(F)}(U=zeros(3, 2), S=zeros(2), V=zeros(2, 2))
            for p in [:U, :S, :V]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = ones(size(Y)...)
                dself, dF, dp = dF_pullback(Ȳ)
                @test dself === NO_FIELDS
                @test dp === DoesNotExist()
                X̄ += dF
            end
            @test X̄.U ≈ ones(3, 2) atol=1e-6
            @test X̄.S ≈ ones(2) atol=1e-6
            @test X̄.V ≈ ones(2, 2) atol=1e-6
        end

        @testset "Helper functions" begin
            X = randn(10, 10)
            Y = randn(10, 10)
            @test ChainRules._mulsubtrans!(copy(X), Y) ≈ Y .* (X - X')
            @test ChainRules._eyesubx!(copy(X)) ≈ I - X
            @test ChainRules._add!(copy(X), Y) ≈ X + Y
        end
    end
    @testset "cholesky" begin
        @testset "the thing" begin
            X = generate_well_conditioned_matrix(10)
            V = generate_well_conditioned_matrix(10)
            F, dX_pullback = rrule(cholesky, X)
            for p in [:U, :L]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = (p === :U ? UpperTriangular : LowerTriangular)(randn(size(Y)))
                (dself, dF, dp) = dF_pullback(Ȳ)
                @test dself === NO_FIELDS
                @test dp === DoesNotExist()

                # NOTE: We're doing Nabla-style testing here and avoiding using the `j′vp`
                # machinery from FiniteDifferences because that isn't set up to respect
                # necessary special properties of the input. In the case of the Cholesky
                # factorization, we need the input to be Hermitian.
                ΔF = unthunk(dF)
                _, dX = dX_pullback(ΔF)
                X̄_ad = dot(unthunk(dX), V)
                X̄_fd = _fdm() do ε
                    dot(Ȳ, getproperty(cholesky(X .+ ε .* V), p))
                end
                @test X̄_ad ≈ X̄_fd rtol=1e-6 atol=1e-6
            end
        end
        @testset "helper functions" begin
            A = randn(5, 5)
            r, d, B2, c = level2partition(A, 4, false)
            R, D, B3, C = level3partition(A, 4, 4, false)
            @test all(r .== R')
            @test all(d .== D)
            @test B2[1] == B3[1]
            @test all(c .== C)

            # Check that level 2 partition with `upper == true` is consistent with `false`
            rᵀ, dᵀ, B2ᵀ, cᵀ = level2partition(transpose(A), 4, true)
            @test r == rᵀ
            @test d == dᵀ
            @test B2' == B2ᵀ
            @test c == cᵀ

            # Check that level 3 partition with `upper == true` is consistent with `false`
            R, D, B3, C = level3partition(A, 2, 4, false)
            Rᵀ, Dᵀ, B3ᵀ, Cᵀ = level3partition(transpose(A), 2, 4, true)
            @test transpose(R) == Rᵀ
            @test transpose(D) == Dᵀ
            @test transpose(B3) == B3ᵀ
            @test transpose(C) == Cᵀ

            A = Matrix(LowerTriangular(randn(10, 10)))
            Ā = Matrix(LowerTriangular(randn(10, 10)))
            # NOTE: BLAS gets angry if we don't materialize the Transpose objects first
            B = Matrix(transpose(A))
            B̄ = Matrix(transpose(Ā))
            @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 1, false)
            @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 3, false)
            @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 5, false)
            @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 10, false)
            @test chol_unblocked_rev(Ā, A, false) ≈ transpose(chol_unblocked_rev(B̄, B, true))

            @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 1, true)
            @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 5, true)
            @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 10, true)
        end
    end
end


text:
"""
Module `MIRTjim` exports the "jiffy image display" method `jim`,
and the two helper methods `caller_name` and `prompt`.
"""
module MIRTjim

    const RealU = Number # Union{Real, Unitful.Length}

    include("jim.jl")
    include("3d.jl")
    include("array.jl")
    include("caller_name.jl")
    include("prompt.jl")

end # module


text:
"""
Provides a definition for a null mutation operator.
"""
type NullMutationDefinition <: MutationDefinition
  stage::AbstractString
end

type NullMutation <: Mutation
  stage::AbstractString
  representation::Representation
end

"""
Composes a null mutation operator from its definition.
"""
function compose!(d::NullMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  NullpMutation(d.stage, sp.stages[d.stage].representation)
end

"""
The null mutation operator accepts a single chromosome as its input and does
absolutely nothing to it. Useful for debugging purposes.

**Parameter:**

* `stage::AbstractString`, the name of the stage that this null operator should
  be applied to. If omitted, this operator will act upon the genotype.
"""
null() = null("")
null(stage::AbstractString) = null(stage)

"""
The null mutation operator accepts an input and does nothing with it.
"""
mutate!(::NullMutation, i::Any) = return i


text:
using DynamicalSystemsBase
using OrdinaryDiffEq, SimpleDiffEq, StaticArrays, LinearAlgebra
using Test

algs = (Vern9(), Tsit5(), SimpleATsit5())

using DynamicalSystemsBase: CDS, DDS, DS
using DynamicalSystemsBase.Systems: hoop, hoop_jac, hiip, hiip_jac
using DynamicalSystemsBase.Systems: loop, loop_jac, liip, liip_jac
using DiffEqBase

println("\nTesting tangent dynamics...")
@testset "tangent space" begin

for alg in algs

u0 = [0, 10.0, 0]
p = [10, 28, 8/3]
u0h = ones(2)
ph = [1.4, 0.3]

FUNCTIONS = [liip, liip_jac, loop, loop_jac, hiip, hiip_jac, hoop, hoop_jac]
INITCOD = [u0, u0h]
PARAMS = [p, ph]

# minimalistic lyapunov
function lyapunov_iip(ds::DS, k)
    D = dimension(ds)
    tode = tangent_integrator(ds, orthonormal(D,k); alg = alg)
    λ = zeros(k)
    for t in 1:1000
        while tode.t < t
            step!(tode)
        end
        # println("K = $K")
        QR = qr(get_deviations(tode))
        Q, R = QR.Q, QR.R
        λ .+= log.(abs.(diag(R)))

        set_deviations!(tode, Matrix(Q))
        u_modified!(tode, true)
    end
    λ = λ/tode.t # woooorks
end
function lyapunov_oop(ds::DS, k)
    D = dimension(ds)
    tode = tangent_integrator(ds, orthonormal(D,k); alg = alg)
    λ = zeros(k)
    ws_idx = SVector{k, Int}(collect(2:k+1))
    for t in 1:1000
        while tode.t < t
            step!(tode)
        end
        # println("K = $K")
        QR = qr(get_deviations(tode))
        Q, R = QR.Q, QR.R
        λ .+= log.(abs.(diag(R)))

        set_deviations!(tode, Q)
        u_modified!(tode, true)
    end
    λ./tode.t # woooorks
end

for i in 1:8
    @testset "$alg combination $i" begin
        sysindx = i < 5 ? 1 : 2
        if i < 5
            if isodd(i)
                ds = CDS(FUNCTIONS[i], INITCOD[sysindx], PARAMS[sysindx])
            else
                ds = CDS(FUNCTIONS[i-1], INITCOD[sysindx], PARAMS[sysindx], FUNCTIONS[i])
            end
        else
            if isodd(i)
                ds = DDS(FUNCTIONS[i], INITCOD[sysindx], PARAMS[sysindx])
            else
                ds = DDS(FUNCTIONS[i-1], INITCOD[sysindx], PARAMS[sysindx], FUNCTIONS[i])
            end
        end

        IIP = isinplace(ds)
        if IIP
            λ = lyapunov_iip(ds, 2)
        else
            λ = lyapunov_oop(ds, 2)
        end

        if i < 5
            @test 0.8 < λ[1] < 0.92
        else
            @test 0.4 < λ[1] < 0.45
        end
    end
end
end

end


text:
"""
# LocalRegistry

Create and maintain local registries for Julia packages. This package
is intended to provide a simple but manual workflow for maintaining
small local registries (private or public) without making any
assumptions about how the registry or the packages are hosted.

Registry creation is done by the function `create_registry`.

Registration of new and updated packages is done by the function `register`.
"""
module LocalRegistry

using RegistryTools: RegistryTools, gitcmd, Compress,
                     check_and_update_registry_files, ReturnStatus, haserror,
                     find_registered_version
using UUIDs: uuid4
using Pkg: Pkg, TOML

export create_registry, register

# Note: The `uuid` keyword argument is intentionally omitted from the
# documentation string since it's not intended for users.
"""
    create_registry(name, repo)
    create_registry(path, repo)

Create a registry with the given `name` or at the local directory
`path`, and with repository URL `repo`. The first argument is
interpreted as a path if it has more than one path component and
otherwise as a name. If a path is given, the last path component is
used as the name of the registry. If a name is given, it is created in
the standard registry location. In both cases the registry path must
not previously exist.

*Keyword arguments*

    create_registry(...; description = nothing, push = false, gitconfig = Dict())

* `description`: Optional description of the purpose of the registry.
* `push`: If `false`, the registry will only be prepared locally. Review the result and `git push` it manually.
* `gitconfig`: Optional configuration parameters for the `git` command.
"""
function create_registry(name_or_path, repo; description = nothing,
                         gitconfig::Dict = Dict(), uuid = nothing,
                         push = false)
    if length(splitpath(name_or_path)) > 1
        path = abspath(expanduser(name_or_path))
    else
        path = joinpath(first(DEPOT_PATH), "registries", name_or_path)
    end
    name = basename(path)
    if isempty(name)
        # If `path` ends with a slash, `basename` becomes empty and
        # the slash needs to be peeled off with `dirname`.
        name = basename(dirname(path))
    end
    if isdir(path) || isfile(path)
        error("$path already exists.")
    end
    mkpath(path)

    # The only reason for the `uuid` function argument is to allow
    # deterministic testing of the package.
    if isnothing(uuid)
        uuid = string(uuid4())
    end

    registry_data = RegistryTools.RegistryData(name, uuid, repo = repo,
                                               description = description)
    RegistryTools.write_registry(joinpath(path, "Registry.toml"), registry_data)

    git = gitcmd(path, gitconfig)
    run(`$git init -q`)
    run(`$git add Registry.toml`)
    run(`$git commit -qm 'Create registry.'`)
    run(`$git remote add origin $repo`)
    if push
        run(`$git push -u origin master`)
    end
    @info "Created registry in directory $(path)"

    return path
end

"""
    register(package, registry)

Register the new `package` in the registry `registry`. The version
number and other information is obtained from the package's
`Project.toml`.

    register(package)

Register a new version of `package`.  The version number is obtained
from the package's `Project.toml`.

Notes:
 * By default this will, in both cases, only update the registry
   locally. Review the result and `git push` it manually. This step
   can be made automatically with the keyword argument `push`
   described below.
 * The package must live in a git working copy, e.g. having been
   cloned by `Pkg.develop`.

`package` can be specified in the following ways:
* By package name. The package must be available in the active `Pkg`
  environment.
* By path. This is distinguished from package name by having more than
  one path component. A path in the current working directory can be
  specified by starting with `"./"`.
* By module. It needs to first be loaded by `using` or `import`.

`registry` can be specified by name or by path in the same way as
`package`. If omitted or `nothing`, it will be automatically located
by `package`.

*Keyword arguments*

    register(package, registry; commit = true, push = false, repo = nothing,
             gitconfig = Dict())

* `commit`: If `false`, only make the changes to the registry but do not commit.
* `push`: If `true`, push the changes to the registry repository automatically. Ignored if `commit` is false.
* `repo`: Specify the package repository explicitly. Otherwise looked up as the `git remote` of the package the first time it is registered.
* `gitconfig`: Optional configuration parameters for the `git` command.
"""
function register(package::Union{Module, AbstractString},
                  registry::Union{Nothing, AbstractString} = nothing;
                  commit = true, push = false, repo = nothing,
                  gitconfig::Dict = Dict())
    # Find and read the `Project.toml` for the package. First look for
    # the alternative `JuliaProject.toml`.
    package_path = find_package_path(package)
    local pkg
    for project_file in Base.project_names
        pkg = Pkg.Types.read_project(joinpath(package_path, project_file))
        if !isnothing(pkg.name)
            break
        end
    end
    if isnothing(pkg.name)
        error("$(package) does not have a Project.toml or JuliaProject.toml file")
    end

    # If the package directory is dirty, a different version could be
    # present in Project.toml.
    if is_dirty(package_path, gitconfig)
        error("Package directory is dirty. Stash or commit files.")
    end

    registry_path = find_registry_path(registry, pkg)
    if is_dirty(registry_path, gitconfig)
        error("Registry directory is dirty. Stash or commit files.")
    end

    # Compute the tree hash for the package and the subdirectory
    # location within the git repository. For normal packages living
    # at the top level of the repository, `subdir` will be the empty
    # string.
    tree_hash, subdir = get_tree_hash(package_path, gitconfig)

    # Check if this exact package has already been registered. Do
    # nothing and don't error in that case.
    if find_registered_version(pkg, registry_path) == tree_hash
        @info "This version has already been registered and is unchanged."
        return
    end

    # Use the `repo` argument or, if this is a new package
    # registration, check for the git remote. If `repo` is `nothing`
    # and this is an existing package, the repository information will
    # not be updated.
    if isnothing(repo)
        if !has_package(registry_path, pkg)
            package_repo = get_remote_repo(package_path, gitconfig)
        else
            package_repo = ""
        end
    else
        package_repo = repo
    end

    @info "Registering package" package_path registry_path package_repo uuid=pkg.uuid version=pkg.version tree_hash subdir
    clean_registry = true

    git = gitcmd(registry_path, gitconfig)
    HEAD = readchomp(`$git rev-parse --verify HEAD`)
    status = ReturnStatus()
    try
        check_and_update_registry_files(pkg, package_repo, tree_hash,
                                        registry_path, String[], status,
                                        subdir = subdir)
        if !haserror(status)
            if commit
                commit_registry(pkg, package_path, package_repo, tree_hash, git)
                if push
                    run(`$git push`)
                end
            end
            clean_registry = false
        end
    finally
        if clean_registry
            run(`$git reset --hard $(HEAD)`)
            run(`$git clean -f -d`)
        end
    end

    # Registration failed. Explain to the user what happened.
    if haserror(status)
        error(explain_registration_error(status))
    end

    return
end

function explain_registration_error(status)
    for triggered_check in status.triggered_checks
        check = triggered_check.id
        data = triggered_check
        if check == :version_exists
            return "Version $(data.version) has already been registered and the content has changed."
        elseif check == :change_package_name
            return "Changing package name is not supported."
        elseif check == :change_package_uuid
            return "Changing package UUID is not allowed."
        elseif check == :package_self_dep
            return "The package depends on itself."
        elseif check == :name_mismatch
            return "Error in (Julia)Project.toml: UUID $(data.uuid) refers to package '$(data.reg_name)' in registry but Project.toml has '$(data.project_name)'."
        elseif check == :wrong_stdlib_uuid
            return "Error in (Julia)Project.toml: UUID $(data.project_uuid) for package $(data.name) should be $(data.stdlib_uuid)"
        elseif check == :package_url_missing
            return "No repo URL provided for a new package."
        elseif check == :unexpected_registration_error
            return "Unexpected registration error."
        end
    end
end

# This does the same thing as `LibGit2.isdirty(LibGit2.GitRepo(path))`
# but also works when `path` is a subdirectory of a git
# repository. Only dirt within the subdirectory is considered.
function is_dirty(path, gitconfig)
    git = gitcmd(path, gitconfig)
    # TODO: There should be no need for the `-u` option but without it
    # a bogus diff is reported in the tests.
    return !isempty(read(`$git diff-index -u HEAD -- .`))
end

# If the package is provided as a module, directly find the package
# path from the loaded code. This works both if the module is loaded
# from the current package environment or found in LOAD_PATH.
function find_package_path(package::Module)
    return dirname(dirname(pathof(package)))
end

# A string argument is either interpreted as a path or as a package
# name, decided by the number of components returned by `splitpath`.
#
# If the package is given by name, it must be available in the current
# package environment as a developed package.
function find_package_path(package_name::AbstractString)
    if length(splitpath(package_name)) > 1
        if !isdir(package_name)
            error("Package path $(package_name) does not exist.")
        end
        return abspath(expanduser(package_name))
    end

    ctx = Pkg.Types.Context()
    if !haskey(ctx.env.project.deps, package_name)
        error("Unknown package $package_name.")
    end
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    if isnothing(pkg_path)
        error("Package must be developed to be registered.")
    elseif !isabspath(pkg_path)
        # If the package is developed with --local, pkg_path is
        # relative to the project path.
        pkg_path = joinpath(dirname(ctx.env.manifest_file), pkg_path)
    end

    # `pkg_path` might be a relative path, in which case it is
    # relative to the directory of Manifest.toml. If `pkg_path`
    # already is an absolute path, this call does not affect it.
    pkg_path = abspath(dirname(ctx.env.manifest_file), pkg_path)

    return pkg_path
end

function find_registry_path(registry::AbstractString, ::Pkg.Types.Project)
    return find_registry_path(registry)
end

function find_registry_path(registry::AbstractString)
    if length(splitpath(registry)) > 1
        return abspath(expanduser(registry))
    end

    all_registries = Pkg.Types.collect_registries()
    matching_registries = filter(r -> r.name == registry, all_registries)
    if isempty(matching_registries)
        error("Registry $(registry) not found.")
    end

    return first(matching_registries).path
end

function find_registry_path(::Nothing, pkg::Pkg.Types.Project)
    all_registries = Pkg.Types.collect_registries()

    matching_registries = filter(all_registries) do reg_spec
        reg_data = Pkg.Types.read_registry(joinpath(reg_spec.path,
                                                    "Registry.toml"))
        haskey(reg_data["packages"], string(pkg.uuid))
    end

    if isempty(matching_registries)
        error("Package $(pkg.name) not found in any registry. Please specify in which registry you want to register it.")
    elseif length(matching_registries) > 1
        error("Package $(pkg.name) is registered in more than one registry, please specify in which you want to register the package.")
    end

    return first(matching_registries).path
end

function has_package(registry_path, pkg::Pkg.Types.Project)
    registry = Pkg.Types.read_registry(joinpath(registry_path, "Registry.toml"))
    return haskey(registry["packages"], string(pkg.uuid))
end

function get_tree_hash(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    subdir = readchomp(`$git rev-parse --show-prefix`)
    tree_hash = readchomp(`$git rev-parse HEAD:$subdir`)
    # Get rid of trailing slash.
    if isempty(basename(subdir))
        subdir = dirname(subdir)
    end
    return tree_hash, subdir
end

function get_remote_repo(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    remote_name = split(readchomp(`$git remote`), '\n')
    length(remote_name) > 1 && error("Repository has multiple remotes.")
    remote_name[1] == "" && error("Repository does not have a remote.")
    return readchomp(`$git remote get-url $(remote_name[1])`)
end

function commit_registry(pkg::Pkg.Types.Project, package_path, package_repo, tree_hash, git)
    @debug("commit changes")
    message = """
    New version: $(pkg.name) v$(pkg.version)

    UUID: $(pkg.uuid)
    Repo: $(package_repo)
    Tree: $(string(tree_hash))
    """
    run(`$git add --all`)
    run(`$git commit -qm $message`)
end

end


text:
using Makie

f(x,y) = x^3*y^3

x = range(-1,1, length=400)
y = range(-1,1, length=400)
z = f.(x, y')

surface(x,y,z)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module ProbabilityDistributionTest

using Test
using ForneyLab
using ForneyLab: isProper, gaussianQuadrature, matches
using LinearAlgebra: Diagonal

@testset "<<" begin
    @test <<(ProbabilityDistribution{Univariate, Gaussian}, ProbabilityDistribution)
    @test <<(ProbabilityDistribution{Univariate, GaussianMeanVariance}, ProbabilityDistribution)
    @test <<(ProbabilityDistribution{Multivariate, GaussianMeanVariance}, ProbabilityDistribution)
    @test !<<(Nothing, ProbabilityDistribution)
end

@testset "sample" begin
    dist = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    samples = sample(dist, 10)
    @test length(samples) == 10
end

@testset "Univariate" begin
    # ProbabilityDistribution should be parameterized on a node type (distribution family)
    gaussian = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    @test isa(gaussian, ProbabilityDistribution{Univariate})
    @test isa(gaussian, ProbabilityDistribution{Univariate, GaussianMeanVariance})
    @test !isa(gaussian, ProbabilityDistribution{Multivariate})
    @test !isa(gaussian, ProbabilityDistribution{MatrixVariate})
    @test gaussian.params[:m] == 0.0
    @test gaussian.params[:v] == 1.0

    # PointMass should be defined as a special family
    point_mass = ProbabilityDistribution(Univariate, PointMass, m=0.0)
    @test isa(point_mass, ProbabilityDistribution{Univariate, PointMass})
    @test point_mass.params == Dict(:m=>0.0)
    @test isProper(point_mass)
    @test mean(point_mass) == 0.0
    @test var(point_mass) == 0.0

    # Remaining DeltaFactors should not be allowed as a distribution family
    @test_throws Exception ProbabilityDistribution(Univariate, Equality)
    @test_throws Exception ProbabilityDistribution(Univariate, Clamp, m=0.0)

    @test ProbabilityDistribution(Univariate, PointMass, m=0.0) != ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    @test ProbabilityDistribution(Univariate, PointMass, m=0.0) != ProbabilityDistribution(Univariate, PointMass, m=1.0)
end

@testset "Multivariate" begin
    gaussian = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0))
    @test isa(gaussian, ProbabilityDistribution{Multivariate})
    @test isa(gaussian, ProbabilityDistribution{Multivariate, GaussianMeanVariance})
    @test !isa(gaussian, ProbabilityDistribution{Univariate})
    @test !isa(gaussian, ProbabilityDistribution{MatrixVariate})
    @test gaussian.params[:m] == [0.0]
    @test gaussian.params[:v] == mat(1.0)

    # PointMass should be defined as a special family
    point_mass = ProbabilityDistribution(Multivariate, PointMass, m=[0.0])
    @test isa(point_mass, ProbabilityDistribution{Multivariate, PointMass})
    @test point_mass.params == Dict(:m=>[0.0])
    @test isProper(point_mass)
    @test mean(point_mass) == [0.0]
    @test var(point_mass) == [0.0]
    @test cov(point_mass) == transpose([0.0])
end

@testset "MatrixVariate" begin
    point_mass = ProbabilityDistribution(MatrixVariate, PointMass, m=transpose([0.0]))
    @test isa(point_mass, ProbabilityDistribution{MatrixVariate, PointMass})
    @test point_mass.params == Dict(:m=>transpose([0.0]))
    @test isProper(point_mass)
    @test mean(point_mass) == transpose([0.0])
end

@testset "convert" begin
    @test convert(ProbabilityDistribution{Multivariate, PointMass}, ProbabilityDistribution(Univariate, PointMass, m=1.0)) == ProbabilityDistribution(Multivariate, PointMass, m=[1.0])
    @test convert(ProbabilityDistribution{MatrixVariate, PointMass}, ProbabilityDistribution(Univariate, PointMass, m=1.0)) == ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0))
    @test convert(ProbabilityDistribution{MatrixVariate, PointMass}, ProbabilityDistribution(Multivariate, PointMass, m=[1.0, 2.0])) == ProbabilityDistribution(MatrixVariate, PointMass, m=reshape([1.0, 2.0], 2, 1))
end

@testset "PointMass ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Univariate, PointMass, m=0.2) == ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>0.2))
    @test_throws Exception ProbabilityDistribution(Multivariate, PointMass, m=0.2)
    @test ProbabilityDistribution(Multivariate, PointMass, m=[0.2]) == ProbabilityDistribution{Multivariate, PointMass}(Dict(:m=>[0.2]))
    @test ProbabilityDistribution(MatrixVariate, PointMass, m=transpose([0.2])) == ProbabilityDistribution{MatrixVariate, PointMass}(Dict(:m=>transpose([0.2])))
    @test ProbabilityDistribution(PointMass, m=0.2) == ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>0.2))
    @test ProbabilityDistribution(Univariate, PointMass) == ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>1.0))
    @test ProbabilityDistribution(PointMass) == ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>1.0))
    @test ProbabilityDistribution(Multivariate, PointMass) == ProbabilityDistribution{Multivariate, PointMass}(Dict(:m=>[1.0]))
    @test ProbabilityDistribution(MatrixVariate, PointMass) == ProbabilityDistribution{MatrixVariate, PointMass}(Dict(:m=>mat(1.0)))
    @test Message(PointMass) == Message{PointMass, Univariate}(ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>1.0)))
    @test Message(Univariate, PointMass) == Message{PointMass, Univariate}(ProbabilityDistribution{Univariate, PointMass}(Dict(:m=>1.0)))
    @test_throws Exception Message(Multivariate, PointMass, m=0.2)
end

@testset "ProbabilityDistribution and Message construction with functions" begin
    f(x) = x
    @test isa(ProbabilityDistribution(Univariate, Function, f=f).params[:f], Function)
    @test isa(ProbabilityDistribution(Univariate, Function, f=f), ProbabilityDistribution{Univariate, Function})
    @test isa(ProbabilityDistribution(Multivariate, Function, f=f), ProbabilityDistribution{Multivariate, Function})
    @test isa(ProbabilityDistribution(MatrixVariate, Function, f=f), ProbabilityDistribution{MatrixVariate, Function})
    @test isa(ProbabilityDistribution(Function, f=f), ProbabilityDistribution{Univariate, Function})
    @test isempty(vague(Function).params)
    @test isa(Message(Univariate, Function, f=()->()).dist.params[:f], Function)
    @test isa(Message(Univariate, Function), Message{Function, Univariate})
    @test isa(Message(Function), Message{Function, Univariate})
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Univariate, PointMass, m=0.0)) == ()
    @test dims(ProbabilityDistribution(Multivariate, PointMass, m=ones(2))) == (2,)
    @test dims(ProbabilityDistribution(MatrixVariate, PointMass, m=eye(2))) == (2,2)
end

@testset "gaussianQuadrature" begin
    @test gaussianQuadrature(x -> 1.0, m=-1.0, v=2.0) == 0.9999996733487053
    @test gaussianQuadrature(x -> x, m=-1.0, v=2.0) == -0.9999996733487055
    @test gaussianQuadrature(x -> x^2, m=-1.0, v=2.0) == 3.0000013419448215
end

end #module


text:
export MixedHomotopy

"""
    MixedHomotopy

A data structure which uses for `evaluate!` and `evaluate_and_jacobian!` a [`CompiledHomotopy`](@ref)
and for `taylor!` an [`InterpretedHomotopy`](@ref).
"""
struct MixedHomotopy{T,ID} <: AbstractHomotopy
    compiled::CompiledHomotopy{ID}
    interpreted::InterpretedHomotopy{T}
end
MixedHomotopy(H::Homotopy) = MixedHomotopy(CompiledHomotopy(H), InterpretedHomotopy(H))

Base.size(H::MixedHomotopy) = size(H.compiled)
ModelKit.variables(H::MixedHomotopy) = variables(H.interpreted)
ModelKit.parameters(H::MixedHomotopy) = parameters(H.interpreted)
ModelKit.variable_groups(H::MixedHomotopy) = variable_groups(H.interpreted)
ModelKit.Homotopy(H::MixedHomotopy) = Homotopy(H.interpreted)

function Base.show(io::IO, H::MixedHomotopy)
    print(io, "Hybrid: ")
    show(io, Homotopy(H))
end

(H::MixedHomotopy)(x, t, p = nothing) = H.interpreted(x, t, p)
ModelKit.evaluate!(u, H::MixedHomotopy, x::AbstractVector, t, p = nothing) =
    evaluate!(u, H.compiled, x, t, p)
ModelKit.evaluate_and_jacobian!(u, U, H::MixedHomotopy, x, t, p = nothing) =
    evaluate_and_jacobian!(u, U, H.compiled, x, t, p)
ModelKit.taylor!(u, v::Val, H::MixedHomotopy, tx, t, p = nothing) =
    taylor!(u, v, H.interpreted, tx, t, p)


text:
"""
    add_obj_fun!(
        ata_model::CCMaximinModel;
        psychometrics = false,
        items_file = "",
        items::Vector{Psychometrics.Item} = Psychometrics.Item[],
        kwargs...
    )

# Description

It adds the objective function as specified in the `settings_file`. It requires the [`start_ata`](#ATA.start_ata) build step.

Computes the IIFs at predefined ability points using the sampled item parameter estimates.
For each triplet of item, ability point, and item parameter sample \$(i, θ_k, r)\$, an IIF\$(θ_k)_{ir}\$ is computed.

The vector of items can be passed either through the argument `items` (`Vector{Psychometrics.Item}`) or through the argument `items_file` (string file name).
In both cases, the items in the vector must have the same order of the items in the item bank and they must contain the R sampled item parameters in the field `parameters.chain`.

# Arguments

- **`ata_model::CCMaximinModel`** : Required.
"""
function add_obj_fun!(
    ata_model::CCMaximinModel;
    psychometrics = false,
    items_file = "",
    items::Vector{Psychometrics.Item} = Psychometrics.Item[],
    kwargs...,
)
    try
        if !psychometrics
            R = ata_model.obj.R
            K = zeros(Int, T)
            IIF = Vector{Array{Float64,3}}(undef, T)
            ICF = Vector{Array{Float64,3}}(undef, T)
            if !isfile("BSPar.jld2")
                push!(
                    ata_model.output.infos,
                    [
                        "danger",
                        "ccmaximin objective requires a jld2 file named \"BSPar.jld2\" with sampled values for the item parameters.\nAlternatively, a vector of `Psychometrics.Item` with samples in chain can be provided.",
                    ],
                )
                return nothing
            end
            JLD2.@load "BSPar.jld2" BSPar
            BSa = Matrix(BSPar[2]) |> x -> x[:, 2:end]
            BSb = Matrix(BSPar[1]) |> x -> x[:, 2:end]
            for t = 1:T
                K[t] = size(ata_model.obj.cores[t].points, 1)
                IIF[t] = zeros(K[t], n_items, R)
                ICF[t] = zeros(K[t], n_items, R)
                for r = 1:R
                    if irt_model == "1PL"
                        df = DataFrames.DataFrame(b = BSb[:, r]) #nqp values in interval\r\n",
                    elseif irt_model == "2PL"
                        df = DataFrames.DataFrame(a = BSa[:, r], b = BSb[:, r]) #nqp values in interval\r\n",
                    elseif irt_model == "3PL"
                        df = DataFrames.DataFrame(
                            a = BSa[:, r],
                            b = BSb[:, r],
                            c = BSc[:, r],
                        )
                    end
                    for k = 1:K[t]
                        IIF[t][k, :, r] = item_info(
                            df,
                            ata_model.obj.cores[t].points[k];
                            model = irt_model,
                            parametrization = irt_parametrization,
                            D = irt_D,
                        ) # K[t] x I x R
                        ICF[t][k, :, r] = item_char(
                            df,
                            ata_model.obj.cores[t].points[k];
                            model = irt_model,
                            parametrization = irt_parametrization,
                            D = irt_D,
                        )[1][
                            :,
                            :,
                            1,
                        ] # K[t] x I x R
                    end
                    ata_model.obj.cores[t].IIF = IIF[t]
                end
            end
            JLD2.@save "opt/IIF_CC.jld2" IIF
            JLD2.@save "opt/ICF_CC.jld2" ICF
        else
            cc_maximin_load_parameters_chain!(
                ata_model;
                items_file = items_file,
                items = items,
                kwargs...,
            )
        end
        success!(ata_model, "IIFs for all item parameters samples computed.")
        # open("opt/Settings.jl", "a") do f
        #     write(f, "K = $K\n\n")
        # end

    catch e
        error!(ata_model, string(sprint(showerror, e)))

    end
    return nothing
end


text:
module Duck
using LinearAlgebra
using FreqTables
export leastsquares

function leastsquares(x, y, Ω; order::T = 1) where T<:Integer
    X = ones((length(x), order))
    for n in 1:order
        X[:, n+1] = x.^n
    end

    if ndims(Ω) == 1
        Ω = diagm(0 => Ω)
    end

    w = inv(Ω)  # Weights defined as 1/σ²
    Xᵀ = transpose(X)
    kω = inv(Xᵀ*w*X)*Xᵀ*w
	coeffs = kω*y
    Σ = kω*Ω*kω'
    coeffs, Σ
end

function leastsquares(x, y; order::T = 1) where T<:Integer
    X = ones((length(x), order+1))
    for n in 1:order
        X[:, n+1] = x.^n
    end
    inv(transpose(X)*X)*transpose(X)*y
end

mutable struct RNG
    #= Simple random number generator =#
    value::Float64
end

function Base.iterate(rng::RNG, state::Int64)
    #= Park-Miller-Carta Pseudo-Random Number Generator =#
    a = 16807
    m = 2147483647
    q = 127773
    r = 2836
    MASK = 123459876
    am = 1/m
    state = state ⊻ MASK
    k::Int64 = round(state/q)
    state = a*(state - k*q) - r*k
    state += if state < 0 m else 0 end
    rng.value = am*state
    state = state ⊻ MASK
    (rng.value, state)
end

function rngarray(N)
    A = zeros(N)
    i = 1
    @inbounds for random in RNG(abs(rand(Int)) % 10492841)
        A[i] = random
        i += 1
        if i ≥ N+1
            break
        end
    end
    A
end

function israndom(randomnums::Array{Int64, 1}, r)
    numsamples = length(randomnums)
    if numsamples ≤ 10r
        return false
    end

    nr = numsamples/r
    # Frequency of randoms
    ht = freqtable(randomnums).array
    χ² = sum((v - nr)^2 for v in ht)/nr
    return abs(χ² - r) / 2*√r
end

end

max_stars_repo_path: exercises/exercise 1/statistics.jl

max_issues_repo_path: exercises/exercise 1/statistics.jl

max_forks_repo_path: exercises/exercise 1/statistics.jl


text:
mutable struct body_wrench_data_t <: LCMType
    timestamp::Int64
    body_name::String
    wrench::SVector{6, Float64}
end

@lcmtypesetup(body_wrench_data_t)


text:
"""
A [`Message`](@ref) attachment.
More details [here](https://discordapp.com/developers/docs/resources/channel#attachment-object).
"""
struct Attachment
    id::Snowflake
    filename::String
    size::Int
    url::String
    proxy_url::String
    height::Union{Int, Missing}
    width::Union{Int, Missing}
end
@boilerplate Attachment :dict :docs :lower :merge


text:
# This example shows how to convert a file with a list of images
# to HDF5 format that Mocha can process
#
# Note this is a simple demo to show how to convert a typical
# dataset into HDF5 file that Mocha could handle. We do not
# perform any pre-processing on the images and simply assume
# all the images are of the same size and with the same
# number of channels.
#
# The format of the file is:
#
# ImageName,Label,TrainTest
# image1.jpg,009,Train
# image2.jpg,018,Train
# ....
# imageK.jpg,053,Test
# ...
# imageN.jpg,009,Test

using DataFrames
using Images
using Color
using HDF5
using Compat

############################################################
# Configuration
############################################################
image_list_file = "example.csv"
data_type       = Float32

############################################################
# Converting dataset
############################################################
datasets = Dict{AbstractString, HDF5File}()
dsets_idx = Dict{AbstractString, Int}()

df = readtable(image_list_file)
for i = 1:size(df,1)
  image_fn = df[i, :ImageName]
  label    = df[i, :Label]
  dset_key = df[i, :TrainTest]

  println(image_fn)

  image = separate(convert(Image{RGB}, imread(image_fn)))
  data  = convert(Array, image)

  if !haskey(datasets, dset_key)
    #---------------
    # create dataset
    #---------------

    # count how many elements in this dataset
    dataset_count = sum(df[:TrainTest] .== dset_key)

    h5 = h5open("$dset_key.hdf5", "w")
    dset_data = d_create(h5, "data", datatype(data_type),
        dataspace(size(data)..., dataset_count))
    dset_label = d_create(h5, "label", datatype(data_type),
        dataspace(1, dataset_count))
    datasets[dset_key] = h5
    dsets_idx[dset_key] = 0
  end

  h5 = datasets[dset_key]
  idx = dsets_idx[dset_key]
  idx += 1
  h5["label"][1,idx] = label
  h5["data"][map(x -> 1:x, size(data))..., idx] = data
end

for h5 in values(datasets)
  close(h5)
end


text:
include("fbase.jl")
include("quadrature.jl")
include("utils.jl")

include("melem_rt0.jl")

using LinearAlgebra

function melemk0P2(verts::Array{Float64, 2})
    #= Compute the A = ∫div(v) div(u) and the B = ∫v u element-matrices when the
        approximation order is k=0. In that case, there are no interior degrees
        of freedom and the amount of dofs is exactly the number of edges =#
    nvtx = size(verts, 2)
    ndof = nvtx
    # 0  : preliminaries
    E  = polygonArea(verts)         # area
    xE = polygonCentroid(verts)     # centroid
    hE = polygonDiameter(verts)     # diameter

    # 1st: compute A; it is easily done thanks to k=0
    DIV = ones(Float64, ndof, ndof) / E

    # 2nd: compute B; somehow a little more difficult as a projection needs to
    # be computed
    gmα = getGradientMonomialsP2(xE, hE)
    ref2d = quadGL2d(3)

    Παβ = zeros(Float64, 5, 5)
    for α=1:5
        for β = 1:5
            Παβ[α,β] = integratePolygon(verts, xE, ref2d,
                x -> dot(gmα[α+1](x), gmα[β+1](x)))
        end
    end
    # RHS: first term: null!
    # RHS: second term
    # initialization of the interior basis functions
    mα = getMonomialsP2(xE, hE)
    Cs = (-1/E)*integratePolygon(verts, xE, ref2d, mα) # average over polygon
    mα = getMonomialsP2(xE, hE, Cs) # 0-averaged monomials
    #
    rhs = zeros(Float64, 5, ndof)
    ref1d = quadGL1d(2, 0., 1.)
    # rhs[α, i] = int_{e_i}{((φ_i⋅n) Φ_α)} where φ_i⋅n = 1
    for i=1:ndof
        ip  = mod_wrap(i+1, ndof)        # return the next index modulo ndof
        e_v = verts[:, ip] - verts[:, i] # tangent to edge
        le  = norm(e_v)                  # length of edge
        for α=1:5
            # integration with quadrature rule
            rhs[α, i] = integrate(ref1d, t -> mα[α+1](verts[:, i] + t*e_v)/le)*le
        end
    end
    Π = Παβ \ rhs # Π is the projection matrix
    # now the Di,α = dof_i(m_α) = ∫(∇mα . n) q^i  matrix
    D = zeros(Float64, ndof, 5)
    for i=1:ndof
       ip  = mod_wrap(i+1, ndof)
       e_v = verts[:, ip] - verts[:, i]
       le  = norm(e_v)
       n_v = e_v / le          # normalized tangent
       n_v = [n_v[2], -n_v[1]] # outside pointing normal vector
       for α = 1:5
           D[i,α] = integrate(ref1d,
            t -> dot(gmα[α+1]( verts[:, i] + t*e_v ), n_v))*le
       end
    end
    # stability constant
    σ    = 1.
    F    = Matrix(1.0I, ndof, ndof) - D*Π
    MASS = Π'*(Παβ*Π) + σ*(F' * F)

    # 3rd:  both parts have been computed
    return DIV, MASS
end

function melemk0P1(verts::Array{Float64, 2})
    #= Compute the A = ∫div(v) div(u) and the B = ∫v u element-matrices when the
        approximation order is k=0. In that case, there are no interior degrees
        of freedom and the amount of dofs is exactly the number of edges =#
    nvtx = size(verts, 2)
    ndof = nvtx
    # 0  : preliminaries
    E  = polygonArea(verts)         # area
    xE = polygonCentroid(verts)     # centroid
    hE = polygonDiameter(verts)     # diameter

    # 1st: compute A; it is easily done thanks to k=0
    DIV = ones(Float64, ndof, ndof) / E

    # 2nd: compute B; somehow a little more difficult as a projection needs to
    # be computed
    gmα = getGradientMonomialsP1(xE, hE)
    ref2d = quadGL2d(1)

    Παβ = zeros(Float64, 2, 2)
    for α=1:2
        for β = 1:2
            # quadrature rule not really necessary as everythin is constant
            # but it mimics the code of melemk0P2
            Παβ[α,β] = integratePolygon(verts, xE, ref2d,
                x -> dot(gmα[α+1](x), gmα[β+1](x)))
        end
    end
    # RHS: first term: null!
    # RHS: second term
    # initialization of the interior basis functions
    mα = getMonomialsP1(xE, hE)
    Cs = (-1/E)*integratePolygon(verts, xE, ref2d, mα) # average over polygon
    mα = getMonomialsP1(xE, hE, Cs) # 0-averaged monomials
    #
    rhs = zeros(Float64, 2, ndof)
    ref1d = quadGL1d(2, 0., 1.)
    # rhs[α, i] = int_{e_i}{((φ_i⋅n) Φ_α)} where φ_i⋅n = 1/le
    for i=1:ndof
        ip  = mod_wrap(i+1, ndof)
        e_v = verts[:, ip] - verts[:, i]
        le  = norm(e_v)
        for α=1:2
            # rhs[α, i] = integrate(ref1d, t -> mα[α+1](verts[:, i] + t*e_v)/le)*le
            rhs[α, i] = mα[α+1]((verts[:, ip] + verts[:, i])/2); # with the midpoint rule
        end
    end
    Π = Παβ \ rhs # Π is the projection matrix
    # now the Di,α = dof_i(m_α) = ∫(∇mα . n) q^i  matrix
    D = zeros(Float64, ndof, 2)
    for i=1:ndof
       ip  = mod_wrap(i+1, ndof)
       e_v = verts[:, ip] - verts[:, i]
       le  = norm(e_v)
       n_v = e_v / le
       n_v = [n_v[2], -n_v[1]]
       for α = 1:2
            D[i,α] = integrate(ref1d,
            t -> dot(gmα[α+1]( verts[:, i] + t*e_v ), n_v))*le
            # D[i,α] = dot(gmα[α+1](0.), n_v)*le # if you wish the midpoint rule, comment out
       end
    end
    # stability constant
    σ    = 1.
    F    = Matrix(1.0I, ndof, ndof) - D*Π
    MASS = Π'*(Παβ*Π) + σ*(F' * F)

    # 3rd:  both parts have been computed
    return DIV, MASS
end


text:
"""
Miscealenous functions for genetic genealogy.
Function to calculate an ethnicity estimate for MyHeritage matches files.
Functions to manipulate phylogenetic trees and to calculate TMRCAs.
"""
module GeneticGenealogy

export ethnicity, readtree, readtrees, writetree, calculateTMRCAs

include("ethnicity.jl")
include("tree.jl")

end # module


text:
using JuliaWebAPI
using Test
using Distributed

include("clnt.jl")

const opts = Base.JLOptions()
const inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
const cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``

function addproc_srvr()
    addprocs(1; exeflags=`$cov_flag $inline_flag`)
end

function rmproc_srvr()
    rmprocs(workers())
end

function time_remotecall()
    addproc_srvr()

    t1 = time()
    for idx in 1:NCALLS
        arg1,arg2,narg1,narg2 = APIARGS[(4*idx-3):(4*idx)]
        @test remotecall_fetch((a1,a2,a3,a4)->(a1*a3 + a2*a4), 2, arg1, arg2, narg1, narg2) == (arg1 * narg1) + (arg2 * narg2)
    end
    t = time() - t1
    println("time for $NCALLS calls with remotecall_fetch: $t secs @ $(t/NCALLS) per call")
    rmproc_srvr()

    println("stopped all workers")
end

# run tests if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runremotecall") && time_remotecall()