text:
using Plots, BenchmarkTools, BioParticleSimulation

f = BioParticleSimulation.FieldSet(100, 100)
t = BioParticleSimulation.CellType(0, 0.3, 0.000001, 0.0007, 100, 2000, 4000, 100, 100, 10, 1)
BioParticleSimulation.addresource(f, 500, 75)
BioParticleSimulation.fillfield(f, 0.005, t) #0.0007
BioParticleSimulation.round(f)
BioParticleSimulation.plotfield(f)

for i = 1:10000

    BioParticleSimulation.round(f)
    BioParticleSimulation.growresource(f, 0.005)
    BioParticleSimulation.round(f)
end

p = BioParticleSimulation.plotfield(f)
BioParticleSimulation.round(f)
BioParticleSimulation.printall(f)

#Graph

f = BioParticleSimulation.FieldSet(100, 100)
BioParticleSimulation.fillfield(f, 0.005)
BioParticleSimulation.addresource(f, 500, 75)
n = Vector{Int}(undef, 0)
e = Vector{Float64}(undef, 0)
anim = @animate for i ∈ 1:600
    BioParticleSimulation.round(f)
    BioParticleSimulation.growresource(f, 0.005)
    BioParticleSimulation.round(f)
    push!(n, length(f.list))
    push!(e, sum(f.rfield))
    plot(BioParticleSimulation.plotfield(f), plot(n, legend = false), plot(e, legend = false),layout = grid(3, 1, heights=[0.6, 0.2, 0.2]), size=(400, 800))
end
gif(anim, "anim_fps15.gif", fps = 15)

n = Vector{Int}(undef, 0)
e = Vector{Float64}(undef, 0)
for i = 1:40000
    BioParticleSimulation.round(f)
    BioParticleSimulation.growresource(f, 0.005)
    BioParticleSimulation.round(f)
    push!(n, length(f.list))
    push!(e, sum(f.rfield))
end
plot(plot(n, legend = false), plot(e, legend = false), layout = grid(2, 1, heights=[0.5, 0.5]), size=(400, 800))
#benchmark

bench = @benchmark for i = 1:2000
    BioParticleSimulation.round(f)
    BioParticleSimulation.growresource(f, 0.005)
end


text:
"""
    broadcast!!(f, dest, As...) -> dest′

A mutate-or-widen version of `dest .= f.(As...)`.
"""
broadcast!!
@inline broadcast!!(f::F, dest, As...) where {F} =
    materialize!!(dest, instantiate(broadcasted(f, As...)))

"""
    broadcast_inplace!!(f, inputoutput, As...) -> inputoutput′

A mutate-or-widen version of `inputoutput .= f.(inputoutput, As...)`.
"""
Extras.broadcast_inplace!!
@inline Extras.broadcast_inplace!!(f::F, inputoutput, As...) where {F} =
    materialize!!(inputoutput, instantiate(broadcasted(f, inputoutput, As...)))

ext: jl


text:
using ManifoldLearning
using Test
import Random

Random.seed!(3483743871)

@testset "ManifoldLearning" begin

    # setup parameters
    k = 12
    d = 2
    X, L = ManifoldLearning.swiss_roll()

    # test NN graph
    D, E = ManifoldLearning.find_nn(X, k)
	@test size(X,2) == size(D,2) && size(D, 1) == k
	@test size(X,2) == size(E,2) && size(E, 1) == k

    # test algorithms
    @testset for algorithm in [Isomap, LEM, LLE, HLLE, LTSA, DiffMap]
        for (k, T) in zip([5, 12], [Float32, Float64])
            # construct KW parameters
            kwargs = [:maxoutdim=>d]
            if algorithm == DiffMap
                push!(kwargs, :t => k)
            else
                push!(kwargs, :k => k)
            end

            # call transformation
            Y = fit(algorithm, convert(Array{T}, X); kwargs...)

            # test results
            @test outdim(Y) == d
            @test size(transform(Y), 2) == size(X, 2)
            @test length(split(sprint(show, Y), '\n')) > 1
            if algorithm !== DiffMap
                @test neighbors(Y) == k
                @test length(eigvals(Y)) == d
            end
            if algorithm === Isomap || algorithm === LEM
                @test length(vertices(Y)) > 1
            end
        end
    end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((2, 4, 1.0),),
    ((2, 4, 0.5),),
    ]
insts["fast"] = [
    ((15, 20, 1.0),),
    ((15, 20, 0.25),),
    ((100, 100, 1.0),),
    ((100, 100, 0.15),),
    ((500, 100, 1.0),),
    ((500, 100, 0.15),),
    ]
insts["various"] = [
    ((250, 1000, 0.1),),
    ((500, 5000, 0.05),),
    ((500, 2000, 1.0),),
    ]
return (LinearOptNative, insts)


text:
import MathProgBase
const MPB = MathProgBase

###############################################################################
# Solver objects
export IpoptSolver
struct IpoptSolver <: MPB.AbstractMathProgSolver
    options::Vector{Tuple} # list of options set in Ipopt on each MPB.optimize! call
end
function IpoptSolver(;kwargs...)
    args = Vector{Tuple}()
    for arg in kwargs
        if isa(arg, Pair)
            push!(args, (arg.first, arg.second))
        else # is a tuple in v0.6
            push!(args, arg)
        end
    end
    IpoptSolver(args)
end

mutable struct IpoptMathProgModel <: MPB.AbstractNonlinearModel
    inner::IpoptProblem
    numvar::Int
    numconstr::Int
    warmstart::Vector{Float64}
    options

    function IpoptMathProgModel(;options...)
        model = new()
        model.options   = options
        model.numvar    = 0
        model.numconstr = 0
        model.warmstart = Float64[]
        model
    end
end
MPB.NonlinearModel(s::IpoptSolver) = IpoptMathProgModel(;s.options...)
MPB.LinearQuadraticModel(s::IpoptSolver) = MPB.NonlinearToLPQPBridge(MPB.NonlinearModel(s))

###############################################################################
# Begin interface implementation

# generic nonlinear interface
function MPB.loadproblem!(m::IpoptMathProgModel, numVar::Integer, numConstr::Integer,
                          x_l, x_u, g_lb, g_ub, sense::Symbol, d::MPB.AbstractNLPEvaluator)

    m.numvar = numVar
    features = MPB.features_available(d)
    has_hessian = (:Hess in features)
    init_feat = [:Grad]
    has_hessian && push!(init_feat, :Hess)
    numConstr > 0 && push!(init_feat, :Jac)

    MPB.initialize(d, init_feat)
    Ihess, Jhess = has_hessian ? MPB.hesslag_structure(d) : (Int[], Int[])
    Ijac, Jjac = numConstr > 0 ? MPB.jac_structure(d) : (Int[], Int[])
    @assert length(Ijac) == length(Jjac)
    @assert length(Ihess) == length(Jhess)
    @assert sense == :Min || sense == :Max

    # Objective callback
    scl = (sense == :Min) ? 1.0 : -1.0
    eval_f_cb(x) = scl * MPB.eval_f(d,x)

    # Objective gradient callback
    function eval_grad_f_cb(x, grad_f)
        MPB.eval_grad_f(d, grad_f, x)
        rmul!(grad_f, scl)
    end

    # Constraint value callback
    eval_g_cb(x, g) = MPB.eval_g(d, g, x)

    # Jacobian callback
    function eval_jac_g_cb(x, mode, rows, cols, values)
        if mode == :Structure
            for i in 1:length(Ijac)
                rows[i] = Ijac[i]
                cols[i] = Jjac[i]
            end
        else
            MPB.eval_jac_g(d, values, x)
        end
    end

    if has_hessian
        # Hessian callback
        function eval_h_cb(x, mode, rows, cols, obj_factor,
            lambda, values)
            if mode == :Structure
                for i in 1:length(Ihess)
                    rows[i] = Ihess[i]
                    cols[i] = Jhess[i]
                end
            else
                obj_factor *= scl
                MPB.eval_hesslag(d, values, x, obj_factor, lambda)
            end
        end
    else
        eval_h_cb = nothing
    end

    m.inner = createProblem(numVar, float(x_l), float(x_u), numConstr,
                            float(g_lb), float(g_ub), length(Ijac), length(Ihess),
                            eval_f_cb, eval_g_cb, eval_grad_f_cb, eval_jac_g_cb,
                            eval_h_cb)
    m.inner.sense = sense

    # Ipopt crashes by default if NaN/Inf values are returned from the
    # evaluation callbacks. This option tells Ipopt to explicitly check for them
    # and return Invalid_Number_Detected instead. This setting may result in a
    # minor performance loss and can be overwritten by specifying
    # check_derivatives_for_naninf="no".
    addOption(m.inner, "check_derivatives_for_naninf", "yes")
    if !has_hessian
        addOption(m.inner, "hessian_approximation", "limited-memory")
    end

end

MPB.getsense(m::IpoptMathProgModel) = m.inner.sense
MPB.numvar(m::IpoptMathProgModel) = m.numvar
MPB.numconstr(m::IpoptMathProgModel) = m.numconstr

MPB.numlinconstr(m::IpoptMathProgModel) = 0

MPB.numquadconstr(m::IpoptMathProgModel) = 0

function MPB.optimize!(m::IpoptMathProgModel)
    copyto!(m.inner.x, m.warmstart) # set warmstart
    for (name,value) in m.options
        sname = string(name)
        if occursin(r"(^resto_)", sname)
            sname = replace(sname, r"(^resto_)" => "resto.")
        end
        addOption(m.inner, sname, value)
    end
    solveProblem(m.inner)
end

function MPB.status(m::IpoptMathProgModel)
    # Map all the possible return codes, as enumerated in
    # Ipopt.ApplicationReturnStatus, to the MPB statuses:
    # :Optimal, :Infeasible, :Unbounded, :UserLimit, and :Error
    stat_sym = ApplicationReturnStatus[m.inner.status]
    if  stat_sym == :Solve_Succeeded ||
        stat_sym == :Solved_To_Acceptable_Level
        return :Optimal
    elseif stat_sym == :Infeasible_Problem_Detected
        return :Infeasible
    elseif stat_sym == :Diverging_Iterates
        return :Unbounded
        # Things that are more likely to be fixable by changing
        # a parameter will be treated as UserLimit, although
        # some are error-like too.
    elseif stat_sym == :User_Requested_Stop ||
        stat_sym == :Maximum_Iterations_Exceeded ||
        stat_sym == :Maximum_CpuTime_Exceeded
        return :UserLimit
    else
        # Default is to not mislead user that it worked
        # Includes:
        #   :Search_Direction_Becomes_Too_Small
        #   :Feasible_Point_Found
        #   :Restoration_Failed
        #   :Error_In_Step_Computation
        #   :Not_Enough_Degrees_Of_Freedom
        #   :Invalid_Problem_Definition
        #   :Invalid_Option
        #   :Invalid_Number_Detected
        #   :Unrecoverable_Exception
        #   :NonIpopt_Exception_Thrown
        #   :Insufficient_Memory
        #   :Internal_Error
        @warn "Ipopt finished with status $stat_sym"
        return :Error
    end

end
MPB.getobjval(m::IpoptMathProgModel) = m.inner.obj_val * (m.inner.sense == :Max ? -1 : +1)
MPB.getsolution(m::IpoptMathProgModel) = m.inner.x

function MPB.getreducedcosts(m::IpoptMathProgModel)
    sense = m.inner.sense
    redcost = m.inner.mult_x_U - m.inner.mult_x_L
    return sense == :Max ? redcost : -redcost
end
function MPB.getconstrduals(m::IpoptMathProgModel)
    v = m.inner.mult_g # return multipliers for all constraints
    return m.inner.sense == :Max ? copy(v) : -v
end

MPB.getrawsolver(m::IpoptMathProgModel) = m.inner
MPB.setwarmstart!(m::IpoptMathProgModel, x) = (m.warmstart = x)


text:
############################################################################################
"Default Configuration for Metropolis sampler"
struct ConfigMetropolis <: AbstractConfiguration
    "Target Acceptance Rate"
    δ::Float64
    "Default size for tuning iterations in each cycle."
    window::ConfigTuningWindow
    function ConfigMetropolis(
        δ::Float64,
        window::ConfigTuningWindow
    )
        @argcheck 0.0 < δ <= 1.0 "Acceptance rate not bounded between 0 and 1"
        return new(
            δ,
            window,
        )
    end
end

