text:
abstract type AbstractSIMP <: TopOptAlgorithm end

@params struct ReuseStatus
    obj::Bool
    constr
end
Base.any(r::ReuseStatus) = r.obj || any(r.constr)
Base.all(r::ReuseStatus) = r.obj && all(r.constr)
function setreuse!(optimizer, reuse::Bool)
    optimizer.obj.f.reuse = reuse
    setproperty!.(optimizer.constr, :reuse, reuse)
    return getreuse(optimizer)
end
function setreuse!(optimizer, reuse::ReuseStatus)
    optimizer.obj.f.reuse = reuse.obj
    setproperty!.(optimizer.constr, :reuse, reuse.constr)
    return getreuse(optimizer)
end
function getreuse(optimizer)
    obj_reuse = optimizer.obj.reuse
    constr_reuse = getproperty.(optimizer.constr, :reuse)
    return ReuseStatus(obj_reuse, constr_reuse)
end
Functions.getfevals(optimizer) = FunctionEvaluations(optimizer)

# MMA wrapper
include("math_optimizers.jl")

# Basic SIMP
include("basic_simp.jl")

## Continuation SIMP
include("continuation_schemes.jl")
include("continuation_simp.jl")

## Adaptive SIMP
include("polynomials.jl")
include("polynomials2.jl")
include("adaptive_simp.jl")


text:
#=
# A triangular reference element.
# Notation follows that used in Nodal Discontinuous Galerkin Methods
# by Hesthaven and Warburton.
#  https://link.springer.com/book/10.1007/978-0-387-72067-8
#
=#

function build_triangle_refel(refel)
    # refel has already been created, but needs quadrature matrices
    (refel.V, refel.Dr, refel.Ds) = triangle_vandermonds(refel, refel.r);
    refel.invV = inv(refel.V);

    (refel.Vg, DVgr, DVgs) = triangle_vandermonds(refel, refel.g);
    #refel.invVg = inv(refel.Vg);

    refel.Q = refel.Vg*refel.invV;
    refel.Qr = DVgr*refel.invV;
    refel.Qs = DVgs*refel.invV;
    refel.Ddr = refel.Dr*refel.invV;
    refel.Dds = refel.Ds*refel.invV;

    ## Surfaces
    nfaces = 3;
    refel.surf_V = fill!(Array{Array{Float64}}(undef, nfaces), []);
    #refel.surf_gradV = fill!(Array{Array{Float64}}(undef, nfaces), []);
    refel.surf_Vg = fill!(Array{Array{Float64}}(undef, nfaces), []);
    #refel.surf_gradVg = fill!(Array{Array{Float64}}(undef, nfaces), []);

    refel.surf_Q = fill!(Array{Array{Float64}}(undef, nfaces), []);
    refel.surf_Qr = fill!(Array{Array{Float64}}(undef, nfaces), []);
    refel.surf_Qs = fill!(Array{Array{Float64}}(undef, nfaces), []);
    refel.surf_Ddr = fill!(Array{Array{Float64}}(undef, nfaces), []);
    refel.surf_Dds = fill!(Array{Array{Float64}}(undef, nfaces), []);

    for fi=1:nfaces
        (refel.surf_V[fi], surf_gradVr, surf_gradVs) = triangle_vandermonds(refel, refel.surf_r[fi]);
        (refel.surf_Vg[fi], surf_gradVgr, surf_gradVgs) = triangle_vandermonds(refel, refel.surf_g[fi]);

        refel.surf_Q[fi] = refel.surf_V[fi] * refel.invV;
        refel.surf_Qr[fi] = surf_gradVgr * refel.invV;
        refel.surf_Qs[fi] = surf_gradVgs * refel.invV;
        refel.surf_Ddr[fi] = surf_gradVr * refel.invV;
        refel.surf_Dds[fi] = surf_gradVs * refel.invV;
    end

    return refel;
end

function triangle_vandermonds(refel, r)
    Np = refel.Np;      # number of nodal points
    Nrp = size(r,1);    # number of quadrature points may be different
    V = zeros(Nrp, Np);
    gradVr = zeros(Nrp, Np);
    gradVs = zeros(Nrp, Np);

    # Transfer (r,s) to (a,b) coordinates
    a = zeros(Nrp,1);
    for ni=1:Nrp
        if r[ni,2] != 1
            a[ni] = 2*(1+r[ni,1])/(1-r[ni,2])-1;
        else
            a[ni] = -1;
        end
    end
    b = r[:,2];

    # build the Vandermonde and gradVandermond matrix
    sk = 1;
    for i=0:refel.N
        h1 = jacobi_polynomial(a,0,0,i);
        dfa = grad_jacobi_polynomial(a, 0, 0, i);
        for j=0:refel.N - i
            h2 = jacobi_polynomial(b,2*i+1,0,j);
            dgb = grad_jacobi_polynomial(b, 2*i+1,0, j);
            V[:,sk] = 1.4142135623730951 .* h1 .* h2 .* (1 .- b).^i;

            # r-derivative
            # d/dr = da/dr d/da + db/dr d/db = (2/(1-s)) d/da = (2/(1-b)) d/da
            dmodedr = dfa.*h2;
            if i>0
                dmodedr = dmodedr.*((0.5 .* (1 .- b)).^(i-1));
            end
            # s-derivative
            # d/ds = ((1+a)/2)/((1-b)/2) d/da + d/db
            dmodeds = dfa.*(h2.*(0.5 .* (1 .+ a)));
            if i>0
                dmodeds = dmodeds.*((0.5 .* (1 .- b)).^(i-1));
            end

            tmp = dgb.*((0.5 .* (1 .- b)).^i);
            if i>0
                tmp = tmp .- 0.5*i .* h2 .* ((0.5 .* (1 .- b)).^(i-1));
            end
            dmodeds = dmodeds + h1.*tmp;

            # Normalize
            dmodedr = 2^(i+0.5) .* dmodedr;
            dmodeds = 2^(i+0.5) .* dmodeds;

            gradVr[:,sk] = dmodedr;
            gradVs[:,sk] = dmodeds;

            sk = sk+1;
        end
    end

    return (V, gradVr, gradVs);
end

ext: jl


text:
using Documenter, Checkpoints

