text:
# Conditionals

for op ∈ (:(==), :(!=), :(<), :(<=), :(>), :(>=) )
    rename = VECTOR_SYMBOLS[op]
    # rename_mask = Symbol(rename, :_mask)
    @eval begin
        # scalar versions handled in floating_point_arithmetic.jl
        # @inline $rename(s1::ScalarTypes, s2::ScalarTypes) = $op(s1,s2)
        @vectordef $rename function Base.$op(v1, v2) where {W,T}
            llvmwrap_bitmask(Val{$(QuoteNode(op))}(), extract_data(v1), extract_data(v2))
        end
        @vectordef $rename function Base.$op(v1, v2::T) where {W,T}
            $rename(extract_data(v1), v2)
        end
        @vectordef $rename function Base.$op(v1::T, v2) where {W,T}
            $rename(extract_data(v1), v2)
        end
    end
end
@inline visfinite(s::ScalarTypes) = isfinite(s)
@inline function visfinite(v1::Vec{W,T}) where {W,T<:FloatingTypes}
    U = uint_type(T)
    em = vbroadcast(Vec{W,U}, exponent_mask(T))
    iv = vreinterpret(Vec{W,U}, v1)
    vnot_equal(vand(iv, em), em)
end
@inline visfinite(v1::SVec) = SVec(visfinite(extract_data(v1)))
@inline Base.isfinite(v1::SVec) = SVec(visfinite(extract_data(v1)))

@inline visinf(s1::ScalarTypes) = isinf(s1)
@inline visinf(v1::Vec{W,T}) where {W,T<:FloatingTypes} = visequal(vabs(v1), vbroadcast(Vec{W,T},Inf))
@inline visinf(v1::SVec) = SVec(visinf(extract_data(v1)))
@inline Base.isinf(v1::SVec) = SVec(visinf(extract_data(v1)))

@inline visnan(s1::ScalarTypes) = isnan(s1)
@inline visnan(v1::Vec{W,T}) where {W,T<:FloatingTypes} = vnot_equal(v1, v1)
@inline visnan(v1::SVec) = SVec(visnan(extract_data(v1)))
@inline Base.isnan(v1::SVec) = SVec(visnan(extract_data(v1)))

@inline vissubnormal(s1::ScalarTypes) = issubnormal(s1)
@inline function vissubnormal(v1::Vec{W,T}) where {W,T<:FloatingTypes}
    U = uint_type(T)
    em = vbroadcast(Vec{W,U}, exponent_mask(T))
    sm = vbroadcast(Vec{W,U}, significand_mask(T))
    iv = vreinterpret(Vec{W,U}, v1)
    vand(visequal(vand(iv, em), vbroadcast(Vec{W,U}, 0)), vnot_equal(vand(iv, sm), vbroadcast(Vec{W,U}, 0)))
end
@inline vissubnormal(v1::SVec) = SVec(vissubnormal(extract_data(v1)))
@inline Base.issubnormal(v1::SVec) = SVec(vissubnormal(extract_data(v1)))

@inline signbit(s1::ScalarTypes) = signbit(s1)
@inline function vsignbit(v1::Vec{W,T}) where {W,T<:FloatingTypes}
    U = uint_type(T)
    sm = vbroadcast(Vec{W,U}, sign_mask(T))
    iv = vreinterpret(Vec{W,U}, v1)
    vnot_equal(vand(iv, sm), vbroadcast(Vec{W,U}, 0))
end
@inline function vsignbit(v1::SVec{W,T}) where {W,T<:FloatingTypes}
    SVec(vsignbit(extract_data(v1)))
end
@inline function Base.signbit(v1::SVec{W,T}) where {W,T<:FloatingTypes}
    SVec(vsignbit(extract_data(v1)))
end

@inline vifelse(c::Bool, x, y) = c ? x : y
# @inline vifelse(c::Bool, x, y) = ifelse(c, x, y)
@generated function vifelse(v1::Vec{W,Bool}, v2::Vec{W,T}, v3::Vec{W,T}) where {W,T}
    btyp = llvmtype(Bool)
    vbtyp = "<$W x $btyp>"
    abtyp = "[$W x $btyp]"
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    atyp = "[$W x $typ]"
    decls = String[]
    instrs = String[]
    push!(instrs, "%cond = trunc $vbtyp %0 to <$W x i1>")
    push!(instrs, "%res = select <$W x i1> %cond, $vtyp %1, $vtyp %2")
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{W,T},
            Tuple{Vec{W,Bool}, Vec{W,T}, Vec{W,T}},
            v1, v2, v3)
    end
end
@inline function vifelse(v1::AbstractSIMDVector{W,Bool},
        v2::AbstractSIMDVector{W,T}, v3::AbstractSIMDVector{W,T}) where {W,T}
    SVec(vifelse(extract_data(v1), extract_data(v2), extract_data(v3)))
end
# @inline function vifelse(v1::VecOrProd{W,Bool},
#         v2::VecOrProd{W,T}, v3::VecOrProd{W,T}) where {W,T}
#     vifelse(extract_data(v1), extract_data(v2), extract_data(v3))
# end
# @inline function Base.ifelse(v1::SVec{W,Bool},
#         v2::SVec{W,T}, v3::SVec{W,T}) where {W,T}
#     SVec(vifelse(extract_data(v1), extract_data(v2), extract_data(v3)))
# end

@generated function vifelse(mask::U, v2::Vec{W,T}, v3::Vec{W,T}) where {W,T,U<:Unsigned}
    @assert 8sizeof(U) >= W
    btyp = llvmtype(Bool)
    vbtyp = "<$W x $btyp>"
    abtyp = "[$W x $btyp]"
    typ = llvmtype(T)
    vtyp = "<$W x $typ>"
    atyp = "[$W x $typ]"
    decls = String[]
    instrs = String[]

    mtyp_input = llvmtype(U)
    mtyp_trunc = "i$W"

    if mtyp_input == mtyp_trunc
        push!(instrs, "%cond = bitcast $mtyp_input %0 to <$W x i1>")
    else
        push!(instrs, "%condtrunc = trunc $mtyp_input %0 to $mtyp_trunc")
        push!(instrs, "%cond = bitcast $mtyp_trunc %condtrunc to <$W x i1>")
    end
    # push!(instrs, "%cond = trunc $vbtyp %0 to <$W x i1>")
    if T <: FloatingTypes && Base.libllvm_version >= v"9"
        push!(instrs, "%res = select fast <$W x i1> %cond, $vtyp %1, $vtyp %2")
    else
        push!(instrs, "%res = select <$W x i1> %cond, $vtyp %1, $vtyp %2")
    end
    push!(instrs, "ret $vtyp %res")
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$W,$T},
            Tuple{$U, Vec{$W,$T}, Vec{$W,$T}},
            mask, v2, v3)
    end
end

@vpromote vifelse 3

@inline vifelse(U::Unsigned, v2::AbstractSIMDVector, v3::AbstractSIMDVector) = SVec(vifelse(U, extract_data(v2), extract_data(v3)))
@inline vifelse(U::Unsigned, v2::Vec{W,T}, s::Union{T,Int}) where {W,T} = vifelse(U, v2, vbroadcast(Vec{W,T}, s))
@inline vifelse(U::Unsigned, v2::AbstractSIMDVector{W,T}, s::Union{T,Int}) where {W,T} = SVec(vifelse(U, extract_data(v2), vbroadcast(Vec{W,T}, s)))
@inline vifelse(U::Unsigned, s::Union{T,Int}, v2::Vec{W,T}) where {W,T} = vifelse(U, vbroadcast(Vec{W,T}, s), v2)
@inline vifelse(U::Unsigned, s::Union{T,Int}, v2::AbstractSIMDVector{W,T}) where {W,T} = SVec(vifelse(U, vbroadcast(Vec{W,T}, s), extract_data(v2)))

@inline function vifelse(f::F, m::Mask{W}, vargs::Vararg{<:Any,N}) where {F<:Function,W,N}
    vifelse(m, f(vargs...), @inbounds(vargs[N]))
end

@vectordef visodd function Base.isodd(v) where {W,T<:Integer}
    visequal(vand(v, one(T)), one(T))
end

# @inline vand(U::Unsigned, vb::Vec{W,Bool}) where {W} = U & vecbool_to_unsigned(vb)
# @inline vand(vb::Vec{W,Bool}, U::Unsigned) where {W} = U & vecbool_to_unsigned(vb)


text:
"""
    NLevelSpace <: HilbertSpace
    NLevelSpace(name::Symbol,levels,GS=levels[1])

Define a [`HilbertSpace`](@ref) for an object consisting of `N` discrete energy
levels. The given `levels` must be an integer specifying the number of levels,
or an iterable collection of levels. The argument `GS` specifies which state
should be treated as ground state and is rewritten using population conservation
during simplification.
See also: [`Transition`](@ref)

Examples:
========
```
julia> ha = NLevelSpace(:a,3)
ℋ(a)

julia> ha = NLevelSpace(:a,(:g,:e))
ℋ(a)
```
"""
struct NLevelSpace{S,L,G} <: ConcreteHilbertSpace
    name::S
    levels::L
    GS::G
    function NLevelSpace{S,L,G}(name::S, levels::L, GS::G) where {S,L,G}
        (GS ∈ levels) || throw(ArgumentError("Ground state $GS not in levels $levels"))
        return new(name, levels, GS)
    end
end
NLevelSpace(name::S, levels::L, GS::G) where {S,L,G} = NLevelSpace{S,L,G}(name, levels, GS)
NLevelSpace(name,N::Int,GS) = NLevelSpace(name,1:N,GS)
NLevelSpace(name,N::Int) = NLevelSpace(name,1:N,1)
NLevelSpace(name,levels) = NLevelSpace(name,levels,levels[1])
Base.:(==)(h1::T,h2::T) where T<:NLevelSpace = (h1.name==h2.name && h1.levels==h2.levels && h1.GS==h2.GS)
Base.hash(n::NLevelSpace, h::UInt) = hash(n.name, hash(n.levels, hash(n.GS, h)))

levels(h::NLevelSpace) = h.levels
levels(h::NLevelSpace,aon) = levels(h)
levels(h::ProductSpace,aon) = levels(h.spaces[aon])
ground_state(h::NLevelSpace) = h.GS
ground_state(h::NLevelSpace,aon) = h.GS
ground_state(h::ProductSpace,aon) = ground_state(h.spaces[aon])

"""
    Transition <: QSym
    Transition(h::NLevelSpace,name::Symbol,i,j)

Fundamental operator defining a transition from level `j` to level `i` on a
[`NLevelSpace`](@ref). The notation corresponds to Dirac notation, i.e. the
above is equivalent to `|i⟩⟨j|`.

Examples
=======
```
julia> ha = NLevelSpace(:a,(:g,:e))
ℋ(a)

julia> σ = Transition(ha,:σ,:g,:e)
σge
```
"""
struct Transition{H,S,I,A,M} <: QSym
    hilbert::H
    name::S
    i::I
    j::I
    aon::A
    metadata::M
    function Transition{H,S,I,A,M}(hilbert::H, name::S, i::I, j::I, aon::A, metadata::M) where {H,S,I,A,M}
        @assert has_hilbert(NLevelSpace,hilbert,aon)
        @assert i∈levels(hilbert,aon) && j∈levels(hilbert,aon)
        new(hilbert, name, i, j, aon, metadata)
    end
end
function Transition(hilbert::H, name::S, i::I, j::I, aon::A; metadata::M=NO_METADATA) where {H,S,I,A,M}
    gs = ground_state(hilbert, aon)
    if isequal(i,j) && isequal(i,gs)
        args = Any[1]
        for k∈levels(hilbert, aon)
            isequal(k,gs) && continue
            push!(args, QMul(-1, [Transition{H,S,I,A,M}(hilbert, name, k, k, aon, metadata)]))
        end
        return QAdd(args)
    else
        return Transition{H,S,I,A,M}(hilbert, name, i, j, aon, metadata)
    end
end
Transition(hilbert::NLevelSpace,name,i,j; metadata=NO_METADATA) = Transition(hilbert,name,i,j,1; metadata)
function Transition(hilbert::ProductSpace,name,i,j; metadata=NO_METADATA)
    inds = findall(x->isa(x,NLevelSpace) || isa(x,ClusterSpace{<:NLevelSpace}),hilbert.spaces)
    if length(inds)==1
        return Transition(hilbert,name,i,j,inds[1]; metadata)
    else
        isempty(inds) && error("Can only create Transition on NLevelSpace! Not included in $(hilbert)")
        length(inds)>1 && error("More than one NLevelSpace in $(hilbert)! Specify on which Hilbert space Transition should be created with Transition(hilbert,name,i,j,acts_on)!")
    end
end
function Transition(hilbert::H,name::S,i::I,j::I,aon::A; metadata::M=NO_METADATA) where {H<:ProductSpace,S,I,A<:Int,M}
    if hilbert.spaces[aon] isa ClusterSpace
        hilbert.spaces[get_i(aon)].op_name[] = name # write original name of operator into hilbertspace
        op = Transition(hilbert.spaces[aon].original_space,name,i,j,1; metadata)
        return _cluster(hilbert, op, aon) # TODO maybe need metadata?
    else
        gs = ground_state(hilbert, aon)
        if isequal(i,j) && isequal(i,gs)
            args = Any[1]
            for k∈levels(hilbert, aon)
                isequal(k,gs) && continue
                push!(args, QMul(-1, [Transition{H,S,I,A,M}(hilbert,name,k,k,aon,metadata)]))
            end
            return QAdd(args)
        else
            return Transition{H,S,I,A,M}(hilbert,name,i,j,aon,metadata)
        end
    end
end
function Transition(h::ClusterSpace{<:NLevelSpace}, name, i, j, aon::Int=1; metadata=NO_METADATA)
    hilbert.spaces[get_i(aon)].op_name[] = name
    op = Transition(hilbert.original_space,name,i,j,aon; metadata)
    return _cluster(h,op,aon) # TODO maybe need metadata?
end

levels(t::Transition,args...) = levels(t.hilbert,args...)
ground_state(t::Transition,args...) = ground_state(t.hilbert,args...)