"""
$(SIGNATURES)
Initialize Metropolis custom configurations.

# Examples
```julia
```

"""
function init(
    ::Type{AbstractConfiguration},
    mcmc::Type{Metropolis},
    objective::Objective,
    proposalconfig::ConfigProposal;
    ## Target Acceptance Rate
    δ=0.234,
    ## MCMC Phase tune variables based on standard HMC best practice target acceptance rate of 0.80
    window= ConfigTuningWindow(
        [1, 5, 1],
        Int.(floor.([75, 50, 25] * (0.80 / 0.234))),
        [Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()],
        )
)
    return ConfigMetropolis(
        δ,
        window,
    )
end

############################################################################################
# Export
export ConfigMetropolis


text:
const BOOK_URL = "https://api-fxpractice.oanda.com/v3/instruments/"

function books(book::String, pair::String, t::DateTime)
    df = Dates.DateFormat("yyyy-mm-ddTHH:MM:00Z")
    dt = floor(t,Dates.Minute(20)) #else "errorMessage" => "Specified time 2017-10-08T01:22:55.230Z is misaligned."
    url = BOOK_URL * pair * "/" * book * "?time=" * Dates.format(dt,df)
    Requests.json(get(url; headers = HEADERS))[book]
end

function books(book::String, pair::String)
    url = BOOK_URL * pair * "/" * book
    Requests.json(get(url; headers = HEADERS))[book]
end

function books(book::String, pair::String, t::String)
    url = BOOK_URL * pair * "/" * book * "?time=" * t
    Requests.json(get(url; headers = HEADERS))[book]
end

function books(book::String, pair::String, t0::DateTime, t1::DateTime)
    out = Array{Any,1}()
    while t0 < t1
        push!(out,books(book,pair,t0))
        t0 += Dates.Minute(20)
        sleep(.1)
    end
    return out
end

function books(book::String, pair::Array{String,1}, t0::DateTime, t1::DateTime)
    out = Dict{String,Array{Any,1}}()
    for i in pair
        out[i] = books(book,i,t0,t1)
    end
    return out
end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
using Documenter, Microbiome

