text:
abstract type BC end

struct NoPBC <:BC end

struct PBC <:BC
    dim::Int64
    Lx::Float64
    Ly::Float64
    Lz::Float64
    xhi::Float64
    xlo::Float64
    yhi::Float64
    ylo::Float64
    zhi::Float64
    zlo::Float64
    function PBC(   ;xhi=0.0,
                    xlo=0.0,
                    yhi=0.0,
                    ylo=0.0,
                    zhi=0.0,
                    zlo=0.0)
        Lx=xhi-xlo
        Ly=yhi-ylo
        Lz=zhi-zlo
        dim=3
        new(dim::Int64,
            Lx::Float64,
            Ly::Float64,
            Lz::Float64,
            xhi::Float64,
            xlo::Float64,
            yhi::Float64,
            ylo::Float64,
            zhi::Float64,
            zlo::Float64

        )
    end
end


text:
# FixMe! Right now the MPI workers are deduced from the DArrays, but if a DArray is distributed on fewer workers that what consistutes the MPI world, then this approach will fail.

type RemoteElementalMatrix
    refs::Matrix{Any}
end

function toback{T<:BlasFloat,S<:StridedMatrix}(A::DArray{T,2,S})
    rs = Array(Any, size(procs(A)))
    @sync for p in eachindex(procs(A))
        ind = A.indexes[p]
        @async rs[p] = remotecall(procs(A)[p]) do
            lA = localpart(A)
            AlA = Elemental.DistMatrix(T)
            zeros!(AlA, size(A)...)
            for j = 1:size(lA,2), i = 1:size(lA, 1)
                queueUpdate(AlA,
                            start(ind[1]) + i - 1,
                            start(ind[2]) + j - 1, lA[i,j])
            end
            processQueues(AlA)
            AlA
        end
    end
    return rs
end

function tofront(r::Base.Matrix)
    tt = Array(Any, length(r))
    for i = 1:length(r)
        tt[i] = remotecall(r[i].where, r[i]) do t
            typeof(fetch(t))
        end
    end

    rType = fetch(tt[1])
    if rType <: ElementalMatrix
        for i = 1:length(r)
            tt[i] = remotecall(r[i].where, r[i]) do t
                v = fetch(t)
                (size(v), eltype(v))
            end
        end
        sizeAndEltype = fetch(tt[1])
        if isa(sizeAndEltype, Exception)
            throw(sizeAndEltype)
        end

        mn = sizeAndEltype[1]
        A = dzeros(sizeAndEltype[2], mn,
                   Int[r[i].where for i in eachindex(r)])

        @sync for p in eachindex(r)
            ind = A.indexes[p]
            rr = r[p]
            @async remotecall_wait(r[p].where) do
                rrr = fetch(rr)
                lA = localpart(A)
                for j = 1:size(lA, 2), i = 1:size(lA, 1)
                    queuePull(rrr,
                              start(ind[1]) + i - 1,
                              start(ind[2]) + j - 1)
                end
                processPullQueue(rrr, lA)
            end
        end
    elseif rType <: Exception
        throw(fetch(r[1]))
    else
        error("Remote type was $rType and is not handled yet")
    end
    return A
end

function (\){T<:BlasFloat,S}(A::DArray{T,2,S}, B::DArray{T,2,S})
    rA = toback(A)
    rB = toback(B)
    pidsAB = union(A.pids, B.pids)
    rvals = Array(Any, length(pidsAB))
    @sync for i = 1:length(pidsAB)
        @async rvals[i] = remotecall_wait(pidsAB[i], rA[i], rB[i]) do t1,t2
            solve!(fetch(t1), fetch(t2))
        end
    end
    return tofront(reshape(rvals, size(procs(B))))
end

function eigvals{T<:BlasFloat}(A::Hermitian{T,DArray{T,2,Array{T,2}}})
    rA = toback(A.data)
    rvals = Array(Any, size(procs(A.data)))
    uplo = A.uplo == 'U' ? UPPER : LOWER
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i].where, rA[i]) do t
            eigvalsHermitian(uplo, fetch(t))
        end
    end
    return tofront(rvals)
end

function svdvals{T<:BlasFloat}(A::DArray{T,2})
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i].where, rA[i]) do t
            svdvals(fetch(t))
        end
    end
    return tofront(rvals)
end

function inv!{T<:BlasFloat}(A::DArray{T,2})
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for j = 1:size(rvals, 2)
        for i = 1:size(rvals, 1)
            @async rvals[i,j] = remotecall_wait(t -> inverse!(fetch(t)), rA[i,j].where, rA[i,j])
        end
    end
    return tofront(rvals)
end

inv{T<:BlasFloat}(A::DArray{T,2}) = inv!(copy(A))

function logdet{T<:BlasFloat}(A::DArray{T,2})
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i].where, rA[i]) do t
            d = safeHPDDeterminant(Elemental.LOWER, fetch(t))
            return log(d.ρ) + d.κ*d.n
        end
    end
    return fetch(rvals[1])
end

function spectralPortrait{T<:BlasReal}(A::DArray{T,2},
                                       realSize::Integer,
                                       imagSize::Integer,
                                       psCtrl::PseudospecCtrl{T}=PseudospecCtrl(T))
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i].where, rA[i]) do t
            spectralPortrait(fetch(t), ElInt(realSize), ElInt(imagSize), psCtrl)[1]
        end
    end
    return tofront(rvals)
end

function spectralPortrait{T<:BlasReal}(A::DArray{Complex{T},2},
                                       realSize::Integer,
                                       imagSize::Integer,
                                       psCtrl::PseudospecCtrl{T}=PseudospecCtrl(T))
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i,j] = remotecall_wait(rA[i].where, rA[i]) do t
            spectralPortrait(fetch(t), ElInt(realSize), ElInt(imagSize), psCtrl)[1]
        end
    end
    return tofront(rvals)
end

function spectralWindow{T<:BlasReal}(A::DArray{T,2},
                                     center::Complex{T},
                                     realWidth::T,
                                     imagWidth::T,
                                     realSize::Integer,
                                     imagSize::Integer,
                                     psCtrl::PseudospecCtrl{T}=PseudospecCtrl(T))
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i].where, rA[i]) do t
            spectralWindow(fetch(t), center, realWidth, imagWidth,
                           ElInt(realSize), ElInt(imagSize), psCtrl)
        end
    end
    return tofront(rvals)
end

function spectralWindow{T<:BlasReal}(A::DArray{Complex{T},2},
                                     center::Complex{T},
                                     realWidth::T,
                                     imagWidth::T,
                                     realSize::Integer,
                                     imagSize::Integer,
                                     psCtrl::PseudospecCtrl{T}=PseudospecCtrl(T))
    rA = toback(A)
    rvals = Array(Any, size(procs(A)))
    @sync for i in eachindex(rvals)
        @async rvals[i] = remotecall_wait(rA[i,j].where, rA[i]) do t
            spectralWindow(fetch(t), center, realWidth, imagWidth,
                           ElInt(realSize), ElInt(imagSize), psCtrl)
        end
    end
    return tofront(rvals)
end

function foxLi{T<:BlasComplex}(::Type{T}, n::Integer, ω::Real)
    sz = tuple(DistributedArrays.defaultdist((n,n), workers())...)
    rvals = Array(Any, sz)
    @sync for j = 1:size(rvals, 2), i = 1:size(rvals, 1)
        @async rvals[i,j] = remotecall_wait(workers()[sub2ind(sz, i, j)]) do
            A = Elemental.DistMatrix(T)
            foxLi!(A, ElInt(n), real(T)(ω))
        end
    end
    return tofront(rvals)
end
foxLi(n::Integer, ω::Real) = foxLi(Complex128, n, ω)

# Andreas: Just saw this one. It is almost identical to the one I wrote above,
# but I don't think that we can return a Elemental array beacause it has to
# live on the MPI cluster and cannot live on the "front end".
# It is assumed that the DArray is distributed over MPI.COMM_WORLD
for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:Complex64, :c),
                    (:Complex128, :z))
    @eval begin
        function convert(::Type{DistSparseMatrix{$elty}}, DA::DistributedArrays.DArray)
            npr, npc = size(procs(DA))
            if npr*npc != MPI.Comm_size(MPI.COMM_WORLD)
                error("Used non MPI.COMM_WORLD DArray for DistSparseMatrix, ",
                      "as procs(DA)=($npr,$npc) is incompatible with ",
                      "MPI.Comm_size(MPI.COMM_WORLD)=$(MPI.Comm_size(MPI.COMM_WORLD))")
            end

            m, n = size(DA)
            A = DistSparseMatrix($elty, m, n)
            @sync begin
                for id in workers()
                    let A = A, DA = DA
                        @async remotecall_fetch(id) do
                            rows, cols = DistributedArrays.localindexes(DA)
                            i,j,v = findnz(DistributedArrays.localpart(DA))
                            gi, gj, gv = (i.+(first(rows)-1), j.+(first(cols)-1), v)
                            numLocal = length(gi)
                            reserve(A,numLocal)
                            for s=1:numLocal
                                queueUpdate(A,gi[s],gj[s],v[s])
                            end
                        end
                    end
                end
            end
            processQueues(A)
            return A
        end
    end
end

ext: jl


text:
using CSV, DataFrames
using DimensionalData
using StanSample
using Statistics, Test

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), DataFrame);

# Define the Stan language model