Base.adjoint(t::Transition) = Transition(t.hilbert,t.name,t.j,t.i,acts_on(t); t.metadata)
Base.isequal(t1::Transition,t2::Transition) = isequal(t1.hilbert, t2.hilbert) && isequal(t1.name,t2.name) && isequal(t1.i,t2.i) && isequal(t1.j,t2.j) && isequal(t1.aon,t2.aon)
Base.hash(t::Transition, h::UInt) = hash(t.hilbert, hash(t.name, hash(t.i, hash(t.j, hash(t.aon, h)))))

"""
    CallableTransition

A [`Transition`](@ref) where no levels have been specified. This type is callable
to allow for easy construction of concrete [`Transition`](@ref) instances.

Examples
========
```
julia> h = NLevelSpace(:atom, (:g,:e))
ℋ(atom)

julia> σ = Transition(h,:σ)
σ

julia> σ(:g,:e)
σge
```
"""
struct CallableTransition{H,S,A,M}
    hilbert::H
    name::S
    aon::A
    metadata::M
end

acts_on(c::CallableTransition) = c.aon

function (c::CallableTransition)(i, j)
    Transition(c.hilbert, c.name, i, j, c.aon; c.metadata)
end

Transition(hilbert, name, aon; metadata=NO_METADATA) = CallableTransition(hilbert, name, aon, metadata)
Transition(hilbert::NLevelSpace, name; metadata=NO_METADATA) = CallableTransition(hilbert, name, 1, metadata)
function Transition(hilbert::ProductSpace,name; metadata=NO_METADATA)
    inds = findall(x->isa(x,NLevelSpace),hilbert.spaces)
    if length(inds)==1
        return CallableTransition(hilbert,name,inds[1],metadata)
    else
        isempty(inds) && error("Can only create Transition on NLevelSpace! Not included in $(hilbert)")
        length(inds)>1 && error("More than one NLevelSpace in $(hilbert)! Specify on which Hilbert space Transition should be created with Transition(hilbert,name,i,j,acts_on)!")
    end
end

# Simplification
function *(a::Transition,b::Transition)
    check_hilbert(a, b)
    aon_a = acts_on(a)
    aon_b = acts_on(b)
    if aon_a == aon_b
        if isequal(a.j, b.i)
            return Transition(a.hilbert, a.name, a.i, b.j, a.aon)
        else
            return 0
        end
    elseif aon_a < aon_b
        return QMul(1, [a,b])
    else
        return QMul(1, [b,a])
    end
end
ismergeable(::Transition,::Transition) = true

ext: jl


text:
include("_dcd_on.jl")
include("_load_lp_1.jl")

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using MCMCBenchmarks, Test, Random