makedocs(
    sitename = "Microbiome.jl",
    pages = [
        "Home" => "index.md",
        "Microbial Abundances" => "abundances.md",
        "Distances & Dissimilarity" => "distances.md",
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/BioJulia/Microbiome.jl.git",
    push_preview=true
)


text:
# a is emax, b is ed50, c is dose
function emax_dose_response_quote(M::Union{Int,Symbol}, T, isvec::NTuple{3,Bool}, track::NTuple{3,Bool}, sp::Bool)# partial::Bool, sp::Bool)
    #    f(a,b,c) =  a*c   / (b + c)
    # ∂f∂a(a,b,c) =    c   / (b + c)
    # ∂f∂b(a,b,c) = -a*c   / (b + c)^2
    # ∂f∂c(a,b,c) =  a * b / (b + c)^2
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    tracka, trackb, trackc = track
    aisvec, bisvec, cisvec = isvec
    need_loop = aisvec | bisvec | cisvec
    partial = tracka | trackb | trackc
    aexpr = aisvec ? :(a[i]) : :a
    bexpr = bisvec ? :(b[i]) : :b
    cexpr = cisvec ? :(c[i]) : :c
    head = quote end
    body = quote
        d = one($T) / ($bexpr + $cexpr)
        ∂f∂a = $cexpr * d
        f = $aexpr * ∂f∂a
    end
    return_f = need_loop ? :vf : :f
    return_expr = partial ? Expr(:tuple, return_f) : return_f
    if partial
        if trackb
            if need_loop && !bisvec
                push!(body.args, :(∂f∂b += -d*f))
            else
                push!(body.args, :(∂f∂b = -d*f))
            end
        end
        if trackc
            if need_loop && !cisvec
                push!(body.args, :(∂f∂c += a*b*d*d))
            else
                push!(body.args, :(∂f∂c = a*b*d*d))
            end
        end
        for (sym, isvec, track) in ( (:∂f∂a,aisvec,tracka),(:∂f∂b,bisvec,trackb),(:∂f∂c,cisvec,trackc) )
            track || continue
            if need_loop && !isvec
                push!(head.args, :($sym = zero($T)))
            elseif isvec
                vsym = Symbol(:v, sym)
                push!(head.args, PaddedMatrices.pointer_vector_expr( vsym, M, T, sp, :sptr ) )
                push!(body.args, :($vsym[i] = $sym))
                push!(return_expr.args, vsym)
            end
            isvec || push!(return_expr.args, sym)
        end
    end
    body = macroexpand(Base, :(@fastmath $body))
    final_ret = sp ? :((sptr, $return_expr)) : return_expr
    if need_loop
        quote
            $head
            @inbounds @simd for i in 1:$M
                $body
            end
            $final_ret
        end
    else
        quote
            $head
            $body
            $final_ret
        end
    end
end

@generated function emax_dose_response(
    a::Union{T,<:AbstractVector{T}},
    b::Union{T,<:AbstractVector{T}},
    c::Union{T,<:AbstractVector{T}},
    ::Val{track} = Val{(false,false,false)}()
) where {T, track}
    (L,isvec) = if a <: AbstractFixedSizeVector
        full_length(a), true
    elseif b <: AbstractFixedSizeVector
        full_length(b), true
    elseif c <: AbstractFixedSizeVector
        full_length(c), true
    elseif a != T
        :(length(a)), true
    elseif b != T
        :(length(b)), true
    elseif c != T
        :(length(c)), true
    else
        1, false
    end
    emax_dose_response_quote(L, T, isvec, track, false)
end
@generated function emax_dose_response(
    sptr::StackPointer,
    a::Union{T,<:AbstractVector{T}},
    b::Union{T,<:AbstractVector{T}},
    c::Union{T,<:AbstractVector{T}},
    ::Val{track} = Val{(false,false,false)}()
) where {T, track}
    (L,isvec) = if a <: AbstractFixedSizeVector
        full_length(a), true
    elseif b <: AbstractFixedSizeVector
        full_length(b), true
    elseif c <: AbstractFixedSizeVector
        full_length(c), true
    elseif a != T
        :(length(a)), true
    elseif b != T
        :(length(b)), true
    elseif c != T
        :(length(c)), true
    else
        1, false
    end
    emax_dose_response_quote(L, T, isvec, track, true)
end

function ITPExpectedValue_quote(
    M::Int, N::Int, T, @nospecialize(track::NTuple{<:Any,Bool}), partial::Bool, sp::Bool = false, P::Int = (M + Wm1) & ~Wm1
)
    Nparamargs = length(track)
    if Nparamargs == 2
        (track_β, track_κ) = track
        track_θ = false
        add_θ = false
    else
        @assert Nparamargs == 3
        (track_β, track_κ, track_θ) = track
        add_θ = true
    end
    #TODO: the first time equals θ, and the last equals θ + β; can skip exponential calculations.
    # M x N output
    # M total times
    # N total β and κs
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W - 1
    if !partial || (!track_β && !track_κ)
        return_expr = track_θ ? :(μ, ProbabilityModels.Reducer{:row}()) : :μ
        return_expr = sp ? :(sp, $return_expr) : return_expr
        return quote
            @inbounds for n ∈ 0:$(N-1)
                βₙ = β[n+1]
                κₙ = κ[n+1]
                $(add_θ ? :(θₙ = θ[n+1]) : nothing)
                @vvectorize $T for m ∈ 1:$M
                    μ[m + $P*n] = $(add_θ ? :(βₙ * (one($T) - exp(-κₙ * τ[m])) + θₙ ) : :(βₙ * (one($T) - exp(-κₙ * τ[m])))  )
                end
            end
            # println("ITP μ")
            # println(μ)
            $return_expr
        end
    end
    if track_β && track_κ
        return_expr = :(μ, StructuredMatrices.BlockDiagonalColumnView(∂β), StructuredMatrices.BlockDiagonalColumnView(∂κ))
        track_θ && push!(return_expr.args, ProbabilityModels.Reducer{:row}())
        return_expr = sp ? :(sp, $return_expr) : return_expr
        return quote
            vone = vbroadcast(Vec{$W,$T},one($T))
            @inbounds for n ∈ 0:$(N-1)
                βₙ = β[n+1]
                nκₙ = -κ[n+1]
                $(add_θ ? :(θₙ = θ[n+1]) : nothing)
                $(macroexpand(ProbabilityModels, quote
                              LoopVectorization.@vvectorize $T for m ∈ 1:$M
                              tₘ = τ[m]
                              ℯκt = exp( nκₙ * tₘ)
                              βₙℯκt = βₙ * ℯκt
                              Omℯκt = vone - ℯκt
                              ∂β[m + $P*n] = Omℯκt
                              μ[m + $P*n] = $(add_θ ? :(βₙ * Omℯκt + θₙ) : :(βₙ - βₙℯκt) )
                              ∂κ[m + $P*n] = βₙℯκt * tₘ
                              end
                end))
            end
            $return_expr
        end
    elseif track_β
        return_expr = :(μ, StructuredMatrices.BlockDiagonalColumnView(∂β))
        track_θ && push!(return_expr.args, ProbabilityModels.Reducer{:row}())
        return_expr = sp ? :(sp, $return_expr) : return_expr

        return quote
            @inbounds for n ∈ 0:$(N-1)
                βₙ = β[n+1]
                κₙ = κ[n+1]
                $(add_θ ? :(θₙ = θ[n+1]) : nothing)
                @vvectorize $T for m ∈ 1:$M
                    tₘ = τ[m]
                    ℯκt = exp(- κₙ * tₘ)
                    Omℯκt = one($T) - ℯκt
                    ∂β[m + $P*n] = Omℯκt
                    μ[m + $P*n] = $(add_θ ? :(βₙ * Omℯκt + θₙ) : :(βₙ * Omℯκt) )
                end
            end
            $return_expr
        end
    else # track_κ
        return_expr = :(μ, StructuredMatrices.BlockDiagonalColumnView(∂κ))
        track_θ && push!(return_expr.args, ProbabilityModels.Reducer{:row}())
        return_expr = sp ? :(sp, $return_expr) : return_expr

        return quote
            @inbounds for n ∈ 0:$(N-1)
                βₙ = β[n+1]
                κₙ = κ[n+1]
                $(add_θ ? :(θₙ = θ[n+1]) : nothing)
                @vvectorize $T for m ∈ 1:$M
                    tₘ = τ[m]
                    ℯκt = exp(- κₙ * tₘ)
                    βₙℯκt = βₙ * ℯκt
                    μ[m + $P*n] = $(add_θ ? :(βₙ - βₙℯκt + θₙ) : :(βₙ - βℯκt) )
                    ∂κ[m + $P*n] = βₙℯκt * tₘ
                end
            end
            $return_expr
        end
    end
end

@generated function ITPExpectedValue(
            τ::PaddedMatrices.AbstractFixedSizeVector{R},
            β::PaddedMatrices.AbstractFixedSizeVector{N,T},
            κ::PaddedMatrices.AbstractFixedSizeVector{N,T}
        ) where {R,N,T}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    quote
        μ = FixedSizeMatrix{$M,$N,$T,$M}(undef)
        $(ITPExpectedValue_quote(M, N, T, (false, false), false, false, M))
    end
end
@generated function ∂ITPExpectedValue(
            τ::PaddedMatrices.AbstractFixedSizeVector{R},
            β::PaddedMatrices.AbstractFixedSizeVector{N,T},
            κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
            ::Val{track}
        ) where {R,N,T,track}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    (track_β, track_κ) = track
    q = quote
        μ = FixedSizeMatrix{$M,$N,$T,$M}(undef)
    end
    if track_β
        push!(q.args, :(∂β = FixedSizeMatrix{$M,$N,$T,$M}(undef)))
    end
    if track_κ
        push!(q.args, :(∂κ = FixedSizeMatrix{$M,$N,$T,$N}(undef)))
    end
    push!(q.args, ITPExpectedValue_quote(M, N, T, track, true,false,M))
    q
end
@generated function ITPExpectedValue(
            τ::PaddedMatrices.AbstractFixedSizeVector{R},
            β::PaddedMatrices.AbstractFixedSizeVector{N,T},
            κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
            θ::PaddedMatrices.AbstractFixedSizeVector{N,T}
        # ) where {R,T,N}
        ) where {R,N,T}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    quote
        μ = FixedSizeMatrix{$M,$N,$T,$M}(undef)
        $(ITPExpectedValue_quote(M, N, T, (false, false, false), false,false,M))
    end
end
@generated function ∂ITPExpectedValue(
            τ::PaddedMatrices.AbstractFixedSizeVector{R},
            β::PaddedMatrices.AbstractFixedSizeVector{N,T},
            κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
            θ::PaddedMatrices.AbstractFixedSizeVector{N,T},
            ::Val{track}
        ) where {R,T,N,track}
        # ) where {R,N,T,track}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    (track_β, track_κ, track_θ) = track
    q = quote
        μ = FixedSizeMatrix{$M,$N,$T,$M}(undef)
    end
    if track_β
        push!(q.args, :(∂β = FixedSizeMatrix{$M,$N,$T,$M}(undef)))
    end
    if track_κ
        push!(q.args, :(∂κ = FixedSizeMatrix{$M,$N,$T,$M}(undef)))
    end
    push!(q.args, ITPExpectedValue_quote(M, N, T, track, true,false,M))
    q
end
@generated function ITPExpectedValue(
    sp::StackPointer,
    τ::PaddedMatrices.AbstractFixedSizeVector{R},
    β::PaddedMatrices.AbstractFixedSizeVector{N,T},
    κ::PaddedMatrices.AbstractFixedSizeVector{N,T}
) where {R,N,T}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    quote
        (sp, μ) = PtrMatrix{$M,$N,$T,$M}(sp)
        $(ITPExpectedValue_quote(M, N, T, (false, false), false, true, M))
    end
end
@generated function ∂ITPExpectedValue(
    sp::StackPointer,
    τ::PaddedMatrices.AbstractFixedSizeVector{R},
    β::PaddedMatrices.AbstractFixedSizeVector{N,T},
    κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
    ::Val{track}
) where {R,N,T,track}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    (track_β, track_κ) = track
    q = quote
        (sp, μ) = PtrMatrix{$M,$N,$T,$M}(sp)
    end
    if track_β
        push!(q.args, :((sp, ∂β) = PtrMatrix{$M,$N,$T,$M}(sp)))
    end
    if track_κ
        push!(q.args, :((sp, ∂κ) = PtrMatrix{$M,$N,$T,$M}(sp)))
    end
    push!(q.args, ITPExpectedValue_quote(M, N, T, track, true, true, M))
    q
end
@generated function ITPExpectedValue(
    sp::StackPointer,
    τ::PaddedMatrices.AbstractFixedSizeVector{R},
    β::PaddedMatrices.AbstractFixedSizeVector{N,T},
    κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
    θ::PaddedMatrices.AbstractFixedSizeVector{N,T}
        # ) where {R,T,N}
) where {R,N,T}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
    quote
        (sp, μ) = PtrMatrix{$M,$N,$T,$M}(sp)
        $(ITPExpectedValue_quote(M, N, T, (false, false, false), false, true, M))
    end
end
@generated function ∂ITPExpectedValue(
    sp::StackPointer,
    τ::PaddedMatrices.AbstractFixedSizeVector{R},
    β::PaddedMatrices.AbstractFixedSizeVector{N,T},
    κ::PaddedMatrices.AbstractFixedSizeVector{N,T},
    θ::PaddedMatrices.AbstractFixedSizeVector{N,T},
    ::Val{track}
) where {R,T,N,track}
    # ) where {R,N,T,track}
    if isa(R, AbstractRange)
        M = length(R)
    else
        M = R
    end
#    @show τ, β, κ, θ, track
    (track_β, track_κ, track_θ) = track
    q = quote
#        @show reinterpret(Int,pointer(sp)) - reinterpret(Int,pointer(ProbabilityModels.STACK_POINTER))
        (sp, μ) = PtrMatrix{$M,$N,$T,$M}(sp)
    end
    if track_β
        push!(q.args, :((sp,∂β) = PtrMatrix{$M,$N,$T,$M}(sp)))
    end
    if track_κ
        push!(q.args, :((sp,∂κ) = PtrMatrix{$M,$N,$T,$M}(sp)))
    end
#    push!(q.args, :(@show reinterpret(Int,pointer(sp)) - reinterpret(Int,pointer(ProbabilityModels.STACK_POINTER))))
    push!(q.args, ITPExpectedValue_quote(M, N, T, track, true, true, M))
    q
end

struct Domains{S} end
Base.@pure Domains(S::NTuple{N,Int}) where {N} = Domains{S}()
Base.@pure Domains(S::Vararg{Int,N}) where {N} = Domains{S}()
Base.getindex(::Domains{S}, i) where {S} = S[i]
Base.length(::Domains{S}) where {S} = length(S)
Base.eltype(::Domains{S}) where {S} = eltype(S)
@inline Base.iterate(::Domains{S}) where {S} = @inbounds (S[1], 2)
@inline function Base.iterate(::Domains{S}, i) where {S}
    i > length(S) && return nothing
    S[i], i+1
end
Base.Array(::Domains{S}) where {S} = [S...]

function HierarchicalCentering_quote(M::Int, @nospecialize(T), μisvec::Bool, σisvec::Bool, (track_y, track_μ, track_σ)::NTuple{3,Bool}, partial, sp::Bool)
    μsym = μisvec ? :(μ[m]) : :μ
    σsym = σisvec ? :(σ[m]) : :σ
    if !partial
        if sp
            return quote
                xout = PtrVector{$M,$T,$M}(pointer(sp,$T))
                @vvectorize $T for m ∈ 1:$M
                    xout[m] = $μsym + $σsym * y[m]
                end
                sp + $(sizeof(T)*M), xout
            end
        else
            return quote
                xout = FixedSizeVector{$M,$T}(undef)
                @vvectorize $T for m ∈ 1:$M
                    xout[m] = $μsym + $σsym * y[m]
                end
                ConstantFixedSizeVector(xout)
            end
        end
    end
    # partial is true
    if sp
        q = quote
            xout = PtrVector{$M,$T,$M}(pointer(sp,$T))
            sp += $(sizeof(T)*M)
        end
        return_expr = Expr(:tuple, :xout )
    else
        q = quote xout = FixedSizeVector{$M,$T}(undef) end
        return_expr = Expr(:tuple, :(ConstantFixedSizeVector(xout)) )
    end
    loop_body = quote xout[m] = $μsym + $σsym * y[m] end
    if track_y
        if σisvec
            push!(return_expr.args, :(Diagonal(σ)) )
        else
            push!(return_expr.args, :(LinearAlgebra.UniformScaling(σ)) )
        end
    end
    if track_μ
        if μisvec
            push!(return_expr.args, :(One()))
        else
            push!(return_expr.args, :(Reducer{true}()))
        end
    end
    if track_σ
        if σisvec
            push!(return_expr.args, :(Diagonal(y)) )
        else
            # push!(return_expr.args, :∂σ)
            # push!(q.args, :(∂σ = zero($T)) )
            # push!(loop_body.args, :(∂σ += y[m]) )
            push!(return_expr.args, :(y) )
        end
    end
    push!(q.args, quote
          @vectorize $T for m ∈ 1:$M
              $loop_body
          end
          $(ProbabilityDistributions.return_expression(return_expr, sp))
      end)
    q
end
"""
This method takes a "Domains" argument, which is a tuple indicating number within each domain.

Depending on whether σ is a scalar or a vector:
∂y is either UniformScaling(σ), or a vector with the same type as y, with replciated σs.

if μ is a vector,
∂μ is Domains{S}()
which, when multiplied by a vector of length = sum(S), reduces the corresponding elements.
Otherwise, is a reducer that sums the vector it is multiplied with.

∂σ is a vector with the same type as y, and...
if σ is a scalar, will dot product on multiplication
if σ is a vector, will do length(σ) mini dot products.
"""
function HierarchicalCentering_quote(
    M::Int, P::Int, @nospecialize(T), μisvec::Bool, σisvec::Bool, @nospecialize(S::NTuple{<:Any,Int}), (track_y, track_μ, track_σ)::NTuple{3,Bool}, sp::Bool
)
    N = length(S)
    @assert sum(S) == M
    @assert μisvec | σisvec
    q = quote end
    if sp
        push!(q.args, :(xout = PtrVector{$M,$T,$M}(pointer(sp,$T))))
        push!(q.args, :(sp += $(VectorizationBase.align(sizeof(T)*M))))
        if track_y && σisvec
            push!(q.args, :(∂y = PtrVector{$M,$T,$M}(pointer(sp,$T))))
            push!(q.args, :(sp += $(VectorizationBase.align(sizeof(T)*M))))
        end
    else
        outtup = Expr(:tuple,)
        if track_y && σisvec
            ∂yexpr = Expr(:tuple,)
        end
    end

    ind = 0
    for (j,s) ∈ enumerate(S)
        for i ∈ 1:s
            ind += 1
            sym = gensym()
            if sp
                if μisvec && σisvec
                    push!(q.args, :($sym = μ[$j] + σ[$j] * y[$ind]))
                    track_y && push!(q.args, :(∂y[$ind] = σ[$j]))
                elseif μisvec
                    push!(q.args, :($sym = μ[$j] + σ * y[$ind]))
                else #if σisvec
                    push!(q.args, :($sym = μ + σ[$j] * y[$ind]))
                    track_y && push!(q.args, :(∂y[$ind] = σ[$j]))
                end
                push!(q.args, :(xout[$ind] = $sym))
            else
                if μisvec && σisvec
                    push!(q.args, :($sym = μ[$j] + σ[$j] * y[$ind]))
                    track_y && push!(∂yexpr.args, :(σ[$j]))
                elseif μisvec
                    push!(q.args, :($sym = μ[$j] + σ * y[$ind]))
                else #if σisvec
                    push!(q.args, :($sym = μ + σ[$j] * y[$ind]))
                    track_y && push!(∂yexpr.args, :(σ[$j]))
                end
                push!(outtup.args, sym)
            end
        end
    end
    if !sp
        for p ∈ M+1:P
            push!(outtup.args, zero(T))
        end
        push!(q.args, :(xout = ConstantFixedSizeVector{$M,$T,$P}($outtup) ))
    end
    return_expr = Expr(:tuple, :xout )
    if track_y
        if σisvec
            if !sp
                for p ∈ M+1:P
                    push!(∂yexpr.args, zero(T))
                end
                push!(q.args, :( ∂y = ConstantFixedSizeVector{$M,$T,$P}($∂yexpr) ))
            end
            push!(return_expr.args, :(LinearAlgebra.Diagonal(∂y)) )
        else
            push!(return_expr.args,  :(LinearAlgebra.UniformScaling(σ)) )
        end
    end
    if track_μ
        push!(return_expr.args, μisvec ? :(Reducer{$S}()) : :(Reducer{true}()) )
    end
    if track_σ
        push!(return_expr.args, σisvec ? :(ReducerWrapper{$S}(y)) :  :y )
    end
    quote
        @fastmath @inbounds begin
            $q
        end
        $(ProbabilityDistributions.return_expression(return_expr,sp))
    end
end

@generated function HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T},
            μ::Union{T, <: AbstractFixedSizeVector{M,T}},
            σ::Union{T, <: AbstractFixedSizeVector{M,T}}
        ) where {M,T}
        # ) where {T,M}

    HierarchicalCentering_quote(M, T, μ <: AbstractFixedSizeVector, σ <: AbstractFixedSizeVector, (false,false,false), false, false)
end

@generated function ∂HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T},
            μ::Union{T, <: AbstractFixedSizeVector{M,T}},
            σ::Union{T, <: AbstractFixedSizeVector{M,T}},
            ::Val{track}
        # ) where {T,M,track}
        ) where {M,T,track}

    HierarchicalCentering_quote(M, T, μ <: AbstractFixedSizeVector, σ <: AbstractFixedSizeVector, track, true, false)
end

@generated function HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::AbstractFixedSizeVector{N,T},
            σ::AbstractFixedSizeVector{N,T},
            ::Domains{S}
        ) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, true, S, (false,false,false),false)
end
@generated function HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::AbstractFixedSizeVector{N,T},
            σ::T,
            ::Domains{S}
        ) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, false, S, (false,false,false),false)
end
@generated function HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::T,
            σ::AbstractFixedSizeVector{N,T},
            ::Domains{S}
        ) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, false, true, S, (false,false,false),false)
end
@generated function ∂HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::AbstractFixedSizeVector{N,T},
            σ::AbstractFixedSizeVector{N,T},
            ::Domains{S}, ::Val{track}
        ) where {M,N,S,T,track,P}
        # ) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, true, S, track,false)
end
@generated function ∂HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::AbstractFixedSizeVector{N,T},
            σ::T,
            ::Domains{S}, ::Val{track}
        ) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, false, S, track,false)
end
@generated function ∂HierarchicalCentering(
            #x::AbstractFixedSizeVector{M,T},
            y::AbstractFixedSizeVector{M,T,P},
            μ::T,
            σ::AbstractFixedSizeVector{N,T},
            ::Domains{S}, ::Val{track}
        ) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, false, true, S, track,false)