makedocs(
    modules=[Checkpoints],
    format=Documenter.HTML(
        assets=["assets/invenia.css"],
        prettyurls=(get(ENV, "CI", nothing) == "true")
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
    repo="https://github.com/invenia/Checkpoints.jl/blob/{commit}{path}#L{line}",
    sitename="Checkpoints.jl",
    authors="Invenia Technical Computing Corporation",
)

deploydocs(;
    repo="github.com/invenia/Checkpoints.jl",
)


text:
# This file contains default definitions for TermInterface methods on Julia
# Builtin Expr type.

istree(x::Type{Expr}) = true
exprhead(e::Expr) = e.head

operation(e::Expr) = expr_operation(e, Val{exprhead(e)}())
arguments(e::Expr) =  expr_arguments(e, Val{exprhead(e)}())

# See https://docs.julialang.org/en/v1/devdocs/ast/
expr_operation(e::Expr, ::Union{Val{:call}, Val{:macrocall}}) = e.args[1]
expr_operation(e::Expr, ::Union{Val{:ref}}) = getindex
expr_operation(e::Expr, ::Val{T}) where {T} = T

expr_arguments(e::Expr, ::Union{Val{:call}, Val{:macrocall}}) = e.args[2:end]
expr_arguments(e::Expr, _) = e.args

function similarterm(x::Type{Expr}, head, args, symtype=nothing; metadata=nothing, exprhead=:call)
    expr_similarterm(head, args, Val{exprhead}())
end

function similarterm(x::Expr, head, args, symtype=nothing; metadata=nothing, exprhead=exprhead(x))
    expr_similarterm(head, args, Val{exprhead}())
end

expr_similarterm(head, args, ::Val{:call}) = Expr(:call, head, args...)
expr_similarterm(head, args, ::Val{:macrocall}) = Expr(:macrocall, head, args...) # discard linenumbernodes?
expr_similarterm(head, args, ::Val{eh}) where {eh} = Expr(eh, args...)


text:
"""
    as_added(model::AbstractModel)

Activates agents at each step in the same order as they have been added to the model.
"""
function as_added(model::AbstractModel)
  agent_ids = [i.id for i in 1:length(model.agents)]
  return sortperm(agent_ids)
end

"""
    random_activation(model::AbstractModel)

Activates agents once per step in a random order.
"""
function random_activation(model::AbstractModel)
  order = shuffle(1:length(model.agents))
end

"""
    partial_activation(model::AbstractModel)

At each step, activates only `activation_prob` number of randomly chosen of individuals with a `activation_prob` probability. `activation_prob` should be a field in the model and between 0 and 1.
"""
function partial_activation(model::AbstractModel)
  agentnum = nagents(model)
  return randsubseq(1:agentnum, model.activation_prob)
end

function return_activation_order(model::AbstractModel)
  order = model.scheduler(model)
end


text:
"""
    cb_backdoor_discrete(X, Y, U, M)

    Implements the backdoor causal bootstrapping procedure for discrete confounders.
"""
function cb_backdoor_discrete(
    X::Union{DataFrame, Matrix},
    Y::Union{DataFrame, Vector},
    U::Matrix,
    M::Int)

    y, u  = unique(Y), unique(U)                  # Unique classes and confounders
    N, K, L = size(X,1), size(y,1), size(u,1)     # N samples, classes, confounders

    # STEP 1: Estimate P(Y,U), P(U), P(Y), and P(Y|U)
    C = [sum((Y.==y[i]) .& (U.==u[j])) for i∈1:K, j∈1:L]    # Counts
    if sum(C .== 0) >= 1
        C = C .+ 1
    end
    Pyu = C ./ sum(C)                                       # P(Y,U)
    Py, Pu = sum(Pyu, dims=2), sum(Pyu, dims=1)'            # P(Y) and P(U)
    Py_u = Pyu ./ repeat(Pu', K, 1)                         # P(Y|U)

    # Compute the number of bootstrap samples to draw per class
    Nₛ = Int.(round.(Py .* M, digits=0)); Nₛ[K] = M - sum(Nₛ[1:K-1])

    # STEP 2: For each unique level of Y
    k = 1; indices = zeros(M); Y′ = zeros(M,1);
    for i ∈ 1:K
        # STEP 3: Resample indicies
        actual = (Y .== y[i])
        expected = map(Uᵢ -> N*Py_u[i, findall(x-> x == Uᵢ, u)[1]], U)
        w = (actual ./ expected) |> vec                                 # Weights
        w = w ./ sum(w)                                                 # Ensure normalization
        indices[k:k+Nₛ[i]-1] = rand(Categorical(w), Nₛ[i])               # New indices
        Y′[k:k+Nₛ[i]-1] .= y[i]                                          # New targets
        k = k+Nₛ[i]
    end
    indices = vec(Int.(indices))
    X′ = X[indices,:]; U′ = U[indices,:]
    return X′, Y′, U′
end

"""
    discrete_causal_bootstrapping(X, Y, Z, M)

    Function that wraps the cb_backdoor_discrete causal bootstrapping procedure to coerce the input and
    output datatypes into formats appropriate for the prediction models.
"""
function discrete_causal_bootstrapping(X, Y, Z, M)
    X′, Y′, Z′ = cb_backdoor_discrete(X, coerce(Y, Count) .- 1, Matrix(Z), M)
    return X′, vec(coerce(Y′, Finite)), Z′
end


text:
include("init.jl")

@testset "Broker API tests" begin

    @test h.helicsIsCoreTypeAvailable("zmq") == 1
    broker1 = h.helicsCreateBroker("zmq", "broker1", "--federates 3 --loglevel 1")
    broker2 = h.helicsBrokerClone(broker1)
    address_string = h.helicsBrokerGetAddress(broker1)
    occursin("tcp://127.0.0.1:23404", address_string)
    occursin("broker1", h.helicsBrokerGetIdentifier(broker1))
    @test h.helicsBrokerIsValid(broker1) == 1
    @test h.helicsBrokerIsConnected(broker1) == 1
    h.helicsBrokerDisconnect(broker1)
    @test h.helicsBrokerIsConnected(broker1) == 0
    h.helicsBrokerDisconnect(broker2)
    h.helicsBrokerFree(broker1)
    h.helicsBrokerFree(broker2)
    h.helicsCloseLibrary()

end

@testset "Core API tests" begin

    core1 = h.helicsCreateCore("test", "core1", "--autobroker")
    core2 = h.helicsCoreClone(core1)
    @test occursin("core1", h.helicsCoreGetIdentifier(core1))

    @test h.helicsCoreIsConnected(core1) == 0

    sourceFilter1 = h.helicsCoreRegisterFilter(core1, h.HELICS_FILTER_TYPE_DELAY, "core1SourceFilter")
    h.helicsFilterAddSourceTarget(sourceFilter1, "ep1")
    destinationFilter1 = h.helicsCoreRegisterFilter(core1, h.HELICS_FILTER_TYPE_DELAY, "core1DestinationFilter")
    h.helicsFilterAddDestinationTarget(destinationFilter1, "ep2")
    cloningFilter1 = h.helicsCoreRegisterCloningFilter(core1, "ep3")
    h.helicsFilterRemoveDeliveryEndpoint(cloningFilter1, "ep3")

    h.helicsCoreSetReadyToInit(core1)
    h.helicsCoreDisconnect(core1)
    h.helicsCoreDisconnect(core2)
    h.helicsCoreFree(core1)
    h.helicsCoreFree(core2)
    h.helicsCloseLibrary()

end

@testset "Misc API tests" begin
    fedInfo1 = h.helicsCreateFederateInfo()
    h.helicsFederateInfoSetCoreInitString(fedInfo1, "-f 1")
    h.helicsFederateInfoSetCoreName(fedInfo1, "core3")
    h.helicsFederateInfoSetCoreType(fedInfo1, 3)
    h.helicsFederateInfoSetCoreTypeFromString(fedInfo1, "zmq")
    h.helicsFederateInfoSetFlagOption(fedInfo1, 1, true)
    h.helicsFederateInfoSetTimeProperty(fedInfo1,
            h.HELICS_PROPERTY_TIME_INPUT_DELAY, 1.0)
    h.helicsFederateInfoSetIntegerProperty(fedInfo1,
            h.HELICS_PROPERTY_INT_LOG_LEVEL, 1)
    h.helicsFederateInfoSetIntegerProperty(fedInfo1,
            h.HELICS_PROPERTY_INT_MAX_ITERATIONS, 100)
    h.helicsFederateInfoSetTimeProperty(fedInfo1,
            h.HELICS_PROPERTY_TIME_OUTPUT_DELAY, 1.0)
    h.helicsFederateInfoSetTimeProperty(fedInfo1,
            h.HELICS_PROPERTY_TIME_PERIOD, 1.0)
    h.helicsFederateInfoSetTimeProperty(fedInfo1,
            h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    h.helicsFederateInfoSetTimeProperty(fedInfo1,
            h.HELICS_PROPERTY_TIME_OFFSET, 0.1)
    h.helicsFederateInfoFree(fedInfo1)

    broker3 = h.helicsCreateBroker("zmq", "broker3", "--federates 1 --loglevel 1")
    fedInfo2 = h.helicsCreateFederateInfo()
    coreInitString = "--federates 1"
    h.helicsFederateInfoSetCoreInitString(fedInfo2, coreInitString)
    h.helicsFederateInfoSetCoreTypeFromString(fedInfo2, "zmq")
    h.helicsFederateInfoSetIntegerProperty(fedInfo2, h.HELICS_PROPERTY_INT_LOG_LEVEL, 1)
    h.helicsFederateInfoSetTimeProperty(fedInfo2, h.HELICS_PROPERTY_TIME_DELTA, 1.0)
    fed1 = h.helicsCreateCombinationFederate("fed1", fedInfo2)
    fed2 = h.helicsFederateClone(fed1)
    fed3 = h.helicsGetFederateByName("fed1")
    h.helicsFederateSetFlagOption(fed2, 1, false)

    h.helicsFederateSetTimeProperty(fed2, h.HELICS_PROPERTY_TIME_INPUT_DELAY, 1.0)
    h.helicsFederateSetIntegerProperty(fed1, h.HELICS_PROPERTY_INT_LOG_LEVEL, 1)
    h.helicsFederateSetIntegerProperty(fed2, h.HELICS_PROPERTY_INT_MAX_ITERATIONS, 100)
    h.helicsFederateSetTimeProperty(fed2, h.HELICS_PROPERTY_TIME_OUTPUT_DELAY, 1.0)
    h.helicsFederateSetTimeProperty(fed2, h.HELICS_PROPERTY_TIME_PERIOD, 0.0)
    h.helicsFederateSetTimeProperty(fed2, h.HELICS_PROPERTY_TIME_DELTA, 1.0)

    fed1CloningFilter = h.helicsFederateRegisterCloningFilter(fed1, "fed1/Ep1")
    fed1DestinationFilter = h.helicsFederateRegisterFilter(fed1, h.HELICS_FILTER_TYPE_DELAY, "fed1DestinationFilter")
    h.helicsFilterAddDestinationTarget(fed1DestinationFilter, "ep2")

    ep1 = h.helicsFederateRegisterEndpoint(fed1, "Ep1", "string")
    ep2 = h.helicsFederateRegisterGlobalEndpoint(fed1, "Ep2", "string")
    pub1 = h.helicsFederateRegisterGlobalPublication(fed1, "pub1", h.HELICS_DATA_TYPE_DOUBLE, "")
    pub2 = h.helicsFederateRegisterGlobalTypePublication(fed1, "pub2", "complex", "")

    sub1 = h.helicsFederateRegisterSubscription(fed1, "pub1")
    sub2 = h.helicsFederateRegisterSubscription(fed1, "pub2")
    h.helicsInputAddTarget(sub2, "Ep2")
    pub3 = h.helicsFederateRegisterPublication(fed1, "pub3", h.HELICS_DATA_TYPE_STRING, "")

    pub1KeyString = h.helicsPublicationGetKey(pub1)
    pub1TypeString = h.helicsPublicationGetType(pub1)
    pub1UnitsString = h.helicsPublicationGetUnits(pub1)
    sub1KeyString = h.helicsSubscriptionGetKey(sub1)
    sub1UnitsString = h.helicsInputGetUnits(sub1)
    @test "pub1" == pub1KeyString
    @test "double" == pub1TypeString
    @test "" == pub1UnitsString
    @test "pub1" == sub1KeyString
    @test "" == sub1UnitsString

    fed1SourceFilter = h.helicsFederateRegisterFilter(fed1,
            h.HELICS_FILTER_TYPE_DELAY, "fed1SourceFilter")
    h.helicsFilterAddSourceTarget(fed1SourceFilter, "Ep2")
    h.helicsFilterAddDestinationTarget(fed1SourceFilter, "fed1/Ep1")
    h.helicsFilterRemoveTarget(fed1SourceFilter, "fed1/Ep1")
    h.helicsFilterAddSourceTarget(fed1SourceFilter, "Ep2")
    h.helicsFilterRemoveTarget(fed1SourceFilter, "Ep2")

    fed1SourceFilterNameString = h.helicsFilterGetName(fed1SourceFilter)
    @test fed1SourceFilterNameString == "fed1/fed1SourceFilter"

    sub3 = h.helicsFederateRegisterSubscription(fed1, "fed1/pub3", "")
    pub4 = h.helicsFederateRegisterTypePublication(fed1, "pub4", "int", "")

    sub4 = h.helicsFederateRegisterSubscription(fed1, "fed1/pub4", "")
    pub5 = h.helicsFederateRegisterGlobalTypePublication(fed1, "pub5", "boolean", "")

    sub5 = h.helicsFederateRegisterSubscription(fed1, "pub5", "")
    pub6 = h.helicsFederateRegisterGlobalPublication(fed1, "pub6", h.HELICS_DATA_TYPE_VECTOR, "")
    sub6 = h.helicsFederateRegisterSubscription(fed1, "pub6", "")
    pub7 = h.helicsFederateRegisterGlobalPublication(fed1, "pub7",
            h.HELICS_DATA_TYPE_NAMED_POINT, "")
    sub7 = h.helicsFederateRegisterSubscription(fed1, "pub7", "")

    h.helicsInputSetDefaultBoolean(sub5, false)
    h.helicsInputSetDefaultComplex(sub2, -9.9 + im * 2.5)
    h.helicsInputSetDefaultDouble(sub1, 3.4)
    h.helicsInputSetDefaultInteger(sub4, 6)
    h.helicsInputSetDefaultNamedPoint(sub7, "hollow", 20.0)
    h.helicsInputSetDefaultString(sub3, "default")
    sub6Default = [ 3.4, 90.9, 4.5 ]
    h.helicsInputSetDefaultVector(sub6, sub6Default)
    h.helicsEndpointSubscribe(ep2, "fed1/pub3")
    h.helicsFederateEnterInitializingModeAsync(fed1)
    rs = h.helicsFederateIsAsyncOperationCompleted(fed1)
    if (rs == 0)
        sleep(0.500)
        rs = h.helicsFederateIsAsyncOperationCompleted(fed1)
        if (rs == 0)
            sleep(.500)
            rs = h.helicsFederateIsAsyncOperationCompleted(fed1)
            if (rs == 0)
                @test true == false
            end
        end
    end
    h.helicsFederateEnterInitializingModeComplete(fed1)
    h.helicsFederateEnterExecutingModeAsync(fed1)
    h.helicsFederateEnterExecutingModeComplete(fed1)

    mesg1 = h.helicsFederateCreateMessageObject(fed1)
    h.helicsMessageSetString(mesg1, "Hello")
    h.helicsMessageSetSource(mesg1, "fed1/Ep1")
    h.helicsMessageSetOriginalSource(mesg1, "fed1/Ep1")
    h.helicsMessageSetDestination(mesg1, "Ep2")
    h.helicsMessageSetOriginalDestination(mesg1, "Ep2")

    h.helicsEndpointSendMessageObject(ep1, mesg1)
    mesg1 = h.helicsFederateCreateMessageObject(fed1)
    h.helicsMessageSetString(mesg1, "There")
    h.helicsMessageSetSource(mesg1, "fed1/Ep1")
    h.helicsMessageSetOriginalSource(mesg1, "fed1/Ep1")
    h.helicsMessageSetDestination(mesg1, "Ep2")
    h.helicsMessageSetOriginalDestination(mesg1, "Ep2")
    h.helicsEndpointSendMessageObject(ep1, mesg1)
    h.helicsEndpointSetDefaultDestination(ep2, "fed1/Ep1")

    ep1NameString = h.helicsEndpointGetName(ep1)
    ep1TypeString = h.helicsEndpointGetType(ep1)

    @test ep1NameString == "fed1/Ep1"
    @test ep1TypeString == "string"

    coreFed1 = h.helicsFederateGetCoreObject(fed1)

    fed1Time = h.helicsFederateGetCurrentTime(fed1)
    @test fed1Time == 0.0
    fed1EndpointCount = h.helicsFederateGetEndpointCount(fed1)
    @test fed1EndpointCount == 2

    fed1NameString = h.helicsFederateGetName(fed1)
    @test fed1NameString == "fed1"

    fed1State = h.helicsFederateGetState(fed1)
    @test fed1State == 2
    fed1PubCount = h.helicsFederateGetPublicationCount(fed1)
    @test fed1PubCount == 7
    fed1SubCount = h.helicsFederateGetInputCount(fed1)
    @test fed1SubCount == 7

    h.helicsPublicationPublishBoolean(pub5, true)
    h.helicsPublicationPublishComplex(pub2, 5.6 + im * -0.67)
    h.helicsPublicationPublishDouble(pub1, 457.234)
    h.helicsPublicationPublishInteger(pub4, 1)
    h.helicsPublicationPublishNamedPoint(pub7, "Blah Blah", 20.0)
    h.helicsPublicationPublishString(pub3, "Mayhem")
    pub6Vector = [ 4.5, 56.5 ]
    h.helicsPublicationPublishVector(pub6, pub6Vector)
    sleep(0.500)
    h.helicsFederateRequestTimeAsync(fed1, 1.0)

    returnTime = h.helicsFederateRequestTimeComplete(fed1)
    @test returnTime == 1.0
    ep2MsgCount = h.helicsEndpointPendingMessages(ep2)
    @test ep2MsgCount == 2
    ep2HasMsg = h.helicsEndpointHasMessage(ep2)
    @test ep2HasMsg == 1

    msg2 = h.helicsEndpointGetMessageObject(ep2)
    @test h.helicsMessageGetTime(msg2) == 1.0
    @test h.helicsMessageGetString(msg2) == "Hello"
    @test h.helicsMessageGetOriginalSource(msg2) == "fed1/Ep1"
    @test h.helicsMessageGetSource(msg2) == "fed1/Ep1"
    @test h.helicsMessageGetDestination(msg2) == "Ep2"
    @test h.helicsMessageGetOriginalDestination(msg2) == "Ep2"

    fed1MsgCount = h.helicsFederatePendingMessages(fed1)
    @test fed1MsgCount == 1

    @test h.helicsFederateHasMessage(fed1) == 1

    msg3 = h.helicsFederateGetMessageObject(fed1)
    @test h.helicsMessageGetTime(msg3) == 1.0
    @test h.helicsMessageGetString(msg3) == "There"
    @test h.helicsMessageGetOriginalSource(msg3) == "fed1/Ep1"
    @test h.helicsMessageGetSource(msg3) == "fed1/Ep1"
    @test h.helicsMessageGetDestination(msg3) == "Ep2"
    @test h.helicsMessageGetOriginalDestination(msg3) == "Ep2"

    sub1Updated = h.helicsInputIsUpdated(sub1)
    @test sub1Updated == 1

    @test h.helicsInputLastUpdateTime(sub2) == 1.0

    @test h.helicsInputGetComplex(sub2) == 5.6 - im * 0.67

    @test h.helicsInputGetDouble(sub1) == 457.234
    @test h.helicsInputGetInteger(sub4) == 1
    sub7PointString, sub7DoubleValue = h.helicsInputGetNamedPoint(sub7)
    @test sub7PointString == "Blah Blah"
    @test sub7DoubleValue == 20.0
    @test h.helicsInputGetBoolean(sub5) == 1
    @test h.helicsInputGetString(sub3) == "Mayhem"

    sub3ValueSize = h.helicsInputGetRawValueSize(sub3)
    @test sub3ValueSize == 6

    if typeof(1) == Int32
        @test_broken h.helicsInputGetVector(sub6) == [4.5, 56.5]
    else
        @test h.helicsInputGetVector(sub6) == [4.5, 56.5]
    end

    h.helicsFederateFinalize(fed1)
    h.helicsFederateFinalize(fed2)
    h.helicsFederateFree(fed1)
    h.helicsFederateFinalize(fed2)
    h.helicsFederateFree(fed2)
    h.helicsFederateInfoFree(fedInfo2)
    h.helicsBrokerDisconnect(broker3)

    h.helicsBrokerFree(broker3)

    h.helicsCleanupLibrary()
    h.helicsCloseLibrary()
end


text:
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.Mesh.Grids
using ClimateMachine.Thermodynamics
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: grav, MSLP
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_test!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    FT = eltype(state)
    param_set = parameter_set(bl)

    z = FT(z)
    _grav::FT = grav(param_set)
    _MSLP::FT = MSLP(param_set)

    # These constants are those used by Stevens et al. (2005)
    qref = FT(9.0e-3)
    q_pt_sfc = PhasePartition(qref)
    Rm_sfc = FT(gas_constant_air(param_set, q_pt_sfc))
    T_sfc = FT(290.4)
    P_sfc = _MSLP

    # Specify moisture profiles
    q_liq = FT(0)
    q_ice = FT(0)

    θ_liq = FT(289.0)
    q_tot = qref

    ugeo = FT(7)
    vgeo = FT(-5.5)
    u, v, w = ugeo, vgeo, FT(0)

    # Pressure
    H = Rm_sfc * T_sfc / _grav
    p = P_sfc * exp(-z / H)

    # Density, Temperature
    ts = PhaseEquil_pθq(param_set, p, θ_liq, q_tot)
    ρ = air_density(ts)

    e_kin = FT(1 / 2) * FT((u^2 + v^2 + w^2))
    e_pot = _grav * z
    E = ρ * total_energy(e_kin, e_pot, ts)

    state.ρ = ρ
    state.ρu = SVector(ρ * u, ρ * v, ρ * w)
    state.energy.ρe = E
    state.moisture.ρq_tot = ρ * q_tot

    return nothing
end

function main()
    @test_throws ArgumentError ClimateMachine.init(dsisable_gpu = true)
    ClimateMachine.init()

    FT = Float64

    # DG polynomial orders
    N = (4, 4)

    # Domain resolution and size
    Δh = FT(40)
    Δv = FT(40)
    resolution = (Δh, Δh, Δv)

    xmax = FT(320)
    ymax = FT(320)
    zmax = FT(400)

    t0 = FT(0)
    timeend = FT(10)
    CFL = FT(0.4)

    ode_solver = ClimateMachine.ExplicitSolverType()
    driver_config = ClimateMachine.AtmosLESConfiguration(
        "Driver test",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_test!,
        solver_type = ode_solver,
    )
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        Courant_number = CFL,
    )

    # Test the courant wrapper
    # by default the CFL should be less than what asked for
    CFL_nondiff = ClimateMachine.DGMethods.courant(
        ClimateMachine.Courant.nondiffusive_courant,
        solver_config,
    )
    @test CFL_nondiff < CFL
    CFL_adv = ClimateMachine.DGMethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config,
    )
    CFL_adv_v = ClimateMachine.DGMethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config;
        direction = VerticalDirection(),
    )
    CFL_adv_h = ClimateMachine.DGMethods.courant(
        ClimateMachine.Courant.advective_courant,
        solver_config;
        direction = HorizontalDirection(),
    )

    # compute known advective Courant number (based on initial conditions)
    ugeo_abs = FT(7)
    vgeo_abs = FT(5.5)
    Δt = solver_config.dt
    ca_h = ugeo_abs * (Δt / Δh) + vgeo_abs * (Δt / Δh)
    # vertical velocity is 0
    caᵥ = FT(0.0)
    @test isapprox(CFL_adv_v, caᵥ, atol = 10 * eps(FT))
    @test isapprox(CFL_adv_h, ca_h, atol = 0.0005)
    @test isapprox(CFL_adv, ca_h, atol = 0.0005)

    cb_test = 0
    result = ClimateMachine.invoke!(solver_config)
    # cb_test should be zero since user_info_callback not specified
    @test cb_test == 0

    result = ClimateMachine.invoke!(
        solver_config,
        user_info_callback = () -> cb_test += 1,
    )
    # cb_test should be greater than one if the user_info_callback got called
    @test cb_test > 0

    # Test that if dt is not adjusted based on final time the CFL is correct
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        Courant_number = CFL,
        timeend_dt_adjust = false,
    )

    CFL_nondiff = ClimateMachine.DGMethods.courant(
        ClimateMachine.Courant.nondiffusive_courant,
        solver_config,
    )
    @test CFL_nondiff ≈ CFL

    # Test that if dt is not adjusted based on final time the CFL is correct
    for fixed_number_of_steps in (-1, 0, 10)
        solver_config = ClimateMachine.SolverConfiguration(
            t0,
            timeend,
            driver_config,
            Courant_number = CFL,
            fixed_number_of_steps = fixed_number_of_steps,
        )

        if fixed_number_of_steps < 0
            @test solver_config.timeend == timeend
        else
            @test solver_config.timeend ==
                  solver_config.dt * fixed_number_of_steps
            @test solver_config.numberofsteps == fixed_number_of_steps
        end
    end