@testset "Gaussian Tests " begin
    Random.seed!(540457)
    path = pathof(MCMCBenchmarks)
    include("Gaussian_Models.jl")
    mu = 0
    sigma = 1
    Nd = 10^3
    Nreps = 1
    ProjDir = @__DIR__
    cd(ProjDir)
    samplers = (
      CmdStanNUTS(CmdStanConfig, ProjDir),
      AHMCNUTS(AHMCGaussian, AHMCconfig),
      DHMCNUTS(sampleDHMC))
    options = (Nsamples=2000,Nadapt=1000,delta=.8,Nd=Nd)
    results = benchmark(samplers, simulateGaussian, Nreps; options...)
    @test results[!,:mu_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ mu atol = .05
    @test results[!,:mu_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ mu atol = .05
    @test results[!,:mu_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ mu atol = .05
    @test results[!,:sigma_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ sigma atol = .05
    @test results[!,:sigma_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ sigma atol = .05
    @test results[!,:sigma_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ sigma atol = .05
    @test results[!,:mu_sampler_rhat][1] ≈ 1 atol = .03
    @test results[!,:sigma_sampler_rhat][1] ≈ 1 atol = .03
    isdir("tmp") && rm("tmp", recursive=true)
end


text:
function ECOM2(arg0::jint, arg1::jint, arg2::jdouble, arg3::ExtendedPVCoordinatesProvider, arg4::jdouble)
    return ECOM2((jint, jint, jdouble, ExtendedPVCoordinatesProvider, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function acceleration(obj::ECOM2, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::ECOM2, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function get_parameters_drivers(obj::ECOM2)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end


text:
export i_softmax_crossentropy, i_relu, i_logsumexp

function (_::PlusEq{typeof(argmax)})(out!, x::AbstractArray)
    out! += argmax(x)
    out!, x
end

function (_::MinusEq{typeof(argmax)})(out!, x::AbstractArray)
    out! -= argmax(x)
    out!, x
end

"""
    i_softmax_crossentropy(x, p, imax, xmax, Z, out)

Softmax-Cross entropy function.
"""
@i function i_softmax_crossentropy(x, p, imax, xmax, Z, out::T) where T
    # subtract maximum
    imax += argmax(x)  # trade off space of xmax to time
    xmax += x[imax]
    # accumulate exp(x) to Z, and finally get logZ
    for i=1:length(x)
        x[i] -= xmax
        Z += Base.exp(x[i])
    end
    @routine begin
        yi ← zero(T)
        logZ ← zero(T)
        logZ += log(Z)
    end
    for i=1:length(x)
        yi += logZ
        yi -= x[i]
        out += yi * p[i]
        yi += x[i]
        yi -= logZ
    end
    ~@routine
end

"""
    i_relu(out!, x)

ReLU in machine learning.
"""
@i function i_relu(out!, x)
    @invcheckoff if (x > 0, ~)
        out! += x
    end
end

"""
    i_logsumexp(logout!, out!, xs!, inds!, x)

Compute `logout! = log(sum(exp(x)))`.

# Arguments

    * `out!`, output,
    * `logout!`, logged output,
    * `xs!`, an empty vector to cache the ascending values (same type as `x`),
    * `inds!`, an empty vector to cache the ascending indices (integer type),
    * `x`, input vector.
"""
@i function i_logsumexp(logout!, out!, xs!, inds!, x::AbstractArray{T}) where T
  	i_ascending!(xs!, inds!, x)
    @routine begin
        mx ← zero(T)
        mx += xs![end]
    end
	@invcheckoff @inbounds for i=1:length(x)
		x[i] -= mx
		out! += exp(x[i])
		x[i] += mx
	end
  	logout! += log(out!)
	logout! += mx
    ~@routine
end


text:
using Printf

using DifferentialEquations

"calculate -i Heff = -i * (H - i/2*∑_i c_i^† c_i)"
get_miHeff(H,J) = -1im*(H - 0.5im*sum(j'*j for j in J))

mutable struct qutip_saver{T}
    saver::T
    outdir::String
    istep::Int
    qutip_saver(saver,outdir) = new{typeof(saver)}(saver,outdir,0)
end

function (f::qutip_saver)(u,t,integrator)
    data = f.saver(u,t,integrator)
    outfi = @sprintf("%s/%08d",f.outdir,f.istep+=1)
    qt.qsave((t,data),outfi)
    data
end

function get_savingcallback(saveat, saveddata=:probabilities, outdir=nothing)
    if saveddata == :probabilities
        savefunc = (u,t,integrator) -> collect(real(diag(u)))
        sv = SavedValues(Float64,Vector{Float64})
    elseif saveddata == :full
        savefunc = (u,t,integrator) -> collect(u) # collect to get a copy
        sv = SavedValues(Float64,Matrix{ComplexF64})
    elseif saveddata isa AbstractVector # e_ops
        savefunc = (u,t,integrator) -> ComplexF64[tr(e*u) for e in saveddata]
        sv = SavedValues(Float64,Vector{ComplexF64})
    else
        raise(ArgumentError("No function implemented for saveddata: $saveddata."))
    end

    if outdir === nothing
        saver = savefunc
    else
        mkdir(outdir)
        saver = qutip_saver(savefunc,outdir)
    end

    SavingCallback(saver, sv, saveat=saveat)
end

function setup(H,J,ρ0,ts,saveddata; backend, assume_herm=false)
    miHeff = get_miHeff(H,J)
    Lrefill = sum(kron.(conj.(J),J))
    if backend == :CUDA
        miHeffc = CuSparseMatrixCSR(miHeff)
        Lrefillc = CuSparseMatrixCSR(Lrefill)
        ρc = CuArray(ρ0)
        tmpρc = copy(ρc)
        tmpρc2 = copy(ρc)
        ps = (miHeffc,Lrefillc,tmpρc,tmpρc2,Val(assume_herm))
        if saveddata isa AbstractVector
            saveddata = CuSparseMatrixCSR.(saveddata)
        end
        return ρc, ts, ps, saveddata
    elseif backend == :CPU
        tmpρ = copy(ρ0)
        ps = (miHeff,Lrefill,tmpρ,tmpρ,Val(assume_herm))
        return ρ0, ts, ps, saveddata
    else
        raise(ArgumentError("Unknown backend $backend."))
    end
end

function mesolve(H,J,ρ0,ts;backend=:CUDA,saveddata=:probabilities,outdir=nothing)
    (ρ0,ts,ps,saveddata) = setup(H,J,ρ0,ts,saveddata; backend=backend)
    tspan = (ts[1],ts[end])
    cb = get_savingcallback(ts, saveddata, outdir)
    sv = cb.affect!.saved_values

    prob = ODEProblem(L!,ρ0,tspan,ps)
    sol = solve(prob,Tsit5(),save_start=false,save_end=false,save_everystep=false,
                callback=cb,reltol=1e-10,abstol=1e-12,maxiters=1e18)
    sol, sv
end

function mesolve_qutip(infile,outdir;backend=:CUDA,saveddata=nothing)
    input, saveddata_qu = load_input_from_qutip(infile)
    saveddata = saveddata === nothing ? saveddata_qu : saveddata
    println("writing output for propagation defined in '$infile' to directory '$outdir/'")
    mesolve(input...;backend=backend,saveddata=saveddata,outdir=outdir)
end


text:
module MendelImpute

    using LinearAlgebra
    using VCFTools
    using SnpArrays
    using BGEN
    using GroupSlices
    using Random
    using ProgressMeter
    using JLD2, FileIO, JLSO
    using Distances
    using LazyArrays
    using VariantCallFormat
    using CSV, DataFrames

    # functions for phasing/imputation/data-compression
    export phase               # main function for imputation and phasing
    export compress_haplotypes # for compressing haplotype panels
    export convert_compressed  # for ultra-compression outputs

    # functions for ancestry inference
    export admixture_global, admixture_local
    export paint, composition
    export thousand_genome_population_to_superpopulation
    export thousand_genome_samples_to_population
    export thousand_genome_samples_to_super_population

    include("compress.jl")
    include("data_structures.jl")
    include("phase.jl")
    include("haplotype_pair.jl")
    include("haplotype_rescreen.jl")
    include("haplotype_thinning.jl")
    include("haplotype_stepscreen.jl")
    include("dynamic_programming.jl")
    include("impute.jl")
    include("breakpoints.jl")
    include("intersect.jl")
    include("ancestry.jl")
    include("ultra_compress.jl")
    include("read.jl")

    # test data directory
    datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)
end # module

max_stars_repo_name: Hua-Zhou/MendelImpute


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

JLLWrappers.@generate_main_file_header("AsmJit")
JLLWrappers.@generate_main_file("AsmJit", UUID("1afa9899-fc99-5d84-9f06-5cc81cdd983b"))
end  # module AsmJit_jll


text:
# NOTE: tests are sensitive to the line number at which statements appear
function f1(x, y)
    # A comment
    return x + y
end

@noinline function throws()
    x = nothing
    error("oops")
end
@inline inlined() = throws()
call_throws() = inlined()

f2(x, y) = x + y

@inline function multilinesig(x::Int,
                              y::String)
    z = x + 1
    return z
end


text:
@testset "Non-nested" begin
    F.set_vars!(F.LOCAL_VARS, [
        "a" => "false",
        "b" => "false"])

    hs = """A{{if a}}B{{elseif b}}C{{else}}D{{end}}E"""

    tokens          = F.find_tokens(hs, F.HTML_TOKENS, F.HTML_1C_TOKENS)
    hblocks, tokens = F.find_all_ocblocks(tokens, F.HTML_OCB)
    qblocks         = F.qualify_html_hblocks(hblocks)

    F.set_vars!(F.LOCAL_VARS, ["a"=>"true","b"=>"false"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ABE")

    F.set_vars!(F.LOCAL_VARS, ["a"=>"false","b"=>"true"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ACE")

    F.set_vars!(F.LOCAL_VARS, ["a"=>"false","b"=>"false"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ADE")
end

@testset "Nested" begin
    F.def_LOCAL_VARS!()
    F.set_vars!(F.LOCAL_VARS, [
        "a" => "false",
        "b" => "false",
        "c" => "false"])

    hs = """A {{if a}} B {{elseif b}} C {{if c}} D {{end}} {{else}} E {{end}} F"""

    tokens          = F.find_tokens(hs, F.HTML_TOKENS, F.HTML_1C_TOKENS)
    hblocks, tokens = F.find_all_ocblocks(tokens, F.HTML_OCB)
    qblocks         = F.qualify_html_hblocks(hblocks)

    F.set_vars!(F.LOCAL_VARS, ["a"=>"true"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ABF")

    F.set_vars!(F.LOCAL_VARS, ["a"=>"false", "b"=>"true", "c"=>"false"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ACF")

    F.set_vars!(F.LOCAL_VARS, ["a"=>"false", "b"=>"true", "c"=>"true"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "ACDF")

    F.set_vars!(F.LOCAL_VARS, ["a"=>"false", "b"=>"false"])
    fhs = F.process_html_qblocks(hs, qblocks)
    @test isapproxstr(fhs, "AEF")
end

@testset "Bad cases" begin
    F.def_LOCAL_VARS!()
    # Lonely End block
    s = """A {{end}}"""
    @test_throws F.HTMLBlockError F.convert_html(s)

    # Inbalanced
    s = """A {{if a}} B {{if b}} C {{else}} {{end}}"""
    @test_throws F.HTMLBlockError F.convert_html(s)

    # Some of the conditions are not bools
    F.set_vars!(F.LOCAL_VARS, [
        "a" => "false",
        "b" => "false",
        "c" => "\"Hello\""])
    s = """A {{if a}} A {{elseif c}} B {{end}}"""
    @test_throws F.HTMLBlockError F.convert_html(s)
end


text:
edit(pkg) =
  isdir(Pkg.dir(pkg)) ?
    run(`atom $(Pkg.dir(pkg))`) :
    error("$pkg not installed")

isuntitled(p) = occursin(r"^(\.\\|\./)?untitled-[\d\w]+(:\d+)?$", p)

appendline(path, line) = line > 0 ? "$path:$line" : path

function realpath′(p)
  try
    ispath(p) ? realpath(p) : p
  catch e
    p
  end
end

function isfile′(p)
  try
    isfile(p)
  catch
    false
  end
end

function basepath(file)
  srcdir = joinpath(Sys.BINDIR,"..","..","base")
  releasedir = joinpath(Sys.BINDIR,"..","share","julia","base")
  normpath(joinpath(isdir(srcdir) ? srcdir : releasedir, file))
end

fullpath(path) =
  (isuntitled(path) || isabspath(path) ? path : basepath(path)) |> realpath′

function pkgpath(path)
  m = match(r"((?:[^/\\]+[/\\]){2}src[/\\].*)$", path)
  m == nothing ? basename(path) : m.captures[1]
end

expandpath(path) =
  isempty(path) ? (path, path) :
    path == "./missing" ? ("<unknown file>", path) :
      isuntitled(path) ? ("untitled", path) :
        !isabspath(path) ? (normpath(joinpath("base", path)), basepath(path)) :
          occursin(joinpath("julia", "stdlib"), path) ?
            begin
              p = last(split(path, joinpath("julia", "stdlib", "")))
              return (normpath(joinpath("stdlib", p)), isfile′(path) ? path :
                                                       normpath(joinpath(basepath(joinpath("..", "stdlib")), p)))
            end :
            (pkgpath(path), path)

function baselink(path, line)
  name, path = expandpath(path)
  name == "<unkown file>" ? span(".fade", "<unknown file>") :
                            link(path, line, Text(appendline(name, line)))
end

using Markdown: MD, HorizontalRule
function md_hlines(md)
  if !isa(md, MD) || !haskey(md.meta, :results) || isempty(md.meta[:results])
      return md
  end
  v = Any[]
  for (n, doc) in enumerate(md.content)
      push!(v, doc)
      n == length(md.content) || push!(v, HorizontalRule())
  end
  return MD(v)
end


text:
macro forward_func(ex, fs)
    T, field = ex.args[1], ex.args[2].value

    T = esc(T)
    fs = isexpr(fs, :tuple) ? map(esc, fs.args) : [esc(fs)]
    :($([:($f(x::$T, args...) = (Base.@_inline_meta; $f(x.$field, args...)))
        for f in fs]...);
    nothing)
end

@forward_func SEAssemblage.occ nthings, nplaces, occupancy, richness, nrecords, occurring, occupied, thingnames
@forward_func SEThings.commatrix nthings, nplaces, thingnames, placenames, occupancy, richness, nrecords, occurring, occupied
@forward_func SEAssemblage.site placenames

#--------------------------------------------------------------------------
# Basic summary functions

occupancy(com::AbstractComMatrix) = occupancy(occurrences(com))
richness(com::AbstractComMatrix) = richness(occurrences(com))
occurring(com::AbstractComMatrix, idx...) = occurring(occurrences(com), idx...)
occupied(com::AbstractComMatrix, idx...) = occupied(occurrences(com), idx...)

const nspecies = nthings
nthings(com::AbstractComMatrix) = size(com.occurrences, 1)

const nsites = nplaces
nplaces(com::AbstractComMatrix) = size(com.occurrences, 2)
nplaces(sd::SiteData) = size(coordinates(sd.site), 1)
nplaces(sd::SELocations) = DataFrames.ncol(sd.sitestats)
nplaces(gr::GridData) = size(gr.indices, 1)
nplaces(pd::PointData) = size(pd.coords, 1)

nrecords(com::AbstractComMatrix) = _nnz(occurrences(com))

const getspecies = thingoccurrences
thingoccurrences(com::AbstractComMatrix, idx) = thingoccurrences(occurrences(com), idx)

const getsite = placeoccurrences
placeoccurrences(com::AbstractComMatrix, idx) = placeoccurrences(occurrences(com), idx)

const speciesnames = thingnames
thingnames(com::AbstractComMatrix) = com.speciesnames

const sitenames = placenames
placenames(com::AbstractComMatrix) = com.sitenames
placenames(sd::SiteData) = sitenames(sd.site)
placenames(sd::SELocations) = collect(sd.sitestats[:sites])

sitetotals(com::AbstractComMatrix) = vec(colsum(com.occurrences))
speciestotals(com::AbstractComMatrix) = vec(rowsum(com.occurrences))

size(com::AbstractComMatrix) = size(occurrences(com))
size(com::AbstractComMatrix, dims...) = size(occurrences(com), dims...)

#----------------------------------------------------------------------------------------------
# Mutating and concatenating

#-------------------------------------------------------------------------------------------------
# show methods

summary(com::AbstractComMatrix) = "$(nsites(com))x$(nspecies(com)) $(typeof(com))"

function createsummaryline(vec::AbstractVector{T}) where T<:AbstractString
    linefunc(vec) = mapreduce(x->x*", ", *, vec[1:(end-1)])*vec[end]
    length(vec) == 1 && return vec[1]
    length(vec) < 6 && return linefunc(vec)
    linefunc(vec[1:3])*"..."*linefunc(vec[(end-1):end])
end

function show(io::IO, com::ComMatrix)
    sp = createsummaryline(speciesnames(com))
    si = createsummaryline(sitenames(com))
    println(io, "ComMatrix with $(nspecies(com)) species in $(nsites(com)) sites\n\nSpecies names:\n$(sp)\n\nSite names:\n$(si)")
end

function show(io::IO, com::Assemblage)
    sp = createsummaryline(speciesnames(com))
    si = createsummaryline(sitenames(com))
    println(io, "Assemblage with $(nspecies(com)) species in $(nsites(com)) sites\n\nSpecies names:\n$(sp)\n\nSite names:\n$(si)")
end

function show(io::IO, sd::SiteData)
    println(io, "Spatial data set with $(nsites(sd)) sites\n\nSite names:\n$(createsummarylines(sitenames(sd)))")
end

#TODO also create render functions for Juno

getindex(com::AbstractComMatrix, inds...) = ComMatrix(getindex(com.occurrences, inds...))

setindex!(com::AbstractComMatrix, X, inds...) = setindex!(com.occurrences, X, inds...)

function getindex(site::S, inds) where S<:SELocations
  S(coordinates(site)[inds,:], site.sitestats[inds,:])
end

function getindex(com::SEAssemblage, ind::Symbol)
    if ind in names(com.site.sitestats)
        return com.site.sitestats[ind]
    elseif ind in names(com.occ.traits)
        return com.occ.traits[ind]
    else
        error("No such name in traits or sitestats")
    end
end


text:
"""
测试正方晶系
"""

using Test
using MARY_fRG.Basics, MARY_fRG.Fermi

@testset "正方晶系" begin
    disp(x, y) = -2 * (cos(x) + cos(y))
    ss = QuadrateSystem{:COM}([disp])
    println(ss)
    println(ss.brillouin.vertex)
    pt1 = Point2D(2., 2.)
    pt2 = ss.kadd(pt1, pt1)
    @test pt2.x == 4 - 2*pi
    pt3 = ss.kadd(pt2, -pt1)
    @test pt3.x == 2
end

@testset "三角晶系" begin
    disp(x, y) = -2*numpy.cos(kxv) - 4*numpy.cos(kyp)*numpy.cos(kxv/2)
    ts = TriangularSystem{:COM}([disp])
    println(ts)
    println(ts.brillouin.vertex)
    pt1 = Point2D(2., 2.)
    pt2 = ts.kadd(pt1, pt1)
    @test pt2.y == 4 - 2pi/sqrt(3)
    @test pt2.x == 4 - 2pi
    pt3 = Point2D(2pi, 0.)
    pt3 = ts.kadd(pt2, pt3)
    @test pt3.y == 4 - 2pi/sqrt(3) - 2pi/sqrt(3)
    @test pt3.x == 4 + 2pi - 4pi
    for idx in 1:1:1000000
        pt1 = Point2D((rand()-0.5)*4π, (rand()-0.5)*4π)
        pt2 = Point2D((rand()-0.5)*4π, (rand()-0.5)*4π)
        pt3 = Fermi.hexagon_kadd(pt1, pt2)
        pt4 = Fermi.hexagon_kadd2(pt1, pt2)
        pt5 = Fermi.hexagon_kadd3(pt1, pt2)
        #println(pt3)
        #println(pt4)
        #if pt3.x != pt4.x || pt3.y != pt4.y
        #    println(pt1)
        #    println(pt2)
        #end
        @test isapprox(pt3.x, pt4.x, atol=1e-12)
        @test isapprox(pt3.y, pt4.y, atol=1e-12)
        @test isapprox(pt3.x, pt5.x, atol=1e-12)
        @test isapprox(pt3.y, pt5.y, atol=1e-12)
    end
    #println(Fermi.hexagon_kadd2(
    #    Point2D(0.8806377184930461, 5.305646990668821),
    #    Point2D(-6.126445471920293, -6.016421208448352)
    #))
end


text:
using Base.Test

excluded = ["runtests.jl"]

if is_windows() && Sys.WORD_SIZE == 32
  push!(excluded, "except.jl")
end

@testset "CxxWrap tests" begin
  @testset "$f" for f in filter(fname -> fname ∉ excluded, readdir())
    include(f)
  end
end


text:
# add one unit to the voltage limit
# SU+
export psp_incrementvlimit

const INCREMENT_VLIMIT = [0x53, 0x55, 0x2B, 0x0D]

"""
    psp_incrementvlimit(io)

Add one unit to the voltage limit.
"""
function psp_incrementvlimit(io_psp::IO)
  write(io_psp, INCREMENT_VLIMIT)
  return nothing
end


text:
using ParserCombinator
import Dates

space = Drop(Star(Space()))

PNonZeroInt64() = Parse(p"-?[1-9][0-9]*", Int64)
PPosInt64() = Parse(p"[1-9][0-9]*", Int64)
PPosZeroInt64() = Parse(p"[0-9][0-9]*", Int64)

@with_pre space begin
    sum = Delayed()

    weekday_short = Alt(map(x -> Pattern(uppercase(x)), Dates.ENGLISH.days_of_week_abbr)...)
    month_short = Alt(map(x -> Pattern(uppercase(x)), Dates.ENGLISH.months_abbr)...)

    brackets = E"(" + space + sum + space + E")"

    rdate_term = Alt()
    rdate_expr = rdate_term | brackets

    neg = Delayed()
    neg.matcher = rdate_expr | (E"-" + neg > -)
    cal_adj = neg + (E"@" + p"[a-zA-Z\\\\\\s\\|]+" + E"[" + Alt(map(Pattern, collect(keys(HOLIDAY_ROUNDING_MAPPINGS)))...) + E"]")[0:1] |> xs -> length(xs) == 1 ? xs[1] : CalendarAdj(map(String, split(xs[2], "|")), xs[1], HOLIDAY_ROUNDING_MAPPINGS[xs[3]])

    # Our decision on the grammar choice for rdates is purely one of history.
    # * will imply 'no roll', so 2*1m == 2m
    # *roll will imply 'roll', so 2*roll(1m) == 1m + 1m
    mult_roll = cal_adj | ((PPosInt64() + space + E"*" + space + E"roll(" + space + sum + space + E")") > (c, rd) -> multiply_roll(rd, c)) | ((E"roll(" + space + sum + space + E")" + space + E"*" + space + PPosInt64()) > (rd, c) -> multiply_roll(rd, c))

    mult_no_roll = Delayed()
    mult_no_roll = mult_roll | ((PPosInt64() + space + E"*" + space + mult_roll) > (c,rd) -> multiply_no_roll(rd, c)) | ((mult_roll + space + E"*" + space + PPosInt64()) > (rd,c) -> multiply_no_roll(rd, c))

    add = E"+" + mult_no_roll
    sub = E"-" + mult_no_roll > -
    sum.matcher = (sub | mult_no_roll) + (add | mult_no_roll)[0:end] |> Base.sum

    entry = sum + Eos()
end

function register_grammar!(term)
    # Handle the spacing correctly
    push!(rdate_term.matchers[2].matchers, term)
end

macro rd_str(arg::String)
    val = parse_one(arg, entry)[1]
    isa(val, RDate) || error("Unable to parse $(arg) as RDate")
    return val
end

function rdate(arg::String)
    val = parse_one(arg, entry)[1]
    isa(val, RDate) || error("Unable to parse $(arg) as RDate")
    return val
end


text:
import Dates

mstosec(t) = t/1000
sectomin(t) = t/60
mintohrs(t) = t/60
hrstodays(t) = t/24
daystoyears(t) = t/365

mstomin(t) = t |> mstosec |> sectomin
mstohrs(t) = t |>  mstomin |> mintohrs
mstodays(t) = t |> mstohrs |> hrstodays
mstoyears(t) = t |> mstodays |> daystoyears

"""
    secondsfloat(t::Dates.Millisecond)::Float64

Convert `t` to seconds.
"""
secondsfloat(t::Dates.Millisecond)::Float64 = t.value |> mstosec

"""
    minutesfloat(t::Dates.Millisecond)::Float64

Convert `t` to minutes.
"""
minutesfloat(t::Dates.Millisecond)::Float64 = t.value |> mstomin

"""
    hoursfloat(t::Dates.Millisecond)::Float64

Convert `t` to hours.
"""
hoursfloat(t::Dates.Millisecond)::Float64 = t.value |> mstohrs

"""
    daysfloat(t::Dates.Millisecond)::Float64

Convert `t` to days.
"""
daysfloat(t::Dates.Millisecond)::Float64 = t.value |> mstodays

"""
    yearsfloat(t::Dates.Millisecond)::Float64

Convert `t` to years.
"""
yearsfloat(t::Dates.Millisecond)::Float64 = t.value |> mstoyears

nothing


text:
@testset "Test LAStopping" begin
  ###############################################################################
  #
  # The Stopping structure eases the implementation of algorithms and the
  # stopping criterion.
  #
  # The following examples illustrate a specialized Stopping for the
  # LinearOperators type.
  # https://github.com/JuliaSmoothOptimizers/LinearOperators.jl
  #
  using LinearAlgebra, SparseArrays, LinearOperators #v.1.0.1

  ##############################################################
  #TEST

  n = 5
  A = rand(n, n)
  xref = 100 * rand(n)
  x0 = zeros(n)
  b = A * xref
  sA = sparse(A)
  opA = LinearOperator(A)

  mLO = LinearSystem(A, b)
  sLO = LLSModel(sA, b)
  opLO = LinearSystem(opA, b)
  meta = StoppingMeta(max_cntrs = init_max_counters_linear_operators())
  mLOstp_meta = LAStopping(mLO, meta, GenericState(x0)) #this is different because of optimality_check
  mLOstp = LAStopping(mLO, GenericState(x0), max_cntrs = init_max_counters_linear_operators())
  sLOstp = LAStopping(sLO, GenericState(x0))
  short_stop = LAStopping(A, b, sparse = true) #note that sparse is true by default
  maxcn = init_max_counters_linear_operators(nprod = 1)
  opLOstp = LAStopping(opLO, GenericState(x0), max_cntrs = maxcn)

  opLOstp2 = LAStopping(opLO, GenericState(x0), optimality_check = linear_system_check)
  @test opLOstp2.meta.max_cntrs == init_max_counters_linear_operators()
  @test opLOstp2.meta.optimality_check == opLOstp.meta.optimality_check

  mLOstp_src = LAStopping(mLO, meta, StopRemoteControl(), GenericState(x0))

  @test typeof(mLOstp.pb.counters) <: LACounters
  @test typeof(sLOstp.pb.counters) <: NLSCounters
  @test typeof(opLOstp.pb.counters) <: LACounters
  @test typeof(short_stop.pb.counters) <: NLSCounters
  @test typeof(short_stop.pb) <: LLSModel

  @test start!(mLOstp) == false
  @test start!(sLOstp) == false
  @test start!(opLOstp) == false

  update_and_stop!(mLOstp, x = xref, res = Float64[])
  @test status(mLOstp) == :Optimal
  @test mLOstp.meta.resources == false
  update_and_stop!(sLOstp, x = xref, res = Float64[])
  @test status(sLOstp) == :Optimal
  @test sLOstp.meta.resources == false
  update_and_stop!(opLOstp, x = xref, res = Float64[])
  @test status(opLOstp) == :Optimal
  @test opLOstp.meta.resources == true

  include("_additional_linear_operator_algebra.jl")

  """
  Randomized block Kaczmarz
  """
  function RandomizedBlockKaczmarz(stp::AbstractStopping; kwargs...)
    b = stp.pb.b
    x0 = stp.current_state.x
    m, n = length(b), length(x0)
    A = try
      stp.pb.A #LSSModel.A does not exist
    catch
      sparse(stp.pb.Arows, stp.pb.Acols, stp.pb.Avals, m, n)
    end

    xk = x0

    OK = start!(stp)

    while !OK
      i = Int(floor(rand() * m) + 1) #rand a number between 1 and m
      Ai = A[i, :]
      xk = Ai == 0 ? x0 : x0 - (dot(Ai, x0) - b[i]) / dot(Ai, Ai) * Ai
      #xk  = Ai == 0 ? x0 : x0 - (Ai' * x0-b[i])/(Ai' * Ai) * Ai

      if xtype(stp.current_state) <: SparseVector && !(typeof(xk) <: SparseVector)
        OK = update_and_stop!(stp, x = sparse(xk), res = sparse(Float64[]))
      elseif xtype(stp.current_state) <: Vector && typeof(xk) <: SparseVector
        OK = update_and_stop!(stp, x = convert(Vector, xk), res = Float64[])
      else
        OK =
          update_and_stop!(stp, x = xk, res = Stopping._init_field(typeof(stp.current_state.res)))
      end
      x0 = xk
    end

    return stp
  end

  m, n = 200, 100 #size of A: m x n
  A = 100 * rand(m, n) #It's a dense matrix :)
  xref = 100 * rand(n)
  b = A * xref

  #Our initial guess
  x0 = zeros(n)

  la_stop = LAStopping(
    A,
    b,
    GenericState(x0),
    max_iter = 150000,
    rtol = 1e-6,
    max_cntrs = init_max_counters_NLS(residual = 150000),
  )
  #Be careful using GenericState(x0) would not work here without forcing convert = true
  #in the update function. As the iterate will be a SparseVector to the contrary of initial guess.
  #Tangi: maybe start! should send a Warning for such problem !?
  sa_stop = LAStopping(
    sparse(A),
    b,
    GenericState(sparse(x0)),
    max_iter = 150000,
    rtol = 1e-6,
    max_cntrs = init_max_counters_NLS(residual = 150000),
  )
  op_stop = LAStopping(
    LinearSystem(LinearOperator(A), b),
    GenericState(x0),
    max_iter = 150000,
    rtol = 1e-6,
    max_cntrs = init_max_counters_linear_operators(nprod = 150000),
  )
  opbis_stop = LAStopping(LinearOperator(A), b)

  try
    @timed RandomizedBlockKaczmarz(la_stop)
    @test status(la_stop) == :Optimal
    @timed RandomizedBlockKaczmarz(sa_stop)
    @test status(sa_stop) == :Optimal
  catch
    @warn "If LSSModel.A does not exist consider [la_stop.pb.Avals[i,j] for (i) in la_stop.pb.Arows, j in la_stop.pb.Acols]"
    #https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/lls_model.jl
  end

  @timed RandomizedBlockKaczmarz(op_stop)
  @test status(op_stop) == :Optimal

  update!(la_stop.current_state, x = xref)
  @test normal_equation_check(la_stop.pb, la_stop.current_state) <= la_stop.meta.atol
  update!(op_stop.current_state, x = xref)
  @test normal_equation_check(op_stop.pb, op_stop.current_state) <= la_stop.meta.atol
end


text:
using Optim, LineSearches
using LinearAlgebra: I, norm
using TimerOutputs
using BCVUMPS: bigleftenv, bigrightenv, ALCtoAC, obs2x2FL, obs2x2FR

"""
    energy(h, bcipeps; χ, tol, maxiter)

return the energy of the `bcipeps` 2-site hamiltonian `h` and calculated via a
BCVUMPS with parameters `χ`, `tol` and `maxiter`.
"""
function energy(h, model::HamiltonianModel, bcipeps::BCIPEPS; χ::Int, tol::Real, maxiter::Int, verbose = false)
    bcipeps = indexperm_symmetrize(bcipeps)  # NOTE: this is not good
    D = getd(bcipeps)^2
    s = gets(bcipeps)
    Ni,Nj = size(bcipeps.bulk)
    ap = [ein"abcdx,ijkly -> aibjckdlxy"(bcipeps.bulk[i], conj(bcipeps.bulk[i])) for i = 1:Ni*Nj]
    ap = [reshape(ap[i], D, D, D, D, s, s) for i = 1:Ni*Nj]
    ap = reshape(ap, Ni, Nj)
    a = [ein"ijklaa -> ijkl"(ap[i]) for i = 1:Ni*Nj]
    a = reshape(a, Ni, Nj)

    folder = "./data/$(model)/"
    mkpath(folder)
    chkp_file = folder*"bcvumps_env_$(Ni)x$(Nj)_D$(D)_chi$(χ).jld2"
    if isfile(chkp_file)
        rt = SquareBCVUMPSRuntime(a, chkp_file, χ; verbose = verbose)
    else
        rt = SquareBCVUMPSRuntime(a, Val(:random), χ; verbose = verbose)
    end

    env = bcvumps(rt; tol=tol, maxiter=maxiter, verbose = verbose)
    save(chkp_file, "env", env)
    e = expectationvalue(h, ap, env)
    return e
end

"""
    expectationvalue(h, ap, env)

return the expectationvalue of a two-site operator `h` with the sites
described by rank-6 tensor `ap` each and an environment described by
a `SquareBCVUMPSRuntime` `env`.
"""
function expectationvalue(h, ap, env::SquareBCVUMPSRuntime)
    M,AL,C,AR,FL,FR = env.M,env.AL,env.C,env.AR,env.FL,env.FR
    Ni,Nj = size(M)
    # Ni,Nj = 1,2
    ap /= norm(ap)
    etol = 0
    for j = 1:Nj, i = 1:Ni
        # ir = i + 1 - Ni * (i == Ni)
        jr = j + 1 - (j==Nj) * Nj
        _, FL = obs2x2FL(AL, M, FL)
        _, FR = obs2x2FR(AR, M, FR)
        e = ein"abc,cde,anm,ef,ml,fgh,lkj,hij,bnodpq,okigrs,pqrs -> "(FL[i,j],AL[i,j],conj(AL[i,j]),C[i,j],conj(C[i,j]),AR[i,jr],conj(AR[i,jr]),FR[i,jr],ap[i,j],ap[i,jr],h)[]
        n = ein"abc,cde,anm,ef,ml,fgh,lkj,hij,bnodpq,okigrs -> pqrs"(FL[i,j],AL[i,j],conj(AL[i,j]),C[i,j],conj(C[i,j]),AR[i,jr],conj(AR[i,jr]),FR[i,jr],ap[i,j],ap[i,jr])
        n = ein"pprr -> "(n)[]
        # @show i,j,e/n
        etol += e/n
    end

    # Zygote.@ignore begin
    #     AC = ALCtoAC(AL,C)
    #     for j = 1:Nj, i = 1:Ni
    #         ir = i + 1 - Ni * (i==Ni)
    #         # irr = i + 2 - Ni * (i + 2 > Ni)
    #         _, BgFL = bigleftenv(AL, M)
    #         _, BgFR = bigrightenv(AR, M)
    #         e2 = ein"dcba,def,aji,fghi,ckgepq,bjhkrs,pqrs -> "(BgFL[i,j],AC[i,j],conj(AC[ir,j]),BgFR[i,j],ap[i,j],ap[ir,j],h)[]
    #         n2 = ein"dcba,def,aji,fghi,ckgepq,bjhkrs -> pqrs"(BgFL[i,j],AC[i,j],conj(AC[ir,j]),BgFR[i,j],ap[i,j],ap[ir,j])
    #         n2 = ein"pprr -> "(n2)[]
    #         @show i,j,e2/n2
    #     end
    # end

    return etol/Ni/Nj
end

"""
    ito12(i,Ni)

checkerboard pattern
```
    │    │
  ──A────B──
    │    │
  ──B────A──
    │    │
```
"""
ito12(i,Ni) = mod(mod(i,Ni) + Ni*(mod(i,Ni)==0) + fld(i,Ni) + 1 - (mod(i,Ni)==0), 2) + 1

buildbcipeps(bulk,Ni,Nj) = reshape([bulk[:,:,:,:,:,ito12(i,Ni)] for i = 1:Ni*Nj], (Ni, Nj))

"""
    init_ipeps(model::HamiltonianModel; D::Int, χ::Int, tol::Real, maxiter::Int)

Initial `bcipeps` and give `key` for use of later optimization. The key include `model`, `D`, `χ`, `tol` and `maxiter`.
The iPEPS is random initial if there isn't any calculation before, otherwise will be load from file `/data/model_D_chi_tol_maxiter.jld2`
"""
function init_ipeps(model::HamiltonianModel; D::Int, χ::Int, tol::Real, maxiter::Int, verbose = true)
    folder = "./data/$(model)/"
    mkpath(folder)
    key = (model, D, χ, tol, maxiter)
    chkp_file = folder*"$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).jld2"
    if isfile(chkp_file)
        bulk = load(chkp_file)["bcipeps"]
        verbose && println("load BCiPEPS from $chkp_file")
    else
        bulk = rand(D,D,D,D,2,2)
        verbose && println("random initial BCiPEPS")
    end
    Ni, Nj = model.Ni, model.Nj
    bcipeps = SquareBCIPEPS(buildbcipeps(bulk,Ni,Nj))
    bcipeps = indexperm_symmetrize(bcipeps)
    return bcipeps, key
end

"""
    optimiseipeps(bcipeps, h; χ, tol, maxiter, optimargs = (), optimmethod = LBFGS(m = 20))

return the tensor `bulk'` that describes an bcipeps that minimises the energy of the
two-site hamiltonian `h`. The minimization is done using `Optim` with default-method
`LBFGS`. Alternative methods can be specified by loading `LineSearches` and
providing `optimmethod`. Other options to optim can be passed with `optimargs`.
The energy is calculated using vumps with key include parameters `χ`, `tol` and `maxiter`.
"""
function optimiseipeps(bcipeps::BCIPEPS{LT}, key; f_tol = 1e-6, verbose= false, optimmethod = LBFGS(m = 20)) where LT
    model, D, χ, tol, maxiter = key
    h = hamiltonian(model)
    Ni, Nj = model.Ni, model.Nj
    to = TimerOutput()
    f(x) = @timeit to "forward" real(energy(h, model, BCIPEPS{LT}(buildbcipeps(x,Ni,Nj)); χ=χ, tol=tol, maxiter=maxiter, verbose=verbose))
    ff(x) = real(energy(h, model, BCIPEPS{LT}(buildbcipeps(x,Ni,Nj)); χ=χ, tol=tol, maxiter=maxiter, verbose=verbose))
    g(x) = @timeit to "backward" Zygote.gradient(ff,x)[1]
    x0 = zeros(D,D,D,D,2,2)
    x0[:,:,:,:,:,1] = bcipeps.bulk[1,1]
    x0[:,:,:,:,:,2] = bcipeps.bulk[2,1]
    res = optimize(f, g,
        x0, optimmethod,inplace = false,
        Optim.Options(f_tol=f_tol,
        extended_trace=true,
        callback=os->writelog(os, key)),
        )
    println(to)
    return res
end

"""
    writelog(os::OptimizationState, key=nothing)

return the optimise infomation of each step, including `time` `iteration` `energy` and `g_norm`, saved in `/data/model_D_chi_tol_maxiter.log`. Save the final `bcipeps` in file `/data/model_D_chi_tol_maxiter.jid2`
"""
function writelog(os::OptimizationState, key=nothing)
    message = "$(round(os.metadata["time"],digits=2))   $(os.iteration)   $(os.value)   $(os.g_norm)\n"

    printstyled(message; bold=true, color=:red)
    flush(stdout)

    model, D, χ, tol, maxiter = key
    if !(key === nothing)
        logfile = open("./data/$(model)/$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).log", "a")
        write(logfile, message)
        close(logfile)
        save("./data/$(model)/$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).jld2", "bcipeps", os.metadata["x"])
    end
    return false
end


text:
__precompile__()

module MXNet

# we put everything in the namespace mx, because there are a lot of
# functions with the same names as built-in utilities like "zeros", etc.
export mx
module mx

using Compat
import Compat.String
import Compat.view

using Formatting

# Functions from base that we can safely extend and that are defined by libmxnet.
import Base: round, ceil, floor, cos, sin, abs, sign, exp, sqrt, exp, log, norm,
             transpose

include("base.jl")
include("context.jl")
include("util.jl")

include("ndarray.jl")
include("random.jl")

include("name.jl")
include("symbolic-node.jl")
include("executor.jl")

include("metric.jl")
include("optimizer.jl")
include("initializer.jl")

include("io.jl")
include("kvstore.jl")

include("callback.jl")
include("model.jl")

include("visualize.jl")

include("nn-factory.jl")

end # mx

end # module MXNet


text:
using Pkg, Pkg.TOML

root = joinpath(@__DIR__, "..")

meta = TOML.parsefile(joinpath(@__DIR__, "Notebooks.toml"))
meta = meta[ARGS[1]]

path = meta["path"]
deps = get(meta, "deps", [])

for d in ["Project.toml", "Manifest.toml", ".gitignore"]
  isfile(joinpath(root, path, d)) && push!(deps, d)
end

mkpath(joinpath(root, "notebooks", path))
for dep in deps
  cp(joinpath(root, path, dep), joinpath(root, "notebooks", path, dep), force = true)
end

pushfirst!(LOAD_PATH, @__DIR__)
Pkg.activate(joinpath(root, "notebooks", path))

using Literate

function init_nb(content)
	content = replace(content, r"#\s*using CuArrays" => "## using CuArrays")
	content = "using Pkg; Pkg.activate(\".\"); Pkg.instantiate();\n\n" * content
	return content
end

scripts = meta["notebook"]
scripts isa String && (scripts = [scripts])

for script in scripts
  Literate.notebook(joinpath(root, path, script),
                    joinpath(root, "notebooks", path),
                    credit = false, preprocess = init_nb)
end

scripts = map(x -> x[1:end - 3] * ".ipynb", scripts)
keep = union(deps, scripts)
files = readdir(joinpath(root, "notebooks", path))

for r in files
  r in keep || rm(joinpath(root, "notebooks", path, r), force = true)
end

max_stars_repo_name: c-p-murphy/model-zoo

max_issues_repo_name: c-p-murphy/model-zoo

max_forks_repo_name: c-p-murphy/model-zoo


text:
using AlgebraOfGraphics, Test
using AlgebraOfGraphics: data,
                         style,
                         Spec,
                         spec,
                         layers,
                         positional,
                         keyword,
                         dims,
                         run_pipeline

using CategoricalArrays: categorical
using NamedDims
using RDatasets: dataset
using Observables: to_value
using AbstractPlotting: default_palettes

@testset "product" begin
    s = dims() * style(1:2, ["a", "b"], color = dims(1))
    ps = pairs(s)
    @test ps[1] == Pair((color = 1,), style(1, "a"))
    @test ps[2] == Pair((color = 2,), style(2, "b"))
end

@testset "contexts" begin
    mpg = dataset("ggplot2", "mpg")
    d = style(:Cyl, :Hwy, color = :Year => categorical)
    s = spec(color = :red, font = 10) + style(markersize = :Year)
    res = data(mpg) * d * s
    @test res[1].options == (color = :red, font = 10)

    idx1 = mpg.Year .== 1999
    idx2 = mpg.Year .== 2008

    styles = pairs.(getproperty.(res, :style))
    @test last(styles[1][1]).value == style(mpg[idx1, :Cyl], mpg[idx1, :Hwy]).value
    @test last(styles[1][2]).value == style(mpg[idx2, :Cyl], mpg[idx2, :Hwy]).value
    @test last(styles[2][1]).value == style(mpg[idx1, :Cyl], mpg[idx1, :Hwy], markersize = mpg[idx1, :Year]).value
    @test last(styles[2][2]).value == style(mpg[idx2, :Cyl], mpg[idx2, :Hwy], markersize = mpg[idx2, :Year]).value

    @test first(styles[1][1]).color isa NamedDimsArray
    @test only(first(styles[1][1]).color) == 1999
    @test dimnames(first(styles[1][1]).color) == (:Year,)

    @test first(styles[1][2]).color isa NamedDimsArray
    @test only(first(styles[1][2]).color) == 2008
    @test dimnames(first(styles[1][2]).color) == (:Year,)

    @test first(styles[2][1]).color isa NamedDimsArray
    @test only(first(styles[2][1]).color) == 1999
    @test dimnames(first(styles[2][1]).color) == (:Year,)

    @test first(styles[2][2]).color isa NamedDimsArray
    @test only(first(styles[2][2]).color) == 2008
    @test dimnames(first(styles[2][2]).color) == (:Year,)

    x = rand(5, 3, 2)
    y = rand(5, 3)
    s = dims(1) * style(x, y, color = dims(2))

    res = pairs(s)
    for (i, r) in enumerate(res)
        group, st = r
        @test group == (; color = mod1(i, 3))
        xsl = x[:, mod1(i, 3), (i > 3) + 1]
        ysl = y[:, mod1(i, 3)]
        @test st.value == style(xsl, ysl).value
    end
end

@testset "compute specs" begin
    wong = default_palettes[:color][]
    t = (x = [1, 2], y = [10, 20], z = [3, 4], c = ["a", "b"])
    d = style(:x, :y, color = :c)
    s = spec(:log) * spec(font = 10) + style(size = :z)
    ds = data(t) * d
    sl = ds * s
    @test length(sl) == 2

    res = run_pipeline(sl)

    @test res[1].options.color == wong[1]
    @test res[2].options.color == wong[2]
    @test dimnames(res[1].pkeys.color) == (:c,)
    @test dimnames(res[2].pkeys.color) == (:c,)
    @test res[1].style.value == style([1], [10]).value
    @test res[2].style.value == style([2], [20]).value

    @test res[3].options.color == wong[1]
    @test res[4].options.color == wong[2]
    @test dimnames(res[3].pkeys.color) == (:c,)
    @test dimnames(res[4].pkeys.color) == (:c,)
    @test res[3].style.value == style([1], [10], size = [3]).value
    @test res[4].style.value == style([2], [20], size = [4]).value
end


text:
# ---
# cover: assets/swirl.gif
# title: Swirl effect using warp operation
# author: Johnny Chen
# date: 2021-07-12
# ---

# In this example, we illustrate how to construct a custom warping map
# and pass it to `warp`. This swirl example comes from
# [the Princeton Computer Graphics course for Image Warping (Fall 2000)](https://www.cs.princeton.edu/courses/archive/fall00/cs426/lectures/warp/warp.pdf)
# and [scikit-image swirl example](https://scikit-image.org/docs/dev/auto_examples/transform/plot_swirl.html).

using ImageTransformations
using OffsetArrays, StaticArrays
using ImageShow, TestImages
using LinearAlgebra

img = imresize(testimage("cameraman"), (256, 256));

# As we've illustrated in [image warping](@ref index_image_warping), a warp operation
# consists of two operations: backward coordinate map `ϕ` and intensity estimator.
# To implement swirl operation, we need to customize the coordinate map `ϕ`.
# A valid coordinate map `q = ϕ(p)` follows the following interface:
#
# ```julia
# # SVector comes from StaticArrays
# ϕ(::SVector{N})::SVector{N} where N
# ```
#
# A cartesian position `(x, y)` can be transfered to/from polar coordinate `(ρ, θ)`
# using formula:
#
# ```julia
# # Cartesian to Polar
# ρ = norm(y-y0, x-x0)
# θ = atan(y/x)
#
# # Polar to Cartesian
# y = y0 + ρ*sin(θ)
# x = x0 + ρ*cos(θ)
# ```
#
# For given input index `p`, a swirl operation enforces more rotations in its polar coordinate using
# `θ̃ = θ + ϕ + s*exp(-ρ/r)`, and returns the cartesian index (x̃, ỹ) from the warped polor coordinate
# (ρ, θ̃). (Here we use the formula from [scikit-image swirl example](https://scikit-image.org/docs/dev/auto_examples/transform/plot_swirl.html)
# to build our version.)

function swirl(rotation, strength, radius)
    x0 = OffsetArrays.center(img)
    r = log(2)*radius/5

    function swirl_map(x::SVector{N}) where N
        xd = x .- x0
        ρ = norm(xd)
        θ = atan(reverse(xd)...)

        ## Note that `x == x0 .+ ρ .* reverse(sincos(θ))`
        ## swirl adds more rotations to θ based on the distance to center point
        θ̃ = θ + rotation + strength * exp(-ρ/r)

        SVector{N}(x0 .+ ρ .* reverse(sincos(θ̃)))
    end

    warp(img, swirl_map, axes(img))
end

# Now let's see how radius argument affects the result

preview = ImageShow.gif([swirl(0, 10, radius) for radius in 10:10:150]; fps=5)

using FileIO #src
save("assets/swirl.gif", preview; fps=5) #src


text:
using ModelingToolkit, Test
using ModelingToolkit: value
using SymbolicUtils: <ₑ
@parameters α β δ
expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
ref = sort([β, δ, α], lt = <ₑ)
sol = sort(Num.(ModelingToolkit.get_variables(expr)), lt = <ₑ)
@test all(x->x isa Num, sol[i] == ref[i] for i in 1:3)
@test all(simplify∘value, sol[i] == ref[i] for i in 1:3)

@parameters γ
s = α => γ
expr = (((1 / β - 1) + δ) / α) ^ (1 / (α - 1))
sol = ModelingToolkit.substitute(expr, s)
new = (((1 / β - 1) + δ) / γ) ^ (1 / (γ - 1))
@test iszero(sol - new)


text:
## Convergence tests on silicon for different lattice constants (higher close the gap)
## see Cancès/Kemlin/Levitt, Convergence analysis of direct minimization and self-consistent iterations, https://arxiv.org/abs/2004.09088
## Some parameter values have been reduced for faster automatic testing

using DFTK
using Plots
using Printf
using LinearAlgebra

# Calculation parameters
kgrid = [1, 1, 1]
Ecut = 15 # 30 in the paper
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]
tol = 1e-10
diagtol = 1e-12
# 100 in the paper; a high value is important to see the divergence of the simple damping algorithms
# Very reduced here for fast automated testing
maxiter = 5

global resids = []
function my_callback(info)
    info.stage == :finalize && return
    global resids
    info.n_iter == 1 && (resids = [])
    err = norm(info.ρout.fourier - info.ρin.fourier)
    println(info.n_iter, " ", err)
    push!(resids, err)
end
my_isconverged = info -> norm(info.ρout.fourier - info.ρin.fourier) < tol
opts = (callback=my_callback, is_converged=my_isconverged, maxiter=maxiter, tol=tol,
        determine_diagtol=info -> diagtol)

global errs = []
global gaps = []
## is (10.26, 11.2, 11.405) in the paper
global as = (10.26, 11.405)
for a in as
    lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
    model = model_LDA(lattice, atoms)
    basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
    res = self_consistent_field(basis; opts...)
    gap = res.eigenvalues[1][5] - res.eigenvalues[1][4]
    errs_anderson = copy(resids)
    self_consistent_field(basis; solver=scf_damping_solver(1), opts...)
    errs_1 = copy(resids)
    self_consistent_field(basis; solver=scf_damping_solver(.5), opts...)
    errs_05 = copy(resids)
    self_consistent_field(basis; solver=scf_damping_solver(.2), opts...)
    errs_02 = copy(resids)
    self_consistent_field(basis; solver=scf_damping_solver(.1), opts...)
    errs_01 = copy(resids)
    global errs
    global gaps
    push!(errs, [errs_anderson, errs_1, errs_05, errs_02, errs_01])
    push!(gaps, gap)
end
# @save "res.jld" as errs gaps

for ia = 1:length(as)
    p = plot(errs[ia][1], label="Anderson", m=true, yaxis=:log, reuse=false)
    plot!(errs[ia][2], label="Damping β=1", m=true, yaxis=:log)
    plot!(errs[ia][3], label="Damping β=0.5", m=true, yaxis=:log)
    plot!(errs[ia][4], label="Damping β=0.2", m=true, yaxis=:log)
    plot!(errs[ia][5], label="Damping β=0.1", m=true, yaxis=:log)
    yaxis!("residual")
    xaxis!("n")
    g = @sprintf "%1.1e" gaps[ia]
    title!("a=$(as[ia])   gap=$g")
    ylims!(1e-10, 10)
    display(p)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ d2ce3d8a-fdc0-11ea-0920-6fd0ca7c09be
using Pkg, DrWatson

# ╔═╡ d2ce824a-fdc0-11ea-1504-a7f6289a0308
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanQuap
	using StatisticalRethinking
end

# ╔═╡ 37a0fd70-fdc0-11ea-3fa2-4968a27d3d35
md"## Clip-05-44s.jl"

# ╔═╡ d2cf0594-fdc0-11ea-37c8-9d589bcc4b7b
md"### snippet 5.44"

# ╔═╡ d2de0dbe-fdc0-11ea-3cc5-dd8dc2417134
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame; delim=';');
	df = filter(row -> row[:age] > 18, df)
	scale!(df, [:height, :weight])
end;

# ╔═╡ d2deae40-fdc0-11ea-0e50-0fc462ae5deb
stan5_8 = "
data{
    int N;
    int male[N];
    vector[N] age;
    vector[N] weight;
    vector[N] height;
    int sex[N];
}
parameters{
    vector[2] a;
    real<lower=0> sigma;
}
model{
    vector[N] mu;
    sigma ~ exponential(1);
    a ~ normal( 178 , 20 );
    for ( i in 1:N ) {
        mu[i] = a[sex[i]];
    }
    height ~ normal( mu , sigma );
}
";

# ╔═╡ d2ea2a4a-fdc0-11ea-2dab-27a0b544fb0e
md"### Define the SampleModel, etc."

# ╔═╡ d2eac996-fdc0-11ea-0ab9-09ab8a1bcdca
begin
	df[!, :sex] = [df[i, :male] == 1 ? 2 : 1 for i in 1:size(df, 1)]
	df_m = filter(row -> row[:sex] == 2, df)
	df_f = filter(row -> row[:sex] == 1, df)
	data = (N = size(df, 1), male = df.male, weight = df.weight,
		height = df.height, age = df.age, sex = df.sex)
	q5_8s, m5_8s, o5_8s = stan_quap("m5.8", stan5_8; data);
	post5_8s_df = read_samples(m5_8s, :dataframe)
	PRECIS(post5_8s_df)
end

# ╔═╡ d2f648e8-fdc0-11ea-11a6-8d47151362d8
if !isnothing(m5_8s)
	part5_8s = Particles(post5_8s_df)
end

# ╔═╡ 6a3e68b2-31a0-11eb-0420-cfa73a6e8f35
begin
	quap5_8s_df = sample(q5_8s)
	PRECIS(quap5_8s_df)
end

# ╔═╡ d30451e0-fdc0-11ea-2c62-9f4cfd82f8fd
if !isnothing(m5_8s)
	plot(title="Densities by sex")
	density!(df_m.height, lab="Male")
	density!(df_f.height, lab="Female")
	vline!([mean(part5_8s[Symbol("a.1")])], lab="Female mean estimate")
	vline!([mean(part5_8s[Symbol("a.2")])], lab="Male mean estimate")
	vline!([mean(quap5_8s_df[:, Symbol("a.2")])], lab="Male (quap) mean estimate")
end

# ╔═╡ d304f5a0-fdc0-11ea-2651-4db63a66bd53
md"## End of clip-05-44s.jl"

# ╔═╡ Cell order:
# ╟─37a0fd70-fdc0-11ea-3fa2-4968a27d3d35
# ╠═d2ce3d8a-fdc0-11ea-0920-6fd0ca7c09be
# ╠═d2ce824a-fdc0-11ea-1504-a7f6289a0308
# ╟─d2cf0594-fdc0-11ea-37c8-9d589bcc4b7b
# ╠═d2de0dbe-fdc0-11ea-3cc5-dd8dc2417134
# ╠═d2deae40-fdc0-11ea-0e50-0fc462ae5deb
# ╟─d2ea2a4a-fdc0-11ea-2dab-27a0b544fb0e
# ╠═d2eac996-fdc0-11ea-0ab9-09ab8a1bcdca
# ╠═d2f648e8-fdc0-11ea-11a6-8d47151362d8
# ╠═6a3e68b2-31a0-11eb-0420-cfa73a6e8f35
# ╠═d30451e0-fdc0-11ea-2c62-9f4cfd82f8fd
# ╟─d304f5a0-fdc0-11ea-2651-4db63a66bd53


text:
####################################################################################
############################### UnitPartials #######################################
####################################################################################

# Tuple with partial derivatives for an UnitDual Number
# Note: Each partial derivative will have a different concrete type because
# they are expected to be Quantities.
# N reflects the length of the tuple
struct UnitPartials{N, TT <: Union{Tuple, NamedTuple}}
    partials::TT
end

# Check if the an UnitPartials type contains a named tuple
isnamed(::Type{UnitPartials{N,TT}}) where {N,TT} = length(TT.parameters) > 0 && TT.parameters[1] isa Tuple

# Constructors
UnitPartials(t::T) where T <: Union{Tuple, NamedTuple} = UnitPartials{length(t), typeof(t)}(t)
UnitPartials(args...) = UnitPartials(args)
UnitPartials(;kwargs...) = UnitPartials(values(kwargs))

# Pretty printing -> includes names for named tuples
function show(io::IO, up::T) where {T <: UnitPartials}
    if isnamed(T)
        names = keys(up.partials)
        vals = values(up.partials)
        pairs = string.(names) .* " = " .* string.(vals)
        print(io, '{', join(pairs, ", "), '}')
    else
        print(io, '{', join(up.partials, ", "), '}')
    end
end

# Forwarding getter methods to the underlying tuple
function getindex(up::UnitPartials, i::Integer)
    @boundscheck i > length(up.partials) && throw(BoundsError(up, i))
    up.partials[i]
end
function getindex(up::T, n::Symbol) where {T <: UnitPartials}
    !isnamed(T) && throw(DomainError(n, "Attempt to access named partial from an UnitPartials object that does not have named partials"))
    getproperty(up.partials, n)
end
function getproperty(up::T, n::Symbol) where {T <: UnitPartials}
    if n == :partials
        getfield(up, :partials)
    else
        !isnamed(T) && throw(DomainError(up, "Attempt to access named partial from an UnitPartials object that does not have named partials"))
        getfield(up.partials, n)
    end
end

# Iterator interface
length(up::UnitPartials) = length(up.partials)
iterate(up::UnitPartials) = iterate(up.partials)
iterate(up::UnitPartials, state) = iterate(up.partials, state)

# Zero UnitPartials (by type or instance)
zero(up::T) where {T <: UnitPartials} = zero(T)

zero(::Type{UnitPartials{N,TT}}) where {N,TT <: Tuple} = UnitPartials((zero(p) for p in TT.parameters)...)

function zero(::Type{UnitPartials{N,TT}}) where {N,TT <: NamedTuple}
    z = Tuple((zero(p) for p in TT.parameters[2].parameters))
    UnitPartials(NamedTupleTools.namedtuple(TT.parameters[1], z))
end

####################################################################################
################################# UnitDual #########################################
####################################################################################

# An UnitDual is composed of a value and a tuple of partial derivatives
struct UnitDual{DT <: Number, PT <: UnitPartials}
    value::DT
    partials::PT
end

# Constructors
UnitDual(val, partials::Tuple) = UnitDual(val, UnitPartials(partials))
UnitDual(val, args...) = UnitDual(val, UnitPartials(args))
UnitDual(val; kwargs...) = UnitDual(val, UnitPartials(; kwargs...))

# Setters and getters
value(x::UnitDual) = x.value
partials(x::UnitDual) = x.partials

# Convenience functions
zero(x::T) where {T <: UnitDual} = zero(T)
zero(::Type{UnitDual{DT, PT}}) where {DT, PT} = UnitDual(zero(DT), zero(PT))

# Pretty printing
function show(io::IO, x::UnitDual)
    print(io, "$(value(x)) ")
    show(io, partials(x))
end

# Iterator interface
length(x::UnitDual) = 1
iterate(x::UnitDual) = (x, nothing)
iterate(x::UnitDual, state::Nothing) = nothing


text:
###############################################################################
#
# Update canopy fluxes
#
###############################################################################
"""
    canopy_fluxes!(
                can::Canopy4RT{FT},
                can_opt::CanopyOpticals{FT},
                can_rad::CanopyRads{FT},
                in_rad::IncomingRadiation{FT},
                soil::SoilOpticals{FT},
                leaves::Array{LeafBios{FT},1},
                wls::WaveLengths{FT},
                rt_con::RTCache{FT}
    ) where {FT<:AbstractFloat}

Computes a variety of integrated fluxes from the spectrally resolved
    computations in the short-wave Canopy RT (e.g. absorbed soil radiation,
    absorbed direct and diffuse PAR by layer (and angles for direct), net
    direct and diffuse energy balance per layer), given
- `can` [`Canopy4RT`](@ref) type struct
- `can_opt` [`CanopyOpticals`](@ref) type struct
- `can_rad` [`CanopyRads`](@ref) type struct
- `in_rad` [`IncomingRadiation`](@ref) type struct
- `soil` [`SoilOpticals`](@ref) type struct
- `leaves` Array of [`LeafBios`](@ref) type struct
- `wls` [`WaveLengths`](@ref) type struct
- `rt_con` [`RTCache`](@ref) type cache
"""
function canopy_fluxes!(
            can::Canopy4RT{FT},
            can_opt::CanopyOpticals{FT},
            can_rad::CanopyRads{FT},
            in_rad::IncomingRadiation{FT},
            soil::SoilOpticals{FT},
            leaves::Array{LeafBios{FT},1},
            wls::WaveLengths{FT},
            rt_con::RTCache{FT}
) where {FT<:AbstractFloat}
    # 1. unpack variables from structures
    @unpack LAI, nLayer, Ω = can;
    @unpack ρ_SW, ε_SW = soil;
    @unpack dWL, dWL_iPAR, iPAR, WL, WL_iPAR = wls;
    cf_con = rt_con.cf_con;

    # 2. compute some useful variables
    tLAI = LAI / nLayer;
    fac  = FT(1e-3);

    # 3. Compute some fluxes, can be done separately if needed
    #    this is absolute fluxes now, for the entire soil
    last_ind_cr            = lastindex(can_rad.E_down,2);
    cf_con.abs_wave       .= view(can_rad.E_down, :, last_ind_cr) .* ε_SW;
    can_rad.RnSoil_diffuse = fac * numerical∫(cf_con.abs_wave, dWL);
    cf_con.abs_wave       .= view(can_opt.Es_, :, last_ind_cr) .* ε_SW;
    can_rad.RnSoil_direct  = fac * numerical∫(cf_con.abs_wave, dWL);
    can_rad.RnSoil         = can_rad.RnSoil_direct + can_rad.RnSoil_diffuse;

    # 4. Normalization factor for leaf direct PAR
    #    weighted sum has to be 1 to conserve net SW direct
    #    Direct PAR is normalized by layer Ps value
    mul!(cf_con.absfs_lidf, adjoint(can_opt.absfs), can.lidf);
    normi       = 1 / mean(cf_con.absfs_lidf);
    cf_con.lPs .= (view(can_opt.Ps, 1:nLayer  ) .+
                   view(can_opt.Ps, 2:nLayer+1)) ./ 2;
    @unpack lPs = cf_con;

    @inbounds for j in 1:nLayer
        if length(leaves)>1
            cf_con.kChlrel .= view(leaves[j].kChlrel, iPAR);
        else
            cf_con.kChlrel .= view(leaves[1].kChlrel, iPAR);
        end

        # for diffuse PAR
        cf_con.E_iPAR .= view(can_rad.netSW_shade, iPAR, j);
        e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_diff);
        cf_con.PAR_diff .*= fac / tLAI;

        # for direct PAR
        cf_con.E_iPAR .= view(can_rad.netSW_sunlit, iPAR, j);
        e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_dir);
        cf_con.PAR_dir .*= fac / tLAI;

        # for leaf absorbed
        cf_con.PAR_diffCab .= cf_con.kChlrel .* cf_con.PAR_diff;
        cf_con.PAR_dirCab  .= cf_con.kChlrel .* cf_con.PAR_dir;

        # Absorbed PAR per leaf for shaded, PAR_diff and PAR_dir changed!
        _dif = numerical∫(cf_con.PAR_diff, dWL_iPAR);
        _dir = numerical∫(cf_con.PAR_dir, dWL_iPAR) * normi;
        can_rad.absPAR_shade[j] = _dif;
        can_rad.absPAR_sun[:,:,j] .= can_opt.absfs .* _dir;

        # absorbed PAR for photosynthesis
        _difCab = numerical∫(cf_con.PAR_diffCab, dWL_iPAR);
        _dirCab = numerical∫(cf_con.PAR_dirCab , dWL_iPAR) * normi;
        can_rad.absPAR_shadeCab[j] = _difCab;
        can_rad.absPAR_sunCab[:,:,j]  .= can_opt.absfs .* _dirCab;
        can_rad.absPAR_sunCab[:,:,j] .+= _difCab;
    end

    # 5. Total PAR
    # TODO considering remove this part, if we are not using it
    #    re-use the PAR_dir and PAR_diff in the rt_con
    cf_con.E_iPAR .= view(in_rad.E_direct, iPAR);
    e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_dir);
    cf_con.E_iPAR .= view(in_rad.E_diffuse, iPAR);
    e2phot!(WL_iPAR, cf_con.E_iPAR, cf_con.PAR_diff);
    can_rad.incomingPAR_direct  = fac * numerical∫(cf_con.PAR_dir , dWL_iPAR);
    can_rad.incomingPAR_diffuse = fac * numerical∫(cf_con.PAR_diff, dWL_iPAR);
    can_rad.incomingPAR         = can_rad.incomingPAR_diffuse +
                                  can_rad.incomingPAR_direct;
    @inbounds for i in 1:nLayer
        cf_con.E_all .= view(can_rad.netSW_shade, :, i);
        can_rad.intNetSW_shade[i]  = numerical∫(cf_con.E_all, dWL) * fac / tLAI;
        cf_con.E_all .= view(can_rad.netSW_sunlit, :, i);
        can_rad.intNetSW_sunlit[i] = numerical∫(cf_con.E_all, dWL) *
                                     fac / tLAI / lPs[i] +
                                     can_rad.intNetSW_shade[i];
    end

    return nothing
end

max_stars_repo_name: climate-machine/Land


text:
# Autogenerated wrapper script for Fontconfig_jll for i686-w64-mingw32
export fc_cache, fc_cat, fc_conflist, fc_list, fc_match, fc_pattern, fc_query, fc_scan, fc_validate, fonts_conf, libfontconfig

using FreeType2_jll
using Bzip2_jll
using Zlib_jll
using Libuuid_jll
using Expat_jll
JLLWrappers.@generate_wrapper_header("Fontconfig")
JLLWrappers.@declare_executable_product(fc_cache)
JLLWrappers.@declare_executable_product(fc_cat)
JLLWrappers.@declare_executable_product(fc_conflist)
JLLWrappers.@declare_executable_product(fc_list)
JLLWrappers.@declare_executable_product(fc_match)
JLLWrappers.@declare_executable_product(fc_pattern)
JLLWrappers.@declare_executable_product(fc_query)
JLLWrappers.@declare_executable_product(fc_scan)
JLLWrappers.@declare_executable_product(fc_validate)
JLLWrappers.@declare_file_product(fonts_conf)
JLLWrappers.@declare_library_product(libfontconfig, "libfontconfig-1.dll")
function __init__()
    JLLWrappers.@generate_init_header(FreeType2_jll, Bzip2_jll, Zlib_jll, Libuuid_jll, Expat_jll)
    JLLWrappers.@init_executable_product(
        fc_cache,
        "bin\\fc-cache.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_cat,
        "bin\\fc-cat.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_conflist,
        "bin\\fc-conflist.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_list,
        "bin\\fc-list.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_match,
        "bin\\fc-match.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_pattern,
        "bin\\fc-pattern.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_query,
        "bin\\fc-query.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_scan,
        "bin\\fc-scan.exe",
    )

    JLLWrappers.@init_executable_product(
        fc_validate,
        "bin\\fc-validate.exe",
    )

    JLLWrappers.@init_file_product(
        fonts_conf,
        "etc\\fonts\\fonts.conf",
    )

    JLLWrappers.@init_library_product(
        libfontconfig,
        "bin\\libfontconfig-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
    ENV["FONTCONFIG_FILE"] = get(ENV, "FONTCONFIG_FILE", fonts_conf)
    ENV["FONTCONFIG_PATH"] = get(ENV, "FONTCONFIG_PATH", dirname(ENV["FONTCONFIG_FILE"]))

end  # __init__()


text:
"""
get the number of photons collected, uses a normal distribution
Arguments:  ebc, Float64, excited buffered calcium molecules
            release_rate, Float64, expected proportion ebc molecules that release a photon
            capture_rate, Float64, bernoulli probability of capturing each individual photon
Returns:    number_of_photons, Int64
"""
function getPhotonsCollected(ebc::Number, release_rate::Number, capture_rate::Number)
  0 <= capture_rate <= 1 || return 0
  released = floor(Int, ebc*release_rate)
  released <= 0 && return 0
  photon_dist = Binomial(released, capture_rate)
  return rand(photon_dist)
end


text:
import Porta
import Test

samples = 360
segments = 30
radius = 0.01
basepoints = convert(Array{Complex}, fill(2im, samples))
fiberactions = [fill(0, samples) fill(2pi, samples)]
q = Porta.ℍ([cos(0); sin(0) .* [sqrt(3)/3; sqrt(3)/3; sqrt(3)/3]])
p = [0.0; 0.0; 0.0]
h = Porta.⭕(basepoints, fiberactions, segments, radius, q, p)
segments2 = Integer(segments ÷ 3)
Test.@test size(h.m) == (samples, segments, segments2, 3)
Test.@test size(h.c) == size(h.m)


text:
# The infite loop was patched with a killswitch after certain number of iterations.
# Now the obtained reconstructions are too good to be real and I need to check what
push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
using SparseInverseProblems
using SuperResModels
using TestCases
using Utils

using PyCall
@pyimport numpy as np

type StaticParameters
    # allowed frequencies
    f_c::Int64
    # spatial max val.
    x_max::Float64
    # ???
    num_x::Int64
    # Weights on the allowed frequencies.
    filter_x::Array{Float64}
end
type DynamicParameters
    # 2K+1 is the number of time samples.
    K::Int64
    # sampling rate
    tau::Float64
    # maxium allowed speed
    v_max::Float64
    # ????
    num_v::Int64
end

### Static parameters and container for use.
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x =  ones(2*f_c+1);
static_parameters = StaticParameters(f_c, x_max, num_x, filter_x)

model_static = SuperResModels.Fourier1d(static_parameters.f_c, static_parameters.x_max, static_parameters.filter_x, static_parameters.num_x)
### Dynamic parameters and container for use.
K = 2
tau = 1.0/(K*1.0)
v_max = 0.05
num_v = 10
dynamic_parameters = DynamicParameters(K, tau, v_max, num_v)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, dynamic_parameters.v_max, dynamic_parameters.tau, dynamic_parameters.K, dynamic_parameters.num_v*dynamic_parameters.K)
### Variabiliy of weights
minWeights = 1
maxWeights = 10
### Location of data folder
dataFolder = "data/1Dsimulations"

dynamic_parameters = DynamicParameters(K, tau, v_max, num_v)
bins = np.load("../"*dataFolder*"/separationDistribBins.npy")
density = np.load("../"*dataFolder*"/separationDistribVal.npy")

### Number of points
@everywhere begin
    model_static = $model_static
    model_dynamic = $model_dynamic
    x_max = $x_max
    v_max = $v_max
    minWeights = $minWeights
    maxWeights = $maxWeights
    K = $K
    tau = $tau
    dataFolder = $dataFolder
    test_case = () -> TestCases.cloud_1d_full(x_max, v_max, minWeights, maxWeights, K, tau, rand(4:10))
    noises_data = linspace(0, 0.1, 5)
    noises_data = noises_data[2:end]
    noises_position = linspace(0, 0.01, 5)
    noises_position = noises_position[2:end]
    bins = $bins
    density = $density
end

num_trials = 1
srand(100000000)
println("### Beggining round ###")
    K = div(length(model_dynamic.times)-1,2)
    tau = maximum(model_dynamic.times)/K
    (thetas, weights) =  Utils.Rejection_sampling(test_case, bins, density, K, tau, model_static.x_max)
    d = dim(model_dynamic)
    results = zeros(1, 7)
    # Obtain the separation of the configuration of particles
    separation = Utils.separation_val(thetas,K,tau,model_static.x_max)
    # And the separationg given by the dynamic norm.
    separation_dyn = Utils.separation_norm(thetas,K,tau,model_static.x_max)
    ## 	Noiseless case
    # Dynamic case
    # We generate the target measure
    target = Utils.generate_target(model_dynamic, thetas, weights)
    # recomputes, adding differents levels of noise, for both dynamic and static.
    for noise_data in noises_data[2]
    	target = Utils.generate_target(model_dynamic,thetas,weights, noise_data, 0.0)
    end
    target = Utils.generate_target(model_dynamic,thetas,weights, noises_data[3], 0.0)
    println("---- Printing setup: ")
    println("thetas:",thetas)
    println("weights:",weights)
    println("target:",target)
    # Dynamic case
    println("#### going for noise_data dynamic ",noises_data[3]," ###")
    (thetas_est, weights_est) =  Utils.run_simulation_target(model_dynamic, thetas, weights, target)

println(" Simulation finished with results: ", thetas_est)


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
Module for constructing self-contained power system objects.
"""
module PowerSystems

#################################################################################
# Exports

export System
export Bus
export Arc
export LoadZones

export PowerSystemType
export Component
export Device
export Branch
export Injection
export ACBranch
export Line
export MonitoredLine
export DCBranch
export HVDCLine
export VSCDCLine
export Transformer2W
export TapTransformer
export PhaseShiftingTransformer

export ThreePartCost
export TwoPartCost

export Generator
export HydroGen
export HydroFix
export HydroDispatch
export HydroStorage
export TechHydro

export RenewableGen
export TechRenewable
export RenewableFix
export RenewableDispatch

export ThermalGen
export TechThermal
export ThermalStandard

export ElectricLoad
export StaticLoad
export PowerLoad
export PowerLoadPF
export FixedAdmittance
export ControllableLoad
export InterruptibleLoad

export Storage
export GenericBattery

export Service
export Reserve
export ProportionalReserve
export StaticReserve
export Transfer

export PTDF
export Ybus
export LODF
export GeneratorCostModel
export BusType

export Forecast
export Deterministic
export Probabilistic
export ScenarioBased

export make_pf
export solve_powerflow!

export BevDemand
export ChargingSegment
export ChargingPlan
export Demand
export Envelope
export FlexibleDemand
export InflexibleDemand
export InterruptibleDemand
export LocatedDemand
export LocatedEnvelope
export MobileDemand
export StationaryInflexibleDemand
export TemporalDemand
export aligntimes
export batterylevels
export chargeamounts
export chargelevels
export chargerates
export consumptionamounts
export consumptionrates
export demands
export durations
export earliestdemands
export envelope
export greedydemands
export latestdemands
export loads
export locateddemand
export locations
export maxchargerates
export populate_BEV_demand
export shortfall
export verify
export verifybattery
export verifylimits

export parse_standard_files
export parse_file
export add_forecasts!
export add_forecast!
export remove_forecast!
export clear_forecasts!
export add_component!
export remove_component!
export remove_components!
export get_component
export get_components
export get_components_by_name
export get_component_forecasts
export get_forecast_initial_times
export get_forecasts
export get_forecasts_horizon
export get_forecasts_initial_time
export get_forecasts_interval
export get_forecasts_resolution
export get_forecast_component_name
export get_forecast_value
export get_horizon
export get_timeseries
export get_data
export iterate_components
export iterate_forecasts
export make_forecasts
export split_forecasts!
export get_name
export to_json

#################################################################################
# Imports

import SparseArrays
import LinearAlgebra: LAPACK.getri!
import LinearAlgebra: LAPACK.getrf!
import LinearAlgebra: BLAS.gemm
import LinearAlgebra
import Dates
import TimeSeries
import DataFrames
import JSON
import JSON2
import CSV
import YAML
import UUIDs
import Base.to_index

import InfrastructureSystems
import InfrastructureSystems: Components, Deterministic, Probabilistic, Forecast,
    ScenarioBased, InfrastructureSystemsType, InfrastructureSystemsInternal,
    FlattenIteratorWrapper, LazyDictFromIterator, DataFormatError, InvalidRange,
    InvalidValue

const IS = InfrastructureSystems

#################################################################################
# Includes

"""
Supertype for all PowerSystems types.
All subtypes must include a InfrastructureSystemsInternal member.
Subtypes should call InfrastructureSystemsInternal() by default, but also must
provide a constructor that allows existing values to be deserialized.
"""
abstract type PowerSystemType <: IS.InfrastructureSystemsType end

abstract type Component <: PowerSystemType end
# supertype for "devices" (bus, line, etc.)
abstract type Device <: Component end
abstract type Injection <: Device end
# supertype for generation technologies (thermal, renewable, etc.)
abstract type TechnicalParams <: PowerSystemType end

include("common.jl")

# Include utilities
include("utils/IO/base_checks.jl")
include("utils/timearray.jl")

# PowerSystems models
include("models/topological_elements.jl")
include("models/branches.jl")
include("models/operational_cost.jl")
#include("models/network.jl")

# Static types
include("models/generation.jl")
include("models/storage.jl")
include("models/loads.jl")
include("models/services.jl")
include("models/demand.jl")

# Include all auto-generated structs.
include("models/generated/includes.jl")
include("models/supplemental_constructors.jl")

# Definitions of PowerSystem
include("base.jl")

#Interfacing with Forecasts
include("forecasts.jl")

#Data Checks
include("utils/IO/system_checks.jl")
include("utils/IO/branchdata_checks.jl")

# network calculations
include("utils/network_calculations/common.jl")
include("utils/network_calculations/ybus_calculations.jl")
include("utils/network_calculations/ptdf_calculations.jl")
include("utils/network_calculations/lodf_calculations.jl")

#PowerFlow
include("utils/power_flow/make_pf.jl")
include("utils/power_flow/power_flow.jl")

# Include Parsing files
include("parsers/common.jl")
include("parsers/enums.jl")
include("parsers/pm_io.jl")
include("parsers/im_io.jl")
include("parsers/standardfiles_parser.jl")
include("parsers/forecast_parser.jl")
include("parsers/power_system_table_data.jl")
include("parsers/pm2ps_parser.jl")

# Better printing
include("utils/print.jl")

include("models/serialization.jl")

# Download test data
include("utils/data.jl")
import .UtilsData: TestData

end # module


text:
using Test
using RealNeuralNetworks.Utils.BoundingBoxes

function create_fake_bounding_box()
    bbox = BoundingBox((1,1,1), (204,204,300))
end

@testset "test bounding box" begin
    bbox = BoundingBox((1,2,3), (2,3,4))
    d = BoundingBoxes.distance_from(bbox, (3,4,5))
    @test d ≈ sqrt(3)
    range = BoundingBoxes.get_unit_range(bbox)
    @test range == (1:2, 2:3, 3:4)
    @test size(bbox) == (2,2,2)

    bbox2 = BoundingBox((2,3,4), (3,4,5))
    bbox3 = union(bbox, bbox2)
    @show BoundingBox((1,2,3), (3,4,5))
    @show bbox3
    @test isequal(bbox3, BoundingBox((1,2,3), (3,4,5)))
    @test bbox3 == BoundingBox((1,2,3), (3,4,5))
end


text:
# Needs to be tested
function runtests_gpu(name)
    return "rm(joinpath(@__DIR__, \"..\", \"Manifest.toml\"), force = true) # Remove local Manifest.toml

test_env = joinpath(@__DIR__, \"..\", \"env\", \"test\")
push!(LOAD_PATH, test_env)
push!(LOAD_PATH, joinpath(@__DIR__, \"..\", \"env\", \"gpu\"))
push!(LOAD_PATH, joinpath(@__DIR__, \"..\"))

using Pkg
Pkg.activate(test_env)
Pkg.instantiate(; verbose = true)

using Test
using $name

@testset \"Test $name\" begin
    @test 1==1
end

"

end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
module BitsTallies

import Base: getindex, +, -

const BITS = 5
const mask = UInt8(2 ^ BITS - 1)
const offsets = [BITS * i for i in 0:25]
const charmasks = [UInt128(1) << x for x in offsets]

"""
The BitsTally uses a single unsigned integer to store counts
for every letter in 'a':'z'. This lets us compute those tallies
without allocating an array, and compare them using simple
integer equality.

The way we do this is to divide up the internal 128-bit integer into 26
blocks. Each block is 5 bits long, except the final block which is only 3 bits
long. We can use each of those blocks to store a 5-bit count of a given
letter. We store the counts for the letter 'a' in the low-order bits, which
leaves the awkward 3-bit block for 'z'. So this will overflow if a word
contains more than 7 'z's. or more than 31 of any other letter.
"""
struct BitsTally
    data::UInt128

    BitsTally(data=0) = new(data)
end

"""
Add a single count for char c to the tally t.
"""
+(t::BitsTally, c::Char) = BitsTally(t.data + charmasks[convert(Int, c - 'a' + 1)])

-(t::BitsTally, c::Char) = BitsTally(t.data - charmasks[convert(Int, c - 'a' + 1)])

function BitsTally(s::String)
    t = BitsTally()
    for c in s
        t += c
    end
    t
end

"""
Strip the appropriate bits out of the storage UInt128 to return the count of
the given letter as a UInt8.
"""
getindex(t::BitsTally, c::Char) = mod((t.data >> offsets[convert(Int, c - 'a' + 1)]) & mask, UInt8)

isanagram(t1::BitsTally, t2::BitsTally) = t1.data == t2.data

"""
Returns true iff t1 is a transaddition of t2. That is, the letters of t2
can be rearranged along with one extra letter to make the letters of t1.
"""
function istransaddition(t1::BitsTally, t2::BitsTally)
    d = t1.data - t2.data
    for m in charmasks
        if d == m
            return true
        end
    end
    false
end

end


text:
export MullerBrown

function Base.getproperty(model::Model, name::Symbol, x)
    if length(String(name)) > 7 && String(name)[end-6:end] == "_kwargs"
        name = Symbol(String(name)[1:end-7])
        field = getfield(model, name)
        x = x == nothing ? Dict([(f, nothing) for f in fieldnames(field)]) : x
        ans = Dict()
        for (key, value) in pairs(x)
            ans[key] = getproperty(field, key, value)
        end
        return ans
    elseif parentmodule(fieldtype(typeof(model), name)) ∉ [Base, Core]
        typename = typeof(getfield(model, name))
        return "$typename"
    else
        return getfield(model, name)
    end
end

@inline get_properties(model::Model, paths::Paths, properties::String, args...; kwargs...) = model(paths, paths.coords, [properties], args...; kwargs...)
@inline get_properties(model::Model, paths::Paths, properties::Array{Any, 1}, args...; kwargs...) = model(paths, paths.coords, properties, args...; kwargs...)
@inline get_properties(model::Model, paths::Paths, properties::String, coords::Coords, args...; kwargs...) = model(paths, coords, [properties], args...; kwargs...)
@inline get_properties(model::Model, paths::Paths, properties::Array{Any, 1}, coords::Coords, args...; kwargs...) = model(paths, coords, properties, args...; kwargs...)

@inline get_kinetics(model::Model, paths::Paths, properties::String, args...; kwargs...) = paths.coords(paths, paths.coords, [properties], args...; kwargs...)
@inline get_kinetics(model::Model, paths::Paths, properties::Array{Any, 1}, args...; kwargs...) = paths.coords(paths, paths.coords, properties, args...; kwargs...)
@inline get_kinetics(model::Model, paths::Paths, properties::String, coords::Coords, args...; kwargs...) = paths.coords(paths, coords, [properties], args...; kwargs...)
@inline get_kinetics(model::Model, paths::Paths, properties::Array{Any, 1}, coords::Coords, args...; kwargs...) = paths.coords(paths, coords, properties, args...; kwargs...)

@inline get_distances(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "displacements", args...; kwargs...)
@inline get_momentums(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "momoentum", args...; kwargs...)
@inline get_kinetic_energies(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "kinetic_energy", args...; kwargs...)
@inline get_kinetic_energy_gradients(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "kinetic_energy_gradient", args...; kwargs...)
@inline get_velocities(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "velocity", args...; kwargs...)
@inline get_accelerations(model::Model, paths::Paths, args...; kwargs...) = get_kinetics(model, paths, "acceleration", args...; kwargs...)

function get_masses(model::Model, paths::Paths, args...; kwargs...)
    get_properties(model, paths, "mass", args...; kwargs...)
end
function get_effective_mass(paths::Paths, args...; kwargs...)
    get_effective_mass(paths.model, args...; kwargs...)
end

""" Calculate potential( energy) """
@inline get_potential(model::Model, paths::Paths, args...; kwargs...) = get_properties(model, paths, "potential", args...; kwargs...)
@inline get_potentials(model::Model, paths::Paths, args...; kwargs...) = get_properties(model, paths, "potentials", args...; kwargs...)
@inline get_gradients(model::Model, paths::Paths, args...; kwargs...) = get_properties(model, paths, "gradients", args...; kwargs...)
@inline get_hessian(model::Model, paths::Paths, args...; kwargs...) = get_properties(model, paths, "hessain", args...; kwargs...)
""" Calculate covariance. It only applies when potential is guassian"""
@inline get_covariance(model::Model, paths::Paths, args...; kwargs...) = get_properties(model, paths, "covariance", args...; kwargs...)

struct MullerBrown <: Model
    A::Array{Float64, 1}
    a::Array{Float64, 1}
    b::Array{Float64, 1}
    c::Array{Float64, 1}
    x0::Array{Float64, 1}
    y0::Array{Float64, 1}
    results::Dict
end

MullerBrown(;
   A=[-2., -1., -1.7, 0.15], a=[-1, -1, -6.5, 0.7], b=[0., 0., 11., 0.6],
   c=[-10, -10, -6.5, 0.7], x0=[1, 0, -0.5, -1], y0=[0, 0.5, 1.5, 1],
   results=Dict()) = MullerBrown(A, a, b, c, x0, y0, results)

@inline (model::MullerBrown)(paths::Paths, coords::Coords, properties::Array{Any, 1}) = (model::MullerBrown)(coords::Coords, properties::Array{Any, 1})
@inline (model::MullerBrown)(coords::Coords, properties::Array{Any, 1}) = (model::MullerBrown)(convert(Cartesian{eltype(coords), 2}, coords), properties::Array{Any, 1})

"""
coords : Cartesian; N x 2
"""
function (model::MullerBrown)(coords::Cartesian{T, 2}, properties::Array{Any, 1}) where {T<:Number}
    N, D = size(coords)
    results = Dict()

    # 4 -> 1x4
    A, a, b, c = model.A[nax, :], model.a[nax, :], model.b[nax, :], model.c[nax, :]
    x0, y0 = model.x0[nax, :], model.y0[nax, :]
    # for i=1:N
    𝑥 = coords[:, 1]
    𝑦 = coords[:, 2]
    # Vk = zeros(N, 4)

    𝑥_x0, 𝑦_y0 = 𝑥 .- x0, 𝑦 .- y0
    Vk = @. A * exp(a*(𝑥 - x0)^2 + b*(𝑥-x0)*(𝑦-y0) + c * (𝑦-y0)^2)
    # end

    # "potential" in properties ? results["potential"] = sum(Vk) / 100 : nothing
    if "potential" in properties
        results[:potential] = sum(Vk, dims=2)
    end

    if "gradients" in properties
        dV = zeros(Float64, 4)

        Fx = dropdims(sum((@. Vk * (2 * a * 𝑥_x0 + b * 𝑦_y0)), dims=2), dims=2)
        Fy = dropdims(sum((@. Vk * (b * 𝑥_x0 + 2 * c * 𝑦_y0)), dims=2), dims=2)
        dV = [Fx, Fy]
        results[:gradients] = dV
    end

    if "hessian" in properties
        nothing
    end

    if "mass" in properties
        results[:mass] = ones(N)
    end

    merge!(model.results, results)

    return results
end

include("./ase.jl")

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
module PkgBenchmarks

using BenchmarkTools
using Base.LibGit2

export benchmark

include("pkg.jl")
include("libgit2.jl")
include("report.jl")
include("travis-ci.jl")

end # module


text:
function collectfails(call::Call)
   vcat(call.parsingfails, collectfails(call.args))
end

function collectfails(ellist::ElementList)
   ellist.parsingfails
end

function collectfails(el)
   Fail[]
end

function collectfails(fail::Fail)
   fail
end

"""
Converts a vector to a vector of the most specific type that all
elements share as common supertype.
"""
function converttomostspecifictype(v::Vector)
   if isempty(v)
      return v
   else
      Vector{mostspecifictype(v)}(v)
   end
end

"""
    mostspecifictype(v)
Returns the most specific supertype for all elements in the vector `v`.
"""
function mostspecifictype(v::Vector)
   mapreduce(typeof, typejoin, v)
end

function parsingcheck(el)
   fails = collectfails(el)
   if !isempty(fails)
      error(fails[1].message)
   end
   nothing
end

function evaluate_checked!(call::Call, communicator::CommunicatoR)
   # check parsing first
   try
      parsingcheck(call)
   catch ex
      return Fail("Parsing in Julia failed.", ex)
   end

   try
      return evaluate!(call, communicator)
   catch ex
      return Fail("Evaluation in Julia failed.", ex)
   end
end

function evaluate!(call::Call, communicator::CommunicatoR)
   # evaluate arguments
   evaluate!(call.args, communicator)

   # the actual function call
   Base.invokelatest(call.fun,
         call.args.positionalelements...; call.args.namedelements...)
end

function evaluate!(list::ElementList, communicator::CommunicatoR)
   # recursively evaluate sub-elements first
   for i in eachindex(list.positionalelements)
      list.positionalelements[i] =
            evaluate!(list.positionalelements[i], communicator)
   end
   for name in list.names
      list.namedelements[name] =
            evaluate!(list.namedelements[name], communicator)
   end

   jltype = get(list.attributes, "JLTYPE", "")
   if !isempty(jltype)
      constructor = findfield(jltype)
      if length(list.names) > 0
         if constructor <: AbstractDict
            return Base.invokelatest(constructor, zip(
                  list.namedelements[:keys],
                  list.namedelements[:values]))
         elseif constructor <: NamedTuple
            return NamedTuple{Tuple(list.names)}(
                  [list.namedelements[name] for name in list.names])
         elseif constructor <: Module
            return Module(Symbol(list.namedelements[:name]))
         elseif constructor <: CircularReference
            error("Circular references cannot be reconstructed")
         else
            return reconstruct_struct(constructor, list, communicator)
         end
      else
         if constructor <: Tuple
            return Base.invokelatest(constructor, list.positionalelements)
         else
            # array type (or AbstractSet)
            return reconstruct_array(constructor, list.positionalelements,
                  list.attributes)
         end
      end
   elseif get(list.attributes, "IS_DF", false)
      return RDataFrame(list)
   elseif isempty(list.namedelements)
      return converttomostspecifictype(list.positionalelements)
   else
      return list
   end
end

function evaluate!(objref::ObjectReference, communicator::CommunicatoR)
   getobj(obj::ImmutableObjectReference) = obj.obj
   getobj(obj) = obj

   try
      obj = communicator.sharedheap[objref.ref].obj
      return getobj(obj)
   catch ex
      error("Object reference " * string(objref.ref, base = 16) *
            " could not be resolved")
   end
end

function evaluate!(item, communicator::CommunicatoR)
   item
end

""" Evaluates Julia code, passed as String, in the Main module """
function maineval(str::AbstractString)
   Main.eval(Meta.parse(str))
end

function mainevalcmd(str::String)
   if isempty(str)
      return nothing
   end
   strippedstr = strip(str)
   ret = include_string(Main, str)
   if strippedstr[end] == ';'
      return nothing
   else
      return ret
   end
end

"""
    mainevallet(str; ...)
Evalutes the string as the body of a let expression with the given
keyword arguments as variables in the let expression.
"""
function mainevallet(str::String; kwargs...)
   # Enclose str with a let block containing the kwargs

   # prevent the creation of a toplevel-Expression in Meta.parse
   str = "begin " * str * " end"

   if !isempty(kwargs)
      e = Expr(:let,
            Expr(:block,
                  [Expr(:(=), key, value) for (key, value) in kwargs]...),
            Expr(:block, Meta.parse(str)))
   else
      e = Expr(:let, Expr(:block), Expr(:block, Meta.parse(str)))
   end
   Main.eval(e)
end

function reconstruct_array(constructor, arr::Vector, attributes)
   if haskey(attributes, "JLDIM")
      return reconstruct_multidimensional_array(constructor,
            arr, attributes["JLDIM"])
   elseif any(isequal(JLUNDEF()), arr)
      ret = Base.invokelatest(constructor, undef, size(arr))
      copy_defined!(ret, arr)
      return ret
   else
      return Base.invokelatest(constructor, arr)
   end
end

function reconstruct_multidimensional_array(constructor, arr::Vector,
         dims::Vector{Int})

   if prod(dims) != length(arr)
      error("Incorrect dimensions of array. Dimensions are" *
            " $dims but number of elements is $(length(arr)).\n" *
            " If the array has been modified, the attribute \"JLDIM\"" *
            " (and perhaps the attribute \"JLTYPE\") must be set accordingly.")
   else
      ret = Base.invokelatest(constructor, undef, dims...)
      copy_defined!(ret, arr)
      return ret
   end
end

function reconstruct_struct(constructor, list::ElementList, communicator::CommunicatoR)
   # normal composite type:
   # forge call to inner constructor (which may be private)
   constructorargs = [list.namedelements[name] for name in list.names]
   obj = Main.eval(Expr(:new, constructor, [:($arg) for arg in constructorargs]...))

   # The reconstructed object may depend on a pre-existing copy and its resources.
   # This is referenced by the "JLREF" attribute.
   jlref::UInt64 = get(list.attributes, "JLREF", UInt64(0))
   if jlref != 0
      sharedfinalizer(communicator, obj, jlref)
   end

   return obj
end