end

@generated function HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T},
    μ::Union{T, <: AbstractFixedSizeVector{M,T}},
    σ::Union{T, <: AbstractFixedSizeVector{M,T}}
) where {M,T}
    # ) where {T,M}
    HierarchicalCentering_quote(M, T, μ <: AbstractFixedSizeVector, σ <: AbstractFixedSizeVector, (false,false,false), false, true)
end

@generated function ∂HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T},
    μ::Union{T, <: AbstractFixedSizeVector{M,T}},
    σ::Union{T, <: AbstractFixedSizeVector{M,T}},
    ::Val{track}
    # ) where {T,M,track}
) where {M,T,track}
    HierarchicalCentering_quote(M, T, μ <: AbstractFixedSizeVector, σ <: AbstractFixedSizeVector, track, true, true)
end

@generated function HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::AbstractFixedSizeVector{N,T},
    σ::AbstractFixedSizeVector{N,T},
    ::Domains{S}
) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, true, S, (false,false,false),true)
end
@generated function HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::AbstractFixedSizeVector{N,T},
    σ::T,
    ::Domains{S}
) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, false, S, (false,false,false),true)
end
@generated function HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::T,
    σ::AbstractFixedSizeVector{N,T},
    ::Domains{S}
) where {M,N,T,P,S}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, false, true, S, (false,false,false),true)
end
@generated function ∂HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::AbstractFixedSizeVector{N,T},
    σ::AbstractFixedSizeVector{N,T},
    ::Domains{S}, ::Val{track}
) where {M,N,S,T,track,P}
    # ) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, true, S, track, true)
end
@generated function ∂HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::AbstractFixedSizeVector{N,T},
    σ::T,
    ::Domains{S}, ::Val{track}
) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, true, false, S, track,true)
end
@generated function ∂HierarchicalCentering(
    sp::StackPointer,
    y::AbstractFixedSizeVector{M,T,P},
    μ::T,
    σ::AbstractFixedSizeVector{N,T},
    ::Domains{S}, ::Val{track}
) where {M,N,T,S,track,P}
    @assert length(S) == N
    HierarchicalCentering_quote(M, P, T, false, true, S, track,true)
end

#@support_stack_pointer HierarchicalCentering
#@support_stack_pointer ∂HierarchicalCentering

struct ReshapeAdjoint{S} end
function ∂vec(a::AbstractMutableFixedSizeArray{S}) where {S}
    vec(a), ReshapeAdjoint{S}()
end
function Base.:*(A::AbstractMutableFixedSizeArray, ::ReshapeAdjoint{S}) where S
    reshape(A, Val{S}())
end

# strip_typeof(x::Symbol) = Symbol(String(x)[8:end-1])


text:
using StanRun, Test

@testset "cmdstan run and results" begin
    # setup environment
    MODELDIR = mktempdir()
    SAMPLEDIR = mktempdir()
    SRC = joinpath(MODELDIR, "test.stan")
    cp(joinpath(@__DIR__, "test.stan"), SRC)

    # run a model
    model = StanModel(SRC)
    @test repr(model) ==
        "Stan model at $(SRC)\n    (CmdStan home: $(StanRun.get_cmdstan_home()))"
    exec_path = StanRun.ensure_executable(model)
    @test isfile(exec_path)
    n_chains = 5
    OUTPUT_BASE = joinpath(SAMPLEDIR, "test")
    time_before = time()
    @test stan_compile(model) ≡ nothing
    chains = stan_sample(model, (N = 100, x = randn(100)), 5; output_base = OUTPUT_BASE)
    for (sample, logfile) in chains
        @test ctime(sample) ≥ time_before
        @test ctime(sample) ≥ time_before
    end
    @test first.(chains) == sort(StanRun.find_samples(OUTPUT_BASE)) ==
        [joinpath(SAMPLEDIR, "test_chain_$(i).csv") for i in 1:n_chains]
end

@testset "unset cmdstan environment" begin
    withenv("JULIA_CMDSTAN_HOME" => nothing) do
        @test_throws ErrorException StanRun.get_cmdstan_home()
    end
end

@testset "model error and message" begin
    model = StanModel(joinpath(@__DIR__, "test_incorrect.stan"))
    try
        stan_compile(model)
    catch e
        @test e isa StanRun.StanModelError
        @test occursin("Variable \"x\" does not exist", e.message)
        io = IOBuffer()
        showerror(io, e)
        e_repr = String(take!(io))
        @test occursin("error when compiling", e_repr)
    end
end


text:
include("./LatticeQCD.jl")
using Plots
using .LatticeQCD
using Random
using Dates
using JLD
#using GFlops

@time demo()


text:
# copied from https://github.com/JuliaLang/julia/pull/30924
# modified to be in compliance with C99: http://port70.net/~nsz/c/c99/n1256.html#6.7.2.2
module CEnum

import Core.Intrinsics.bitcast
export @cenum

function namemap end
function name_value_pairs end

"""
    Cenum{T<:Integer}
The abstract supertype of all enumerated types defined with [`@cenum`](@ref).
"""
abstract type Cenum{T<:Integer} end

basetype(::Type{<:Cenum{T}}) where {T<:Integer} = T

(::Type{T})(x::Cenum{T2}) where {T<:Integer,T2<:Integer} = T(bitcast(T2, x))::T
Base.cconvert(::Type{T}, x::Cenum{T2}) where {T<:Integer,T2<:Integer} = T(x)
Base.write(io::IO, x::Cenum{T}) where {T<:Integer} = write(io, T(x))
Base.read(io::IO, ::Type{T}) where {T<:Cenum} = T(read(io, basetype(T)))

Base.isless(x::T, y::T) where {T<:Cenum} = isless(basetype(T)(x), basetype(T)(y))

Base.Symbol(x::Cenum)::Symbol = get(namemap(typeof(x)), Integer(x), :UnknownMember)

Base.print(io::IO, x::Cenum) = print(io, Symbol(x))

function Base.show(io::IO, x::Cenum)
    sym = Symbol(x)
    if !get(io, :compact, false)
        from = get(io, :module, Main)
        def = typeof(x).name.module
        if from === nothing || !Base.isvisible(sym, def, from)
            show(io, def)
            print(io, ".")
        end
    end
    print(io, sym)
end

function Base.show(io::IO, ::MIME"text/plain", x::Cenum)
    print(io, x, "::")
    show(IOContext(io, :compact => true), typeof(x))
    print(io, " = ")
    show(io, Integer(x))
end

function Base.show(io::IO, ::MIME"text/plain", t::Type{<:Cenum})
    print(io, "Cenum ")
    Base.show_datatype(io, t)
    print(io, ":")
    for (s, i) in name_value_pairs(t)
        print(io, "\n", Symbol(s), " = ")
        show(io, Integer(i))
    end
end

# give Cenum types scalar behavior in broadcasting
Base.broadcastable(x::Cenum) = Ref(x)

@noinline enum_argument_error(typename, x) = throw(ArgumentError(string("input value out of range for Cenum $(typename): $x")))

macro cenum(T, syms...)
    if isempty(syms)
        throw(ArgumentError("no arguments given for Cenum $T"))
    end
    basetype = Int32
    typename = T
    if isa(T, Expr) && T.head == :(::) && length(T.args) == 2 && isa(T.args[1], Symbol)
        typename = T.args[1]
        basetype = Core.eval(__module__, T.args[2])
        if !isa(basetype, DataType) || !(basetype <: Integer) || !isbitstype(basetype)
            throw(ArgumentError("invalid base type for Cenum $typename, $T=::$basetype; base type must be an integer primitive type"))
        end
    elseif !isa(T, Symbol)
        throw(ArgumentError("invalid type expression for Cenum $T"))
    end
    seen = Set{Symbol}()
    name_values = Tuple{Symbol,basetype}[]
    namemap = Dict{basetype,Symbol}()
    lo = hi = 0
    i = zero(basetype)

    if length(syms) == 1 && syms[1] isa Expr && syms[1].head == :block
        syms = syms[1].args
    end
    for s in syms
        s isa LineNumberNode && continue
        if isa(s, Symbol)
            if i == typemin(basetype) && !isempty(name_values)
                throw(ArgumentError("overflow in value \"$s\" of Cenum $typename"))
            end
        elseif isa(s, Expr) &&
               (s.head == :(=) || s.head == :kw) &&
               length(s.args) == 2 && isa(s.args[1], Symbol)
            i = Core.eval(__module__, s.args[2]) # allow exprs, e.g. uint128"1"
            if !isa(i, Integer)
                throw(ArgumentError("invalid value for Cenum $typename, $s; values must be integers"))
            end
            i = convert(basetype, i)
            s = s.args[1]
        else
            throw(ArgumentError(string("invalid argument for Cenum ", typename, ": ", s)))
        end
        if !Base.isidentifier(s)
            throw(ArgumentError("invalid name for Cenum $typename; \"$s\" is not a valid identifier"))
        end
        haskey(namemap, i) || (namemap[i] = s;)
        if s in seen
            throw(ArgumentError("name \"$s\" in Cenum $typename is not unique"))
        end
        push!(seen, s)
        push!(name_values, (s,i))
        if length(name_values) == 1
            lo = hi = i
        else
            lo = min(lo, i)
            hi = max(hi, i)
        end
        i += oneunit(i)
    end
    blk = quote
        # enum definition
        Base.@__doc__(primitive type $(esc(typename)) <: Cenum{$(basetype)} $(sizeof(basetype) * 8) end)
        function $(esc(typename))(x::Integer)
            x ≤ typemax(x) || enum_argument_error($(Expr(:quote, typename)), x)
            return bitcast($(esc(typename)), convert($(basetype), x))
        end
        CEnum.namemap(::Type{$(esc(typename))}) = $(esc(namemap))
        CEnum.name_value_pairs(::Type{$(esc(typename))}) = $(esc(name_values))
        Base.typemin(x::Type{$(esc(typename))}) = $(esc(typename))($lo)
        Base.typemax(x::Type{$(esc(typename))}) = $(esc(typename))($hi)
        let insts = (Any[ $(esc(typename))(v[2]) for v in $name_values ]...,)
            Base.instances(::Type{$(esc(typename))}) = insts
        end
    end
    if isa(typename, Symbol)
        for (sym, i) in name_values
            push!(blk.args, :(const $(esc(sym)) = $(esc(typename))($i)))
        end
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end

include("cenumoperators.jl")

end # module

max_stars_repo_name: MiguelMarcelino/PyEnum

max_issues_repo_name: MiguelMarcelino/PyEnum

max_forks_repo_name: MiguelMarcelino/PyEnum


text:
module TestAsNumber

using Test
using ArbFloats

af0 = zero(ArbFloat)
af1 = one(ArbFloat)
af2 = af1 + af1
afhalf = af1 / 2
afNaN = ArbFloat{precision(ArbFloat)}(NaN)
afInf = ArbFloat(Inf)

@test af1 <  af2
@test af1 <= af2
@test af2 >  af1
@test af2 >= af1

@test min(af1,af2) != af2
@test max(af1,af2) == af2

@test  asin(log(acosh(cosh(exp(sin(afhalf)))))) > afhalf
@test  upperbound(asin(log(atanh(tanh(exp(sin(afhalf))))))) >= afhalf
@test  lowerbound(asin(log(atanh(tanh(exp(sin(afhalf))))))) <= afhalf

@test midpoint(afhalf) > midpoint(prevfloat(afhalf))
@test midpoint(afhalf) < midpoint(nextfloat(afhalf))

end # module


text:
module LiteQTL

#put all external pacakges here.
using DelimitedFiles
using LinearAlgebra
using Base.Threads
using CUDA
using CSV
using DocStringExtensions
using Statistics
using Distributions
using Dates

#put all your source file here.
include("common_func.jl")
include("data_io.jl")
export get_geno_data, get_pheno_data
include("util.jl")
include("cpu.jl")
include("gpu.jl")
include("scan.jl")
export scan

end # module


text:
"""
Robust initialization method for model parameters in Hamiltonian samplers.
"""
struct SampleFromUniform <: AbstractSampler end
struct SampleFromPrior <: AbstractSampler end

getspace(::Union{SampleFromPrior, SampleFromUniform}) = ()

# Initializations.
init(dist, ::SampleFromPrior) = rand(dist)
init(dist, ::SampleFromUniform) = istransformable(dist) ? inittrans(dist) : rand(dist)

init(dist, ::SampleFromPrior, n::Int) = rand(dist, n)
function init(dist, ::SampleFromUniform, n::Int)
    return istransformable(dist) ? inittrans(dist, n) : rand(dist, n)