end

main()


text:
# src/core.jl
# Core Abstract Definitions

import Base: numerator,
             denominator,
             divgcd,
             iszero,
             isone,
             isinteger,
             iseven,
             isodd,
             isinf,
             isnan,
             decompose

export missing_api, AbstractRatio, AbstractRational, build_ratio, numden

"""```
missing_api(name, args...; kwargs...)
```
Throw an `ArgumentError` if the function with name `name` is not defined on
the given `args` and `kwargs`.
"""
function missing_api(name::Union{Symbol,AbstractString}, args...; kwargs...)
    throw(ArgumentError("`$name` is not defined for `$args, $kwargs`"))
end

"""```
AbstractRatio{T <: Integer} <: Real
```
Abstraction of generic ratios. Numerator and denominator need not be in canonical form.
"""
abstract type AbstractRatio{T<:Integer} <: Real end

"""```
AbstractRational{T} <: AbstractRatio{T}
```
Abstraction of generic rational numbers. Numerator and denominator must always be in canonical form.
"""
abstract type AbstractRational{T} <: AbstractRatio{T} end

"""```
build_ratio(R::Type{<:AbstractRatio}, num, den, args...; kwargs...)
```
Replacement for `R` ratio constructor which accepts two integers: `num` and `den`.
"""
function build_ratio(R::Type{<:AbstractRatio}, num::Integer, den::Integer, args...; kwargs...)::R
    missing_api("build_ratio", R, num, den, args...; kwargs...)
end

"""```
build_ratio(R::Type{<:AbstractRatio}, r::AbstractRatio, args...; kwargs...)
```
Replacement for `R` ratio constructor which accepts a preexisting ratio: `ratio`.
"""
function build_ratio(R::Type{<:AbstractRatio}, ratio::AbstractRatio, args...; kwargs...)::R
    missing_api("build_ratio", R, ratio, args...; kwargs...)
end

"""```
numerator(x::AbstractRatio)
```
Return the numerator of the ratio `x`.
"""
function numerator(x::AbstractRatio)
    missing_api("numerator", x)
end

"""```
denominator(x::AbstractRatio)
```
Return the denominator of the ratio `x`.
"""
function denominator(x::AbstractRatio)
    missing_api("denominator", x)
end

"""```
numden(x)
    === numerator(x), denominator(x)
```
Return the numerator and denominator of the ratio `x`.
"""
numden(x) = numerator(x), denominator(x)

"""```
crossterms(x, y; product = *)
    === product(numerator(x), denominator(y)), product(numerator(y), denominator(x))
```
Return the cross terms in the construction of the common denominator.
"""
function crossterms(x, y; product = *)
    return product(numerator(x), denominator(y)), product(numerator(y), denominator(x))
end

"""```
resolve(R::Type{<:AbstractRational}, x)
```
Resolve a ratio `x` into a rational number in canonical form with type `R`.
"""
function resolve(R::Type{<:AbstractRational}, x)::R
    num, den = numden(x)
    return build_ratio(R, (signbit(den) ? divgcd(-num, -den) : divgcd(num, den))...)
end

"""```
iszero(x::AbstractRatio)
    === iszero(numerator(x)) && !iszero(denominator(x))
```
Check if the ratio `x` is the zero ratio.
"""
iszero(x::AbstractRatio) = iszero(numerator(x)) && !iszero(denominator(x))

"""```
iszero(x::AbstractRational)
    === iszero(numerator(x))
```
Check if the rational `x` is the zero rational.
"""
iszero(x::AbstractRational) = iszero(numerator(x))

"""```
isone(x::AbstractRatio)
    === numerator(x) == denominator(x)
```
Check if the ratio `x` is the one ratio.
"""
isone(x::AbstractRatio) = numerator(x) == denominator(x)

"""```
isone(x::AbstractRational)
    === isone(numerator(x)) & isone(denominator(x))
```
Check if the rational `x` is the one rational.
"""
isone(x::AbstractRational) = isone(numerator(x)) & isone(denominator(x))

"""```
isinteger(x::AbstractRatio)
    === iszero(rem(numden(x)...))
```
Check if the ratio `x` is an integer.
"""
isinteger(x::AbstractRatio) = iszero(rem(numden(x)...))

"""```
isinteger(x::AbstractRational)
    === isone(denominator(x))
```
Check if the rational `x` is an integer.
"""
isinteger(x::AbstractRational) = isone(denominator(x))

"""```
iseven(x::AbstractRatio)
```
Check if the ratio `x` is an even integer.
"""
function iseven(x::AbstractRatio)
    d, r = divrem(numden(x)...)
    return iszero(r) && iseven(d)
end

"""```
iseven(x::AbstractRational)
    === isinteger(x) && iseven(numerator(x))
```
Check if the rational `x` is an even integer.
"""
function iseven(x::AbstractRational)
    return isinteger(x) && iseven(numerator(x))
end

"""```
isodd(x::AbstractRatio)
```
Check if the ratio `x` is an odd integer.
"""
function isodd(x::AbstractRatio)
    d, r = divrem(numden(x)...)
    return iszero(r) && isodd(d)
end

"""```
isodd(x::AbstractRational)
    === isinteger(x) && isodd(numerator(x))
```
Check if the rational `x` is an odd integer.
"""
function isodd(x::AbstractRational)
    return isinteger(x) && isodd(numerator(x))
end

"""```
isinf(x::AbstractRatio)
```
Check if the ratio `x` is infinite.
"""
function isinf(x::AbstractRatio)
    return iszero(denominator(x)) & !iszero(numerator(x))
end

"""```
isinf(::AbstractRational)
    === false
```
There are no infinite rationals.
"""
isinf(::AbstractRational) = false

"""```
isnan(x::AbstractRatio)
    === iszero(numerator(x)) & iszero(denominator(x))
```
Check if the ratio `x` is not a number.
"""
isnan(x::AbstractRatio) = all(iszero, numden(x))

"""```
isnan(::AbstractRational)
    === false
```
There are no rationals which are not numbers.
"""
isnan(::AbstractRational) = false

"""```
decompose(x::AbstractRatio)
    === numerator(x), zero(x), denominator(x)
```
Decompose the ratio `x`.
"""
decompose(x::AbstractRatio) = numerator(x), zero(x), denominator(x)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

using Documenter, FEMBeam

makedocs(modules = [FEMBeam],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "FEMBeam.jl",
         pages = ["index.md"])


text:
const _GW_WORKER_FILE_NAME = "gitworker.toml"
const _GW_WORKER_FILE_SYSROOT_KEY = "sys_root"
const _GW_WORKER_FILE_REMOTE_URL_KEY = "remote_url"

_worker_file(gw::GitWorker) = joinpath(worker_dir(gw), _GW_WORKER_FILE_NAME)

function _write_worker_file(gw::GitWorker)
    wdict = Dict{String, String}()
    wdict[_GW_WORKER_FILE_SYSROOT_KEY] = sys_root(gw)
    wdict[_GW_WORKER_FILE_REMOTE_URL_KEY] = remote_url(gw)
    wfile = _worker_file(gw)
    return _write_toml(wfile, wdict)
end

function _find_worker_file(path)
    @show path
    if isdir(path) # Base
        wfile = joinpath(path, _GW_WORKER_FILE_NAME)
        isfile(wfile) && return wfile
    end

    dir = dirname(path)
    (dir == path) && return "" # Base
    return _find_worker_file(dir)
end

function _gw_from_toml(gwfile::String)
    gdat = _read_toml(gwfile)
    sys_root = get(gdat, _GW_WORKER_FILE_SYSROOT_KEY, "")
    remote_url = get(gdat, _GW_WORKER_FILE_REMOTE_URL_KEY, "")
    (isempty(sys_root) || isempty(remote_url)) && return nothing
    return GitWorker(;sys_root, remote_url)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
#=
Building graphs from CIF files using PyCall to the pymatgen package.
=#

using PyCall
using ChemistryFeaturization
using Serialization

# options for decay of bond weights with distance...
inverse_square(x) = x^-2.0
exp_decay(x) = exp(-x)

# a few helper fcns for readability...
"Return the index of a given site in the structure."
site_index(site) = convert(UInt16, get(site, 2)) + 1