stan10_4 = "
data{
    int N;
    int N_actors;
    int pulled_left[N];
    int prosoc_left[N];
    int condition[N];
    int actor[N];
}
parameters{
    vector[N_actors] a;
    real bp;
    real bpC;
}
model{
    vector[N] p;
    bpC ~ normal( 0 , 10 );
    bp ~ normal( 0 , 10 );
    a ~ normal( 0 , 10 );
    for ( i in 1:504 ) {
        p[i] = a[actor[i]] + (bp + bpC * condition[i]) * prosoc_left[i];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
";

data10_4 = (N = size(df, 1), N_actors = length(unique(df.actor)),
    actor = df.actor, pulled_left = df.pulled_left,
    prosoc_left = df.prosoc_left, condition = df.condition);

# Sample using cmdstan

m10_4s = SampleModel("m10.4s", stan10_4)
rc10_4s = stan_sample(m10_4s; data=data10_4);

# Result rethinking

rethinking = "
      mean   sd  5.5% 94.5% n_eff Rhat
bp    0.84 0.26  0.43  1.26  2271    1
bpC  -0.13 0.29 -0.59  0.34  2949    1

a[1] -0.74 0.27 -1.16 -0.31  3310    1
a[2] 10.88 5.20  4.57 20.73  1634    1
a[3] -1.05 0.28 -1.52 -0.59  4206    1
a[4] -1.05 0.28 -1.50 -0.60  4133    1
a[5] -0.75 0.27 -1.18 -0.32  4049    1
a[6]  0.22 0.27 -0.22  0.65  3877    1
a[7]  1.81 0.39  1.22  2.48  3807    1
";

# Update sections

@testset ":dimarrays tests" begin
    if success(rc10_4s)
        da = read_samples(m10_4s, :dimarrays)
        da1 = da[param=At(Symbol("a.1"))]

        # Other manipulations

        @test Tables.istable(da) == true

        # All of parameters
        dar = reshape(Array(da), 4000, 9);
        @test size(dar) == (4000, 9)

        # Check :param axis names
        @test dims(da, :param).val == vcat([Symbol("a.$i") for i in 1:7], :bp, :bpC)

        # Test combining vector param 'a'
        ma = matrix(da, "a");
        rma = reshape(ma, 4000, size(ma, 3))
        @test mean(rma, dims=1) ≈ [-0.7 10.9 -1 -1 -0.7 0.2 1.8] atol=0.7
    end
end

@testset ":dimarray tests" begin
    if success(rc10_4s)
        da = read_samples(m10_4s, :dimarray)
        da |> display
        da1 = da[param=At(Symbol("a.1"))]

        # Other manipulations

        @test Tables.istable(da) == true

        # Check :param axis names
        @test dims(da, :param).val == vcat([Symbol("a.$i") for i in 1:7], :bp, :bpC)

        # Test combining vector param 'a'
        ma = matrix(da, "a");
        @test mean(ma, dims=1) ≈ [-0.7 10.9 -1 -1 -0.7 0.2 1.8] atol=0.7
    end
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
using JuMP
using LightGraphs
using MetaGraphs

using GraphJuMP

using Base.Test
using Cbc

@testset "GraphJuMP" begin
  @testset "Integration tests" begin
    @testset "Edge flow formulation" begin
      @testset "Case 1: basic" begin
        g = CompleteDiGraph(4)
        mg = MetaDiGraph(g)

        m = Model(solver=CbcSolver())
        GraphJuMP.GraphModelExt{Int64, Float64}(mg)
        setgraph!(m, mg)
      end
    end
  end
end


text:
#=These are the initial XB Code States for the I5 code,
initial_xb_code_states[1] is a 1 3 chip array which represent the shift
register values initial_xb_code_states[3][4] represents the 4th shift register
of the GPS Signal with PRN numver 3  =#
const INITIAL_XB_CODE_STATES = [                #sat PRN number
    [0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0],    #01
    [1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1],    #02
    [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],    #03
    [1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0],    #04
    [1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1],    #05
    [0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0],    #06
    [1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1],    #07
    [1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0],    #08
    [1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1],    #09
    [0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0],    #10
    [0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0],    #11
    [1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1],    #12
    [0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0],    #13
    [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1],    #14
    [0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0],    #15
    [0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1],    #16
    [0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],    #17
    [1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0],    #18
    [1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1],    #19
    [0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1],    #20
    [0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],    #21
    [1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1],    #22
    [1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0],    #23
    [1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0],    #24
    [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1],    #25
    [1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0],    #26
    [0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0],    #27
    [0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0],    #28
    [0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1],    #29
    [1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1],    #30
    [0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0],    #31
    [0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1],    #32
    [1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1],    #33
    [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1],    #34
    [1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0],    #35
    [1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0],    #36
    [0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0]     #37
]

"""
$(SIGNATURES)

Takes the status of the registers as an `integer` and returns them as an array.
# Examples
```julia-repl
julia> reshape(8190)
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
```
"""
function reshape(integer)
    a = zeros(13)
    for i = 1:13
        b = (integer >> (i-1)) & 1
        a[14-i] = b
    end
    a
end

"""
$(SIGNATURES)

Takes the status of the registers as an Int `registers`, and an array of register `indices`
to calculate and return the new register values and the register output.

```julia-repl
julia> reshape(8190)
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
julia> output, registers = shift_register(8910,  [9, 10, 12, 13])
julia> reshape(registers)
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
julia> output == 1
true
```
"""
function shift_register(register, indices)
    update = 0
    for i in indices
        update = update ⊻ ((register >> (13 - i)) & 1)
    end
    (register & 1), (register >> 1) + 2^12 * update
end

"""
$(SIGNATURES)

Calculate the gps L5 PRN `satellite_code` for the initial XB register states
`initial_xb_code_states`.
```julia-repl
julia> initial_states_PRN_num_1_I = [0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0]
julia> prn_code_sat_1_I_signal = gen_l5_code(initial_states_PRN_num_1_I)
```
"""
function gen_l5_code(initial_xb_code_states)
    XA = 8191 # int with 3 leading zeros and then 13*1
    XB = initial_xb_code_states' * [4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1]
    satellite_code = zeros(Int8, 10230)
    XA_indices = [9, 10, 12, 13]
    XB_indices = [1, 3, 4, 6, 7, 8, 12, 13]
    for i = 1:10230
        output_xa, XA = shift_register(XA, XA_indices)
        output_xb, XB = shift_register(XB, XB_indices)
        satellite_code[i] =  2 * (output_xa ⊻ output_xb) - 1
        if (i == 8190)
            XA = 8191
        end
    end
    return satellite_code
end

"""
$(SIGNATURES)

Generate 10 periods of the PRN L5 code, with `initial_xb_code_states`, each ⊻ with one bit
of the 10bit Neuman-Hofman sequence 0000110101.
"""
function add_neuman_hofman_code(l5_code, neuman_hofman_code)
    vec(l5_code .* Vector{Int8}(neuman_hofman_code)')
end

function get_neuman_hofman_code()
    [1, 1, 1, 1, -1, -1, 1, -1, 1, -1]
end

const GPS_L5_CODES = mapreduce(
    sat -> add_neuman_hofman_code(
        gen_l5_code(INITIAL_XB_CODE_STATES[sat]), get_neuman_hofman_code()
    ),
    hcat,
    1:37
)::Array{Int8, 2}

"""
$(SIGNATURES)

Get codes of type GPSL5 as a Matrix where each column represents a PRN.
```julia-repl
julia> get_code(GPSL5)
```
"""
function get_codes(::Type{GPSL5})
    GPS_L5_CODES
end

"""
$(SIGNATURES)

Get code length of GNSS system GPSL5.
```julia-repl
julia> get_code_length(GPSL5)
```
"""
@inline function get_code_length(::Type{GPSL5})
    10230
end

"""
$(SIGNATURES)

Get secondary code length of GNSS system GPSL5.
```julia-repl
julia> get_secondary_code_length(GPSL5)
```
"""
@inline function get_secondary_code_length(::Type{GPSL5})
    10
end

"""
$(SIGNATURES)

Get center frequency of GNSS system GPSL5.
```julia-repl
julia> get_center_frequency(GPSL5)
```
"""
@inline function get_center_frequency(::Type{GPSL5})
    1_176_450_000Hz
end

"""
$(SIGNATURES)

Get code frequency of GNSS system GPSL5.
```julia-repl
julia> get_code_frequency(GPSL5)
```
"""
@inline function get_code_frequency(::Type{GPSL5})
    10_230_000Hz
end

"""
$(SIGNATURES)

Get data frequency of GNSS system GPSL5.
```julia-repl
julia> get_data_frequency(GPSL5)
```
"""
@inline function get_data_frequency(::Type{GPSL5})
    100Hz
end

"""
$(SIGNATURES)

Get code of GNSS system GPSL5 at phase `phase` of PRN `prn`.
The phase will not be wrapped by the code length. The phase has to smaller than the code
length incl. secondary code and must be an integer.
```julia-repl
julia> get_code_unsafe(GPSL5, 10, 1)
```
"""
Base.@propagate_inbounds function get_code_unsafe(::Type{GPSL5}, phase::Int, prn::Int)
    GPS_L5_CODES[1 + phase, prn]
end


text:
using Foo
using Test

# write your own tests here
@test 1 == 1
Foo.greet()


text:
using Documenter, Bio

custom_deps() = run(`pip install --user pygments mkdocs mkdocs-biojulia`)
makedocs()
deploydocs(
           repo = "github.com/BioJulia/Bio.jl.git",
           julia = "0.5",
           osname = "linux",
           deps = custom_deps
           )


text:
"""
This function returns all the subtypes at the leaf nodes (i.e. which do not have any further subtypes) of a given type.
"""
function leaftypes(t::Any, level=1, leaves = Vector{Any}())
    for s in InteractiveUtils.subtypes(t)
        if length(InteractiveUtils.subtypes(s)) == 0
            push!(leaves, s)
        else
            leaftypes(s, level+1, leaves)
        end
    end
    return leaves
end

"""
This function gets the capacity market forward years for the simulation.
"""
function get_capacity_forward_years(sim::Union{AgentSimulation, AgentSimulationData})
    forward_years = 1
    if get_markets(sim)[:Capacity]
        capacity_mkt_params = DataFrames.DataFrame(CSV.File(joinpath(get_data_dir(get_case(sim)), "markets_data", "capacity_mkt_param.csv")))
        forward_years = capacity_mkt_params.forward_years[1]
    end
    return forward_years
end

"""
This function returns the size of the project (in MW)
based on the project type and size type.
"""
function size_in_MW(investor_dir::String,
                    projecttype::String,
                    sizetype::String)
    filename = joinpath(investor_dir, "sizedict.csv")
    sizedata = read_data(filename)

    size = sizedata[findfirst(x-> x == projecttype, sizedata.type),
                findfirst(x -> x == sizetype, names(sizedata))]
    return size
end

"""
This function returns the size of the project (in MW)
if already specified in MW.
"""
function size_in_MW(investor_dir::String,
                    projecttype::String,
                    size::Float64)
    return size
end

"""
This function returns the size of the PSY Device.
"""
get_device_size(device::PSY.ThermalStandard) = PSY.get_active_power_limits(device)[:max]
get_device_size(device::PSY.RenewableDispatch) = PSY.get_rating(device)
get_device_size(device::PSY.HydroEnergyReservoir) = PSY.get_active_power_limits(device)[:max]
get_device_size(device::PSY.GenericBattery) = PSY.get_rating(device)

"""
This function returns line rating for PSY Lines.
"""
get_line_rating(line::PSY.Line) = PSY.get_rate(line)
get_line_rating(line::PSY.HVDCLine) = line.active_power_limits_from[:max]

"""
This function removes leap days from the data if year is not leap year.
"""
function remove_leap_day!(df:: DataFrames.DataFrame, start_year::Int64)
    if start_year % 4 != 0
        filter!(row -> row["Month"] != 2 || row["Day"] != 29, df)
    end

    return
end

"""
This function finds the start and end years
for reading price data.
"""
function find_price_years(construction_year::Int64,
                          end_life_year::Int64,
                          iteration_year::Int64,
                          yearly_horizon::Int64)

    price_years = (start_year = max((construction_year -
                                    iteration_year + 1),
                                    1),
                   end_year = min((end_life_year -
                                  iteration_year + 1),
                                  (yearly_horizon)))

    return price_years
end

"""
This function finds the start and end years
for updating project profit and npv data.
"""
function find_update_years(construction_year::Int64,
                          end_life_year::Int64,
                          iteration_year::Int64,
                          yearly_horizon::Int64)

    update_years = (start_year = max(construction_year,
                                     iteration_year),
                    end_year = min(end_life_year,
                                   yearly_horizon + iteration_year - 1))

    return update_years
end

"""
This function returns the number of hours modeled in the market clearing problem.
"""
function get_num_hours(system::MarketClearingProblem{Z, T}) where {Z, T}
    num_hours = T

    return num_hours
end

"""
This function returns an O&M cost array filled with the fixed
annual O&M cost values from update year to the end life year.
"""
function create_OMcost_array(finance::Finance,
                            scenario_name::String,
                            project_life_end::Int64,
                            update_start_year::Int64,
                            iteration_year::Int64)

    OM_cost_array = zeros(size(get_scenario_profit(finance)[scenario_name][iteration_year], 2))

    for i in update_start_year:project_life_end
        OM_cost_array[i] = get_fixed_OM_cost(finance)
    end

    return OM_cost_array

end

"""
This function return a queue cost array filled with the queue cost values
corresponding to each year spent in the queue.
"""
function create_queue_cost_array(size_profit::Int64,
         decision_year::Int64,
         queue_time::Int64,
         queue_cost::Vector{Float64})

    queue_cost_array = zeros(size_profit)

    for i in decision_year:(decision_year + queue_time - 1)
        queue_cost_array[i] = queue_cost[i - decision_year + 1]
    end

    return queue_cost_array

end
"""
This function does nothing if the project is not an Option project.
"""
function update_lifecycle!(project::P,
                           iteration_year::Int64,
                           simulation_years::Int64) where P <: Project{<: BuildPhase}
    return
end

"""
This function updates the decision, construction and endlife years of Option projects.
"""
function update_lifecycle!(project::P,
                           iteration_year::Int64,
                           simulation_years::Int64) where P <: Project{Option}

        finance_data = get_finance_data(project)

        queue_time = length(get_queue_cost(finance_data))

        if queue_time + iteration_year <  simulation_years
            set_decision_year!(project, (iteration_year + 1))
            set_construction_year!(project, (get_decision_year(project) +
                                            queue_time +
                                            get_lag_time(finance_data)))

            end_life_year = get_construction_year(project) + get_life_time(finance_data)  - 1
            set_end_life_year!(project, end_life_year)

            set_effective_investment_cost!(project, get_investment_cost(finance_data)[queue_time + iteration_year + 1])

            products = get_products(project)
            finance_data.realized_profit = AxisArrays.AxisArray(zeros(length(products), end_life_year),
                                                    AxisArrays.Axis{:prod}(get_name.(products)),
                                                    AxisArrays.Axis{:year}(1:1:end_life_year))

            finance_data.annual_cashflow = zeros(end_life_year)

            scenario_profit = get_scenario_profit(finance_data)

            for scenario_name in keys(scenario_profit)
                rows = size(scenario_profit[scenario_name][iteration_year], 1)
                cols = get_end_life_year(project)

                set_expected_npv!(finance_data, zeros(cols))
                set_expected_utility!(finance_data, zeros(cols))

                finance_data.scenario_profit[scenario_name] = [AxisArrays.AxisArray(zeros(length(products), cols),
                                                     AxisArrays.Axis{:prod}(get_name.(products)),
                                                     AxisArrays.Axis{:year}(collect(1:1:cols)))
                                                     for y in 1:cols]
                set_scenario_npv!(finance_data, scenario_name, zeros(cols))
                set_scenario_utility!(finance_data, scenario_name, zeros(cols))
            end
        end

    return
end

"""
This function updates the total installed capacity each year.
"""
function update_installed_cap!(installed_capacity::Vector{Float64},
                               announced_projects::Vector{<: Project{<: BuildPhase}},
                               iteration_year::Int64,
                               simulation_years::Int64)

    for year in iteration_year:simulation_years
        installed_capacity[year] = 0.0
        for project in announced_projects
            if (get_construction_year(project) <= year) && (get_end_life_year(project) >= year)
                installed_capacity[year] += get_maxcap(project)
            end
        end
    end
    return installed_capacity
end

"""
This function does nothing if project is not in Existing phase.
"""
function extrapolate_profits!(project::Project{<:BuildPhase}, final_year::Int64)
    return
end

"""
This function extrapolates the profits of Existing projects beyond the simulation horizon.
"""
function extrapolate_profits!(project::Project{Existing}, final_year::Int64)
    end_life_year = get_end_life_year(project)
    finance_data = get_finance_data(project)
    if end_life_year > final_year
        for y in (final_year + 1):end_life_year
            set_annual_cashflow!(finance_data, y, get_annual_cashflow(finance_data)[final_year])
            for product in get_products(project)
                name = get_name(product)
                set_realized_profit!(finance_data, name, y, get_realized_profit(finance_data)[name, final_year])
            end
        end
    end
    return
end


text:
using DataFrames
using CSV
using Statistics
using StatsBase: sample
using Revise
using EvoTrees
using BenchmarkTools

# prepare a dataset
features = rand(100_000, 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X,1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# train model
params1 = EvoTreeRegressor(
    loss=:linear, metric=:mae,
    nrounds=10,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=0.5, nbins=32)

# with Int64 set 100K rows: 654.218 ms (583199 allocations: 488.97 MiB)
# with BitSet 𝑖 100K rows: 563.297 ms (576642 allocations: 475.41 MiB)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n = 2)
@btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train = predict(model, X_train)
@btime pred_train = predict($model, $X_train)
mean(abs.(pred_train .- Y_train))

# train model
params1 = EvoTreeRegressor(
    loss=:logistic, metric=:logloss,
    nrounds=10,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=1.0, colsample=1.0, nbins=64)
@time model = grow_gbtree(X_train, Y_train, params1, X_eval = X_eval, Y_eval = Y_eval, print_every_n=2)
@time pred_train = predict(model, X_train)


text:
using .NavAbilityCaesarExtensions
using NavAbilitySDK

function addDataEntry!(dfg::NavAbilityDFG, label::Symbol, entry::BlobStoreEntry)
  throw("Not yet implemented.")
end

function addDataBlob!(blobstore::NavAbilityBlobStore, entry::BlobStoreEntry, blob::Vector{UInt8})
  throw("Not yet implemented.")
end

function getDataEntry(dfg::NavAbilityDFG, label::Symbol, key::Symbol)
  throw("Not yet implemented.")
end

function getDataBlob(blobstore::NavAbilityBlobStore, entry::BlobStoreEntry) where T
  throw("Not yet implemented.")
end


text:
function f1()
	s = Vector{String}(100000)

	for i = 1:length(s)
		s[i] = string(i)
	end
end

function f2()
	s = Vector{String}(10000)

	for i = 1:length(s)
		s[i] = string(i)
	end
end

function test()
	for i = 1:500
		f1()
		f2()
	end
end

test()
@time test()

Profile.init(n=10^7, delay = 0.00001)
@profile test()
Profile.print(format=:flat)

# Wyniki:
#
# delay   | f1 count | f2 count | f1 %
# -------------------------------------
# 0       | N/A      | N/A      | 90.00
# -------------------------------------
# 0.5     | 6        | 1        | 85.71
# 0.3     | 6        | 1        | 85.71
# 0.1     | 32       | 5        | 86.48
# 0.1     | 33       | 2        | 94.28
# 0.05    | 64       | 7        | 90.14
# 0.05    | 65       | 6        | 91.54
# 0.01    | 377      | 30       | 92.62
# 0.01    | 329      | 32       | 91.13
# 0.005   | 671      | 59       | 91.91
# 0.005   | 671      | 61       | 91.66
# 0.001   | 3554     | 318      | 91.78
# 0.001   | 3613     | 290      | 92.56
# 0.0001  | 61467    | 5249     | 92.13
# 0.0001  | 54516    | 4809     | 91.89
# 0.00001 | 208978   | 18508    | 91.86


text:
######################################################################
# isogenies.jl: Isogenies between elliptic curves
######################################################################

#Only odd-degree isogenies are implemented; create a special class for degree 2 isogenies

export Isogeny

export degree, rational_fractions, areisomorphic, compose, squarefree_kernel, isomorphism_coefficients, subgroup, dual, isomorphism, find_isomorphism, multiplication_isogeny

######################################################################
# Type definitions
######################################################################

"""
Concrete type for isogenies between Weierstrass or Montgomery elliptic curves.
Only odd degree isogenies are currently implemented.
The kernel polynomial of an l-isogeny is a monic, squarefree polynomial of degree (l-1)/2.
In the rational fractions, the usual factor y in the ordinate is dropped in order to use only univariate polynomials.
"""
#TODO: make this a Nemo.Map
mutable struct Isogeny{T}
	domain::EllipticCurve{T}
	degree::Integer
	kernel::PolyElem{T}
	u::T
	r::T
	s::T
	t::T
	image::EllipticCurve{T}
	x::Union{Frac{U}, Nothing} where U <: PolyElem{T}
	y::Union{Frac{U}, Nothing} where U <: PolyElem{T}
end

base_ring(phi::Isogeny) = base_ring(domain(phi))

domain(phi::Isogeny) = phi.domain

image(phi::Isogeny) = phi.image

kernel(phi::Isogeny) = phi.kernel

degree(phi::Isogeny) = phi.degree

"""
Get the u, r, s, t coefficients of the isomophism coming after the normalized part of the isogeny
"""
isomorphism_coefficients(phi::Isogeny) = (phi.u, phi.r, phi.s, phi.t)

"""
Get the rational fractions defining the isogeny. They are computed at the first call.
"""
function rational_fractions(phi::Isogeny)
	try
		phix = get(phi.rationalx)
		phiy = get(phi.rationaly)
		return (phix, phiy)
	catch
		fx, fy = _compute_rational_fractions(phi)
		phi.x = fx
		phi.y = fy
		return fx, fy
	end
end

"""
Decide if two isogenies differ by postcomposition by an isomorphism.
"""
areisomorphic(phi::Isogeny, psi::Isogeny) = domain(phi) == domain(psi) && degree(phi) == degree(psi) && kernel(phi) == kernel(psi)

"""
Decide if two isogenies are the same.
"""
==(phi::Isogeny, psi::Isogeny) =  areisomorphic(phi, psi) && isomorphism_coefficients(phi) == isomorphism_coefficients(psi)

"""
Shows a description of an isogeny between elliptic curves.
"""
function show(io::IO, phi::Isogeny)
	print("Isogeny of degree ")
	show(phi.degree)
	print(" between ")
	show(io, phi.domain)
	print(" and ")
	show(io, phi.image)
	print(" with kernel ")
	show(io, phi.kernel)
end

"""
Decide if an isogeny has valid formulas.
"""
function isvalid(phi::Isogeny)
	E1 = domain(phi)
	E2 = image(phi)
	a11, a21, a31, a41, a61 = a_invariants(E1) #we have a1 = a3 = 0
	a12, a22, a32, a42, a62 = a_invariants(E2) #we have a1 = a3 = 0
	phix, phiy = rational_fractions(phi)
	X = gen(base_ring(parent(phix)))
	lhs = phiy^2 * (X^3 + a21 * X^2 + a41 * X + a61)
	rhs = (phix^3 + a22 * phix^2 + a42 * phix + a62)
	test = lhs - rhs
	return test==0
end

######################################################################
# Rational fractions and evaluation
######################################################################

"""
Get the evaluation of an isogeny on a point.
"""
function (phi::Isogeny)(P::EllipticPoint)
	Q = normalized(P)
	x1 = Q.X
	y1 = Q.Y
	phix, phiy = rational_fractions(phi)
	try
		x2 = phix(x1)
		y2 = y1 * phiy(x1)
		return EllipticPoint(x2, y2, one(base_ring(Q)), image(phi))
	catch
		return infinity(image(phi))
	end
end

#DOES NOT WORK if u or t != 0...
"""
Compute the rational functions defining an isogeny
"""
function _compute_rational_fractions(phi::Isogeny)
	phixn, phiyn = _normalized_rational_fractions(phi)
	u, r, s, t = isomorphism_coefficients(phi)
	return (u^2 * phixn + r, u^3 * phiyn + s * u^2 * phixn + t)
end

function _normalized_rational_fractions(phi::Isogeny)
	l = degree(phi)
	E = domain(phi)
	K = base_ring(E)
	b2, b4, b6, b8 = b_invariants(E)
	a1, a2, a3, a4, a6 = a_invariants(E)

	if l==1 #phi is an isomorphism
		A, x = PolynomialRing(K, "x")
		Ffield = FractionField(A)
		return (Ffield(x), Ffield(1))

	else #remember l is odd

		n = div(l - 1, 2)
		kpoly = kernel(phi)
		s1 = - coeff(kpoly, n-1)
		kpolysquare = kpoly^2
		dkpoly = derivative(kpoly)
		d2kpoly = derivative(dkpoly)

		FX = parent(kpoly)
		psi2 = FX(divisionpolynomial(E, -1))
		X = gen(FX)
		N = psi2 * (dkpoly^2 - d2kpoly * kpoly) - (6 * X^2 + b2 * X + b4) * dkpoly * kpoly + (l * X - 2 * s1) * kpolysquare

		if K(2) != 0
			dN = derivative(N)
			omega = dN * kpoly - 2 * N * dkpoly    #+ (1 // K(2)) * ((a1 * N + a3 * kpolysquare) * kpoly)(x): we force to have a1 = a3 = 0
			return (N // kpolysquare , omega // (kpoly * kpolysquare))

		else #K has characteristic 2
			throw(NotImplementedError("Rational fractions not yet implemented in characteristic 2"))
		end
	end
end

######################################################################
# Composition
######################################################################

function compose(phi1::Isogeny, phi2::Isogeny)
	E1 = domain(phi1)
	E2 = image(phi1)
	E3 = image(phi2)
	@assert E2 == domain(phi2)

	phi1x, phi1y = rational_fractions(phi1)
	phi2x, phi2y = rational_fractions(phi2)

	resx = phi2x(phi1x)
	resy = phi1y * phi2y(phi1x)

	return Isogeny(E1, E3, resx, resy)
end

######################################################################
# Defining isogenies with various inputs
######################################################################

"""
Build an isogeny given its domain, kernel, degree, and image
"""
function Isogeny(E::AbstractWeierstrass{T}, poly::PolyElem{T}, Eprime::AbstractWeierstrass{T}) where T
	phi1 = Isogeny(E, poly)
	postisomorphism = find_isomorphism(image(phi1), Eprime)
	phi = compose(phi1, postisomorphism)
	return phi
end

function Isogeny(E1::AbstractWeierstrass{T}, E2::AbstractWeierstrass{T}, phix::Frac{U}, phiy::Frac{U}) where U<:PolyElem{T} where T
	K = sqrt(denominator(phix))
	F = base_ring(K)
	normalizedphix, normalizedphiy = rational_fractions(Isogeny(E1, K))
	usquare = phix//normalizedphix
	ucube = phiy//normalizedphiy

	#convert to base field
	num = numerator(usquare)
	den = denominator(usquare)
	@assert den == 1
	@assert degree(num) == 0
	usquare = coeff(num, 0)
	num = numerator(ucube)
	den = denominator(ucube)
	@assert den == 1
	@assert degree(num) == 0
	ucube = coeff(num, 0)

	u = ucube // usquare
	testphiy = u^3 * normalizedphiy
	if testphiy == - phiy
		u = -u
	end
	@assert degree(numerator(phix)) == 2 * degree(K) + 1
	res = Isogeny(E1, 2 * degree(K) + 1, K, u, F(0), F(0), F(0), E2, phix, phiy)
	@assert isvalid(res)
	return res
end

"""
Build an isogeny given its domain and its kernel polynomial.
"""
function Isogeny(E::AbstractWeierstrass{T}, poly::PolyElem{T}) where T

	CurveType = curvetype(E)

	K = base_ring(E)
	a1, a2, a3, a4, a6 = a_invariants(E)
	b2, b4, b6, b8 = b_invariants(E)
	n = Nemo.degree(poly)

	coeff(poly, n) == 1 || throw(ArgumentError("kernel polynomial must be monic"))

	s1, s2, s3 = zero(K), zero(K), zero(K)

	n > 0 && (s1 = - coeff(poly, n - 1))
	n > 1 && (s2 = coeff(poly, n - 2))
	n > 2 && (s3 = - coeff(poly, n - 3))
	t = 6 * (s1^2 - 2 * s2) + b2 * s1 + n * b4
	w = 10 * (s1^3 - 3 * s1 * s2 + 3 * s3) + 2 * b2 * (s1^2 - 2 * s2) + 3 * b4 * s1 + n * b6
	E1 = CurveType(a1, a2, a3, a4 - 5 * t, a6 - b2 * t - 7 * w)

	return Isogeny(E, 2 * n + 1, poly, K(1), K(0), K(0), K(0), E1, nothing, nothing)
end

######################################################################
# BMSS algorithm
######################################################################

include("bmss.jl")

######################################################################
# Isogenies between Weierstrass curves
######################################################################

"""
Get the polynomial associated to an l-torsion rational point, l being an odd prime.

The input is not checked.
"""
function subgroup(Q::EllipticPoint{T}, l::Nemo.Integer) where T<:Nemo.FieldElem
	normalize!(Q)
	K = base_ring(Q)
	R, x = PolynomialRing(K, "x")
	poly = Nemo.one(R)
	point = Q
	for k = 1 : div(l-1, 2)
		poly *= (x - point.X)
		point += Q
	end
	return poly
end

"""
Build an isogeny given an elliptic curve in short Weierstrass form, the j-invariant of the targetted elliptic curve, and the degree, assuming it exists.
Returns a normalized isogeny.
If the keyword argument USE_ATKIN is set to false (default), we use classical modular equations. If USE_ATKIN is set to true, we use Atkin modular equations.
"""

function Isogeny(E::ShortWeierstrass{T}, degree::Nemo.Integer, jprime::T; USE_ATKIN=false) where T
	K = base_ring(E)
	j = j_invariant(E)
	KX, X = PolynomialRing(K, "X")
	KXY, Y = PolynomialRing(KX, "Y")
	if USE_ATKIN
		print("using Atkin polys")
		A_ell_j = ClassPolynomials.AtkinModularPolynomial(degree, X, j)
		A_ell_jprime = ClassPolynomials.AtkinModularPolynomial(degree, X, jprime)
		poly = gcd(A_ell_j, A_ell_jprime)
		test, f = any_root(poly)
		if !test
			throw(ArgumentError("The curves are not linked by a rational isogeny of this degree"))
		end
		return Isogeny(E, degree, jprime, f)
	else
		Phi_l = ClassPolynomials.ClassicalModularPolynomial(degree, X, Y)
		#polynomial in Y over polynomials in X
		l = K(degree)

		derx = derivative(Phi_l)(j)(jprime)
		dery = derivative(Phi_l)(jprime)(j) #works because Phi_l is symmetric

		J = - K(18) // l * (E.b // E.a) * (derx // dery) * j
		jj = jprime * (jprime - 1728)

		aprime = - l^4 * J^2 // (48 * jj)
		bprime = - l^6 * J^3 // (864 * jprime * jj)
		Eprime = ShortWeierstrass(aprime, bprime)
		poly = kernelpoly(E, Eprime, degree)
		return Isogeny(E, degree, poly, K(1), K(0), K(0), K(0), Eprime, nothing, nothing)
	end
end

"""
Another element in the base field may be given as argument. In this case we use the Atkin modular A_\ell equations instead of the classical ones,
and f must satisfy A_\ell(j, f) == A_\ell(jprime, f) == 0. The curve must be given in Short Weierstrass form.
"""
function Isogeny(E::ShortWeierstrass{T}, degree::Nemo.Integer, jprime::T, f::T) where T
	j = j_invariant(E)
	_, _, _, a, b = a_invariants(E)
	K = base_ring(E)
	@assert ClassPolynomials.AtkinModularPolynomial(degree, f, j) == 0
	@assert ClassPolynomials.AtkinModularPolynomial(degree, f, jprime) == 0 #sanity checks
	KX, X = PolynomialRing(K, "X")
	dF = derivative(ClassPolynomials.AtkinModularPolynomial(degree, X, j))(f)
	dF *= f
	dF2 = derivative(ClassPolynomials.AtkinModularPolynomial(degree, X, jprime))(f)
	dF2 *= f
	poly = derivative(ClassPolynomials.AtkinModularPolynomial(degree, f, X))
	dJ = j * poly(j)
	dJ2 = degree * jprime * poly(jprime)
	jj = jprime // (jprime - 1728)
	lambda = (dF2 * dJ * b) // (dJ2 * dF * a)
	aprime = - K(27)//K(4) * degree^4 * lambda^2 * jj
    bprime = - K(27)//K(4) * degree^6 * lambda^3 * jj
	Eprime = ShortWeierstrass(aprime, bprime)
	poly = kernelpoly(E, Eprime, degree)

	return Isogeny(E, degree, poly, K(1), K(0), K(0), K(0), Eprime, nothing, nothing)
end

"""
Build an isogeny given an odd integer l and a rational torsion point of this order. The input is not checked.
"""

function Isogeny(E::AbstractWeierstrass{T}, Q::EllipticPoint{T}, l::Nemo.Integer) where T<:Nemo.FieldElem
	poly = subgroup(Q, l)
	return Isogeny(E, poly)
end

"""
Build an isogeny of given degree between two curves, assuming a normalized separable isogeny of this degree exists.
"""

function Isogeny(E1::AbstractWeierstrass{T}, degree::Nemo.Integer, E2::AbstractWeierstrass{T}) where T<:Nemo.FieldElem
	poly = kernelpoly(E1, E2, degree)
	return Isogeny(E1, poly)
end

"""
Compute the scalar multiplication by m on a weierstrass curve as a normalized isogeny.
"""
function multiplication_isogeny(E::AbstractWeierstrass, m::Nemo.Integer)
	@assert isodd(m)
	K = base_ring(E)
	psim = divisionpolynomial(E, m)
	n = degree(psim)
	phi1 = Isogeny(E, 1 // coeff(psim, n) * psim)
	E2 = image(phi1)
	postisomorphism = isomorphism(E2, 1 // K(m))
	@assert image(postisomorphism) == E
	return compose(phi1, postisomorphism)
end

######################################################################
# Dual isogeny
######################################################################

function dual(phi::Isogeny)
	K = base_ring(phi)
	E1 = domain(phi)
	E2 = image(phi)
	d = degree(phi)
	E1_0 = image(isomorphism(E1, K(d)))
	phi_hat = Isogeny(E2, d, E1_0) #use BMSS
	#de-normalize !
	postisomorphism = isomorphism(image(phi_hat), 1//K(d))
	phi_hat = compose(phi_hat, postisomorphism)
	@assert isvalid(phi_hat)
	@assert image(phi_hat) == E1
	return phi_hat
end

function isomorphism(E::AbstractWeierstrass{T}, u::T) where T
	a1, a2, a3, a4, a6 = a_invariants(E) #we have a1 = a3 = 0
	a2prime = a2 * u^2
	a4prime = a4 * u^4
	a6prime = a6 * u^6
	K = base_ring(E)
	R, X = PolynomialRing(base_ring(E), "X")
	phi = Isogeny(E, 1, R(1), u, K(0), K(0), K(0), EllipticCurve(a2prime, a4prime, a6prime), nothing, nothing)
	return phi
end

function find_isomorphism(E1::AbstractWeierstrass{T}, E2::AbstractWeierstrass{T}) where T
	a11, a21, a31, a41, a61 = a_invariants(E1)
	a12, a22, a32, a42, a62 = a_invariants(E2)
	u = sqrt(a22 // a21)
	@assert a42 = a41 * u^4
	@assert a62 = a61 * u^6
	return isomophism(E1, u)
end

######################################################################
# Operations
######################################################################

function +(phi::Isogeny, psi::Isogeny)
	E = domain(phi)
	@assert E == domain(psi)
	Eprime = image(phi)
	@assert Eprime == image(psi)
	phix, phiy = rational_fractions(phi)
	psix, psiy = rational_fractions(phi)
	throw(NotImplementedError("Use the addition laws on the image curve"))
	#return Isogeny(E, Eprime, phix + psix, phiy + psiy)
end

function -(phi::Isogeny)
	pix, phiy = rational_fractions(phi)
	return Isogeny(domain(phi), image(phi), phix, -phiy)
end

function -(phi::Isogeny, psi::Isogeny)
	return phi + (-psi)
end


text:
using BEAST
using CompScienceMeshes

fn = joinpath(@__DIR__,"assets","thick_cladding.msh")
G01 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma01")
G02 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma02")
G12 = CompScienceMeshes.read_gmsh_mesh(fn, physical="Gamma12")

G = CompScienceMeshes.read_gmsh_mesh(fn)

@assert numcells(G01) + numcells(G02) + numcells(G12) == numcells(G)

# Build region boundaries with normals pointing inward:
G1 = weld(-G01, -G12)
@assert CompScienceMeshes.isoriented(G1)

G2 = weld(-G02, G12)
@assert CompScienceMeshes.isoriented(G2)

E1 = CompScienceMeshes.interior(G1)
E2 = CompScienceMeshes.interior(G2)
E = CompScienceMeshes.interior(G)

Nd = BEAST.nedelec(G, E)

Nd1 = BEAST.nedelec(G1, E1)
Nd2 = BEAST.nedelec(G2, E2)

RT1 = n × Nd1
RT2 = n × Nd2

A1 = CompScienceMeshes.embedding(E1,E)
A2 = CompScienceMeshes.embedding(E2,E)


text:
using LinearAlgebra
using Statistics
using StatsBase

# N-dimensional box
struct BoxN
    N::Integer # dimension of the box
    mins::Vector{Float64} # vector of lower bounds on each dimension
    maxs::Vector{Float64} # vector of upper bounds
    points::Vector{Vector{Float64}} # points inside the box
end

"""
    subdivide_box(b::BoxN)

Subdivide an N-dimensional box `b` into 2^N subboxes with half the side lengths of `b`.

The points in `b` are assigned to their respective subboxes. Any point on the border of multiple boxes is put into the box whose mean point is closest to it. Returns an array of non empty boxes.
"""
function subdivide_box(b::BoxN)
    halves = (b.mins+b.maxs)/2
    boxes = Dict{BitVector,BoxN}()
    border_points = []
    for p in b.points
        if any(p .== halves) # does p lie on the border of multiple boxes?
            push!(border_points, p)
        else
            key = p .> halves
            if !haskey(boxes, key)
                mins = (.!key) .* b.mins + key .* halves
                maxs = (.!key) .* halves + key .* b.maxs
                boxes[key] = BoxN(b.N, mins, maxs, [p])
            else
                push!(boxes[key].points, p)
            end
        end
    end
    if isempty(boxes)
        return []
    end
    centers = [(k, mean(v.points)) for (k,v) in boxes]
    # assign each border point to the box whose center is closest to it
    for p in border_points
        i = argmin(Dict([k=> norm(v-p) for (k,v) in centers]))
        push!(boxes[i].points, p)
    end
    return collect(values(boxes))
end

"""
    box_counting_dimension(points, maxiter=10, tol=0.001)

Approximate the box counting dimension of `points` defined by
`` d = \\lim_{ε\\longrightarrow0}\\frac{\\ln N_ε(E)}{-\\ln ε}. ``

 This is done by taking the smallest box containing all points and subdividing it up to `maxiter` times or until d seems to have converged (with tolerance `tol`).
 Returns vectors `d`, `N` and `ε`, containing respectively the iterates of d, the number of boxes N, and the box length ε.
"""
function box_counting_dimension(points::Vector{Vector{T}}, maxiter=10, tol=0.001) where T <: Real
    if isempty(points)
        error("points must be non empty")
    end
    n = length(points[1])
    mins = Vector(undef, n)
    maxs = Vector(undef, n)
    for i in 1:n
        tmp = [points[j][i] for j in 1:length(points)]
        mins[i] = minimum(tmp)
        maxs[i] = maximum(tmp)
    end
    boxes = [BoxN(n, mins, maxs, points)]
    ε = geomean(abs.(maxs-mins)) # use geomtric mean for side length
    N = [1]
    d = Vector{Float64}()
    push!(d, log(N[1])/(0*log(2)-log(ε)))
    i = 2
    while i <= maxiter && !isempty(boxes)
        tmp = map(subdivide_box, boxes)
        boxes = isempty(tmp) ? [] : foldl(append!, tmp)
        push!(N, length(boxes))
        push!(d, log(N[i])/((i-1)*log(2)-log(ε)))
        if abs(d[i]-d[i-1]) < tol
            break # exit if d seems to have converged. This is not entirely foolproof but it will avoid unnecessary iterations when `points` has too few elements. Set tol=0.0 to disable
        end
        i += 1
    end
    return d, N, ε*2.0 .^ -(0:(length(N)-1))
end

function box_counting_dimension(points::Vector{T}, maxiter=10, tol=0.001) where T<:Real
    return box_counting_dimension(map(x->[x], points), maxiter, tol)
end

"""
    box_counting_dimension(A, maxiter=10, tol=0.001)

Approximate the box counting dimension for a 2D array of bools.

This is meant for use with images and 2D arrays by selecting which pixels you want with a predicate.
"""
function box_counting_dimension(A::BitArray{2}, maxiter=10, tol=0.001)
    return box_counting_dimension(array_to_points(A), maxiter, tol)
end

"""
    array_to_points(A)

Convert a 2D Array of bools to a list of 2D points corresponding to indices where `A` is `true`.
"""
function array_to_points(A::BitArray{2})
    (m,n) = size(A)
    points::Vector{Vector{Float64}} = []
    for j in 1:n, i in 1:m
        if A[i,j]
            push!(points, [(i-0.5)/m, (j-0.5)/n])
        end
    end
    return points
end

"""
    boundary(A)

Compute the boundary a 2D region.

The region is represented by a 2D BitArray where entries equal to `true` denote points in the region.
"""
function boundary(A::BitArray{2})
    (m,n) = size(A)
    B = falses(m,n)
    for j in 1:n, i in 1:m
        sum = 0
        if i > 1 sum += A[i-1,j]; end
        if i < m sum += A[i+1,j]; end
        if j > 1 sum += A[i,j-1]; end
        if j < n sum += A[i,j+1]; end
        if sum < 4 && sum > 1
            B[i,j] = true
        end
    end
    return B
end

"""
    correlation_dimension(x, maxiter=10)

Approximate the correlation dimension of `x` using `maxiter` iterations.

Returns vectors `C` and `r` containing the iterates of the correlation function ``C(r)`` and the radii ``r``.
"""
function correlation_dimension(x::Vector{Vector{T}}, maxiter::Int=10) where T <: Real
    if isempty(x)
        error("x must be non empty")
    end
    r = maximum(norm.(x)) # quick and dirty max radius
    N = length(x)
    Cnums = zeros(Int, maxiter)
    C = zeros(maxiter)
    for i in 1:N-1, j in i+1:N
        d = norm(x[i]-x[j])
        for k in 1:maxiter
            if d > r*2.0^(-k)
                break
            else
                Cnums[k] += 1
            end
        end
    end
    C = float(Cnums)/((N*(N-1))/2)
    return C, r * 2.0 .^ -(1:maxiter)
end

"""
    loglog_regression(y,x)

Compute the OLS regression of `log.(y)` by `log.(x)`.

This only returns the coefficients, so you may want to use the `GLM` package for more information.
"""
function loglog_regression(y::Vector,x::Vector)
    X = ones(length(x),2)
    X[:,1] = log.(x)
    return X\log.(y)
end

max_stars_repo_name: csimal/chaos-and-fractals

max_issues_repo_name: csimal/chaos-and-fractals

max_forks_repo_name: csimal/chaos-and-fractals


text:
using Memoize
"""
    Statevector

Represents a pure state.
"""
struct Statevector <: AbstractArray{ComplexF64,1}
    state::Vector{ComplexF64}
    vec::Vector{ComplexF64}
    perm::Vector{Int}
    N::Int

    @doc """
        Statevector(v::Vector{ComplexF64})

    Create a `Statevector` object representing a pure quantum state.
    """
    function Statevector(v::Vector{ComplexF64})
        ispow2(length(v)) || error("Statevector must have length which is a power of 2")
        N = intlog2(length(v))
        new(v, zeros(ComplexF64, 2^N), collect(1:2^N), N)
    end
    @doc """
    Statevector(v::Vector{ComplexF64})

    Create a `Statevector` object representing a pure quantum state with 'N' qubits.
    """
    function Statevector(N::Int)
        new(zeros(ComplexF64, 2^N), zeros(ComplexF64, 2^N), collect(1:2^N), N)
    end
end

Base.setindex!(ψ::Statevector, v, i::Int64) = setindex!(ψ.state, v, i)
Base.getindex(ψ::Statevector, i::Int) = getindex(ψ.state, i)
Base.size(ψ::Statevector) = size(ψ.state)
Base.length(ψ::Statevector) = length(ψ.state)
Base.IndexStyle(ψ::Type{Statevector}) = IndexLinear()

"""
    rdm(N, iwire, ψ, χ, d=2)

Compute the reduced density matrix ``tr_B[|ψ⟩⟨χ|]``, where the trace runs over
the subsystem complementary to the qubits specified by `iwire`.
"""
function rdm(N::Integer, iwire::NTuple{M,<:Integer}, ψ::Statevector, χ::Statevector, d::Int=2) where {M}
    M ≥ 1 || error("Need at least one wire to act on.")
    M ≤ N || error("Number of gate wires cannot be larger than total number of wires.")
    length(unique(iwire)) == M || error("Wire indices must be unique.")
    minimum(iwire) ≥ 1 || error("Wire index cannot be smaller than 1.")
    maximum(iwire) ≤ N || error("Wire index cannot be larger than total number of wires.")
    length(ψ) == d^N || error("Input vector 'ψ' has wrong length.")
    length(χ) == d^N || error("Input vector 'χ' has wrong length.")

    # convert to array
    iwire = collect(iwire)
    # complementary wires
    iwcompl = setdiff(1:N, iwire)
    @assert length(iwire) + length(iwcompl) == N

    ρ = zeros(eltype(ψ), d^M, d^M)

    # Note: following the ordering convention of `kron` here, i.e.,
    # last qubit corresponds to fastest varying index
    strides = [d^(j-1) for j in 1:N]
    wstrides = strides[iwire]
    cstrides = strides[iwcompl]

    if d == 2
        iw = BitArray{1}(undef, M)
        jw = BitArray{1}(undef, M)
        kw = BitArray{1}(undef, N - M)
        # TODO: optimize memory access pattern
        for k in 0:2^(N-M)-1
            binary_digits!(kw, k)
            koffset = dot(kw, cstrides)
            for i in 1:2^M
                binary_digits!(iw, i - 1)
                for j in 1:2^M
                    binary_digits!(jw, j - 1)
                    rowind = koffset + dot(iw, wstrides) + 1
                    colind = koffset + dot(jw, wstrides) + 1
                    ρ[i, j] += ψ[rowind] * conj(χ[colind])
                end
            end
        end
    else
        error("d = $d not supported yet.")
    end

    return ρ
end


text:
__precompile__()
module QBO

    using HttpCommon
    using Requests
    using JSON
    using AbstractTrees
    using FixedPointDecimals
    using Missings
    using Base.Meta

    abstract type QBObject; end

    export Line, Purchase

    # Missing => May have a value, but we don't know what it is
    # Void => Definitely no value there
    const Maybe{T} = Union{T, Missing}
    const Optional{T} = Union{T, Missing, Void}

    const Decimal = FixedDecimal{Int64, 2}

    include("utils.jl")
    include("refs.jl")
    include("mangling.jl")
    include("types.jl")
    include("auth.jl")
    include("request.jl")
    include("api.jl")
    include("query.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/OpenBLAS_jll.jl
module OpenBLAS_jll

using Libdl, CompilerSupportLibraries_jll, Base.BinaryPlatforms

const PATH_list = String[]
const LIBPATH_list = String[]

export libopenblas

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libopenblas_handle = C_NULL
libopenblas_path = ""

if arch(HostPlatform()) in ("x86_64", "ppc64le")
    const libsuffix = "64_"
else
    const libsuffix = ""
end

if Sys.iswindows()
    const libopenblas = "libopenblas$(libsuffix).dll"
elseif Sys.isapple()
    const libopenblas = "@rpath/libopenblas$(libsuffix).dylib"
else
    const libopenblas = "libopenblas$(libsuffix).so"
end

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    global libopenblas_handle = dlopen(libopenblas)
    global libopenblas_path = dlpath(libopenblas_handle)
end

is_available() = true

end  # module OpenBLAS_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
"""
```
LinkedCells
```

`firstatom` is a vector that contains for each `index i3D = index3D(i,j,k)` of a box, which
is the first atom of the list of atoms in that box.

`nextatom` is a vector that contains the index of the next atom of that box given the index
of the previous atom (starting with `firstatom`).

"""
struct LinkedCells
    firstatom::Vector{Int}
    nextatom::Vector{Int}
end

LinkedCells(n) = LinkedCells(
    zeros(Int, n), # firstatom (actual required size is nc1*nc2*nc*3
    zeros(Int, n), # nextatom
)


text:
using Documenter, DynamicGrids, CUDAKernels

CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing

makedocs(
    modules = [DynamicGrids],
    sitename = "DynamicGrids.jl",
    checkdocs = :all,
    strict = true,
    format = Documenter.HTML(
        prettyurls = CI,
    ),
    pages = [
        "DynamicGrids" => "index.md",
        # "Examples" => "examples.md",
    ],
)

deploydocs(
    repo = "github.com/cesaraustralia/DynamicGrids.jl.git",
)


text:
module TestXGBoost

using MLJBase
using Test
import MLJModels
import XGBoost
using CategoricalArrays
using MLJModels.XGBoost_
using Random: seed!

@test_logs (:warn, r"Changing ") XGBoostClassifier(objective="wrong")
@test_logs (:warn, r"Changing ") XGBoostCount(objective="wrong")
@test_logs (:warn, r"Your") XGBoostRegressor(objective="binary:logistic")

## REGRESSOR

seed!(0)
plain_regressor = XGBoostRegressor()
n,m = 10^3, 5 ;
features = rand(n,m);
weights = rand(-1:1,m);
labels = features * weights;
features = MLJBase.table(features)
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 0, features, labels);
rpred = predict(plain_regressor, fitresultR, features);
info_dict(XGBoostRegressor)

plain_regressor.objective = "gamma"
labels = abs.(labels)
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 0, features, labels);
rpred = predict(plain_regressor, fitresultR, features);

importances = reportR.feature_importances

## COUNT

count_regressor = XGBoostCount(num_round=10)
using Random: seed!
using Distributions

seed!(0)

X = randn(100, 3) .* randn(3)'
Xtable = table(X)

α = 0.1
β = [-0.3, 0.2, -0.1]
λ = exp.(α .+ X * β)
y = [rand(Poisson(λᵢ)) for λᵢ ∈ λ]

fitresultC, cacheC, reportC = MLJBase.fit(count_regressor, 0, Xtable, y);
cpred = predict(count_regressor, fitresultC, Xtable);
info_dict(XGBoostCount)

importances = reportC.feature_importances

## CLASSIFIER

plain_classifier = XGBoostClassifier(num_round=100, seed=0)

# test binary case:
N=2
seed!(0)
X = (x1=rand(1000), x2=rand(1000), x3=rand(1000))
ycat = map(X.x1) do x
    string(mod(round(Int, 10*x), N))
end |> categorical
y = identity.(ycat) # make plain Vector with categ. elements
train, test = partition(eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(plain_classifier, 0,
                                            selectrows(X, train), y[train];)
yhat = mode.(predict(plain_classifier, fitresult, selectrows(X, test)))
misclassification_rate = sum(yhat .!= y[test])/length(test)
@test misclassification_rate < 0.01

importances = report.feature_importances

# Multiclass{10} case:
N=10
seed!(0)
X = (x1=rand(1000), x2=rand(1000), x3=rand(1000))
ycat = map(X.x1) do x
    string(mod(round(Int, 10*x), N))
end |> categorical
y = identity.(ycat) # make plain Vector with categ. elements

train, test = partition(eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(plain_classifier, 0,
                                            selectrows(X, train), y[train];)
yhat = mode.(predict(plain_classifier, fitresult, selectrows(X, test)))
misclassification_rate = sum(yhat .!= y[test])/length(test)
@test misclassification_rate < 0.01

# check target pool preserved:
X = (x1=rand(400), x2=rand(400), x3=rand(400))
ycat = vcat(fill(:x, 100), fill(:y, 100), fill(:z, 200)) |>categorical
y = identity.(ycat)
train, test = partition(eachindex(y), 0.5)
@test length(unique(y[train])) == 2
@test length(unique(y[test])) == 1
fitresult, cache, report = MLJBase.fit(plain_classifier, 0,
                                            selectrows(X, train), y[train];)
yhat = predict_mode(plain_classifier, fitresult, selectrows(X, test))
@test Set(MLJBase.classes(yhat[1])) == Set(MLJBase.classes(y[train][1]))

info_dict(XGBoostClassifier)

end
true


text:
struct ProductNode{T,C} <: AbstractProductNode
    data::T
    metadata::C

    function ProductNode{T,C}(data::T, metadata::C) where {T, C}
        @assert length(data) >= 1 && all(x -> nobs(x) == nobs(data[1]), data)
        new(data, metadata)
    end
end

ProductNode(data::T) where {T} = ProductNode{T, Nothing}(data, nothing)
ProductNode(data::T, metadata::C) where {T, C} = ProductNode{T, C}(data, metadata)

@deprecate TreeNode(data) ProductNode(data)
@deprecate TreeNode(data, metadata) ProductNode(data, metadata)

Flux.@functor ProductNode

mapdata(f, x::ProductNode) = ProductNode(map(i -> mapdata(f, i), x.data), x.metadata)

Base.getindex(x::ProductNode, i::Symbol) = x.data[i]
Base.keys(x::ProductNode) = keys(x.data)

Base.ndims(x::AbstractProductNode) = Colon()
LearnBase.nobs(a::AbstractProductNode) = nobs(a.data[1], ObsDim.Last)
LearnBase.nobs(a::AbstractProductNode, ::Type{ObsDim.Last}) = nobs(a)

function reduce(::typeof(catobs), as::Vector{T}) where {T <: ProductNode}
    data = _cattrees([x.data for x in as])
    metadata = reduce(catobs, [a.metadata for a in as])
    ProductNode(data, metadata)
end

Base.getindex(x::ProductNode, i::VecOrRange) = ProductNode(subset(x.data, i), subset(x.metadata, i))

Base.hash(e::ProductNode{T,C}, h::UInt) where {T,C} = hash((T, C, e.data, e.metadata), h)
Base.:(==)(e1::ProductNode{T,C}, e2::ProductNode{T,C}) where {T,C} = e1.data == e2.data && e1.metadata == e2.metadata


text:
using MultimodeCavity
const mm = MultimodeCavity

N = 100
mass = 0.5

k = 1.
delta = -1
eta = 1
kappa = 2
U0 = 1

particles = mm.classical.Particles(N, mass)
mode1 = mm.classical.CavityMode(k, delta, eta, kappa, U0)
mode2 = mm.classical.CavityMode(k, delta, eta, kappa, U0)
system = mm.classical.MultimodeSystem(particles, [mode1, mode2])


text:
# Parameters:
#	N::Int
#	- the problem size
#	K::Int
#	- the problem is homogenous at N-K
# Return Values
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1} (Young's Orthogonal Representations)
#	- YOR[n][p][k] is Young's Orthogonal Representation for the Adjacent Transposition (K, K + 1) for the pth Partition of n that is needed for the bandlimited functionality
#	- length(YOR[n]) = 0 for n = 1:(N - K - 1)
#	- if p < ZFI[n], length(YOR[n][p] = 1) and YOR[n][p][1,1] contains the dimension of the full Young's Orthogonal Representation
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output1 of partition_tree_bl()
#	ZFI::Array{Int, 1}
#	- output2 for paritions_bl()
function yor_bl(N::Int, K::Int)
	P, ZFI, WI = partitions_bl(N, K)
	PT = partition_tree_bl(N, K, P, ZFI, WI)
	YS = ys_symbols_bl(N, K, P, ZFI, PT)
	YOR = Array(Array{Array{SparseMatrixCSC, 1}, 1}, N)
	for n = 1:(N - K - 1)
		YOR[n] = Array(Array{SparseMatrixCSC, 1}, 0)
	end
	for n = (N - K):N
		YSn = YS[n]
		Pn = P[n]
		Pn_L = length(Pn)
		YORn = Array(Array{SparseMatrixCSC, 1}, Pn_L)
		for p = 1:ZFI[n]
			YORnp_sparse = spzeros(1, 1)
			YORnp_sparse[1,1] = length(YSn[p])
			YORnp = Array(SparseMatrixCSC, 1)
			YORnp[1] = YORnp_sparse
			YORn[p] = YORnp
		end
		for p = (ZFI[n] + 1):Pn_L
			YSnp = YSn[p]
			Pnp = Pn[p]
			R = length(Pnp)
			YORn[p] = yor_p(n, R, YSnp)
		end
		YOR[n] = YORn
	end
	return YOR, PT, ZFI
end


text:
using UncNLPrograms
using Test, CUTEst, LinearAlgebra, Printf
import NLPModels as jo

# Tests translation of CUTE problems against CUTEst
for nlp ∈ values(UncNLPrograms.TestSet)
    @printf "\nTesting %s to CUTEst.jl Model at default iterate:\n" nlp.name
    qt_nlp = CUTEstModel(nlp.name)
    qt_x0  = qt_nlp.meta.x0
    x0 = nlp.x0

    # set test points in a random proximity of the default itterate, both +/=
    t1 = x0 - rand(nlp.n)
    t2 = x0 - rand(nlp.n)
    t3 = x0 - rand(nlp.n)

    t4 = x0 + rand(nlp.n)
    t5 = x0 + rand(nlp.n)
    t6 = x0 + rand(nlp.n)
    try
        # test similarity of default iterates
        @test obj(nlp, x0) ≈ jo.obj(qt_nlp, qt_x0)
        @test grad(nlp, x0) ≈ jo.grad(qt_nlp, qt_x0)
        @test hessAD(nlp, x0) ≈ jo.hess(qt_nlp, qt_x0)

        # test at 6 randomly generated points
        @test obj(nlp, t1) ≈ jo.obj(qt_nlp, t1)
        @test grad(nlp, t1) ≈ jo.grad(qt_nlp, t1)
        @test hessAD(nlp, t1) ≈ jo.hess(qt_nlp, t1)

        @test obj(nlp, t2) ≈ jo.obj(qt_nlp, t2)
        @test grad(nlp, t2) ≈ jo.grad(qt_nlp, t2)
        @test hessAD(nlp, t2) ≈ jo.hess(qt_nlp, t2)
        @test obj(nlp, t3) ≈ jo.obj(qt_nlp, t3)
        @test grad(nlp, t3) ≈ jo.grad(qt_nlp, t3)
        @test hessAD(nlp, t3) ≈ jo.hess(qt_nlp, t3)

        @test obj(nlp, t4) ≈ jo.obj(qt_nlp, t4)
        @test grad(nlp, t4) ≈ jo.grad(qt_nlp, t4)
        @test hessAD(nlp, t4) ≈ jo.hess(qt_nlp, t4)

        @test obj(nlp, t5) ≈ jo.obj(qt_nlp, t5)
        @test grad(nlp, t5) ≈ jo.grad(qt_nlp, t5)
        @test hessAD(nlp, t5) ≈ jo.hess(qt_nlp, t5)

        @test obj(nlp, t6) ≈ jo.obj(qt_nlp, t6)
        @test grad(nlp, t6) ≈ jo.grad(qt_nlp, t6)
        @test hessAD(nlp, t6) ≈ jo.hess(qt_nlp, t6)
    catch e
        println(e)
    end

    finalize(qt_nlp)
end


text:
using MathOptInterface
using ParametricOptInterface
using BenchmarkTools
const MOI = MathOptInterface
const POI = ParametricOptInterface
import Pkg

function moi_add_variables(N::Int)
    model = MOI.Utilities.Model{Float64}()
    MOI.add_variables(model, N)
    return nothing
end

function poi_add_variables(N::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    MOI.add_variables(model, N)
    return nothing
end

function poi_add_parameters(N::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    MOI.add_constrained_variable.(model, POI.Parameter.(ones(N)));
    return nothing
end

function poi_add_parameters_and_variables(N::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    MOI.add_variables(model, N/2)
    MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2))))
    return nothing
end

function poi_add_parameters_and_variables_alternating(N::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    for i in 1:Int(N/2)
        MOI.add_variable(model)
        MOI.add_constrained_variable(model, POI.Parameter(1))
    end
    return nothing
end

function moi_add_saf_ctr(N::Int, M::Int)
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.(1.0, x),
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_saf_ctr(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.(1.0, x),
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_ctr(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm.(1.0, x); MOI.ScalarAffineTerm.(1.0, y)],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_ctr_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm.(1.0, x); MOI.ScalarAffineTerm.(1.0, y)],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    MOI.set.(model, POI.ParameterValue(), y, 0.5)
    POI.update_parameters!(model)
    return nothing
end

function moi_add_sqf_variables_ctr(N::Int, M::Int)
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, x),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_sqf_variables_ctr(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, x),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_sqf_variables_parameters_ctr(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_sqf_variables_parameters_ctr_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    MOI.set.(model, POI.ParameterValue(), y, 0.5)
    POI.update_parameters!(model)
    return nothing
end

function poi_add_sqf_parameters_parameters_ctr(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, y, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    return nothing
end

function poi_add_sqf_parameters_parameters_ctr_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.add_constraint(
                model,
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, y, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                ),
                MOI.GreaterThan(1.0),
            )
    end
    MOI.set.(model, POI.ParameterValue(), y, 0.5)
    POI.update_parameters!(model)
    return nothing
end

function moi_add_saf_obj(N::Int, M::Int)
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.set(
            model,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.(1.0, x),
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_saf_obj(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                MOI.ScalarAffineFunction(
                    MOI.ScalarAffineTerm.(1.0, x),
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_obj(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm.(1.0, x); MOI.ScalarAffineTerm.(1.0, y)],
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_obj_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                MOI.ScalarAffineFunction(
                    [MOI.ScalarAffineTerm.(1.0, x); MOI.ScalarAffineTerm.(1.0, y)],
                    0.0,
                )
            )
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), y, 0.5)
        POI.update_parameters!(model)
    end
    return nothing
end

function moi_add_sqf_variables_obj(N::Int, M::Int)
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, x),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_sqf_variables_obj(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N)
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, x),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_sqf_variables_parameters_obj(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_sqf_variables_parameters_obj_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, x, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), y, 0.5)
        POI.update_parameters!(model)
    end
    return nothing
end

function poi_add_sqf_parameters_parameters_obj(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, y, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    return nothing
end

function poi_add_sqf_parameters_parameters_obj_parameter_update(N::Int, M::Int)
    model = POI.Optimizer(MOI.Utilities.Model{Float64}())
    x = MOI.add_variables(model, N/2)
    y = first.(MOI.add_constrained_variable.(model, POI.Parameter.(ones(Int(N/2)))))
    for _ in 1:M
        MOI.set(
                model,
                MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
                MOI.ScalarQuadraticFunction(
                    MOI.ScalarQuadraticTerm.(1.0, y, y),
                    MOI.ScalarAffineTerm{Float64}[],
                    0.0,
                )
            )
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), y, 0.5)
        POI.update_parameters!(model)
    end
    return nothing
end

function run_benchmarks(N::Int, M::Int)
    println("Pkg status:")
    Pkg.status()
    println("")
    GC.gc()
    println("variables on a MOIU.Model.")
    @btime moi_add_variables($N)
    GC.gc()
    println("variables on a POI.Optimizer.")
    @btime poi_add_variables($N)
    GC.gc()
    println("parameters on a POI.Optimizer.")
    @btime poi_add_parameters($N)
    GC.gc()
    println("parameters and variables on a POI.Optimizer.")
    @btime poi_add_parameters_and_variables($N)
    GC.gc()
    println("alternating parameters and variables on a POI.Optimizer.")
    @btime poi_add_parameters_and_variables_alternating($N)
    GC.gc()
    println("SAF constraint with variables on a MOIU.Model.")
    @btime moi_add_saf_ctr($N, $M)
    GC.gc()
    println("SAF constraint with variables on a POI.Optimizer.")
    @btime poi_add_saf_ctr($N, $M)
    GC.gc()
    println("SAF constraint with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_ctr($N, $M)
    GC.gc()
    println("SQF constraint with variables on a MOIU.Model{Float64}.")
    @btime moi_add_sqf_variables_ctr($N, $M)
    GC.gc()
    println("SQF constraint with variables on a POI.Optimizer.")
    @btime poi_add_sqf_variables_ctr($N, $M)
    GC.gc()
    println("SQF constraint with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_ctr($N, $M)
    GC.gc()
    println("SQF constraint with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_ctr($N, $M)
    GC.gc()
    println("SAF objective with variables on a MOIU.Model.")
    @btime moi_add_saf_obj($N, $M)
    GC.gc()
    println("SAF objective with variables on a POI.Optimizer.")
    @btime poi_add_saf_obj($N, $M)
    GC.gc()
    println("SAF objective with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_obj($N, $M)
    GC.gc()
    println("SQF objective with variables on a MOIU.Model.")
    @btime moi_add_sqf_variables_obj($N, $M)
    GC.gc()
    println("SQF objective with variables on a POI.Optimizer.")
    @btime poi_add_sqf_variables_obj($N, $M)
    GC.gc()
    println("SQF objective with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_obj($N, $M)
    GC.gc()
    println("SQF objective with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_obj($N, $M)
    GC.gc()
    println("Update parameters in SAF constraint with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SAF objective with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_obj_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF constraint with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF constraint with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF objective with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_obj_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF objective with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_obj_parameter_update($N, $M)
    return nothing
end

N = 10_000
M = 100
run_benchmarks(N, M)


text:
include("Dataset.jl")

import .Dataset

root = mkpath("build")
Dataset.clone_gem5(root)
gem5b = Dataset.build_gem5_x86(root)
Dataset.exec_gem5opt_command(
    gem5b, Dataset.base_cacheconfig, joinpath(root,"gem5out"),
    joinpath(root,"gem5","tests","test-progs","hello","bin","x86","linux","hello"),
    "", 0)

Dataset.clone_mibench(root)
mibenchb = Dataset.build_mibench_x86(root)

max_stars_repo_name: emcores/CacheConfigurations

max_issues_repo_name: emcores/CacheConfigurations

max_forks_repo_name: emcores/CacheConfigurations


text:
# Autogenerated wrapper script for Geant4_jll for aarch64-apple-darwin
export libG4Analysis, libG4CLHEP, libG4Digits, libG4ErrorPropagation, libG4Event, libG4FR, libG4Geometry, libG4Global, libG4Graphics, libG4Intercoms, libG4Interfaces, libG4Material, libG4Mocren, libG4Modeling, libG4ParModels, libG4Particles, libG4Persistency, libG4PhysicsLists, libG4Processes, libG4RayTracer, libG4Readout, libG4Run, libG4Track, libG4Tracking, libG4Tree, libG4VRML, libG4VisHepRep, libG4VisXXX, libG4Zlib, libG4visManagement

using Expat_jll
using Xorg_libXmu_jll
using Libglvnd_jll
JLLWrappers.@generate_wrapper_header("Geant4")
JLLWrappers.@declare_library_product(libG4Analysis, "@rpath/libG4analysis.dylib")
JLLWrappers.@declare_library_product(libG4CLHEP, "@rpath/libG4clhep.dylib")
JLLWrappers.@declare_library_product(libG4Digits, "@rpath/libG4digits_hits.dylib")
JLLWrappers.@declare_library_product(libG4ErrorPropagation, "@rpath/libG4error_propagation.dylib")
JLLWrappers.@declare_library_product(libG4Event, "@rpath/libG4event.dylib")
JLLWrappers.@declare_library_product(libG4FR, "@rpath/libG4FR.dylib")
JLLWrappers.@declare_library_product(libG4Geometry, "@rpath/libG4geometry.dylib")
JLLWrappers.@declare_library_product(libG4Global, "@rpath/libG4global.dylib")
JLLWrappers.@declare_library_product(libG4Graphics, "@rpath/libG4graphics_reps.dylib")
JLLWrappers.@declare_library_product(libG4Intercoms, "@rpath/libG4intercoms.dylib")
JLLWrappers.@declare_library_product(libG4Interfaces, "@rpath/libG4interfaces.dylib")
JLLWrappers.@declare_library_product(libG4Material, "@rpath/libG4materials.dylib")
JLLWrappers.@declare_library_product(libG4Mocren, "@rpath/libG4GMocren.dylib")
JLLWrappers.@declare_library_product(libG4Modeling, "@rpath/libG4modeling.dylib")
JLLWrappers.@declare_library_product(libG4ParModels, "@rpath/libG4parmodels.dylib")
JLLWrappers.@declare_library_product(libG4Particles, "@rpath/libG4particles.dylib")
JLLWrappers.@declare_library_product(libG4Persistency, "@rpath/libG4persistency.dylib")
JLLWrappers.@declare_library_product(libG4PhysicsLists, "@rpath/libG4physicslists.dylib")
JLLWrappers.@declare_library_product(libG4Processes, "@rpath/libG4processes.dylib")
JLLWrappers.@declare_library_product(libG4RayTracer, "@rpath/libG4RayTracer.dylib")
JLLWrappers.@declare_library_product(libG4Readout, "@rpath/libG4readout.dylib")
JLLWrappers.@declare_library_product(libG4Run, "@rpath/libG4run.dylib")
JLLWrappers.@declare_library_product(libG4Track, "@rpath/libG4track.dylib")
JLLWrappers.@declare_library_product(libG4Tracking, "@rpath/libG4tracking.dylib")
JLLWrappers.@declare_library_product(libG4Tree, "@rpath/libG4Tree.dylib")
JLLWrappers.@declare_library_product(libG4VRML, "@rpath/libG4VRML.dylib")
JLLWrappers.@declare_library_product(libG4VisHepRep, "@rpath/libG4visHepRep.dylib")
JLLWrappers.@declare_library_product(libG4VisXXX, "@rpath/libG4visXXX.dylib")
JLLWrappers.@declare_library_product(libG4Zlib, "@rpath/libG4zlib.dylib")
JLLWrappers.@declare_library_product(libG4visManagement, "@rpath/libG4vis_management.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Expat_jll, Xorg_libXmu_jll, Libglvnd_jll)
    JLLWrappers.@init_library_product(
        libG4Analysis,
        "lib/libG4analysis.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4CLHEP,
        "lib/libG4clhep.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Digits,
        "lib/libG4digits_hits.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4ErrorPropagation,
        "lib/libG4error_propagation.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Event,
        "lib/libG4event.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4FR,
        "lib/libG4FR.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Geometry,
        "lib/libG4geometry.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Global,
        "lib/libG4global.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Graphics,
        "lib/libG4graphics_reps.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Intercoms,
        "lib/libG4intercoms.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Interfaces,
        "lib/libG4interfaces.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Material,
        "lib/libG4materials.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Mocren,
        "lib/libG4GMocren.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Modeling,
        "lib/libG4modeling.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4ParModels,
        "lib/libG4parmodels.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Particles,
        "lib/libG4particles.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Persistency,
        "lib/libG4persistency.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4PhysicsLists,
        "lib/libG4physicslists.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Processes,
        "lib/libG4processes.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4RayTracer,
        "lib/libG4RayTracer.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Readout,
        "lib/libG4readout.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Run,
        "lib/libG4run.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Track,
        "lib/libG4track.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Tracking,
        "lib/libG4tracking.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Tree,
        "lib/libG4Tree.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4VRML,
        "lib/libG4VRML.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4VisHepRep,
        "lib/libG4visHepRep.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4VisXXX,
        "lib/libG4visXXX.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4Zlib,
        "lib/libG4zlib.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libG4visManagement,
        "lib/libG4vis_management.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Plots
function gr_hist_2d(xx::Array{Float64,2} ; nbins=20)
      gr()
      Plots.histogram2d( xx[1,:] , xx[2,:] , nbins=nbins )
end

function plot_heatmap_2d_function( f::Any , xlim , ylim , xres , yres )
    xygrids = xygrid(xlim,ylim,xres,yres)
    zgrid   = NaN*ones(xres,yres)
    for zi=1:length(zgrid); zgrid[zi] = f_ec([xygrids[1][zi];xygrids[2][zi]] ); end
    using Plots
    gr()
    Plots.heatmap(zgrid)
end


text:
###############################################################################
#
# Initialize the RT module
#
###############################################################################
"""
    initialize_rt_module(FT; nLayer::Int = 20, LAI::Number = FT(3))

Initialize the RT module and return the sturctures, given
- `FT` Floating number type
- `nLayer` Number of canopy layers
- `LAI` Leaf area index

This function initializes and returns
- `angles` [`SolarAngles`](@ref)
- `can` [`Canopy4RT`](@ref)
- `can_opt` [`CanopyOpticals`](@ref)
- `can_rad` [`CanopyRads`](@ref)
- `in_rad` [`IncomingRadiation`](@ref)
- `leaves` Array{[`LeafBios`](@ref),1}
- `rt_con` [`RTCache`](@ref)
- `rt_dim` [`RTDimensions`](@ref)
- `soil` [`SoilOpticals`](@ref)
- `wls` [`WaveLengths`](@ref)
"""
function initialize_rt_module(FT; nLayer::Int = 20, LAI::Number = FT(3))
    # 1. create wls and can, which are rt_dims independent
    can = create_canopy_rt(FT, nLayer=nLayer, LAI=LAI);
    wls = create_wave_length(FT);

    # 2. create rt_dims from wls and can
    rt_dim = create_rt_dims(can, wls);

    # 3. create can_rad, can_opt, and etc from rt_dim and wls
    can_rad = create_canopy_rads(FT, rt_dim);
    can_opt = create_canopy_opticals(FT, rt_dim);
    in_rad  = create_incoming_radiation(wls);
    soil    = SoilOpticals{FT}(wls);
    angles  = SolarAngles{FT}();
    rt_con  = create_rt_cache(FT, rt_dim);

    # Create an array of standard leaves
    leaves = [create_leaf_bios(FT, rt_dim) for i in 1:nLayer];
    for i in 1:nLayer
        fluspect!(leaves[i], wls);
    end

    # Four Different steps to compute Short-Wave RT
    canopy_geometry!(can, angles, can_opt, rt_con);
    canopy_matrices!(leaves, can_opt);
    short_wave!(can, can_opt, can_rad, in_rad, soil, rt_con);
    canopy_fluxes!(can, can_opt, can_rad, in_rad, soil, leaves, wls, rt_con);
    SIF_fluxes!(leaves, can_opt, can_rad, can, soil, wls, rt_con, rt_dim);

    # # Compute Long Wave (Last term is LW incoming in W m^-2)
    thermal_fluxes!(leaves, can_opt, can_rad, can, soil, [FT(400.0)], wls);

    return angles, can, can_opt, can_rad, in_rad, leaves, rt_con, rt_dim, soil, wls
end

max_stars_repo_name: climate-machine/Land