end

"""
    has_eval_num(spl::AbstractSampler)

Check whether `spl` has a field called `eval_num` in its state variables or not.
"""
has_eval_num(spl::SampleFromUniform) = false
has_eval_num(spl::SampleFromPrior) = false
has_eval_num(spl::AbstractSampler) = :eval_num in fieldnames(typeof(spl.state))

"""
An abstract type that mutable sampler state structs inherit from.
"""
abstract type AbstractSamplerState end

"""
    Sampler{T}

Generic interface for implementing inference algorithms.
An implementation of an algorithm should include the following:

1. A type specifying the algorithm and its parameters, derived from InferenceAlgorithm
2. A method of `sample` function that produces results of inference, which is where actual inference happens.

DynamicPPL translates models to chunks that call the modelling functions at specified points.
The dispatch is based on the value of a `sampler` variable.
To include a new inference algorithm implements the requirements mentioned above in a separate file,
then include that file at the end of this one.
"""
mutable struct Sampler{T, S<:AbstractSamplerState} <: AbstractSampler
    alg      ::  T
    info     ::  Dict{Symbol, Any} # sampler infomation
    selector ::  Selector
    state    ::  S
end
Sampler(alg) = Sampler(alg, Selector())
Sampler(alg, model::Model) = Sampler(alg, model, Selector())
Sampler(alg, model::Model, s::Selector) = Sampler(alg, model, s)

# AbstractMCMC interface for SampleFromUniform and SampleFromPrior

function AbstractMCMC.step!(
    rng::Random.AbstractRNG,
    model::Model,
    sampler::Union{SampleFromUniform,SampleFromPrior},
    ::Integer,
    transition;
    kwargs...
)
    vi = VarInfo()
    model(vi, sampler)
    return vi
end


text:
using SFML

server = TcpListener()
set_blocking(server, false)

listen(server, 53000)
println("Listening on port 53000")

client = SocketTCP()
while accept(server, client) != SOCKET_DONE
end

packet = Packet()
println(receive(client, packet))

println(read_string(packet))
println(read_string(packet))


text:
using ITensors, Test, Random

using ITensors: nsite, set_nsite!

@testset "Basic DMRG" begin
  @testset "Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1", N)

    os = OpSum()
    for j in 1:(N - 1)
      add!(os, "Sz", j, "Sz", j + 1)
      add!(os, 0.5, "S+", j, "S-", j + 1)
      add!(os, 0.5, "S-", j, "S+", j + 1)
    end
    H = MPO(os, sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps, 10, 20, 40)
    mindim!(sweeps, 1, 10)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps, 1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy, psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "QN-conserving Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1", N; conserve_qns=true)

    os = OpSum()
    for j in 1:(N - 1)
      os += "Sz", j, "Sz", j + 1
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
    end
    H = MPO(os, sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites, state; linkdims=4)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps, 10, 20, 40)
    mindim!(sweeps, 1, 10)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps, 1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy, psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "QN-conserving Spin-one Heisenberg with disk caching" begin
    N = 10
    sites = siteinds("S=1", N; conserve_qns=true)

    os = OpSum()
    for j in 1:(N - 1)
      os += "Sz", j, "Sz", j + 1
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
    end
    H = MPO(os, sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites, state; linkdims=4)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps, 10, 20, 40)
    mindim!(sweeps, 1, 10)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps, 1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy, psi = dmrg(H, psi, sweeps; outputlevel=0, write_when_maxdim_exceeds=15)
    @test energy < -12.0
  end

  @testset "ProjMPO with disk caching" begin
    N = 10
    sites = siteinds("S=1", N; conserve_qns=true)

    os = OpSum()
    for j in 1:(N - 1)
      os += "Sz", j, "Sz", j + 1
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
    end
    H = MPO(os, sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites, state; linkdims=4)
    PH = ProjMPO(H)

    PHc = copy(PH)

    n = 4
    orthogonalize!(psi, n)
    position!(PH, psi, n)
    PHdisk = ITensors.disk(PH)

    @test length(PH) == N
    @test length(PHdisk) == N
    @test ITensors.site_range(PH) == n:(n + 1)
    @test eltype(PH) == Float64
    @test size(PH) == (3^2 * 4^2, 3^2 * 4^2)
    @test PH.lpos == n - 1
    @test PH.rpos == n + 2
    @test PHc.lpos == 0
    @test PHc.rpos == N + 1
    @test rproj(PH) ≈ rproj(PHdisk)
    @test PHdisk.LR isa ITensors.DiskVector{ITensor}
    @test PHdisk.LR[PHdisk.rpos] ≈ PHdisk.Rcache
    position!(PH, psi, N)
    @test PH.lpos == N - 1
  end

  @testset "ProjMPO: nsite" begin
    N = 10
    sites = siteinds("S=1", N)

    os1 = OpSum()
    for j in 1:(N - 1)
      os1 += 0.5, "S+", j, "S-", j + 1
      os1 += 0.5, "S-", j, "S+", j + 1
    end
    os2 = OpSum()
    for j in 1:(N - 1)
      os2 += "Sz", j, "Sz", j + 1
    end
    H1 = MPO(os1, sites)
    H2 = MPO(os2, sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites, state; linkdims=4)
    PH1 = ProjMPO(H1)
    PH = ProjMPOSum([H1, H2])
    PH1c = copy(PH1)
    PHc = copy(PH)
    @test nsite(PH1) == 2
    @test nsite(PH) == 2
    @test nsite(PH1c) == 2
    @test nsite(PHc) == 2

    set_nsite!(PH1, 3)
    @test nsite(PH1) == 3
    @test nsite(PH1c) == 2
    @test nsite(PHc) == 2

    set_nsite!(PH, 4)
    @test nsite(PH) == 4
    @test nsite(PH1c) == 2
    @test nsite(PHc) == 2
  end

  @testset "Transverse field Ising" begin
    N = 32
    sites = siteinds("S=1/2", N)
    Random.seed!(432)
    psi0 = randomMPS(sites)

    os = OpSum()
    for j in 1:N
      j < N && add!(os, -1.0, "Z", j, "Z", j + 1)
      add!(os, -1.0, "X", j)
    end
    H = MPO(os, sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps, 10, 20)
    cutoff!(sweeps, 1E-12)
    noise!(sweeps, 1E-10)
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 1.0 - 1.0 / sin(π / (4 * N + 2))
    @test abs((energy - energy_exact) / energy_exact) < 1e-4
  end

  @testset "Compact Sweeps syntax" begin
    N = 32
    sites = siteinds("S=1/2", N)
    Random.seed!(432)
    psi0 = randomMPS(sites)

    function ising(N; h=1.0)
      os = OpSum()
      for j in 1:N
        j < N && (os -= ("Z", j, "Z", j + 1))
        os -= h, "X", j
      end
      return os
    end

    h = 1.0
    H = MPO(ising(N; h=h), sites)
    energy, psi = dmrg(
      H, psi0; nsweeps=5, maxdim=[10, 20], cutoff=1e-12, noise=1e-10, outputlevel=0
    )

    energy_exact = 1.0 - 1.0 / sin(π / (4 * N + 2))
    @test abs((energy - energy_exact) / energy_exact) < 1e-4
  end

  @testset "Transverse field Ising, conserve Sz parity" begin
    N = 32
    sites = siteinds("S=1/2", N; conserve_szparity=true)
    Random.seed!(432)

    state = [isodd(j) ? "↑" : "↓" for j in 1:N]
    psi0 = randomMPS(sites, state)

    os = OpSum()
    for j in 1:N
      j < N && add!(os, -1.0, "X", j, "X", j + 1)
      add!(os, -1.0, "Z", j)
    end
    H = MPO(os, sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps, 10, 20)
    cutoff!(sweeps, 1E-12)
    noise!(sweeps, 1E-10)
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 1.0 - 1.0 / sin(π / (4 * N + 2))
    @test abs((energy - energy_exact) / energy_exact) < 1e-4
  end

  @testset "DMRGObserver" begin

    # Test that basic constructors work
    observer = DMRGObserver()
    observer = DMRGObserver(; minsweeps=2, energy_tol=1E-4)

    # Test in a DMRG calculation
    N = 10
    sites = siteinds("S=1/2", N)
    Random.seed!(42)
    psi0 = randomMPS(sites)

    os = OpSum()
    for j in 1:(N - 1)
      os += -1, "Sz", j, "Sz", j + 1
    end
    for j in 1:N
      os += -0.2, "Sx", j
    end
    H = MPO(os, sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps, 10)
    cutoff!(sweeps, 1E-12)

    observer = DMRGObserver(["Sz", "Sx"], sites)

    E, psi = dmrg(H, psi0, sweeps; observer=observer, outputlevel=0)
    @test length(measurements(observer)["Sz"]) == 3
    @test length(measurements(observer)["Sx"]) == 3
    @test all(length.(measurements(observer)["Sz"]) .== N)
    @test all(length.(measurements(observer)["Sx"]) .== N)
    @test length(energies(observer)) == 3
    @test length(truncerrors(observer)) == 3
    @test energies(observer)[end] == E
    @test all(truncerrors(observer) .< 1E-9)

    orthogonalize!(psi, 1)
    m = scalar(dag(psi[1]) * noprime(op(sites, "Sz", 1) * psi[1]))
    @test measurements(observer)["Sz"][end][1] ≈ m
  end

  @testset "Sum of MPOs (ProjMPOSum)" begin
    N = 10
    sites = siteinds("S=1", N)

    osZ = OpSum()
    for j in 1:(N - 1)
      osZ += "Sz", j, "Sz", j + 1
    end
    HZ = MPO(osZ, sites)

    osXY = OpSum()
    for j in 1:(N - 1)
      osXY += 0.5, "S+", j, "S-", j + 1
      osXY += 0.5, "S-", j, "S+", j + 1
    end
    HXY = MPO(osXY, sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps, 10, 20, 40)
    mindim!(sweeps, 1, 10, 10)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps, 1E-10)
    energy, psi = dmrg([HZ, HXY], psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "Excited-state DMRG" begin
    N = 10
    weight = 15.0

    sites = siteinds("S=1", N)
    sites[1] = Index(2, "S=1/2,n=1,Site")
    sites[N] = Index(2, "S=1/2,n=$N,Site")

    os = OpSum()
    for j in 1:(N - 1)
      os += "Sz", j, "Sz", j + 1
      os += 0.5, "S+", j, "S-", j + 1
      os += 0.5, "S-", j, "S+", j + 1
    end
    H = MPO(os, sites)

    psi0i = randomMPS(sites; linkdims=10)

    sweeps = Sweeps(4)
    maxdim!(sweeps, 10, 20, 100, 100)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps, 1E-10)

    energy0, psi0 = dmrg(H, psi0i, sweeps; outputlevel=0)
    @test energy0 < -11.5

    psi1i = randomMPS(sites; linkdims=10)
    energy1, psi1 = dmrg(H, [psi0], psi1i, sweeps; outputlevel=0, weight=weight)

    @test energy1 > energy0
    @test energy1 < -11.1

    @test inner(psi1, psi0) < 1E-5
  end

  @testset "Fermionic Hamiltonian" begin
    N = 10
    t1 = 1.0
    t2 = 0.5
    V = 0.2
    s = siteinds("Fermion", N; conserve_qns=true)

    state = fill(1, N)
    state[1] = 2
    state[3] = 2
    state[5] = 2
    state[7] = 2
    psi0 = productMPS(s, state)

    os = OpSum()
    for j in 1:(N - 1)
      os += -t1, "Cdag", j, "C", j + 1
      os += -t1, "Cdag", j + 1, "C", j
      os += V, "N", j, "N", j + 1
    end
    for j in 1:(N - 2)
      os += -t2, "Cdag", j, "C", j + 2
      os += -t2, "Cdag", j + 2, "C", j
    end
    H = MPO(os, s)

    sweeps = Sweeps(5)
    maxdim!(sweeps, 10, 20, 100, 100, 200)
    cutoff!(sweeps, 1E-8)
    noise!(sweeps, 1E-10)

    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
    @test (-6.5 < energy < -6.4)
  end

  @testset "Hubbard model" begin
    N = 10
    Npart = 8
    t1 = 1.0
    U = 1.0
    V1 = 0.5
    sites = siteinds("Electron", N; conserve_qns=true)
    os = OpSum()
    for i in 1:N
      os += (U, "Nupdn", i)
    end
    for b in 1:(N - 1)
      os += -t1, "Cdagup", b, "Cup", b + 1
      os += -t1, "Cdagup", b + 1, "Cup", b
      os += -t1, "Cdagdn", b, "Cdn", b + 1
      os += -t1, "Cdagdn", b + 1, "Cdn", b
      os += V1, "Ntot", b, "Ntot", b + 1
    end
    H = MPO(os, sites)
    sweeps = Sweeps(6)
    maxdim!(sweeps, 50, 100, 200, 400, 800, 800)
    cutoff!(sweeps, 1E-10)
    state = ["Up", "Dn", "Dn", "Up", "Emp", "Up", "Up", "Emp", "Dn", "Dn"]
    psi0 = randomMPS(sites, state; linkdims=10)
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
    @test (-8.02 < energy < -8.01)
  end

  @testset "Input Without Ortho Center or Not at 1" begin
    N = 6
    sites = siteinds("S=1", N)

    os = OpSum()
    for j in 1:(N - 1)
      add!(os, "Sz", j, "Sz", j + 1)
      add!(os, 0.5, "S+", j, "S-", j + 1)
      add!(os, 0.5, "S-", j, "S+", j + 1)
    end
    H = MPO(os, sites)

    sweeps = Sweeps(1)
    maxdim!(sweeps, 10)
    cutoff!(sweeps, 1E-11)

    psi0 = randomMPS(sites; linkdims=4)

    # Test that input works with wrong ortho center:
    orthogonalize!(psi0, 5)
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)

    # Test that input works with no ortho center:
    for j in 1:N
      psi0[j] = randomITensor(inds(psi0[j]))
    end
    energy, psi = dmrg(H, psi0, sweeps; outputlevel=0)
  end