"Return the distance associated with a site in a neighbor list."
site_distance(site) = convert(Float64, get(site, 1))

# these next two functions return the information for the first species in a list – there should only be one because otherwise the structure would be disordered and we probably shouldn't be building a graph...(or maybe we add some fancy functionality later to do superpositions of species?)

"Return atomic number associated with a site."
site_atno(site) = [e.Z for e in site.species.elements][1]

"Return atomic symbol associated with a site."
site_element(site) = [e.symbol for e in site.species.elements][1]

"""
    are_equidistant(site1, site2)

Check if site1 and site2 are equidistant to within tolerance atol, in angstroms (for cutting off neighbor lists consistently).

Note that this only works if site1 and site2 are from a neighbor list from the same central atom.
"""
function are_equidistant(site1, site2, atol=1e-4)
    isapprox(site_distance(site1), site_distance(site2), atol=atol)
end

# TODO: figure out best/idiomatic way to pass through the keyword arguments, surely the copy/paste is not it
"""
Function to build graph from a file storing a crystal structure (currently supports anything ase.io.read can read in). Returns an AtomGraph object.

# Arguments
- `struc`: pymatgen Structure object
- `use_voronoi::bool`: if true, use Voronoi method for neighbor lists, if false use cutoff method

    (The rest of these parameters are only used if use_voronoi is false)

- `radius::Float=8.0`: cutoff radius for atoms to be considered neighbors (in angstroms)
- `max_num_nbr::Integer=12`: maximum number of neighbors to include (even if more fall within cutoff radius)
- `dist_decay_func`: function (e.g. inverse_square or exp_decay) to determine falloff of graph edge weights with neighbor distance
"""

# TODO: featurize here
function build_graph(file_path::String; use_voronoi=false, radius=8.0, max_num_nbr=12, dist_decay_func=inverse_square, normalize=true)
    s = pyimport("pymatgen.core.structure")
    aseio = pyimport("ase.io")
    pmgase = pyimport("pymatgen.io.ase")
    atoms_object = aseio.read(file_path)
    aa = pmgase.AseAtomsAdaptor()
    struc = aa.get_structure(atoms_object)

    num_atoms = size(struc)[1]

    # list of atom symbols
    atom_ids = [site_element(s) for s in struc]

    if use_voronoi
        weight_mat = weights_voronoi(struc)
    else
        weight_mat = weights_cutoff(struc; radius=radius, max_num_nbr=max_num_nbr, dist_decay_func=dist_decay_func)
    end

    # normalize weights
    if normalize
        weight_mat = weight_mat ./ maximum(weight_mat)
    end

    g = SimpleWeightedGraph{Int32}(Float32.(weight_mat))
    return AtomGraph(g, atom_ids)
end

"""
Build graph using neighbors from faces of Voronoi polyedra and weights from areas. Based on the approach from https://github.com/ulissigroup/uncertainty_benchmarking/blob/aabb407807e35b5fd6ad06b14b440609ae09e6ef/BNN/data_pyro.py#L268
"""
function weights_voronoi(struc)
    num_atoms = size(struc)[1]
    sa = pyimport_conda("pymatgen.analysis.structure_analyzer", "pymatgen", "conda-forge")
    vc = sa.VoronoiConnectivity(struc)
    conn = vc.connectivity_array
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    # loop over central atoms
    for atom_ind in 1:size(conn)[1]
        # loop over neighbor atoms
        for nb_ind in 1:size(conn)[2]
            # loop over each possible PBC image for chosen image
            for image_ind in 1:size(conn)[3]
                # only add as neighbor if atom is not current center one AND there is connectivity to image
                if (atom_ind != image_ind) & (conn[atom_ind, nb_ind, image_ind] != 0)
                    weight_mat[atom_ind, nb_ind] += conn[atom_ind, nb_ind, image_ind]/maximum(conn[atom_ind, :, :])
                end
            end
        end
    end

    # average across diagonal (because neighborness isn't strictly symmetric in the way we're defining it here)
    weight_mat = 0.5.* (weight_mat .+ weight_mat')
end

"""
Build graph using neighbor number cutoff method adapted from original CGCNN. Note that `max_num_nbr` is a "soft" max, in that if there are more of the same distance as the last, all of those will be added.

# TODO
- option to cut off by nearest, next-nearest, etc. by DISTANCE rather than NUMBER of neighbors
"""
function weights_cutoff(struc; radius=8.0, max_num_nbr=12, dist_decay_func=inverse_square)
    #= find neighbors, requires a cutoff radius
    returns a NxM Array of PyObject PeriodicSite
    ...except when it returns a list of N of lists of length M...
    each PeriodicSite is (site, distance, index, image)
    N is num sites in crystal, M is num neighbors
    =#
    num_atoms = size(struc)[1]
    all_nbrs = struc.get_all_neighbors(radius, include_index=true)

    # sort by distance
    # returns list of length N of lists of length M
    if length(size(all_nbrs)) == 2
        all_nbrs = [sort(all_nbrs[i,:], lt=(x,y)->isless(site_distance(x), site_distance(y))) for i in 1:num_atoms]
    elseif length(size(all_nbrs)) == 1
        all_nbrs = [sort(all_nbrs[i][:], lt=(x,y)->isless(site_distance(x), site_distance(y))) for i in 1:num_atoms]
    end

    # iterate through each list of neighbors (corresponding to neighbors of a given atom) to find bonds (eventually, graph edges)
    weight_mat = zeros(Float32, num_atoms, num_atoms)
    for atom_ind in 1:num_atoms
        this_atom = get(struc, atom_ind-1)
        atom_nbs = all_nbrs[atom_ind]
        # iterate over each neighbor...
        for nb_num in 1:size(all_nbrs[atom_ind])[1]
            nb = atom_nbs[nb_num]
            global nb_ind = site_index(nb)
            # if we're under the max, add it for sure
            if nb_num < max_num_nbr
                weight_mat[atom_ind, nb_ind] = weight_mat[atom_ind, nb_ind] + dist_decay_func(site_distance(nb))
            # if we're at/above the max, add if distance is the same
            else
                # check we're not on the last one
                if nb_ind < size(atom_nbs)[1] - 1
                    next_nb = atom_nbs[nb_ind + 1]
                    # add another bond if it's the exact same distance to the next neighbor in the list
                    if are_equidistant(nb, next_nb)
                        weight_mat[atom_ind, nb_ind] = weight_mat[atom_ind, nb_ind] + dist_decay_func(site_distance(nb))
                    end
                end
            end
        end
    end

    # average across diagonal (because neighborness isn't strictly symmetric in the way we're defining it here)
    weight_mat = 0.5.* (weight_mat .+ weight_mat')
end

"""
Function to build and serialize to file a batch of CIFs, optionally featurizing them as well. Saved .jls files will have the same names as the .cif ones but with the extensions modified.

# Arguments
- `cif_folder::String`: path to folder containing CIF files
- `output_folder::String`: path to folder where .jls files containing AtomGraph objects should be saved (will be created if it doesn't exist already)

If you want to featurize the graphs, at least the vector of `AtomFeat` objects describing the featurization procedure is required, and optionally the Dict mapping from elemental symbols to feature vectors (if it is not provided, it will be generated).

Other optional arguments are the optional arguments to `build_graph`: `use_voronoi`, `radius`, `max_num_nbr`, `dist_decay_func`, `normalize`

This function does not return anything.
    TODO: decide if there should be an option to return the graphs
"""
function build_graphs_batch(input_folder::String, output_folder::String, featurization=AtomFeat[]; atom_featurevecs=Dict{String, Vector{Float32}}(), use_voronoi=false, radius=8.0, max_num_nbr=12, dist_decay_func=inverse_square, normalize=true)
    # check if input folder exists and contains things, if not throw error
    file_list = readdir(input_folder, join=true)
    if length(file_list)==0
        error("No files in input directory!")
    end

    # check if output folder exists, if not create it
    if !isdir(output_folder)
        mkdir(output_folder)
        @info "Output path provided did not exist, creating folder there."
    end

    # check if there is enough information to actually featurize
    featurize = length(featurization)>0
    atom_featurevecs = (length(atom_featurevecs)>0 & featurize) ? atom_featurevecs : make_feature_vectors(featurization)[1]

    if (length(atom_featurevecs)>0) & (length(featurization)==0)
        @warn "You have supplied only feature vectors but no featurization scheme, so graphs will be built but not featurized."
    end

    # loop over CIFs and build graphs from all files that we can...
    for file in file_list
        id = split(splitpath(file)[end], ".")[1]
        local ag
        try
            ag = build_graph(file; use_voronoi=use_voronoi, radius=radius, max_num_nbr=max_num_nbr, dist_decay_func=dist_decay_func, normalize=normalize)
        catch
            return id
        end
        if featurize
            add_features!(ag, atom_featurevecs, featurization)
        end
        graph_path = joinpath(output_folder, string(id, ".jls"))
        serialize(graph_path, ag)
    end
end

# alternate call signature where featurization is generated
function build_graphs_batch(cif_folder::String, output_folder::String, feature_names::Vector{Symbol}; nbins::Vector{<:Integer}=default_nbins*ones(Int64, size(feature_names,1)), logspaced=false)
    atom_featurevecs, featurization = make_feature_vectors(build_atom_feats(feature_names; nbins=nbins, logspaced=logspaced))
    build_graphs_batch(cif_folder, output_folder, featurization; atom_featurevecs = atom_featurevecs)
end


text:
module gargamel

export legendrepoly, derlegendre

"""
    legendrepoly(x ::Float64,k)

    Evaluate all Legendre polynomials up to degree k at x
"""
function legendrepoly(x ::Float64,k)

    P = zeros(k+1)
    P[1] = 1.0
    if k>0
       P[2] = x
    end
    for i = 2:k
            P[i+1]=((2.0*i-1)*x*P[i+1-1]-(i-1.0)*P[i+1-2])/Float64(i)
    end
    return P
end

"""
    legendrepoly(x ::Float64,k)

    Evaluate all derivatives of Legendre polynomials up to degree k at x
"""
function derlegendre(x ::Float64,k)

    P = legendrepoly(x,k)
    PX = zeros(k+1)
    PX[1] = 0.0
    if k > 1
        PX[2] = 1.0
    end
    for i=2:k
        PX[i+1] = PX[i+1-2]+(2*i-1.0)*P[i+1-1]
    end
    return PX
end
# Write your package code here.

end


text:
using HealthMLBase
using Documenter

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

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

deploydocs(;
    repo="github.com/JuliaHealth/HealthMLBase.jl",
)


text:
####################################################################################################
function displayIteration(i, funceval, residual, itlinear = 0)
	(i == 0) && println("\n Newton Iterations \n   Iterations      Func-count      f(x)      Linear-Iterations\n")
	if length(itlinear)==1
		@printf("%9d %16d %14.4e %9d\n", i, funceval, residual, itlinear);
	else
		if itlinear isa Tuple{Int64,Int64}
			@printf("%9d %16d %14.4e (%6d, %6d)\n", i, funceval, residual, itlinear[1], itlinear[2]);
		else
			# used for nested linear solves
			@printf("%9d %16d %14.4e  ", i, funceval, residual);println(itlinear);
		end
	end
end
####################################################################################################
function computeEigenvalues(contparams::ContinuationPar, n_unstable::Int64, J)
	# the next line is to ensure we compute enough eigenvalues to probe stability
	nev_ = max(n_unstable + 6, contparams.nev)
	eiginfo = contparams.newtonOptions.eigsolver(J, nev_)
	isstable, n_unstable, n_imag = is_stable(contparams, eiginfo[1])
	return eiginfo, isstable, n_unstable, n_imag
end

function computeEigenvalues(iter::PALCIterable, state::PALCStateVariables)
	n_unstable = state.n_unstable[2]
	nev_ = max(n_unstable + 5, iter.contParams.nev)
	J = iter.J(getx(state), set(iter.par, iter.param_lens, getp(state)))
	eiginfo = iter.contParams.newtonOptions.eigsolver(J, nev_)
	isstable, n_unstable, n_imag = is_stable(iter.contParams, eiginfo[1])
	return eiginfo, isstable, n_unstable, n_imag
end

function computeEigenvalues!(iter::PALCIterable, state::PALCStateVariables)
	# we compute the eigen-elements
	n_unstable = state.n_unstable[2]
	nev_ = max(n_unstable + 5, iter.contParams.nev)
	J = iter.J(getx(state), set(iter.par, iter.param_lens, getp(state)))
	eigvals, eigvecs, flag, it_number = iter.contParams.newtonOptions.eigsolver(J, nev_)
	# we update the state
	isstable, n_unstable, n_imag = is_stable(iter.contParams, eigvals)
	updatestability!(state, n_unstable, n_imag)
	if isnothing(state.eigvals) == false
		state.eigvals = eigvals
	end
	if iter.contParams.saveEigenvectors
		state.eigvecs = eigvecs
	end
	return it_number
end

####################################################################################################
"""
	finiteDifferences(F, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences
"""
function finiteDifferences(F, x::AbstractVector; δ = 1e-9)
	N = length(x)
	Nf = length(F(x))
	J = zeros(eltype(x), Nf, N)
	finiteDifferences!(F, J, x; δ = δ)
	return J
end

"""
	finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences, update J
"""
function finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)
	f = F(x)
	x1 = copy(x)
	for i ∈ eachindex(x)
		x1[i] += δ
		J[:, i] .= (F(x1) .- F(x)) / δ
		x1[i] -= δ
	end
	return J
end

####################################################################################################
"""
Save solution / data in JLD2 file
- `filename` is for example "example.jld2"
- `sol` is the solution
- `p` is the parameter
- `i` is the index of the solution to be saved
"""
function saveToFile(filename, sol, p, i::Int64, br::ContResult, contParam::ContinuationPar)
	try
		# create a group in the JLD format
		jldopen(filename*".jld2", "a+") do file
			mygroup = JLD2.Group(file, "sol-$i")
			mygroup["sol"] = sol
			mygroup["param"] = p
		end

		jldopen(filename*"-branch.jld2", "w") do file
			file["branch"] = br
			file["contParam"] = contParam
		end
	catch
		@warn "Could not write in jld2 file"
	end
end
####################################################################################################
using BlockArrays, SparseArrays
"""
Function waiting to be accepted to BlockArrays.jl
"""
function blockToSparse(J::AbstractBlockArray)
	nl, nc = size(J.blocks)
	# form the first line of blocks
	res = J[Block(1,1)]
	for j in 2:nc
		res = hcat(res, J[Block(1,j)])
	end
	# continue with the other lines
	for i in 2:nl
		line = J[Block(i,1)]
		for j in 2:nc
			line = hcat(line, J[Block(i,j)])
		end
		res = vcat(res,line)
	end
	return res
end
####################################################################################################
# this function extracts the indices of the blocks composing the matrix A. We assume that the blocks have the same sparsity
function getBlocks(A::SparseMatrixCSC, N, M)
	I,J,K = findnz(A)
	out = [Vector{Int}() for i in 1:M+1, j in 1:M+1];
	for k in eachindex(I)
		m, l = div(I[k]-1, N), div(J[k]-1, N)
		push!(out[1+m,1+l], k)
	end
	res = [length(m) for m in out]
	out
end
####################################################################################################
# this trick is extracted from KrylovKit. It allows for the Jacobian to be specified as a matrix (sparse / dense) or as a function.
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)

# apply!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
# apply!(y, f, x) = copyto!(y, f(x))

# empty eigenvectors to save memory
# _empty(a::AbstractVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()
# _empty(a::AbstractMatrix{T}, ::Type{U}=T) where {T,U} = similar(a, (0,0))


text:
#
# Wrappers for implementing SCIP constraint handlers in Julia.
#
# Please study the corresponding SCIP documentation first, to become familiar
# with basic concepts and terms: https://scip.zib.de/doc-6.0.2/html/CONS.php
#
# The basic idea is that you create a new subtype of `AbstractConstraintHandler`
# to store the constraint handler data and implement the fundamental callbacks
# by adding methods to the functions `check`, `enforce_lp_sol`,
# `enforce_pseudo_sol` and `lock`.
#
# If your constraint handler uses constraint objects (`needs_constraints`), you
# create a subtype of the parametrized `AbstractConstraint`.
#
#
# Current limitations:
# - We now use fixed values for some properties:
#   - SEPAPRIORITY = 0
#   - SEPAFREQ = -1
#   - DELAYSEPA = FALSE
#   - PROPFREQ = -1
#   - DELAYPROP = -1
#   - PROP_TIMING = SCIP_PROPTIMING_BEFORELP
#   - PRESOLTIMING = SCIP_PRESOLTIMING_MEDIUM
#   - MAXPREROUNDS = -1
# - We don't support these optional methods: CONSHDLRCOPY, CONSFREE, CONSINIT,
#   CONSEXIT, CONSINITPRE, CONSEXITPRE, CONSINITSOL, CONSEXITSOL, CONSDELETE,
#   CONSTRANS, CONSINITLP, CONSSEPALP, CONSSEPASOL, CONSENFORELAX, CONSPROP,
#   CONSPRESOL, CONSRESPROP, CONSACTIVE, CONSDEACTIVE, CONSENABLE, CONSDISABLE,
#   CONSDELVARS, CONSPRINT, CONSCOPY, CONSPARSE, CONSGETVARS, CONSGETNVARS,
#   CONSGETDIVEBDCHGS
# - We don't support linear or nonlinear constraint upgrading.
#

#
# Abstract Supertypes:
#

"""
    AbstractConstraintHandler

Abstract supertype for objects that store all user data belonging to the
constraint handler. This object is also used for method dispatch with the
callback functions.

From SCIP's point-of-view, this objects corresponds to the SCIP_CONSHDLRDATA,
but its memory is managed by Julia's GC.

It's recommended to store a reference to your instance of `ManagedSCIP` or
`SCIP.Optimizer` here, so that you can use it within your callback methods.
"""
abstract type AbstractConstraintHandler end

"""
    AbstractConstraint{Handler}

Abstract supertype for objects that store all user data belonging to an
individual constraint. It's parameterized by the type of constraint handler.

From SCIP's point-of-view, this objects corresponds to the SCIP_CONSDATA,
but its memory is managed by Julia's GC.
"""
abstract type AbstractConstraint{Handler} end

#
# Fundamental Callbacks for Julia.
#

"""
    check(
        constraint_handler::CH,
        constraints::Array{Ptr{SCIP_CONS}},
        sol::Ptr{SCIP_SOL},
        checkintegrality::Bool,
        checklprows::Bool,
        printreason::Bool,
        completely::Bool
    )::SCIP_RESULT

Check whether the solution candidate given by `sol` satisfies all constraints in
`constraints`.

Use the functions `user_constraint` and `sol_values` to access your
constraint-specific user data and solution values, respectively.

Acceptable result values are:
* `SCIP_FEASIBLE`: The solution candidate satisfies all constraints.
* `SCIP_INFEASIBLE`: The solution candidate violates at least one constraint.

There is nothing else to do for the user in terms of dealing with the violation.
"""
function check end

"""
    enforce_lp_sol(
        constraint_handler::CH,
        constraints::Array{Ptr{SCIP_CONS}},
        nusefulconss::Cint,
        solinfeasible::Bool
    )::SCIP_RESULT

Enforce the current solution for the LP relaxation. That is, check all given
constraints for violation and deal with it in some way (see below).

Use the functions `user_constraint` and `sol_values` to access your
constraint-specific user data and solution values, respectively.

Acceptable result values are:
* `SCIP_FEASIBLE`: The solution candidate satisfies all constraints.
* `SCIP_CUTOFF`: The current subproblem is infeasible.
* `SCIP_CONSADDED`: Added a constraint that resolves the infeasibility.
* `SCIP_REDUCEDDOM`: Reduced the domain of a variable.
* `SCIP_SEPARATED`: Added a cutting plane.
* `SCIP_BRANCHED`: Performed a branching.
* `SCIP_INFEASIBLE`: The solution candidate violates at least one constraint.
"""
function enforce_lp_sol end

"""
    enforce_pseudo_sol(
        constraint_handler::CH,
        constraints::Array{Ptr{SCIP_CONS}},
        nusefulconss::Cint,
        solinfeasible::Bool
    )::SCIP_RESULT

Enforce the current pseudo solution (the LP relaxation was not solved). That is,
check all given constraints for violation and deal with it in some way (see
below).

Use the functions `user_constraint` and `sol_values` to access your
constraint-specific user data and solution values, respectively.

Acceptable result values are:
* `SCIP_FEASIBLE`: The solution candidate satisfies all constraints.
* `SCIP_CUTOFF`: The current subproblem is infeasible.
* `SCIP_CONSADDED`: Added a constraint that resolves the infeasibility.
* `SCIP_REDUCEDDOM`: Reduced the domain of a variable.
* `SCIP_SEPARATED`: Added a cutting plane.
* `SCIP_BRANCHED`: Performed a branching.
* `SCIP_SOLVELP`: Force the solving of the LP relaxation.
* `SCIP_INFEASIBLE`: The solution candidate violates at least one constraint.
"""
function enforce_pseudo_sol end

"""
    lock(
        constraint_handler::CH,
        constraint::Ptr{SCIP_CONS},
        locktype::SCIP_LOCKTYPE,
        nlockspos::Cint,
        nlocksneg::Cint
    )

Define the variable locks that the given constraint implies. For each related
variable, call `SCIPaddVarLocksType` to let SCIP know whether rounding the value
up or down might lead to constraint violation.
"""
function lock end

#
# Fundamental callback functions with SCIP's C API.
#
# There is only a single, generic implementation for each of these, which are
# passed to all user-defined SCIP constraint handlers, but they will call the
# user's method in the function body.
#

"""
Generic `check` function, matching the signature from SCIP's C API.
"""
function _conscheck(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR},
                    conss::Ptr{Ptr{SCIP_CONS}}, nconss::Cint,
                    sol::Ptr{SCIP_SOL}, checkintegrality::SCIP_Bool,
                    checklprows::SCIP_Bool, printreason::SCIP_Bool,
                    completely::SCIP_Bool, result::Ptr{SCIP_RESULT})
    # get Julia object out of constraint handler data
    conshdlrdata::Ptr{SCIP_CONSHDLRDATA} = SCIPconshdlrGetData(conshdlr)
    constraint_handler = unsafe_pointer_to_objref(conshdlrdata)

    # get Julia array from C pointer
    constraints = unsafe_wrap(Array{Ptr{SCIP_CONS}}, conss, nconss)

    # call user method via dispatch
    res = check(constraint_handler, constraints, sol,
                convert(Bool, checkintegrality),
                convert(Bool, checklprows),
                convert(Bool, printreason),
                convert(Bool, completely))
    unsafe_store!(result, res)

    return SCIP_OKAY
end

"""
Generic `enfolp` function, matching the signature from SCIP's C API.
"""
function _consenfolp(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR},
                     conss::Ptr{Ptr{SCIP_CONS}}, nconss::Cint,
                     nusefulconss::Cint, solinfeasible::SCIP_Bool,
                     result::Ptr{SCIP_RESULT})
    # get Julia object out of constraint handler data
    conshdlrdata::Ptr{SCIP_CONSHDLRDATA} = SCIPconshdlrGetData(conshdlr)
    constraint_handler = unsafe_pointer_to_objref(conshdlrdata)

    # get Julia array from C pointer
    constraints = unsafe_wrap(Array{Ptr{SCIP_CONS}}, conss, nconss)

    # call user method via dispatch
    res = enforce_lp_sol(constraint_handler, constraints, nusefulconss,
                         convert(Bool, solinfeasible))
    unsafe_store!(result, res)

    return SCIP_OKAY
end

"""
Generic `enfops` function, matching the signature from SCIP's C API.
"""
function _consenfops(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR},
                     conss::Ptr{Ptr{SCIP_CONS}}, nconss::Cint,
                     nusefulconss::Cint, solinfeasible::SCIP_Bool,
                     objinfeasible::SCIP_Bool, result::Ptr{SCIP_RESULT})
    # get Julia object out of constraint handler data
    conshdlrdata::Ptr{SCIP_CONSHDLRDATA} = SCIPconshdlrGetData(conshdlr)
    constraint_handler = unsafe_pointer_to_objref(conshdlrdata)

    # get Julia array from C pointer
    constraints = unsafe_wrap(Array{Ptr{SCIP_CONS}}, conss, nconss)

    # call user method via dispatch
    res = enforce_pseudo_sol(constraint_handler, constraints, nusefulconss,
                             convert(Bool, solinfeasible),
                             convert(Bool, objinfeasible))
    unsafe_store!(result, res)

    return SCIP_OKAY
end

"""
Generic `lock` function, matching the signature from SCIP's C API.
"""
function _conslock(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR},
                   cons::Ptr{SCIP_CONS}, locktype::SCIP_LOCKTYPE,
                   nlockspos::Cint, nlocksneg::Cint)
    # get Julia object out of constraint handler data
    conshdlrdata::Ptr{SCIP_CONSHDLRDATA} = SCIPconshdlrGetData(conshdlr)
    constraint_handler = unsafe_pointer_to_objref(conshdlrdata)

    # call user method via dispatch
    lock(constraint_handler, cons, locktype, nlockspos, nlocksneg)

    return SCIP_OKAY
end

#
# Additional callback functions for memory management.
#
# The implementation should work for all constraint handlers defined in Julia,
# so there is no method for the user to implement.
#

"""
Generic `free` function, matching the signature from SCIP's C API.
"""
function _consfree(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR})
    # Here, we should free the constraint handler data. But because this is an
    # object created and owned by Julia, we will let GC do it.
    # Instead, we will just set the pointer to NULL, so that SCIP will think
    # that it is taken care of.
    SCIPconshdlrSetData(conshdlr, C_NULL)

    return SCIP_OKAY
end

"""
Generic `delete` function, matching the signature from SCIP's C API.
"""
function _consdelete(scip::Ptr{SCIP_}, conshdlr::Ptr{SCIP_CONSHDLR},
                     cons::Ptr{SCIP_CONS}, consdata::Ptr{Ptr{SCIP_CONSDATA}})
    # Here, we should free the constraint data. But because this is an object
    # created and owned by Julia, we will let GC do it.
    # Instead, we will just set the pointer to NULL, so that SCIP will think
    # that it is taken care of.
    unsafe_store!(consdata, C_NULL)

    return SCIP_OKAY
end

#
# Adding constraint handlers and constraints to ManagedSCIP.
#