end

nothing


text:
include("swap.jl")
include("scal.jl")
include("copy.jl")
include("axpy.jl")
include("dot.jl")
include("dotu.jl")
include("dotc.jl")
include("nrm2.jl")
include("asum.jl")
include("sum.jl")
include("iamax.jl")
include("iamin.jl")
include("had.jl")


text:
# Quadmath conversions

@inline Quadmath.Float128(x::Double64) = Quadmath.Float128(x.hi) + x.lo
@inline Quadmath.Float128(x::Double32) = Quadmath.Float128(x.hi) + x.lo
@inline Quadmath.Float128(x::Double16) = Quadmath.Float128(x.hi) + Float32(x.lo)

@inline Quadmath.ComplexF128(x::ComplexDF64) = Quadmath.ComplexF128(Quadmath.Float128(real(x)), Quadmath.Float128(imag(x)))
@inline Quadmath.ComplexF128(x::ComplexDF32) = Quadmath.ComplexF128(Quadmath.Float128(real(x)), Quadmath.Float128(imag(x)))
@inline Quadmath.ComplexF128(x::ComplexDF16) = Quadmath.ComplexF128(Quadmath.Float128(real(x)), Quadmath.Float128(imag(x)))

@inline function Double64(x::Quadmath.Float128)
    hi = Float64(x)
    lo = Float64(x - hi)
    return Double64(hi,lo)
end
@inline function Double32(x::Quadmath.Float128)
    hi = Float32(x)
    lo = Float32(x - hi)
    return Double32(hi,lo)
end
@inline function Double16(x::Quadmath.Float128)
    hi = Float16(x)
    lo = Float16(x - hi)
    return Double16(hi,lo)
end

@inline ComplexDF64(x::Quadmath.ComplexF128) =
    ComplexDF64(Double64(real(x)), Double64(imag(x)))

@inline ComplexDF32(x::Quadmath.ComplexF128) =
    ComplexDF32(Double32(real(x)), Double32(imag(x)))

@inline ComplexDF16(x::Quadmath.ComplexF128) =
    ComplexDF16(Double16(real(x)), Double16(imag(x)))

@inline Double64Float128(fn::F, x::Double64) where {F<:Function} =
    Double64(fn(Quadmath.Float128(x)))

@inline Double32Float128(fn::F, x::Double32) where {F<:Function} =
    Double32(fn(Quadmath.Float128(x)))

@inline Double16Float128(fn::F, x::Double16) where {F<:Function} =
    Double16(fn(Quadmath.Float128(x)))

@inline ComplexDF64ComplexF128(fn::F, x::ComplexDF64) where {F<:Function} =
    ComplexDF64(fn(Quadmath.ComplexF128(x)))

@inline ComplexDF32ComplexF128(fn::F, x::ComplexDF32) where {F<:Function} =
    ComplexDF32(fn(Quadmath.ComplexF128(x)))

@inline ComplexDF16ComplexF128(fn::F, x::ComplexDF32) where {F<:Function} =
    ComplexDF16(fn(Quadmath.ComplexF128(x)))

convert(::Type{Quadmath.Float128}, x::Double64) = Quadmath.Float128(x)
convert(::Type{Quadmath.Float128}, x::Double32) = Quadmath.Float128(x)
convert(::Type{Quadmath.Float128}, x::Double16) = Quadmath.Float128(x)

convert(::Type{Double64}, x::Quadmath.Float128) = Double64(x)
convert(::Type{Double32}, x::Quadmath.Float128) = Double32(x)
convert(::Type{Double16}, x::Quadmath.Float128) = Double16(x)

promote_rule(::Type{Double64}, ::Type{Quadmath.Float128}) = Double64
promote_rule(::Type{Double32}, ::Type{Quadmath.Float128}) = Double32
promote_rule(::Type{Double16}, ::Type{Quadmath.Float128}) = Double16

# Complex

@inline Quadmath.ComplexF128(x::Double64) = Quadmath.ComplexF128(Quadmath.Float128(x))
@inline Quadmath.ComplexF128(x::Double32) = Quadmath.ComplexF128(Quadmath.Float128(x))
@inline Quadmath.ComplexF128(x::Double16) = Quadmath.ComplexF128(Quadmath.Float128(x))

@inline Double64(x::Quadmath.ComplexF128) = Double64(real(x))
@inline Double32(x::Quadmath.ComplexF128) = Double32(real(x))
@inline Double16(x::Quadmath.ComplexF128) = Double16(real(x))

convert(::Type{Quadmath.ComplexF128}, x::ComplexDF64) = Quadmath.ComplexF128(x)
convert(::Type{Quadmath.ComplexF128}, x::ComplexDF32) = Quadmath.ComplexF128(x)
convert(::Type{Quadmath.ComplexF128}, x::ComplexDF16) = Quadmath.ComplexF128(x)

convert(::Type{ComplexDF64}, x::Quadmath.ComplexF128) = ComplexDF64(x)
convert(::Type{ComplexDF32}, x::Quadmath.ComplexF128) = ComplexDF32(x)
convert(::Type{ComplexDF16}, x::Quadmath.ComplexF128) = ComplexDF16(x)

promote_rule(::Type{ComplexDF64}, ::Type{Quadmath.ComplexF128}) = ComplexDF64
promote_rule(::Type{ComplexDF32}, ::Type{Quadmath.ComplexF128}) = ComplexDF32
promote_rule(::Type{ComplexDF16}, ::Type{Quadmath.ComplexF128}) = ComplexDF16


text:
#Structure containg design matrix and related info for each submodel"
mutable struct UmDesign
  name::Symbol
  formula::FormulaTerm
  link::Link
  data::DataFrame
  coefnames::Array{String}
  mat::Array{Float64,2}
  idx::Union{UnitRange{Int64}, Nothing}
end

#Check formulas for validity
function check_formula(formula::FormulaTerm, data::DataFrame)
  vars = varnames(formula).rhs
  if isnothing(vars) return nothing end
  vars = String.(vars)
  covs = names(data)
  for x in vars
    if x ∉ covs error(string(x," not found in data")) end
  end
  return nothing
end

#Add dummy column for response variable to data
add_resp = function(x::DataFrame, f::FormulaTerm)
  lhs = f.lhs.sym
  if lhs in names(x) return x end
  dc = deepcopy(x)
  dc[!,lhs] = zeros(nrow(dc))
  return dc
end

#Workaround for predicting from newdata with missing cat variable levels

#Get number of unique levels for a categorical variable
function n_levels(x::AbstractArray)
  typeof(x)<:CategoricalArray ? length(levels(x)) : 0
end

#Repeat levels of cat variable to certain length
function rep_levels(x::AbstractArray, n::Int)
  if typeof(x)<:CategoricalArray
    reps = Int(ceil(n / length(levels(x))))
    out = categorical(repeat(levels(x), reps))
    levels!(out, levels(x))
    return out[1:n]
  end
  return repeat([0], n)
end

#Create data frame containing all categorical levels to bind to original data
function add_levels(df::DataFrame)
  nlev = [n_levels(col) for col = eachcol(df)]
  maxlev = max(nlev...)
  ndf = DataFrame([rep_levels(col, maxlev) for col = eachcol(df)], :auto)
  rename!(ndf, names(df))
  return ndf
end

#Build and apply schema
function get_schema(formula::FormulaTerm, data::DataFrame)
  check_formula(formula, data)
  dat = add_resp(data, formula)
  dat_aug = [dat; add_levels(dat)]
  sch = schema(formula, dat_aug)
  apply_schema(formula, sch, StatisticalModel)
end

#Outer constructor for UmDesign objects
function UmDesign(name::Symbol, formula::FormulaTerm, link::Link,
                  data::DataFrame)

  asch = get_schema(formula, data)
  cn = coefnames(asch.rhs)
  cn = cn isa String ? [cn] : cn
  mat = modelcols(asch.rhs, data)
  UmDesign(name, formula, link, data, cn, mat, nothing)
end

#Add indexes to map combined coefs back to submodels
function add_idx!(dm::Array{UmDesign})
  idx = 1
  for x in dm
    np = length(x.coefnames)
    x.idx = idx:(idx+np-1)
    idx += np
  end
  return nothing
  #return idx - 1
end

#Get total parameters
function get_np(dm::Array{UmDesign})
  sum(map(x -> length(x.coefnames), dm))
end

#Transform linear predictor back to response scale
function transform(dm::UmDesign, coefs::Array)
  lp = dm.mat * coefs[dm.idx]
  invlink(lp, dm.link)
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
using IterativeSolvers, KrylovKit, Arpack
# In this file, we regroud a way to provide eigen solvers

abstract type AbstractEigenSolver end
abstract type AbstractMFEigenSolver <: AbstractEigenSolver end
abstract type AbstractFloquetSolver <: AbstractEigenSolver end

# the following function returns the n-th eigenvectors computed by an eigen solver. This function is necessary given the different return types each eigensolver has
geteigenvector(eigsolve::ES, vecs, n::Int) where {ES <: AbstractEigenSolver} = vecs[:, n]
geteigenvector(eigsolve::ES, vecs, I::Array{Int64,1}) where {ES <: AbstractEigenSolver} = vecs[:, I]
####################################################################################################
# Solvers for default \ operator (backslash)
####################################################################################################
"""
The struct `Default` is used to  provide the backslash operator to our Package
"""
@with_kw struct DefaultEig{Tby} <: AbstractEigenSolver
	which::Tby = real		# how do we sort the computed eigenvalues
end

function (l::DefaultEig)(J, nev::Int64)
	# I put Array so we can call it on small sparse matrices
	F = eigen(Array(J))
	I = sortperm(F.values, by = l.which, rev = true)
	nev2 = min(nev, length(I))
	return Complex.(F.values[I[1:nev2]]), F.vectors[:, I[1:nev2]], true, 1
end

# case of sparse matrices or matrix free method
@with_kw struct EigArpack{T, Tby, Tw} <: AbstractEigenSolver
	sigma::T = nothing
	which::Symbol = :LR
	by::Tby = real			# how do we sort the computed eigenvalues.
	kwargs::Tw = nothing
end

EigArpack(sigma = nothing, which = :LR; kwargs...) = EigArpack(sigma, which, real, kwargs)