"""
    include_conshdlr(
        mscip::ManagedSCIP,
        ch::CH;
        name::String,
        description::String,
        enforce_priority::Int,
        check_priority::Int,
        eager_frequency::Int,
        needs_constraints::Bool
    )

Include a user defined constraint handler `ch` to the SCIP instance `mscip`.

All parameters have default values that can be set as keyword arguments.
In particular, note the boolean `needs_constraints`:
* If set to `true`, then the callbacks are only called for the constraints that
  were added explicitly using `add_constraint`.
* If set to `false`, the callback functions will always be called, even if no
  corresponding constraint was added. It probably makes sense to set
  `misc/allowdualreds` to `FALSE` in this case.
"""
function include_conshdlr(mscip::ManagedSCIP, ch::CH;
                          name="", description="", enforce_priority=-15,
                          check_priority=-7000000, eager_frequency=100,
                          needs_constraints=true) where CH <: AbstractConstraintHandler
    # Get C function pointers from Julia functions
    _enfolp = @cfunction(_consenfolp, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Ptr{SCIP_CONS}}, Cint, Cint, SCIP_Bool, Ptr{SCIP_RESULT}))
    _enfops = @cfunction(_consenfops, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Ptr{SCIP_CONS}}, Cint, Cint, SCIP_Bool, SCIP_Bool, Ptr{SCIP_RESULT}))
    _check = @cfunction(_conscheck, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Ptr{SCIP_CONS}}, Cint, Ptr{SCIP_SOL}, SCIP_Bool, SCIP_Bool, SCIP_Bool, SCIP_Bool, Ptr{SCIP_RESULT}))
    _lock = @cfunction(_conslock, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_CONS}, SCIP_LOCKTYPE, Cint, Cint))

    # Store pointer to SCIP structure (for future C API calls)
    conshdlr__ = Ref{Ptr{SCIP_CONSHDLR}}(C_NULL)

    # Hand over Julia object as constraint handler data:
    conshdlrdata_ = pointer_from_objref(ch)

    # Try to create unique name, or else SCIP will complain!
    if name == ""
        name = "__ch__$(length(mscip.conshdlrs))"
    end

    # Register constraint handler with SCIP instance.
    @SC SCIPincludeConshdlrBasic(mscip, conshdlr__, name, description,
                                 enforce_priority, check_priority,
                                 eager_frequency, needs_constraints,
                                 _enfolp, _enfops, _check, _lock,
                                 conshdlrdata_)

    # Sanity checks
    @assert conshdlr__[] != C_NULL

    # Set additional callbacks.
    @SC SCIPsetConshdlrFree(
        mscip, conshdlr__[],
        @cfunction(_consfree, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR})))
    @SC SCIPsetConshdlrDelete(
        mscip, conshdlr__[],
        @cfunction(_consdelete, SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_CONSDATA}})))

    # Register constraint handler (for GC-protection and mapping).
    mscip.conshdlrs[ch] = conshdlr__[]
end

"""
    add_constraint(
        mscip::ManagedSCIP,
        ch::CH,
        c::C;
        initial=true,
        separate=true,
        enforce=true,
        check=true,
        propagate=true,
        _local=false,
        modifiable=false,
        dynamic=false,
        removable=false,
        stickingatnode=false
    )::ConsRef

Add constraint `c` belonging to user-defined constraint handler `ch` to model.
Returns constraint reference.

All keyword arguments are passed to the `SCIPcreateCons` call.
"""
function add_constraint(mscip::ManagedSCIP, ch::CH, c::C;
                        initial=true, separate=true, enforce=true, check=true,
                        propagate=true, _local=false, modifiable=false,
                        dynamic=false, removable=false, stickingatnode=false) where {CH <:AbstractConstraintHandler, C <: AbstractConstraint{CH}}
    # Find matching SCIP constraint handler plugin.
    conshdlr_::Ptr{SCIP_CONSHDLR} = get(mscip.conshdlrs, ch, C_NULL)
    conshdlr_ != C_NULL || error("No matching constraint handler registered!")

    # Hand over Julia object as constraint data:
    consdata_ = pointer_from_objref(c)

    # Create SCIP constraint (and attach constraint data).
    cons__ = Ref{Ptr{SCIP_CONS}}(C_NULL)
    @SC SCIPcreateCons(mscip, cons__, "", conshdlr_, consdata_,
                       initial, separate, enforce, check, propagate,
                       _local, modifiable, dynamic, removable, stickingatnode)

    # Sanity check.
    @assert cons__[] != C_NULL

    # Register constraint data (for GC-protection).
    mscip.conshdlrconss[c] = cons__[]

    # Add constraint to problem.
    @SC SCIPaddCons(mscip, cons__[])

    # Register constraint and return reference.
    return store_cons!(mscip, cons__)
end