function (l::EigArpack)(J, nev::Int64)
	if J isa AbstractMatrix
		λ, ϕ, ncv = Arpack.eigs(J; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	else
		N = length(l.kwargs[:v0])
		T = eltype(l.kwargs[:v0])
		Jmap = LinearMap{T}(J, N, N; ismutating = false)
		λ, ϕ, ncv, = Arpack.eigs(Jmap; nev = nev, which = l.which, sigma = l.sigma, l.kwargs...)
	end
	I = sortperm(λ, by = l.by, rev = true)
	ncv < nev && @warn "$ncv eigenvalues have converged using Arpack.eigs, you requested $nev"
	return λ[I], ϕ[:, I], true, 1
end
####################################################################################################
# Solvers for KrylovKit
####################################################################################################
@with_kw struct EigKrylovKit{T, vectype} <: AbstractMFEigenSolver
	dim::Int64 = KrylovDefaults.krylovdim	# Krylov Dimension
	tol::T = 1e-4							# tolerance for solver
	restart::Int64 = 200					# number of restarts
	maxiter::Int64 = KrylovDefaults.maxiter
	verbose::Int = 0
	which::Symbol = :LR
	issymmetric::Bool = false				# if the linear map is symmetric, only meaningful if T<:Real
	ishermitian::Bool = false 				# if the linear map is hermitian
	x₀::vectype = nothing					# example of vector in case of a matrix-free operator
end

function (l::EigKrylovKit{T, vectype})(J, nev::Int64) where {T, vectype}
	if J isa AbstractMatrix && isnothing(l.x₀)
		vals, vec, info = KrylovKit.eigsolve(J, nev, l.which;  verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, tol = l.tol, issymmetric = l.issymmetric, ishermitian = l.ishermitian)
	else
		vals, vec, info = KrylovKit.eigsolve(J, l.x₀, nev, l.which;  verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, tol = l.tol, issymmetric = l.issymmetric, ishermitian = l.ishermitian)
	end
	info.converged == 0 && (@warn "KrylovKit.eigsolve solver did not converge")
	return vals, vec, true, info.numops
end

geteigenvector(eigsolve::EigKrylovKit{T, vectype}, vecs, n::Int) where {T, vectype} = vecs[n]
geteigenvector(eigsolve::EigKrylovKit{T, vectype}, vecs, I::Array{Int64,1}) where {T, vectype} = vecs[I]


text:
@testset "Parameters" begin
    # Swarm size
    n = 3
    # Test ranges
    r1 = range(String, :r1; values=["a", "b", "c"])
    r2 = range(Int, :r2; lower=1, upper=3, scale=exp10)
    r3 = range(Float64, :r3; lower=1, upper=Inf, origin=4, unit=1)
    r4 = range(Float32, :r4; lower=-Inf, upper=Inf, origin=0, unit=1)
    rs = (r1, r2, r3, r4)
    # Test distribution types
    Ds = (Uniform, Gamma, Normal)
    # Manually fitted distributions for test ranges
    d1 = Dirichlet(ones(3))
    d2 = Uniform(1, 3)
    d3 = truncated(Gamma(16, 0.25), 1, Inf)
    d4 = Normal(0, 1)
    ds = (d1, d2, d3, d4)
    # Test range representation lengths
    lengths = (3, 1, 1, 1)
    # Test range's corresponding indices in internal representation state.X
    indices = (1:3, 4, 5, 6)
    # Initial hyperparameter representation
    rng = StableRNG(1234)
    X1 = rand(rng, d1, n)'
    X2 = rand(rng, d2, n)
    X3 = rand(rng, d3, n)
    X4 = rand(rng, d4, n)
    Xs = (X1, X2, X3, X4)

    @testset "Initializer" begin
        for (r, D) in zip(rs[2:end], Ds)
            @test PSO._initializer(MLJTuning.boundedness(r)) === D
        end
    end

    @testset "Initialize with distribution types" begin
        rng = StableRNG(1234)
        PSO._initialize(rng, r1, n)
        for (r, D, l, X) in zip(rs[2:end], Ds, lengths[2:end], Xs[2:end])
            r̂, l̂, X̂ = PSO._initialize(rng, r, D, n)[[1,3,4]]
            @test r̂ === r
            @test l̂ == l
            @test X̂ ≈ X
        end
    end

    @testset "Initialize with distributions" begin
        rng = StableRNG(1234)
        for (r, d, l, X) in zip(rs, ds, lengths, Xs)
            r̂, l̂, X̂ = PSO._initialize(rng, r, d, n)[[1,3,4]]
            @test r̂ === r
            @test l̂ == l
            @test X̂ ≈ X
        end
    end

    @testset "Range Indices" begin
        @test PSO._to_indices(lengths) == indices
    end

    @testset "Unsupported distributions" begin
        rng = StableRNG(1234)
        @test_throws ArgumentError PSO._initialize(rng, r1, Uniform, n)
        @test_throws ArgumentError PSO._initialize(rng, r1, Dirichlet, n)
        @test_throws ArgumentError PSO._initialize(rng, r1, d2, n)
        @test_throws ArgumentError PSO._initialize(rng, r1, Dirichlet(ones(4)), n)
        @test_throws ArgumentError PSO._initialize(rng, r2, Dirichlet, n)
        @test_throws ArgumentError PSO._initialize(rng, r2, d1, n)
    end

    @testset "Initialize one range" begin
        rng = StableRNG(1234)
        for (r, l, i, X) in zip(rs, lengths, indices, Xs)
            state = PSO.initialize(rng, r, n)
            @test state.ranges == (r,)
            @test state.indices == (l == 1 ? 1 : 1:l,)
            @test state.X ≈ X
        end
    end

    @testset "Initialize multiple ranges" begin
        rng = StableRNG(1234)
        ranges = [r1, (r2, Uniform), (r3, d3), r4]
        state = PSO.initialize(rng, ranges, n)
        @test state.ranges == rs
        @test state.indices == indices
        @test state.X ≈ hcat(Xs...)
    end

    @testset "Retrieve parameters" begin
        rng = StableRNG(1234)
        ranges = [r1, (r2, Uniform), (r3, d3), r4]
        state = PSO.initialize(rng, ranges, n)
        PSO.retrieve!(rng, state)
        @test state.parameters == (
            ["a", "a", "c"],
            [553, 250, 375],
            [3.9372495283243105, 3.6569395920512977, 3.6354556967115146],
            [-0.8067647f0, 0.4209916f0, 0.6736019f0]
        )
    end
end


text:
"""
sice用のモジュール
"""
module SiceKinematics

# export x1, x2, x3, x4
# export J1, J2, J3, J4
# export J1_dot, J2_dot, J3_dot, J4_dot

export calc_x, calc_J, calc_dJ

const l1 = 1.0
const l2 = 1.0
const l3 = 1.0
const l4 = 1.0

function x1(q)
    [
        l1*cos(q[1])
        l1*sin(q[1])
    ]
end

function x2(q)
    [
        l1*cos(q[1]) + l2*cos(q[1] + q[2])
        l1*sin(q[1]) + l2*sin(q[1] + q[2])
    ]
end

function x3(q)
    [
        l1*cos(q[1]) + l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3])
        l1*sin(q[1]) + l2*sin(q[1] + q[2]) + l3*sin(q[1] + q[2] + q[3])
    ]
end

function x4(q)
    [
        l1*cos(q[1]) + l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3]) + l4*cos(q[1] + q[2] + q[3] + q[4])
        l1*sin(q[1]) + l2*sin(q[1] + q[2]) + l3*sin(q[1] + q[2] + q[3]) + l4*sin(q[1] + q[2] + q[3] + q[4])
    ]
end

function J1(q::Vector{T}) where T
    [
        -l1*sin(q[1]) 0 0 0
        l1*cos(q[1]) 0 0 0
    ]
end

function J2(q::Vector{T}) where T
    [
        -l1*sin(q[1]) - l2*sin(q[1] + q[2]) -l2*sin(q[1] + q[2]) 0 0
        l1*cos(q[1]) + l2*cos(q[1] + q[2]) l2*cos(q[1] + q[2]) 0 0
    ]
end

function J3(q::Vector{T}) where T
    [
        -l1*sin(q[1]) - l2*sin(q[1] + q[2]) - l3*sin(q[1] + q[2] + q[3]) -l2*sin(q[1] + q[2]) - l3*sin(q[1] + q[2] + q[3]) -l3*sin(q[1] + q[2] + q[3]) 0
        l1*cos(q[1]) + l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3]) l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3]) l3*cos(q[1] + q[2] + q[3]) 0
    ]
end

function J4(q::Vector{T}) where T
    [
        -l1*sin(q[1]) - l2*sin(q[1] + q[2]) - l3*sin(q[1] + q[2] + q[3]) - l4*sin(q[1] + q[2] + q[3] + q[4]) -l2*sin(q[1] + q[2]) - l3*sin(q[1] + q[2] + q[3]) - l4*sin(q[1] + q[2] + q[3] + q[4]) -l3*sin(q[1] + q[2] + q[3]) - l4*sin(q[1] + q[2] + q[3] + q[4]) -l4*sin(q[1] + q[2] + q[3] + q[4])
        l1*cos(q[1]) + l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3]) + l4*cos(q[1] + q[2] + q[3] + q[4]) l2*cos(q[1] + q[2]) + l3*cos(q[1] + q[2] + q[3]) + l4*cos(q[1] + q[2] + q[3] + q[4]) l3*cos(q[1] + q[2] + q[3]) + l4*cos(q[1] + q[2] + q[3] + q[4]) l4*cos(q[1] + q[2] + q[3] + q[4])
    ]
end

function J1_dot(q::Vector{T}, dq::Vector{T}) where T
    [
        -dq[1]*l1*cos(q[1]) 0 0 0
        -dq[1]*l1*sin(q[1]) 0 0 0
    ]
end

function J2_dot(q::Vector{T}, dq::Vector{T}) where T
    [
        -dq[1]*l1*cos(q[1]) - l2*(dq[1] + dq[2])*cos(q[1] + q[2]) -l2*(dq[1] + dq[2])*cos(q[1] + q[2]) 0 0
        -dq[1]*l1*sin(q[1]) - l2*(dq[1] + dq[2])*sin(q[1] + q[2]) -l2*(dq[1] + dq[2])*sin(q[1] + q[2]) 0 0
    ]
end

function J3_dot(q::Vector{T}, dq::Vector{T}) where T
    [
        -dq[1]*l1*cos(q[1]) - l2*(dq[1] + dq[2])*cos(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) -l2*(dq[1] + dq[2])*cos(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) -l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) 0
        -dq[1]*l1*sin(q[1]) - l2*(dq[1] + dq[2])*sin(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) -l2*(dq[1] + dq[2])*sin(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) -l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) 0
    ]
end

function J4_dot(q::Vector{T}, dq::Vector{T}) where T
    [
        -dq[1]*l1*cos(q[1]) - l2*(dq[1] + dq[2])*cos(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*cos(q[1] + q[2] + q[3] + q[4]) -l2*(dq[1] + dq[2])*cos(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*cos(q[1] + q[2] + q[3] + q[4]) -l3*(dq[1] + dq[2] + dq[3])*cos(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*cos(q[1] + q[2] + q[3] + q[4]) -l4*(dq[1] + dq[2] + dq[3] + dq[4])*cos(q[1] + q[2] + q[3] + q[4])
        -dq[1]*l1*sin(q[1]) - l2*(dq[1] + dq[2])*sin(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*sin(q[1] + q[2] + q[3] + q[4]) -l2*(dq[1] + dq[2])*sin(q[1] + q[2]) - l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*sin(q[1] + q[2] + q[3] + q[4]) -l3*(dq[1] + dq[2] + dq[3])*sin(q[1] + q[2] + q[3]) - l4*(dq[1] + dq[2] + dq[3] + dq[4])*sin(q[1] + q[2] + q[3] + q[4]) -l4*(dq[1] + dq[2] + dq[3] + dq[4])*sin(q[1] + q[2] + q[3] + q[4])
    ]
end

"""位置を全部計算"""
function calc_x(q::Vector{T}) where T
    x1(q), x2(q), x3(q), x4(q)
end

"""ヤコビ行列を全部計算"""
function calc_J(q::Vector{T}) where T
    return J1(q), J2(q), J3(q), J4(q)
end

"""ヤコビ行列の時間微分を全部計算"""
function calc_dJ(q::Vector{T}, dq::Vector{T}) where T
    J1_dot(q, dq), J2_dot(q, dq), J3_dot(q, dq), J4_dot(q, dq)
end

end


text:
using BenchmarkTools
s = read("../data/input1.txt", String)
arr = parse.(Int, split(s))

part1(arr) = count(i -> arr[i+1] > arr[i], 1:length(arr)-1)

@show part1(arr)
@btime part1($arr)

part2(arr) = count(i -> arr[i+3] > arr[i], 1:length(arr)-3)

@show part2(arr)
@btime part2($arr)

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
"""
    Mlm(B, varB, sigma, data, weights, targetType, lambda)

Type for storing the results of an mlm model fit.

"""
mutable struct Mlm

    # Coefficient estimates
    B::Array{Float64,2}
    # Coefficient variance estimates
    varB::Array{Float64,2}
    # Estimated variance of errors
    sigma::Array{Float64,2}

    # Response and predictor matrices
    data::RawData

    # Column weights for `Y`, or `nothing`
    weights
    # String indicating target type to shrink the variance toward, or `nothing`
    targetType
    # Estimated shrinkage coefficient
    lambda::Float64
end

"""
    mlm_fit(data, weights, targetType)

Matrix linear model using least squares method. Optionally incorporates
shrinkage of the variance of the errors.

# Arguments

- data = RawData object
- weights = `nothing`
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D".
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries

# Value

An Mlm object

"""
function mlm_fit(data::RawData, weights::Nothing, targetType)

    # Calculate and store transpose(X)*X
    XTX = transpose(get_X(data))*get_X(data)
    # Calculate and store transpose(Z)*Z
    ZTZ = transpose(get_Z(data))*get_Z(data)

    # Estimate MLM coefficients
    B = calc_coeffs(get_X(data), get_Y(data), get_Z(data), XTX, ZTZ)

    # Calculate residuals
    resid = calc_resid(get_X(data), get_Y(data), get_Z(data), B)

    # Estimate variance of errors, optionally with variance shrinkage
    sigma, lambda = calc_sigma(resid, targetType)

    # Estimate variance of coefficient estimates
    varB = calc_var(get_X(data), get_Z(data), XTX, ZTZ, sigma)

    # Return Mlm object with estimates for coefficients, variances, and sigma
    return Mlm(B, varB, sigma, data, weights, targetType, lambda)
end

"""
    mlm_fit(data, weights, targetType)

Matrix linear model using column weighted least squares method. Optionally
incorporates shrinkage of the variance of the errors.

# Arguments

- data = RawData object
- weights = 1d array of floats to use as column weights for `Y`. Must be the
  same length as the number of columns of `Y`.
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D".
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries

# Value

An Mlm object

"""
function mlm_fit(data::RawData, weights::Array{Float64,1}, targetType)

    # Throw an error if the length of weights is incorrect
    if (length(weights) != size(get_Y(data),2))
        error("Weights must be same length as number of columns in Y.")
    end

    # Diagonal matrix with weights on diagonal.
    W = diagm(weights)

    # Calculate and store WZ product
    WZ = W*get_Z(data)

    # Calculate and store transpose(X)*X
    XTX = transpose(get_X(data))*get_X(data)
    # Calculate and store transpose(Z)*W*Z
    ZTWZ = transpose(get_Z(data))*WZ

    # Estimate MLM coefficients
    B = calc_coeffs(get_X(data), get_Y(data), WZ, XTX, ZTWZ)

    # Calculate residuals
    resid = calc_resid(get_X(data), get_Y(data), get_Z(data), B)

    # Estimate variance of errors, optionally with variance shrinkage
    sigma, lambda = calc_sigma(resid, targetType)

    # Estimate variance of coefficient estimates
    varB = calc_var(get_X(data), WZ, XTX, ZTWZ, sigma)

    # Return Mlm object with estimates for coefficients, variances, and sigma
    return Mlm(B, varB, sigma, data, weights, targetType, lambda)
end

"""
    mlm(data; hasXIntercept, hasZIntercept, weights, targetType)

Matrix linear model using least squares method. Column weighted least squares
and shrinkage of the variance of the errors are options.

# Arguments

- data = RawData object

# Keyword arguments

- hasXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- hasZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- weights = 1d array of floats to use as column weights for `Y`, or `nothing`.
  If the former, must be the same length as the number of columns of `Y`.
  Defaults to `nothing`.
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D". Defaults to `nothing`.
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries

# Value

An Mlm object

"""
function mlm(data::RawData; hasXIntercept::Bool=true, hasZIntercept::Bool=true,
             weights=nothing, targetType=nothing, kwargs...)

    #= Deprecation:check for previous version keyword arguments
    ---------------------------------------------------------------------------------------------=#
    if haskey(kwargs, :isXIntercept)
        @warn "Keyword arguments `isXIntercept` and `isZIntercept` are deprecated, use
                    `hasXIntercept` and `hasZIntercept` instead."

        hasXIntercept = values(kwargs).isXIntercept
    end

    if haskey(kwargs, :isZIntercept)
        @warn "Keyword arguments `isXIntercept` and `isZIntercept` are deprecated, use
                    `hasXIntercept` and `hasZIntercept` instead."

        hasZIntercept = values(kwargs).isZIntercept
    end
    #-----------------------------------------------------------------------------------------------

    # Add X and Z intercepts if necessary
    if hasXIntercept==true && data.predictors.hasXIntercept==false
        data.predictors.X = add_intercept(data.predictors.X)
        data.predictors.hasXIntercept = true
        data.p = data.p + 1
    end
    if hasZIntercept==true && data.predictors.hasZIntercept==false
        data.predictors.Z = add_intercept(data.predictors.Z)
        data.predictors.hasZIntercept = true
        data.q = data.q + 1
    end

    # Remove X and Z intercepts in new predictors if necessary
    if hasXIntercept==false && data.predictors.hasXIntercept==true
        data.predictors.X = remove_intercept(data.predictors.X)
        data.predictors.hasXIntercept = false
        data.p = data.p - 1
    end
    if hasZIntercept==false && data.predictors.hasZIntercept==true
        data.predictors.Z = remove_intercept(data.predictors.Z)
        data.predictors.hasZIntercept = false
        data.q = data.q - 1
    end

    if (typeof(targetType) != Nothing) & !(targetType in ["A", "B", "C", "D"])
        println("Unrecognizable targetType will be ignored and no variance shrinkage will be performed.")
        targetType = nothing
    end

    # Run matrix linear models
    return mlm_fit(data, weights, targetType)
end

"""
    t_stat(MLM, isMainEff)

Calculates t-statistics of an Mlm object

# Arguments

- MLM = Mlm object
- isMainEff = boolean flag indicating whether or not to include t-statistics
  for the main effects

# Value

2d array of floats

"""
function t_stat(MLM::Mlm, isMainEff::Bool=false)

    # Cases when not including main effects
    if isMainEff== false
        # If X and Z intercepts were both included
        if (MLM.data.predictors.hasXIntercept==true) &&
           (MLM.data.predictors.hasZIntercept==true)
            return MLM.B[2:end, 2:end]./sqrt.(MLM.varB[2:end, 2:end])

        # If only X intercept was included
        elseif MLM.data.predictors.hasXIntercept==true
            return MLM.B[2:end, :]./sqrt.(MLM.varB[2:end, :])

        # If only Z intercept was included
        elseif  MLM.data.predictors.hasZIntercept==true
            return MLM.B[:, 2:end]./sqrt.(MLM.varB[:, 2:end])

        end
    end

    # Case when including main effects
    return MLM.B ./ sqrt.(MLM.varB)
end


text:
# Warp Shuffle (B.14)

# TODO: does not work on sub-word (ie. Int16) or non-word divisible sized types

# TODO: these functions should dispatch based on the actual warp size
const ws = Int32(32)

# TODO: this functionality should throw <sm_30

# primitive intrinsics

# "two packed values specifying a mask for logically splitting warps into sub-segments
# and an upper bound for clamping the source lane index"
@inline pack(width, mask) = (convert(UInt32, ws - width) << 8) | convert(UInt32, mask)

# NOTE: CUDA C disagrees with PTX on how shuffles are called
for (name, mode, mask, offset) in (("_up",   :up,   UInt32(0x00), src->src),
                                   ("_down", :down, UInt32(0x1f), src->src),
                                   ("_xor",  :bfly, UInt32(0x1f), src->src),
                                   ("",      :idx,  UInt32(0x1f), src->:($src-1)))
    fname = Symbol("shfl$name")
    @eval export $fname

    # newer hardware/CUDA versions use synchronizing intrinsics, which take an extra
    # mask argument indicating which threads in the lane should be synchronized
    let intrinsic = "llvm.nvvm.shfl.sync.$mode.i32"
        fname_sync = Symbol("$(fname)_sync")
        __fname_sync = Symbol("__$(fname)_sync")
        @eval begin
            export $fname_sync

            # HACK: recurse_value_invocation and friends split the first argument of a call,
            #       so swap mask and val for these tools to works.
            @inline $fname_sync(mask, val, src, width=$ws) =
                $__fname_sync(val, mask, src, width)
            @inline $__fname_sync(val::UInt32, mask, src, width) =
                ccall($intrinsic, llvmcall, UInt32,
                      (UInt32, UInt32, UInt32, UInt32),
                      mask, val, $(offset(:src)), pack(width, $mask))
        end
    end

    # FIXME: if we only support CUDA 9+, should we not always use the sync version?
    #@inline $fname(val::UInt32, src, width=$ws, mask::UInt32=0xffffffff) =
    #    $fname_sync(mask, val, src, width)
    let intrinsic = "llvm.nvvm.shfl.$mode.i32"
        @eval begin
            @inline $fname(val::UInt32, src, width=$ws) =
                ccall($intrinsic, llvmcall, UInt32,
                      (UInt32, UInt32, UInt32),
                      val, $(offset(:src)), pack(width, $mask))
        end
    end
end

# wide and aggregate intrinsics

for name in ["_up", "_down", "_xor", ""]
    fname = Symbol("shfl$name")
    @eval @inline $fname(src, args...) = recurse_value_invocation($fname, src, args...)

    fname_sync = Symbol("__$(fname)_sync")
    @eval @inline $fname_sync(src, args...) = recurse_value_invocation($fname_sync, src, args...)
end

# documentation

@doc """
    shfl(val, lane::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a directly indexed lane `lane`. The argument `threadmask` for selecting
which threads to synchronize is only available on recent hardware, and defaults to all
threads in the warp.
""" shfl

@doc """
    shfl_up(val, delta::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane with lower ID relative to caller. The argument `threadmask` for
selecting which threads to synchronize is only available on recent hardware, and defaults to
all threads in the warp.
""" shfl_up

@doc """
    shfl_down(val, delta::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane with higher ID relative to caller. The argument `threadmask` for
selecting which threads to synchronize is only available on recent hardware, and defaults to
all threads in the warp.
""" shfl_down

@doc """
    shfl_xor(val, lanemask::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane based on bitwise XOR of own lane ID with `lanemask`. The
argument `threadmask` for selecting which threads to synchronize is only available on recent
hardware, and defaults to all threads in the warp.
""" shfl_xor

@doc """
    shfl_sync(threadmask::UInt32, val, lane::Integer, width::Integer=32)

Shuffle a value from a directly indexed lane `lane`, and synchronize threads according to
`threadmask`.
""" shfl_sync

@doc """
    shfl_up_sync(threadmask::UInt32, val, delta::Integer, width::Integer=32)

Shuffle a value from a lane with lower ID relative to caller, and synchronize threads
according to `threadmask`.
""" shfl_up_sync

@doc """
    shfl_down_sync(threadmask::UInt32, val, delta::Integer, width::Integer=32)

Shuffle a value from a lane with higher ID relative to caller, and synchronize threads
according to `threadmask`.
""" shfl_down_sync

@doc """
    shfl_xor_sync(threadmask::UInt32, val, mask::Integer, width::Integer=32)

Shuffle a value from a lane based on bitwise XOR of own lane ID with `mask`, and synchronize
threads according to `threadmask`.
""" shfl_xor_sync


text:
module ErrorFlow

export @catch
export @chain
export Error

struct Error
    value
    hints
    stacktrace
end

function Base.showerror(io::IO, err::Error)
    println(io, "Error:")
    Base.showerror(io, err.value)
    println(io)
    println(io)
    println(io, "Hints:")
    for (i, hint) in pairs(err.hints)
        print(io, i, ") ")
        println(io, hint)
    end
    Base.show_backtrace(
        io, err.stacktrace
    )
end

function append_hints(err::Error, hints)
    Error(err.value, [err.hints; hints], err.stacktrace)
end

macro chain(code, err)
    chain_macro(code, err)
end

macro var"catch"(code, hint=nothing)
    code = esc(code)
    if (hint === nothing)
        hints = []
    else
        hints = esc(:([$hint]))
    end
    quote
        try
            $code
        catch err
            st = $stacktrace($catch_backtrace())
            $Error(err, $hints, st)
        end
    end
end

function chain_macro(code, hint=nothing)
    code = esc(code)
    if (hint === nothing)
        hints = []
    else
        hints = esc(:([$hint]))
    end
    quote
        result = try
            $code
        catch err
            st = $stacktrace($catch_backtrace())
            ce = $Error(err, $hints, st)
            return ce
        end
        if result isa $Error
            return $append_hints(result, $hints)
        elseif result isa Exception
            st = $stacktrace($catch_backtrace())
            return $Error(result, $hints, st)
        else
            result
        end
    end
end

end


text:
module Funds
include("math.jl")
include("data.jl")
end