text:
mutable struct MemoryPoolHandle <: SEALObject
  handle::Ptr{Cvoid}

  function MemoryPoolHandle(handle::Ptr{Cvoid})
    object = new(handle)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::MemoryPoolHandle)
  if isallocated(object)
    @check_return_value ccall((:MemoryPoolHandle_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

function alloc_byte_count(handle::MemoryPoolHandle)
  count = Ref{UInt64}(0)
  retval = ccall((:MemoryPoolHandle_AllocByteCount, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 handle, count)
  @check_return_value retval
  return Int(count[])
end

function memory_manager_get_pool()
  handleref = Ref{Ptr{Cvoid}}(C_NULL)
  retval = ccall((:MemoryManager_GetPool2, libsealc), Clong,
                 (Ref{Ptr{Cvoid}},),
                 handleref)
  @check_return_value retval
  return MemoryPoolHandle(handleref[])
end


text:
# This file holds the definition of the Personnel type. This type is used to
#   define members of a military organisation, in this case, the Belgian Army.

# The Personnel type requires no additional types.
requiredTypes = [ ]

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

export Personnel
type Personnel
    # Information about the personnel member.
    persData::Dict{Symbol, Any}

    function Personnel( id::String = "-----" )
        newPers = new()
        newPers.persData = Dict{Symbol, Any}( :id => id )

        return newPers
    end  # Personnel( id )

    function Personnel( idKey::Symbol, id::String = "-----" )
        newPers = new()
        newPers.persData = Dict{Symbol, Any}( idKey => id )

        return newPers
    end  # Personnel( idKey, id )

    function Personnel( idKey::String, id::String )
        newPers = new()
        newPers.persData = Dict{Symbol, Any}( Symbol( idKey ) => id )

        return newPers
    end  # Personnel( idKey, id )
end  # type Personnel


text:
# This file was generated, do not modify it. # hide
y1 = predict(m1, X);
y2 = predict(m2, X);
X_judge = MLJ.table(hcat(y1, y2))
yhat = predict(m_judge, X_judge)


text:
function xavier_uniform!(t::Tensor, gain::Real=0)
    k = sqrt(6 / sum(size(t)))
    t.values = (rand(size(t)...) * 2k .- k) * gain
end

function zero!(t::Tensor)
    fill!(t.values, 0)
end

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
using IHT, SnpArrays, MendelBase, CSV, DataFrames

srand(2018)

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

function test_data()
	#dataset with 2 SNP and 6 people. The SNP matrix is 6x3 (with column of intercept)
	x = SnpData("test")
	y = CSV.read("test.fam", delim = ' ', header = false)
	y = convert(Array{Float64,1}, y[:, 6])
	k = 2
	v = IHTVariables(x, y, k)
	return (x, y, k, v)
end

function gwas1_data()
	#dataset with 10000 SNP and 2200 people. The SNP matrix is 2200x10001 (with column of intercept)
	x = SnpData("gwas 1 data")
	y = CSV.read("gwas 1 data_kevin.fam", delim = ',', header = false) # same file, comma separated
	y = convert(Array{Float64,1}, y[:, 6])
	k = 10
	v = IHTVariables(x, y, k)
	return (x, y, k, v)
end

@testset "initilize IHTVariables" begin
	(x, y, k, v) = test_data()

	#k must be an integer, but could be larger than
	@test_throws(ArgumentError, IHTVariables(x, y, 0))
	@test_throws(ArgumentError, IHTVariables(x, y, -1))
	@test_throws(MethodError, IHTVariables(x, y, 1.1))
	@test_throws(MethodError, IHTVariables(x, y, NaN))
	@test_throws(MethodError, IHTVariables(x, y, missing))
	@test_throws(MethodError, IHTVariables(x, y, Inf))

	#Different types of inputs for IHTVariables(x, y, k) is
	@test typeof(v) == IHTVariable
	@test typeof(x) == SnpData || typeof(x) <: SnpArray

	@test size(v.b)    == (3,)
	@test size(v.b0)   == (3,)
	@test size(v.xb)   == (6,)
	@test size(v.xb0)  == (6,)
	@test size(v.xk)   == (6, 2)
	@test size(v.gk)   == (2,)
	@test size(v.xgk)  == (6,)
	@test size(v.idx)  == (3,)
	@test size(v.idx0) == (3,)
	@test size(v.r)	   == (6,)
	@test size(v.df)   == (3,)

	@test typeof(v.b)    == Array{Float64, 1}
	@test typeof(v.b0)   == Array{Float64, 1}
	@test typeof(v.xb)   == Array{Float64, 1}
	@test typeof(v.xb0)  == Array{Float64, 1}
	@test typeof(v.xk)   == Array{Float64, 2}
	@test typeof(v.gk)   == Array{Float64, 1}
	@test typeof(v.xgk)  == Array{Float64, 1}
	@test typeof(v.idx)  == BitArray{1}
	@test typeof(v.idx0) == BitArray{1}
	@test typeof(v.r)	 == Array{Float64, 1}
	@test typeof(v.df)   == Array{Float64, 1}
end

@testset "_iht_indices" begin
	(x, y, k, v) = gwas1_data()

	# if v.idx is zero vector (i.e. first iteration of L0_reg), running _iht_indices should
	# set v.idx = 1 for the k largest terms in v.df
	v.df[1:10001] .= rand(10001)
	p = sortperm(v.df, rev = true)
	top_k_index = p[1:10]
	IHT._iht_indices(v, k)

	@test all(v.idx[top_k_index] .== 1)

	# if v.idx is not zero vector, then _iht_indices should find non-0 entries of b and
	# set v.idx = 1 for those entries
	v.b[1:k] .= rand(10)
	shuffle!(v.b)
	p = sortperm(v.b, rev = true)
	top_k_index = p[1:10]
	IHT._iht_indices(v, k)

	@test all(v.idx[top_k_index] .== 1)
end

@testset "project_k!" begin
    x = rand(100000)
    k = 100
    p = sortperm(x, rev = true)
    top_k_index = p[1:k]
	last_k_index = p[k+1:end]
	IHT.project_k!(x, k)

	@test all(x[top_k_index] .!= 0.0)
	@test all(x[last_k_index] .== 0.0)
end

@testset "compute_ω!" begin

end

@testset "use_A2_as_minor_allele" begin
	(x, y, k, v) = test_data()
    result = use_A2_as_minor_allele(x.snpmatrix)
    answer = [[0.0 1.0]; [1.0 1.0]; [2.0 0.0]; [1.0 2.0]; [2.0 1.0]; [2.0 2.0]]

    @test all(result .== answer)
end

@testset "_iht_backtrack" begin
    (x, y, k, v) = gwas1_data()
    μ, ω = 1.0, 1.0
    @test IHT._iht_backtrack(v, ω, μ) == false

    μ, ω = 0.2, 0.3
    @test IHT._iht_backtrack(v, ω, μ) == true

    μ, ω = 0.8, 0.792
    @test IHT._iht_backtrack(v, ω, μ) == false

    μ, ω = 0.5, 0.2
    @test IHT._iht_backtrack(v, ω, μ) == false

    μ, ω = 0.98, 1.0
    @test IHT._iht_backtrack(v, ω, μ) == true
end

@testset "_iht_gradstep" begin
    (x, y, k, v) = test_data()
    v.b .= rand(3)
    v.df .= rand(3)
    b = copy(v.b)
    df = copy(v.df)
    k = 2
    μ = 0.9

    IHT._iht_gradstep(v, μ, k)
    @test v.b[1] == 0.0 # because first entry is smallest, it should be set to 0
    @test v.b[2] == (b + μ*df)[2]
    @test v.b[3] == (b + μ*df)[3]
end


text:
using Documenter, BioSymbols

makedocs(
    format = Documenter.HTML(),
    sitename = "BioSymbols",
    pages = [
        "Home" => "index.md",
        "Nucleic Acids" => "nucleicacids.md",
        "Amino Acids" => "aminoacids.md",
        "Sequences" => "sequences.md",
        #"References" => "references.md",
        "Library" => [
            "Public" => "lib/public.md",
            hide("Internals" => "lib/internals.md", [
                "lib/internals/conversion-tables.md",
                "lib/internals/nucleic-acid-encoding.md"
            ])
        ]
    ],
    authors = "Ben J. Ward, The BioJulia Organisation and other contributors."
)

deploydocs(repo = "github.com/BioJulia/BioSymbols.jl.git")


text:
"""
`module  JAC.CloseCoupling`
    ... a submodel of JAC that contains methods for constructing (close-coupling) scattering states including the
        proper in-coming or outgoing boundary conditions.
"""
module CloseCoupling

    using Printf, ..AngularMomentum, ..Basics, ..Continuum, ..Defaults, ..ManyElectron, ..Radial

    export closeDummy

end # module


text:
mutable struct LazyOrderedDict{K,V} <: AbstractDict{K,V}
    # function to generate (key => value) for every element in list
    fun
    getindex_guess
    list
    data::Union{Nothing,OrderedDict{K,V}}
end

#LazyOrderedDict(load!) = LazyOrderedDict(load!,nothing)

Base.length(ld::LazyOrderedDict) = length(ld.list)

function load!(ld::LazyOrderedDict)
    if ld.data == nothing
        ld.data = OrderedDict(ld.fun.(ld.list)...)
    end
end

function Base.keys(ld::LazyOrderedDict)
    load!(ld)
    return collect(keys(ld.data))
end

function Base.getindex(ld::LazyOrderedDict{K,V},id::K) where {K,V}
    if ld.getindex_guess != nothing
        guess = ld.getindex_guess(ld.list,ld.fun,id)
        if guess != nothing
            return guess
        end
    end
    load!(ld)
    return ld.data[id]
end

#=
function Base.iterate(ld::LazyOrderedDict, state = keys(ld))
    if length(state) == 0
        return nothing
    end

    return (state[1] => ld[popfirst!(state)], state)
end
=#
function Base.iterate(ld::LazyOrderedDict, state = copy(ld.list))
    if length(state) == 0
        return nothing
    end

    return (ld.fun(popfirst!(state)), state)
end

for (item_color,default) in (
    (:title_color, :error_color),
    (:catalog_color, :info_color),
    (:item_color, :info_color),
    (:asset_color, :info_color),
)

    item_color_str = String(item_color)
    item_str = split(item_color_str,"_")[1]
    default_str = String(default)

    @eval begin
        $item_color = Ref(Symbol(load_preference(STAC,$(item_color_str), Base.$default())))

        """
        STAC.set_$($item_color_str)(color::Symbol)

Set the $($item_str) color. The default color is `Base.$($default_str)()`.
"""
        function $(Symbol(:set_,item_color))(color::Symbol)
            @set_preferences!($(item_color_str) => String(color))
            $item_color[] = color
        end

    end
end

# print unless one argument is nothing
function _print(args...; kwargs...)
    if !any(isnothing,args)
        print(args...; kwargs...)
    end
end

function _printstyled(args...; kwargs...)
    if !any(isnothing,args)
        printstyled(args...; kwargs...)
    end
end

_assets(data) = OrderedDict((String(k) => Asset(v) for (k,v) in get(data,:assets,[])))

function _show_assets(io,item)
    ident = "  "

    if length(item.assets) > 0
        println(io,"Assets:")
        for (id,asset) in item.assets
            print(io,ident," * ")
            printstyled(io, id, color=asset_color[])

            _print(io,": ",title(asset),"")
            printstyled(io, "\n")
        end
    end

end


text:
# Functions for preprocessing PDDL constructs into canonical forms

"Preprocess domain by instantiating and regularizing logical formulae."
function preprocess(domain::Domain, problem::Union{Problem,Nothing}=nothing;
                    options...)
    domain = copy(domain)
    # Get object type declarations
    objtypes = get(options, :objtypes, nothing)
    if !isnothing(problem) && isnothing(objtypes)
        objtypes = get_obj_clauses(problem)
    end
    # Unpack flags
    regularize = get(options, :regularize, true)
    instantiate = get(options, :instantiate, !isnothing(objtypes))
    # Remove universal quantifiers in axioms
    if instantiate && !isnothing(objtypes)
        domain.axioms = [deuniversalize(ax, objtypes) for ax in domain.axioms]
    end
    # Regularize all axiom bodies to conjunctions of literals
    if regularize
        domain.axioms = regularize_clauses(domain.axioms)
    end
    # Preprocess action and event definitions
    domain.actions = Dict{Symbol,Action}(
        k => preprocess(v, objtypes; options...) for (k, v) in domain.actions)
    domain.events = [preprocess(e, objtypes; options...) for e in domain.events]
    return domain
end

function preprocess(act::Action, objtypes=nothing; options...)
    regularize = get(options, :regularize, true)
    instantiate = get(options, :instantiate, !isnothing(objtypes))
    precond, effect = act.precond, act.effect
    if instantiate && !isnothing(objtypes)
        precond = deuniversalize(precond, objtypes)
        effect = deuniversalize(effect, objtypes)
    end
    if regularize
        precond = to_dnf(precond)
    end
    return Action(act.name, act.args, act.types, precond, effect)
end

function preprocess(event::Event, objtypes=nothing; options...)
    regularize = get(options, :regularize, true)
    instantiate = get(options, :instantiate, !isnothing(objtypes))
    precond, effect = event.precond, event.effect
    if instantiate && !isnothing(objtypes)
        precond = deuniversalize(precond, objtypes)
        effect = deuniversalize(effect, objtypes)
    end
    if regularize
        precond = to_dnf(precond)
    end
    return Event(event.name, precond, effect)
end


text:
# A version of the conjugate gradient method for linesearch.

"""A version of the conjugate gradient method for linesearch.
CG(A, b, ϵa, ϵr, itmax; quad) solves the symmetric linear system 'A * x = b'
If quad = true, the values of the quadratic model are computed
A can be positive definite or not.
"""
function CGlin(A, b, ϵa::Float64=1e-8, ϵr::Float64=1e-6, itmax::Int=0, ε::Float64=1e-6; quad::Bool=false)
    n = size(b, 1) # size of the problem
    (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size")
    @info(loggerCGlin, @sprintf("CGlin: system of %d equations in %d variables", n, n))

    x = zeros(n) # initial estimation x = 0
    xNorm = 0.0
    xNorms = [xNorm] # Values of ‖x‖
    r = -b # initial residual r = Ax-b = -b
    d = b # first descent direction
    rNorm = norm(r, 2)
    ϵ = ϵa + ϵr * rNorm

    iter = 0
    itmax == 0 && (itmax = 2 * n)

    if quad
        q = 0.0
        qvalues = [q] # values of the quadratic model
        @info(loggerCGlin, @sprintf("%5s %10s %10s\n", "Iter", "‖r‖", "q"))
        @info(loggerCGlin, @sprintf("    %d    %8.1e    %8.1e", iter, rNorm, q))
    end

    solved = rNorm ≤ ϵ
    tired = iter ≥ itmax

    while ! (solved || tired)
        iter += 1
        Ad = A * d
        dAd = dot(d, Ad)

        # if the model is not convexe, the algorithm stops
        if dAd ≤ ε * dot(d, d)
            @debug(loggerCGlin, @sprintf("non positive curvature dAd = %8.1e", dAd))
            iter == 1 && return b
            return x
        end

        α = rNorm^2 / dAd # step for x estimation
        x = x + α * d # new estimation
        xNorm = norm(x, 2)
        push!(xNorms, xNorm)
        roldNorm = rNorm
        r = r + α * Ad # new residual
        rNorm = norm(r, 2)

        if quad
            q = -dot(b, x) + 0.5 * dot(x, A * x)
            push!(qvalues, q)
            @info(loggerCGlin, @sprintf("    %d    %8.1e    %8.1e", iter, rNorm, q))
        end

        solved = rNorm ≤ ϵ
        tired = iter ≥ itmax

        (solved || tired) && continue
        β = rNorm^2 / roldNorm^2 # step for the next descent direction
        d = -r + β * d # new descent direction

        end
        return x
    end


text:
module HelloWorldApp

using Matte

const title = "Matte Example :: Hello World"

function ui()
    sidebar_layout(
        side_panel(
            h1("Hello World!"),
            text_input("my_input", "Input some text")
        ),
        main_panel(
            h1("You wrote:"),
            p(text_output("my_output"))
        )
    )
end

module Server

function my_output(my_input)
    my_input
end

end

end


text:
"""
    tdvp1sitesweep!(dt, mps, mpo, env [; verbose])

Performs a single sweep (left to right and back) of the 1site TDVP on
`mps` where the `mpo` is the Hamiltonian and environment matrices are
stored in `env`. This process symplectically evolves the `mps` in time
for the times step `dt`.

See arXiv:1408.5056 for details of the algorithm and explanations

"""
function tdvp1sitesweep!(dt::Float64,
                         mps::MPS{T},
                         mpo::MPO{T},
                         env::Vector{Array{T, 3}};
                         verbose::Bool=false) where {T<:Number}

    !(T <: Complex) && error("TDVP only accpets complex MPSs. Convert first!")

    lx = length(mps)

    center_at!(mps, 1)

    A = mps.As[1]
    for l = 1:lx-1
        # forward evolution of mps at site l
        A, info = exponentiate(v->_applymps1site(v, env[l], env[l+2], mpo.Ws[l]),
                               -im*dt, A, ishermitian=true)

        if verbose
            e = dot(A, _applymps1site(A, env[l], env[l+2], mpo.Ws[l]))
            println("Sweep L2R: mps site $l -> energy $e")
        end

        Q, Λ = qr(reshape(A, size(A,1)*size(A,2), size(A,3)))
        mps.As[l] = reshape(Matrix(Q), size(A))
        env[l+1] = _mpsupdateleft(env[l], mps.As[l], mpo.Ws[l])

        # backward evolution of Λ
        Λ, info = exponentiate(v->_applymps0site(v, env[l+1], env[l+2]),
                               +im*dt, Λ)

        if verbose
            e = dot(Λ, _applymps0site(Λ, env[l+1], env[l+2]))
            println("Sweep L2R: Λ between  site $l and $(l+1) -> energy $e")
        end

        @tensor A[l,o,r] := Λ[l,m] * mps.As[l+1][m,o,r]
    end

    l = lx
    A, info = exponentiate(v->_applymps1site(v, env[l], env[l+2], mpo.Ws[l]),
                           -im*dt, A, ishermitian=true)

    if verbose
        e = dot(A, _applymps1site(A, env[l], env[l+2], mpo.Ws[l]))
        println("Sweep L2R: mps site $l -> energy $e")
    end

    for l = lx-1:-1:1
        Q, Λ = qr(transpose(reshape(A, size(A, 1), size(A,2)*size(A,3))))
        Λ = Matrix(transpose(Λ))

        mps.As[l+1] = reshape(transpose(Matrix(Q)), size(A))
        env[l+2] = _mpsupdateright(env[l+3], mps.As[l+1], mpo.Ws[l+1])

        # backward evolution of Λ
        Λ, info = exponentiate(v->_applymps0site(v, env[l+1], env[l+2]),
                               +im*dt, Λ)
        if verbose
            e = dot(Λ, _applymps0site(Λ, env[l+1], env[l+2]))
            println("Sweep L2R: Λ between  site $l and $(l+1) -> energy $e")
        end

        @tensor A[l,o,r] := mps.As[l][l,o,m] * Λ[m,r]

        # forward evolution of mps at site l
        A, info = exponentiate(v->_applymps1site(v, env[l], env[l+2], mpo.Ws[l]),
                               -im*dt, A; ishermitian = true)

        if verbose
            e = dot(A, _applymps1site(A, env[l], env[l+2], mpo.Ws[l]))
            println("Sweep R2L: mps site $l -> energy $e")
        end
    end
    mps.As[1] = A
    mps.center = 1
    mps
end

"""
    tdvp2sitesweep!(dt, mps, mpo, env [; maxdim, tol, verbose])

Performs a single sweep (left to right and back) of the 2site TDVP on
`mps` where the `mpo` is the Hamiltonian and environment As are
stored in `env`. This process is not a result of direct formulation of
tangent spaces, because there are now a hierarchy of manifold of MPS
with different system sizes, but it is a direct suggestion from the
one-site version.

See arXiv:1408.5056 for details of the algorithm and explanations/
justification of the two-site proposed algorithm.

"""
function tdvp2sitesweep!(dt::Float64,
                         mps::MPState{Array{T, 3}},
                         mpo::MPOperator{Array{T, 4}},
                         env::Vector{Array{T, 3}};
                         maxdim::Int=200,
                         tol::Float64=1.e-14,
                         verbose::Bool=false) where {T<:Number}

    !(T <: Complex) && error("TDVP only accpets complex MPSs. Convert first!")
    lx = length(mps)

    center_at!(mps, 1)

    A = mps.As[1]
    for l = 1:lx-2
        # forward evolution of mps at site l and site l+1
        @tensor AA[-1,-2,-3,-4] := A[-1,-2,1] * mps.As[l+1][1,-3,-4]
        AA, info = exponentiate(v->_applymps2site(v, env[l], env[l+3],
                                                  mpo.Ws[l], mpo.Ws[l+1]),
                                -im*dt, AA, ishermitian=true)

        if verbose
            e = dot(AA, _applymps2site(AA, env[l], env[l+3],
                                       mpo.Ws[l], mpo.Ws[l+1]))
            println("Sweep L2R: mps site $l, $(l+1) -> energy $e")
        end

        u,s,v = svdtrunc(reshape(AA, prod(size(AA)[1:2]), prod(size(AA)[3:4])),
                         maxdim=maxdim, tol=tol)
        mps.As[l] = reshape(Matrix(u), (size(AA)[1:2]..., size(s,1)))
        env[l+1] = _mpsupdateleft(env[l], mps.As[l], mpo.Ws[l])

        # backward evolution of Λ at site l+1
        Λ = reshape(s*v, (size(s,1), size(AA)[3:4]...))
        Λ, info = exponentiate(v->_applymps1site(v, env[l+1], env[l+3], mpo.Ws[l+1]),
                               +im*dt, Λ)

        if verbose
            e = dot(Λ, _applymps1site(Λ, env[l+1], env[l+3], mpo.Ws[l+1]))
            println("Sweep L2R: Λ on site $(l+1) -> energy $e")
        end

        A = Λ
        #@tensor A[l,o,r] := Λ[l,m] * mps.As[l+1][m,o,r]
    end

    l = lx-1
    @tensor AA[-1,-2,-3,-4] := A[-1,-2,1] * mps.As[l+1][1,-3,-4]
    AA, info = exponentiate(v->_applymps2site(v, env[l], env[l+3],
                                              mpo.Ws[l], mpo.Ws[l+1]),
                            -im*dt, AA, ishermitian=true)

    if verbose
        e = dot(AA, _applymps2site(AA, env[l], env[l+3],
                                   mpo.Ws[l], mpo.Ws[l+1]))
        println("Sweep L2R: mps site $l, $(l+1) -> energy $e")
    end

    for l = lx-1:-1:2

        u, s, v = svdtrunc(reshape(AA, prod(size(AA)[1:2]), prod(size(AA)[3:4])),
                           maxdim=maxdim, tol=tol)
        mps.As[l+1] = reshape(v, (size(s,1),size(AA)[3:4]...))
        env[l+2] = _mpsupdateright(env[l+3], mps.As[l+1], mpo.Ws[l+1])

        # backward evolution of Λ at site l
        Λ = reshape(u*s, (size(AA)[1:2]...,size(s,2)))
        Λ, info = exponentiate(v->_applymps1site(v, env[l], env[l+2], mpo.Ws[l]),
                               +im*dt, Λ)
        if verbose
            e = dot(Λ, _applymps1site(Λ, env[l], env[l+2], mpo.Ws[l]))
            println("Sweep R2L: Λ on site $l -> energy $e")
        end

        A = Λ
        #@tensor A[l,o,r] := mps.As[l-1][l,o,m] * Λ[m,r]
        @tensor AA[-1,-2,-3,-4] := mps.As[l-1][-1,-2,1] * A[1,-3,-4]

        # forward evolution of mps at site l
        AA, info = exponentiate(v->_applymps2site(v, env[l-1], env[l+2],
                                                  mpo.Ws[l-1], mpo.Ws[l]),
                                -im*dt, AA; ishermitian = true)

        if verbose
            e = dot(AA, _applymps2site(AA, env[l-1], env[l+2],
                                       mpo.Ws[l-1], mpo.Ws[l]))
            println("Sweep R2L: mps site $(l-1), $(l) -> energy $e")
        end
    end
    l=1
    u, s, v = svdtrunc(reshape(AA, prod(size(AA)[1:2]), prod(size(AA)[3:4])),
                       maxdim=maxdim, tol=tol)
    mps.As[l+1] = reshape(v, (size(s,1),size(AA)[3:4]...))
    env[l+2] = _mpsupdateright(env[l+3], mps.As[l+1], mpo.Ws[l+1])

    Λ = reshape(u*s, (size(AA)[1:2]..., size(s,2)))
    ### (IMPORTANT) QUESTION? do I need a final backward evolution on site 1?
    # Λ, info = exponentiate(v->_applymps1site(v, env[l], env[l+2], mpo.Ws[l]),
    #                        +im*dt, Λ)

    mps.As[1] = Λ
    mps.center = 1

    mps
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
#=
Created on 22/01/2021 14:26:51
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Example from the main documentation about interphases.
Illustrates iterator over the squares for integer values.
=#

struct Squares
    count::Int
end

Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

for i in Squares(7)
    println(i)
end

25 in Squares(10)

sum(Squares(1803))

Base.eltype(::Type{Squares}) = Int

Base.length(S::Squares) = S.count

collect(Squares(4))

Base.sum(S::Squares) = (n = S.count; return n*(n+1)*(2n+1)÷6)

sum(Squares(1803))


text:
@test_throws VulkanError Instance(["nonexistent_layer"], [])
@test_throws VulkanError Instance([], ["nonexistent_extension"])
@test unwrap_error(create_instance([], ["nonexistent_extension"])).code == ERROR_EXTENSION_NOT_PRESENT

const instance = Instance(INSTANCE_LAYERS, INSTANCE_EXTENSIONS;
    application_info = ApplicationInfo(v"0.0.1", v"0.0.1", API_VERSION; application_name="Test", engine_name="Experimental engine"))
@info "Instance created: $instance."

if WITH_DEBUG[]
    debug_messenger = DebugUtilsMessengerEXT(instance, debug_callback_c)
end

const device = let pdevices = unwrap(enumerate_physical_devices(instance))
    if isempty(pdevices)
        error("No physical devices available for testing.")
    end

    @info "$(length(pdevices)) physical device(s) found."

    pdevice = first(pdevices)

    @info "Selected $(get_physical_device_properties(pdevice))"
    @test_throws ErrorException("No queue with the desired capabilities could be found.") find_queue_family(pdevice, typemax(QueueFlag))
    @test_throws ErrorException("Invalid physical device features: no_feature") PhysicalDeviceFeatures(:no_feature)

    Device(
        pdevice,
        [DeviceQueueCreateInfo(find_queue_family(pdevice, QUEUE_GRAPHICS_BIT & QUEUE_COMPUTE_BIT), [1.0])],
        [], DEVICE_EXTENSIONS; enabled_features = ENABLED_FEATURES
    )
end

@info "Logical device created: $device."

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Lehmann

eta = 1e-3          # noise level

println("Generate DLR ...")
d = DLRGrid(Euv = 1.0, β = 1000.0, isFermi = false) # Initialize DLR object
# β is the inverse temperature, Euv is the ultraviolt energy cutoff of the Green's function. By definition, Λ = Euv*β

# Generate noisy data
tau_i = collect(LinRange(0, d.β, 100))

println("Prepare the Green's function sample ...")
G_i = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k
# G_i = true_G_tau(tau_i, beta)

G_i_noisy = G_i .+ eta .* (rand(length(G_i)) .- 0.5)

# Fit DLR coeffs from noisy data

G_x = tau2dlr(d, G_i_noisy, tau_i)

# Evaluate DLR

G_i_fit = dlr2tau(d, G_x, tau_i)

diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
println("Maximum difference: ", diff(G_i_noisy, G_i_fit))
# for i in 1:length(G_i_fit)
#     println("$i   $(G_i_noisy[i])         $(G_i_fit[i])")
# end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random, OrdinaryDiffEq

a = 0.2
b = 0.2
c = 5.7
amu = 0.02
ω1 = 1. + amu
ω2 = 1. - amu
kl = 201
kmin = 0
kmax = 0.05

config = @strdict ω1 ω2 a b c kmin kl kmax

function produce_croesslers_spectrum(config)

@unpack ω1, ω2, a, b, c, kmin, kl, kmax = config

y_init = [0.0, 0.2, 0., 0.11, 0.19, 0.1]

croesslers = Systems.coupled_roessler(y_init; ω1, ω2, a, b, c)

diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9, maxiters = Inf)

Δt = 0.01
T = 1000
t = 0:Δt:T

krange = range(kmin, kmax; length = kl)

λs = zeros(length(krange), 4)
sync_error= zeros(length(krange))

for (i, k) in enumerate(krange)
    @show i, k
    set_parameter!(croesslers, 6, k)
    set_parameter!(croesslers, 7, k)
    λs[i, :] .= lyapunovspectrum(croesslers, 100000, 4; Ttr = 5000, diffeq...)

    # computation of phase difference -----
    T = 20000 # length of time series
    Δt = 0.1
    X = trajectory(croesslers, T; Ttr = 5000, Δt, diffeq...)  # which initial conditions are used?
    t = 0:Δt:T
    x1vec = X[:, 1]
    x2vec = X[:, 2]
    x4vec = X[:, 4]
    x5vec = X[:, 5]

    theta1 = 0.
    theta1_old = atan.( x2vec[1] , x1vec[1] )
    theta2 = 0.
    theta2_old = atan.( x5vec[1] , x4vec[1] )

    nts = length(t)
    for its = 2:nts
        theta1_new = atan.( x2vec[its] , x1vec[its]  )    # is decreasing = clockwise rotation
        delta_theta1 = mod( theta1_new  - theta1_old  + pi, 2 .* pi) - pi

        theta1 = theta1 + delta_theta1
        theta1_old = theta1_new

        theta2_new = atan.( x5vec[its] , x4vec[its]  )    # is decreasing = clockwise rotation
        delta_theta2 = mod( theta2_new  - theta2_old  + pi, 2 .* pi) - pi
        theta2 = theta2 + delta_theta2
        theta2_old = theta2_new
    end

    sync_error[i]  = abs(theta1 - theta2) / (t[end] - t[1])
end

ret = @strdict λs sync_error
ret = merge(config, ret)
return ret
end

file, path = produce_or_load(
    datadir("croesslers"), config, produce_croesslers_spectrum;
    prefix = "croesslers_spectrum", tag = false,
)

@unpack λs, sync_error = file

# %% Plot
fig = figure(figsize = (0.66figx, 1.2figy))
gs = matplotlib.gridspec.GridSpec(2, 1, height_ratios=[2,1])
ax1 = plt.subplot(gs[1])
ax2 = plt.subplot(gs[2])

ax1.plot(krange, λs)
# ax.set_xlabel("\$k\$")
ax1.set_xticklabels([])
ax1.set_ylabel("\$\\lambda_i\$")
ax1.set_xlim(kmin,kmax)
ax1.set_ylim(-0.07,0.08)

ax2.plot(krange, sync_error)
ax2.set_xlabel("\$k\$"; labelpad = -20)
ax2.set_ylabel("\$\\Delta \\Omega\$")
ax2.set_xlim(krange[1],krange[end])

add_identifiers!(fig)

fig.tight_layout(;pad = 0.33)
wsave(plotsdir("9", "coupled_roessler_spectrum"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
"""
Creates a GridLayout with all pairs contained in `content`. Each pair consists
of an iterable with row and column spans, and a content object. Each content
object is then placed in the GridLayout at the span from its pair.

Example:

grid!(
    [1, 1] => obj1,
    [1, 2] => obj2,
    [2, :] => obj3,
)
"""
function grid!(content::Vararg{Pair}; kwargs...)
    g = GridLayout(; kwargs...)
    for ((rows, cols), element) in content
        g[rows, cols] = element
    end
    g
end

"""
Creates a single-row GridLayout with all elements contained in `content` placed
from left to right.
"""
function hbox!(content::Vararg; kwargs...)
    ncols = length(content)
    g = GridLayout(1, ncols; kwargs...)
    for (i, element) in enumerate(content)
        g[1, i] = element
    end
    g
end

"""
Creates a single-column GridLayout with all elements contained in `content` placed
from top to bottom.
"""
function vbox!(content::Vararg; kwargs...)
    nrows = length(content)
    g = GridLayout(nrows, 1; kwargs...)
    for (i, element) in enumerate(content)
        g[i, 1] = element
    end
    g
end

"""
Creates a GridLayout in Outside align mode whose bounding box is tied to the pixelarea of `scene` minus
a padding of `padding` pixels.
"""
function grid!(scene::Scene, padding::Real, args...; kwargs...)
    grid!(args...; bbox = lift(x -> BBox(x), pixelarea(scene)), alignmode = Outside(padding), kwargs...)
end

"""
Creates a single-row GridLayout in Outside align mode whose bounding box is tied to the pixelarea
of `scene` minus a padding of `padding` pixels.
"""
function hbox!(scene::Scene, padding::Real, args...; kwargs...)
    hbox!(args...; bbox = lift(x -> BBox(x), pixelarea(scene)), alignmode = Outside(padding), kwargs...)
end

"""
Creates a single-column GridLayout in Outside align mode whose bounding box is tied to the pixelarea
of `scene` minus a padding of `padding` pixels.
"""
function vbox!(scene::Scene, padding::Real, args...; kwargs...)
    vbox!(args...; bbox = lift(x -> BBox(x), pixelarea(scene)), alignmode = Outside(padding), kwargs...)
end

"""
Creates a GridLayout filled with matrix-like content. The size of the grid will
be the size of the matrix.
"""
function grid!(content::AbstractMatrix; kwargs...)
    nrows, ncols = size(content)
    g = GridLayout(nrows, ncols; kwargs...)
    for i in 1:nrows, j in 1:ncols
        g[i, j] = content[i, j]
    end
    g
end

function layoutscene(padding = 30; kwargs...)
    scene = Scene(; camera = campixel!, kwargs...)
    gl = GridLayout(scene, alignmode = Outside(padding))
    scene, gl
end

function layoutscene(nrows::Int, ncols::Int, padding = 30; kwargs...)
    scene = Scene(; camera = campixel!, kwargs...)
    gl = GridLayout(scene, nrows, ncols, alignmode = Outside(padding))
    scene, gl
end


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

JLLWrappers.@generate_main_file_header("CUDNN")
JLLWrappers.@generate_main_file("CUDNN", UUID("62b44479-cb7b-5706-934f-f13b2eb2e645"))
end  # module CUDNN_jll


text:
get_scalefactor(lb, ub) = max(maximum(abs.(lb)), maximum(abs.(ub)))
get_scalefactor(model::MetNet) = get_scalefactor(model.lb, model.ub)