text:
using DiffEqBase

using Reexport
@reexport using DiffEqBase

include("twoscales_pure_ab.jl")

struct HOODEFunction{iip,N} <: DiffEqBase.AbstractODEFunction{iip}
    f::Any
end
(fct::HOODEFunction{false,3})(u, p, t) = fct.f(u, p, t)
(fct::HOODEFunction{false,2})(u, p, t) = fct.f(u, p)
(fct::HOODEFunction{false,1})(u, p, t) = fct.f(u)
function (fct::HOODEFunction{true,4})(u, p, t)
    du = zero(u)
    fct.f(du, u, p, t)
    return du
end

"""
    HOODEProblem(f, u0, tspan, p, A, epsilon, B)

The HOODE problem is :

```math
\\frac{du}{dt} = ( \\frac{1}{\\epsilon} A + B)  u + f(u,p,t)
```

- The initial condition is ``u(tspan[1]) = u0``.
- The solution ``u(t)`` will be computed for ``tspan[1] ≤ t ≤ tspan[2]``
- Constant parameters to be supplied as the second argument of ``f``.
- Periodic Matrix of the problem.
- epsilon of the problem.
- Matrix of linear problem to get the exact solution

"""
struct HOODEProblem{T} <: DiffEqBase.DEProblem
    f::HOODEFunction
    u0::Vector{T}
    tspan::Tuple{T,T}
    p::Any
    A::AbstractMatrix
    epsilon::T
    B::Union{Matrix,Missing}
    function HOODEProblem(
        f,
        u0::Vector{T},
        tspan::Tuple{T,T},
        p,
        A::AbstractMatrix,
        epsilon::T,
        B::Union{Matrix,Missing},
    ) where {T}
        fct = if hasmethod(f, (Vector{T}, Vector{T}, Any, T))
            HOODEFunction{true,4}(f)
        elseif hasmethod(f, (Vector{T}, Any, T))
            HOODEFunction{false,3}(f)
        elseif hasmethod(f, (Vector{T}, Any))
            HOODEFunction{false,2}(f)
        elseif hasmethod(f, (Vector{T},))
            HOODEFunction{false,1}(f)
        else
            println("err !!!!!")
        end
        return new{T}(fct, u0, tspan, p, A, epsilon, B)
    end # end of function
end # end of struct
function HOODEProblem(f, u0, tspan, p, A, epsilon)
    return HOODEProblem(f, u0, tspan, p, A, epsilon, missing)
end
struct HOODEInterpolation{T} <: DiffEqBase.AbstractDiffEqInterpolation
    t::Vector{T}
    u_caret::Vector{Array{Complex{T},2}}
    parphi::PreparePhi
    order::Any
end
(interp::HOODEInterpolation)(t, idxs, deriv, p, continuity) = interp(t)
function (interp::HOODEInterpolation)(t)
    return _getresult(
        interp.t,
        interp.u_caret,
        t,
        interp.parphi,
        interp.t[1],
        interp.t[end],
        interp.order,
    )
    #    return _getresult(interp.u_caret, t,
    #    interp.parphi,
    #    interp.t[1], interp.t[end],
    #    interp.order)
end
(interp::HOODEInterpolation)(vt::Vector{Float64}) = interp.(vt)
(interp::HOODEInterpolation)(vt::Vector{BigFloat}) = interp.(vt)
if typeof(DiffEqBase.AbstractTimeseriesSolution{Float64,Float64,Float64}) == DataType
    abstract type AbstractHOODESolution{T,N} <: DiffEqBase.AbstractTimeseriesSolution{T,N,N} end
else
    abstract type AbstractHOODESolution{T,N} <: DiffEqBase.AbstractTimeseriesSolution{T,N} end
end
struct HOODESolution{T} <: AbstractHOODESolution{T,T}
    u::Vector{Vector{T}}
    u_tr::Union{Vector{Vector{T}},Nothing}
    t::Vector{T}
    dense::Bool
    order::Integer
    par_u0::PrepareU0
    p_coef::CoefExpAB
    prob::HOODEProblem{T}
    retcode::Any
    interp::Union{HOODEInterpolation,Nothing}
    tslocation::Any
    absprec::Any
    relprec::Any
end
function HOODESolution(retcode::Symbol)
    return HOODESolution(
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        retcode,
        nothing,
        nothing,
        nothing,
    )
end
(sol::HOODESolution)(t) = sol.dense ? sol.interp(t) : undef

abstract type AbstractHOODEAlgorithm end

struct HOODETwoScalesAB <: AbstractHOODEAlgorithm end
struct HOODEETDRK2 <: AbstractHOODEAlgorithm end
struct HOODEETDRK3 <: AbstractHOODEAlgorithm end
struct HOODEETDRK4 <: AbstractHOODEAlgorithm end

# function DiffEqBase.build_solution{T}(prob::HOODEProblem{T},
#     sol::Vector{Vector{T}},
#     t::Vector{T},
#     fftsol::Vector{Array{T,2}}) where T<:AbstractFloat
#     return HOODESolution(sol, t, fftsol)
# end
# function DiffEqBase.build_solution{T}(prob::HOODEProblem{T},
#     sol::Vector{Vector{T}},
#     t::Vector{T}) where T<:AbstractFloat
#     return HOODESolution(sol, t, undef)
# end

"""
    function DiffEqBase.solve(prob::HOODEProblem{T};
    nb_tau::Integer=32,
    order::Integer=4,
    order_prep::Integer=order+2,
    dense::Bool=true,
    nb_t::Integer=100,
    getprecision::Bool=dense,
    verbose=100,
    par_u0::Union{PrepareU0,Missing}=missing,
    p_coef::Union{CoefExpAB,Missing}=missing
    ) where T<:AbstractFloat

specific interface solver for Highly oscillatory problems, that an ODE of this form:
```math
\\frac{\\delta u(t)}{\\delta t} = \\frac{1}{\\varepsilon} A + F(u(t), t)
```
where ``u \\in \\R^n`` and  ``0 < \\varepsilon < 1``
``A`` must be a **periodic matrix** i.e. ``e^{t A} = e^{(t+\\pi) A}`` for any ``t \\in \\R``

## Argument :
- `prob::HOODEProblem{T}` : The problem to solve

## Keywords :
- `nb_tau::Integer=32` : number of values of FFT transform, must be power of twoscales_pure_ab
- `order::Integer=4` : order of Adams-Bashforth method, and also of the interpolatation
- `order_prep::Integer=order+2` : order of the preparation
- `dense::Bool=true` : if true it is possible to compute solution at any time of the interval
- `nb_t::Integer=100` : number of period slices
- `getprecision::Bool=dense` : compute the absolute and relative precision
- `par_u0::Union{PrepareU0,Missing}=missing` : preparation data for u0
- `p_coef::Union{CoefExpAB,Missing}=missing` : coefficients for Adams-Bashforth method

## Examples :
"""
function DiffEqBase.solve(prob::HOODEProblem{T}; kwargs...) where {T<:AbstractFloat}
    return DiffEqBase.solve(prob, HOODETwoScalesAB(); kwargs...)
end
function DiffEqBase.solve(
    prob::HOODEProblem{T},
    alg::HOODETwoScalesAB;
    nb_tau::Integer = 32,
    order::Integer = 4,
    order_prep::Integer = order + 2,
    dense::Bool = true,
    nb_t::Integer = 100,
    getprecision::Bool = dense,
    verbose = 100,
    par_u0::Union{PrepareU0,Missing} = missing,
    p_coef::Union{CoefExpAB,Missing} = missing,
) where {T<:AbstractFloat}
    retcode = :Success
    nb_tau = prevpow(2, nb_tau)
    traceln(
        100,
        "solve function prob=$prob,\n nb_tau=$nb_tau, order=$order, order_prep=$order_prep, dense=$dense,\n nb_t=$nb_t, getprecision=$getprecision, verbose=$verbose";
        verbose = verbose,
    )
    if getprecision
        s1 = DiffEqBase.solve(
            prob;
            nb_tau = nb_tau,
            order = order,
            order_prep = order_prep,
            dense = dense,
            nb_t = nb_t,
            getprecision = false,
            verbose = verbose - 10,
            par_u0 = par_u0,
        )
        n_nb_t = Integer(floor(1.1373 * nb_t + 1))
        s2 = DiffEqBase.solve(
            prob;
            nb_tau = nb_tau,
            order = order,
            order_prep = order_prep,
            dense = dense,
            nb_t = n_nb_t,
            getprecision = false,
            verbose = verbose - 10,
            par_u0 = s1.par_u0,
        )
        absprec = norm(s1.u[end] - s2.u[end])
        relprec = absprec / max(norm(s1.u[end]), norm(s2.u[end]))
        # if dense
        #     for i = 1:10
        #         t = rand(T) * (prob.tspan[2]-prob.tspan[1])
        #             + prob.tspan[1]
        #         a, b = s1(t), s2(t)
        #         ap = norm(a-b)
        #         rp = ap/max(norm(a),norm(b))
        #         absprec = max(absprec,ap)
        #         relprec = max(relprec,rp)
        #     end
        # end
        return HOODESolution(
            s1.u,
            s1.u_tr,
            s1.t,
            dense,
            order,
            s1.par_u0,
            s1.p_coef,
            prob,
            retcode,
            s1.interp,
            0,
            Float64(absprec),
            Float64(relprec),
        )
    end

    par_u0 = if ismissing(par_u0)
        parphi = PreparePhi(
            prob.epsilon,
            nb_tau,
            prob.A,
            prob.f,
            prob.B;
            t_0 = prob.tspan[1],
            paramfct = prob.p,
        )
        PrepareU0(parphi, order_prep, prob.u0)
    else
        par_u0
    end
    pargen = PrepareTwoScalesPureAB(
        nb_t,
        prob.tspan[2],
        order,
        par_u0,
        p_coef = p_coef,
        verbose = verbose,
    )
    if dense
        u_mat, _, u_caret = twoscales_pure_ab(pargen, res_fft = true, verbose = verbose)
        t = collect(0:nb_t) * (prob.tspan[2] - prob.tspan[1]) / nb_t .+ prob.tspan[1]
        interp = HOODEInterpolation{T}(t, u_caret, par_u0.parphi, order)
        u = reshape(mapslices(x -> [x], u_mat, dims = 1), size(u_mat, 2))
        u_tr = reshape(mapslices(x -> [x], transpose(u_mat), dims = 1), size(u_mat, 1))
    else
        u = Array{Array{T,1},1}(undef, 2)
        u[1] = copy(prob.u0)
        u[end] = twoscales_pure_ab(pargen, verbose = verbose, only_end = true)
        t = [prob.tspan[1], prob.tspan[2]]
        u_tr = nothing
        interp = nothing
    end
    return HOODESolution(
        u,
        u_tr,
        t,
        dense,
        order,
        par_u0,
        pargen.p_coef,
        prob,
        retcode,
        interp,
        0,
        nothing,
        nothing,
    )
end


text:
"""
    cl6(x::Float64)::Float64

Returns the value of the Clausen function ``\\operatorname{Cl}_6(x)``
for a real angle ``x`` of type `Float64`.  This function is defined as

```math
\\operatorname{Cl}_6(x) = \\Im[\\operatorname{Li}_6(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\sin(kx)}{k^6}
```

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> cl6(1.0)
0.855629273183937
```
"""
function cl6(x::Float64)::Float64
    zeta3 = 1.2020569031595943
    pi28 = pi*pi/8.0

    (x, sgn) = range_reduce_even(x)

    if x == zero(x) || x == pi
        return zero(x)
    end

    if x < 0.5*pi
        P = (1.0369277551433699e+00, -2.087195444107175e-01,
             2.0652251045312954e-02, -1.383438138256840e-04)
        Q = (1.0000000000000000e+00, -8.0784096827362542e-03,
             5.8074568862993102e-06, -5.1960620033050114e-10)
        y = x*x
        y2 = y*y
        p = P[1] + y * P[2] + y2 * (P[3] + y * P[4])
        q = Q[1] + y * Q[2] + y2 * (Q[3] + y * Q[4])
        sgn*x*(p/q - 1/120*y2*log(x))
    else
        P = (7.9544504578027050e-01, -1.9255025309738589e-01,
             1.5805208288846591e-02, -5.4175380521534706e-04,
             6.7577493541009068e-06)
        Q = (1.0000000000000000e+00, -7.0798422394109274e-02,
             7.1744189715634762e-04,  3.9098747334347093e-06,
             3.5669441618295266e-08,  2.5315391843409925e-10)
        y = pi - x
        z = y*y - pi28
        z2 = z*z
        z4 = z2*z2
        p = P[1] + z * P[2] + z2 * (P[3] + z * P[4]) + z4 * P[5]
        q = Q[1] + z * Q[2] + z2 * (Q[3] + z * Q[4]) + z4 * (Q[5] + z * Q[6])
        sgn*y*p/q
    end
end

ext: jl


text:
using Documenter, GmshTools

makedocs(
    doctest=false,
    modules = [GmshTools],
    sitename = "GmshTools",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Index" => "index.md",
        "Element Types" => "element_types.md",
    ],
)

deploydocs(
    repo = "github.com/shipengcheng1230/GmshTools.jl.git",
    target = "build",
)


text:
module MCPolicyPlotting
using Plots

using Curiosity

actions = [1,2,3]

function load_policy(policy_name)
    Curiosity.LearnedPolicy("/home/matthewmcleod/Documents/Masters/curiosity/src/data/MC_learned_policies/$(policy_name).bson")
end

function get_q(LP,state_space)
    q = zeros(length(actions),length(state_space),length(state_space))
    for a in actions
        for (i,x) in enumerate(state_space)
            for (j,y) in enumerate(state_space)
                q[a,i,j] = LP([x,y])[a]
            end
        end
    end
    return q
end

function plot_hms(qs,state_space)
    hms = []
    for a in actions
        hm = heatmap(state_space,state_space,qs[a,:,:],xlabel="position", ylabel="normed velocity", title = "Action $(a)",figsize=(1000,1000))
        push!(hms,hm)
    end
    p = plot(hms...)
    return p
end

function load_and_plot_hms(policy_name)
    LP = load_policy(policy_name)
    state_space = collect(0:0.01:1)
    qs = get_q(LP,state_space)
    p = plot_hms(qs, state_space)
    display(p)
    savefig("./plotting/plots/tmp/$(policy_name).png")
end

#
# LP = load_policy("Wall")
# state_space = collect(0:0.01:1)
#
# ps = []
# hms = []
# for a in actions
#     # p = plot(state_space,state_space,q[a,:,:],st=:surface, xlabel="position", ylabel="normed velocity", camera=(-30,30), title = "Action $(a)")
#     # push!(ps,p)
#     hm = heatmap(state_space,state_space,q[a,:,:],xlabel="position", ylabel="normed velocity", title = "Action $(a)",figsize=(1000,1000))
#     push!(hms,hm)
# end
#
# plot(hms..., figsize=(1000,1000))
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using DEMC
using Random
using LinearAlgebra
using Distributions
Random.seed!(31953150)

# generate data
nobs = 100
npar = 3
Σ = ones(npar, npar)
for i = 1:npar
    for iprime = 1:i
        if i != iprime
            Σ[i, iprime] = 0.25
            Σ[iprime, i] = 0.25
        end
    end
end

DataDistribution = MvNormal(zeros(npar), Σ)
X = zeros(nobs, npar+1)
X[:,1] .= 1.
for iobs = 1:nobs
    X[iobs, 2:end] .= rand(DataDistribution)
end

β = rand(npar+1)
y = X*β + randn(nobs)

# log objective function - for correct std errors need to use optimal weights
log_obj(b) = -0.5 * sum((y.-X*b).^2)

# set up of DEMCz chain
Npar = length(β)
ndim = Npar
blockindex = [1:Npar] # parameter blocks: here choose all parameters to be updated simultaenously
Nblocks = length(blockindex)
eps_scale = 1e-5*ones(Npar) # scale of random error around DE update
γ = 2.38 # scale of DE update, 2.38 is the "optimal" number for a normal distribution
N = 3
K = 10
Z = randn((10*ndim, ndim))

# Number of iterations in Chain
Ngen = 2500
mc, Z = DEMC.demcz_sample(log_obj, Z, N, K, Ngen, Nblocks, blockindex, eps_scale, γ; verbose =false)

# drop first half of chain
keep = Int(round(Ngen/2))+1:Ngen
Ngen_burned = length(keep)
chain_burned = mc.chain[:,:,keep]
chainflat = DEMC.flatten_chain(chain_burned, N, Ngen_burned, Npar)'
bhat = mean(chainflat,dims=1)[:]
println("\n estimates: ", bhat, "\n dist to true: ", bhat - β)
# covariance of estimates
b, Σb = DEMC.mean_cov_chain(chain_burned, N, Ngen_burned, Npar)
# correct std errors for incorrect weighting matrix

# OLS
bols = (X'*X)\(X'*y)
ui = y .- X*bols
s2hat = (1/(nobs-npar)) * sum(ui.^2)
Σols = inv(X'X)*s2hat


text:
module HelloWorldServer

function allsources(dir)
  contents = String[]
  for (root, dirs, files) in walkdir(dir)
    push!.(Ref(contents), joinpath.(root, filter(x -> splitext(x)[2] == ".jl", files)))
  end
  return contents
end

srcdir = @__DIR__
for file in allsources(srcdir * "/models")
  println(file)
  include(file)
end

for file in allsources(srcdir * "/controllers")
  println(file)
  include(file)
end
end


text:
module TestActivate

using Test
using IonBase
using IonBase.ActivateCmd

include(joinpath(pkgdir(IonBase), "test", "utils.jl"))

function read_julia_version(julia::String)
    VersionNumber(readchomp(`$(IonBase.dot_ion("bin", julia)) --version`)[15:end])
end

@testset "test real installation" begin
    with_test_ion() do
        IonBase.activate("1.5.3")
        @test read_julia_version("julia") == v"1.5.3"
        IonBase.activate("stable")
        @test read_julia_version("julia") >= v"1.5.3"
        IonBase.activate("nightly")
        @test read_julia_version("julia") >= v"1.6-DEV"
    end
end

end # TestActivate


text:
x = 5
x1 = 0
y2 = 10
y = readline()
while (x1 < y) || (x1 == y)
    if x1 < 3
        println(2)
    elseif x1 == 3
        println(20)
    elseif x1 == 4
        println(200)
    else
        println(2000)
    end

    if (!(x1 > 3) && (y2 == 10))
        println(10)
    else
        println(100)
    end

    println((x1 > 3) && (x1 < y2))
    println(x1)
    x1 = x1 + 1
end


text:
# [WEBSITE]

supported_args(::GLVisualizeBackend) = merge_with_base_supported([
    # :annotations,
    # :background_color_legend, :background_color_inside, :background_color_outside,
    # :foreground_color_grid, :foreground_color_legend, :foreground_color_title,
    # :foreground_color_axis, :foreground_color_border, :foreground_color_guide, :foreground_color_text,
    # :label,
    # :linecolor, :linestyle, :linewidth, :linealpha,
    # :markershape, :markercolor, :markersize, :markeralpha,
    # :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
    # :fillrange, :fillcolor, :fillalpha,
    # :bins, :bar_width, :bar_edges, :bar_position,
    # :title, :title_location, :titlefont,
    # :window_title,
    # :guide, :lims, :ticks, :scale, :flip, :rotation,
    # :tickfont, :guidefont, :legendfont,
    # :grid, :legend, :colorbar,
    # :marker_z, :levels,
    # :ribbon, :quiver, :arrow,
    # :orientation,
    # :overwrite_figure,
    # :polar,
    # :normalize, :weights,
    # :contours, :aspect_ratio,
    # :match_dimensions,
    # :clims,
    # :inset_subplots,
  ])
supported_types(::GLVisualizeBackend) = [:surface, :scatter, :scatter3d, :path, :path3d]
supported_styles(::GLVisualizeBackend) = [:auto, :solid]
supported_markers(::GLVisualizeBackend) = [:none, :auto, :circle]
supported_scales(::GLVisualizeBackend) = [:identity]
is_subplot_supported(::GLVisualizeBackend) = false

# --------------------------------------------------------------------------------------

function _initialize_backend(::GLVisualizeBackend; kw...)
    @eval begin
        import GLVisualize, GeometryTypes, GLAbstraction, GLWindow
        import GeometryTypes: Point2f0, Point3f0
        export GLVisualize
    end
end

# ---------------------------------------------------------------------------

function _create_backend_figure(plt::Plot{GLVisualizeBackend})
    # init a window
    window = GLVisualize.glscreen()
    @async GLVisualize.renderloop(window)
    window
end

function gl_relative_size(plt::Plot{GLVisualizeBackend}, msize::Number)
    winsz = min(plt[:size]...)
    Float32(msize / winsz)
end

function gl_marker(shape::Symbol, msize::Number, _3d::Bool)
    GeometryTypes.HyperSphere((_3d ? Point3f0 : Point2f0)(0), msize)
end

function gl_display(plt::Plot{GLVisualizeBackend})
    window = plt.o
    for sp in plt.subplots
        # TODO: setup subplot

        for series in series_list(sp)
            d = series.d
            st = d[:seriestype]
            x, y = map(Float32, d[:x]), map(Float32, d[:y])
            msize = gl_relative_size(plt, d[:markersize])

            viz = if st == :surface
                ismatrix(x) || (x = repmat(x', length(y), 1))
                ismatrix(y) || (y = repmat(y, 1, length(x)))
                z = transpose_z(d, map(Float32, d[:z].surf), false)
                viz = GLVisualize.visualize((x, y, z), :surface)
                GLVisualize.view(viz, window, camera = :perspective)

            else
                points = if is3d(st)
                    z = map(Float32, d[:z])
                    Point3f0[Point3f0(xi,yi,zi) for (xi,yi,zi) in zip(x, y, z)]
                else
                    Point2f0[Point2f0(xi,yi) for (xi,yi) in zip(x, y)]
                end

                # markers?
                if st in (:scatter, :scatter3d) || d[:markershape] != :none
                    marker = gl_marker(d[:markershape], msize, is3d(st))
                    viz = GLVisualize.visualize((marker, points))
                    GLVisualize.view(viz, window, camera = :perspective)

                    # TODO: might need to switch to these forms later?
                    # GLVisualize.visualize((marker ,(x, y, z)))
                    #GLVisualize.visualize((marker , map(Point3f0, zip(x, y, z),
                    # billboard=true
                    #))
                end

                # paths
                lw = d[:linewidth]
                if !(st in (:scatter, :scatter3d)) && lw > 0
                    viz = GLVisualize.visualize(points, :lines) #, color=colors, model=rotation)
                    GLVisualize.view(viz, window, camera = :perspective)
                end
            end
        end
    end
    # GLAbstraction.center!(window)

    # TODO: render one frame at a time?  (no renderloop)
    # GLWindow.render_frame(window)
end

# ----------------------------------------------------------------

# function _writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{GLVisualizeBackend})
#     # TODO: write a png to io
# end

function _display(plt::Plot{GLVisualizeBackend})
    gl_display(plt)
end


text:
# auxiliary function to compute equal time correlator for u1-dm symmetry
function compute_eq_corr_u1_dm(
    file    :: HDF5.File,
    Λ       :: Float64
    ;
    verbose :: Bool = true
    )       :: Float64

    # load diagonal correlations
    m, χxx = read_χ(file, Λ, "xx", verbose = verbose)
    m, χzz = read_χ(file, Λ, "zz", verbose = verbose)

    # compute equal-time on-site correlation
    eq_corr = 0.0

    for i in 1 : length(m) - 1
        eq_corr += 2.0 * (m[i + 1] - m[i]) * (χxx[1, i + 1] + χxx[1, i])
        eq_corr += 1.0 * (m[i + 1] - m[i]) * (χzz[1, i + 1] + χzz[1, i])
    end

    # normalization from Fourier transformation
    eq_corr /= 2.0 * pi

    return eq_corr
end


text:
using .Photosynthesis

export AbstractFvCBCAssim, BallBerryPotentialCAssim, BallBerryCAssim, EmaxCAssim

export KooijmanWaterPotentialPhotosynthesis

abstract type AbstractFvCBCAssim <: AbstractCAssim end

@mix @flattenable @columns struct MixinFvCB{P,V,MG}
    vars::V        | false | Photosynthesis.EmaxVars() | _ | _ | _ | _
    photoparams::P | true  | nothing | _         | _           | _ | _
    SLA::MG        | true  | 24.0    | m^2*kg^-1 | (5.0, 30.0) | _ | "Specific leaf Area. Ferns: 17.4, Forbs: 26.2, Graminoids: 24.0, Shrubs: 9.10, Trees: 8.30"
end

"""
    BallBerryCAssim(vars, photoparams, SLA)

FCVB photosyntyhesis with Ball-Berry stomatal conductance.

Requires Photosynthesis.jl.

$(FIELDDOCTABLE)
"""
@MixinFvCB struct BallBerryCAssim{} <: AbstractFvCBCAssim end

@default BallBerryCAssim begin
    vars        | Photosynthesis.EmaxVars()
    photoparams | Photosynthesis.FvCBEnergyBalance(
                      photosynthesis=FvCBPhotosynthesis(
                          stomatal_conductance=BallBerryStomatalConductance(
                              soilmethod = PotentialSoilMethod()
                          ),
                          flux=Flux(),
                         ))
end

"""
    BallBerryCAssim(vars, photoparams, SLA)

FCVB photosyntyhesis with Ball-Berry stomatal conductance,
and a soil water potential model.

Requires Photosynthesis.jl.

$(FIELDDOCTABLE)
"""
@MixinFvCB struct BallBerryPotentialCAssim{} <: AbstractFvCBCAssim end

@default BallBerryPotentialCAssim begin
    vars        | Photosynthesis.EmaxVars()
    photoparams | Photosynthesis.FvCBEnergyBalance(
                      photosynthesis=FvCBPhotosynthesis(
                          stomatal_conductance=BallBerryStomatalConductance(
                              soilmethod = PotentialSoilMethod()
                          ),
                          flux=PotentialModifiedFlux(),
                      ))

end

photosynthesis(f::AbstractFvCBCAssim, o, u) = f.vars.aleaf * f.SLA * w_V(o)

apply_environment!(a::AbstractFvCBCAssim, o, u, shootenv, rootenv) = begin
    v = a.vars

    v.tair = airtemperature(shootenv)
    v.windspeed = windspeed(shootenv)
    v.rh = relhumidity(shootenv)
    v.rnet = radiation(shootenv)
    v.par = radiation(shootenv) * parconv
    v.vpd = vapour_pressure_deficit(v.tair, v.rh)
    # v.soilmoist = mean_soilwatercontent(rootenv)
    # swp = mean_soilwaterpotential(rootenv.microclimate, depth(o), rootenv.t)
    swp = if typeof(rootenv) <: MicroclimControl
        soilwaterpotential(rootenv)
    else
        layermax(soilwaterpotential(rootenv.microclimate), rootenv)
    end
    v.swp = swp
    set_swp!(o, v.swp)

    # This runs energy balance which contains photosynthesis
    # Later in assimilation we can just use the result
    enbal!(v, assimilation_pars(o).photoparams)
    set_soilcorrection!(o, v.fsoil)

    update_temp!(o, v.tleaf)
end

"""
    KooijmanSLAPhotosynthesis(vars, k_C_binding, k_O_binding, K_C, K_O, J_L_K, j_L_Amax, j_C_Amax, j_O_Amax)

Koojman photosynthesis formulation modified by soil water potential.

Requires Photosynthesis.jl. Untested and experimental.

$(FIELDDOCTABLE)
"""
@MixinKooijmanPhoto struct KooijmanWaterPotentialPhotosynthesis{PL} <: AbstractKooijmanPhoto
    potential_modifier::PL | Photosynthesis.ZhouPotentialDependence() | _ | _ | _ | _ | "Modify photosynthesis with a water potential model"
end

photosynthesis(f::KooijmanWaterPotentialPhotosynthesis, o, u) = begin
    va = f.vars
    mass_area_coef = w_V(o) * f.SLA

    j1_l = half_saturation(f.j_L_Amax, f.J_L_K, va.J_L_F) * mass_area_coef / 2

    # Modify CO2 and O2 intake by water availability to simulate stomatal closure.
    potentialcorrection = Photosynthesis.non_stomatal_potential_dependence(f.potential_modifier, va.soilwaterpotential)
    j1_c = half_saturation(f.j_C_Amax, f.K_C, va.X_C) * mass_area_coef * potentialcorrection * tempcorrection(o)
    j1_o = half_saturation(f.j_O_Amax, f.K_O, va.X_O) * mass_area_coef * potentialcorrection * tempcorrection(o)

    # photorespiration.
    bound_o = j1_o/f.k_O_binding # mol/mol
    bound_c = j1_c/f.k_C_binding # mol/mol

    # C flux
    j_c_intake = (j1_c - j1_o)

    j1_co = j1_c + j1_o
    co_l = j1_co/j1_l - j1_co/(j1_l + j1_co)

    j_c_intake / (1 + bound_c + bound_o + co_l)
end


text:
"""
This module provides tools to perform numerical integration on the sphere
for that, it provides different quadratures:

1. intSBazant21: Bazant quadrature 21 points
2. SlanSphere10k: Quadrature Slan 10000 points
3. intPpal: Quadrature with just 3 points in the ppal dir
"""
module SphereQuadratures

include("./Types.jl")
export SphereQuadraturervec

include("./Utils.jl")
export readquadrature

"""
This function computes the integral over the unit sphere of f(r1,r2,r3)
using the SphereQuadrature quad.

f: function of 3 director cosines f(r1,r2,r3)
quad: SphereQuadrature with fields r1v,r2v,r3v,wv
"""
function computeSphereQuad(f::Function, quad::SphereQuadraturervec)
    r1Vec = quad.r1Vec
    r2Vec = quad.r2Vec
    r3Vec = quad.r3Vec
    wVec = quad.wVec
    wsum = sum(wVec)
    Sintegral = (0.)*f(r1Vec[1], r1Vec[2], r1Vec[3])
    for i = 1 : length(r1Vec)
        r1 = r1Vec[i]
        r2 = r2Vec[i]
        r3 = r3Vec[i]
        w = wVec[i]/wsum
        Sintegral += f(r1,r2,r3)*w
    end
    return Sintegral
end

# Bazant Quadrature
BazantSphere21 = SphereQuadraturervec(readquadrature("BazantSphere21.txt", " ")...)

intSBazant21(f) = computeSphereQuad(f, BazantSphere21)

# Slan Quadrature 10k points

SlanSphere10k = SphereQuadraturervec(readquadrature("SlanSphere10k.txt", ";")...)

intSlan10k(f) = computeSphereQuad(f, SlanSphere10k)

# Ppal Quadrature

PpalQuadrature = SphereQuadraturervec(readquadrature("ppaldir.txt", " ")...)

intPpal(f) = computeSphereQuad(f, PpalQuadrature)

export computeSphereQuad, intSBazant21, intSlan10k, intPpal

export BazantSphere21, SlanSphere10k, PpalQuadrature

end #module


text:
#Julia Winston Plotting

using Winston
# fix the random seed so we have reproducible results
srand(111)
# generate a plot
pl = plot(cumsum(rand(100) .- 0.5), "g", cumsum(rand(100) .- 0.5), "b")
# display the plot
display(pl)

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter


text:
"ZDT1 test function"
function ZDT1(x::Vector{Float64})
    f1 = x[1]
    g = 1 + (9/(length(x)-1)) * sum(x[2:end])
    h = 1-sqrt(f1/g)
    return [f1, g*h]
end

function ZDT1(sol::Solution)
    return ZDT1(sol.x)
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
using InvertibleNetworks, NNlib, LinearAlgebra, Test

# Data
nx = 28
ny = 28
n_in = 8
batchsize = 3
n_hidden = 8
k = 3   # kernel size
s = 1   # stride
p = 1   # padding

###################################################################################################
# Test layer invertibility

function test_inv(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    print("\nInvertibility test hyperbolic layer with action=$(action)\n")
    # Data
    X0 = randn(Float32, nx, ny, n_in, batchsize)
    X1 = randn(Float32, nx, ny, n_in, batchsize)

    # Layer
    H1 = HyperbolicLayer(n_in, k, s, p; action=action, α=.2f0, n_hidden=n_hidden)

    Y0, Y1 = H1.forward(X0, X1)
    X0_, X1_ = H1.inverse(Y0, Y1)

    # Test invertibility
    @test isapprox(norm(X0 - X0_)/norm(X0_), 0f0; atol=1f-6)
    @test isapprox(norm(X1 - X1_)/norm(X1_), 0f0; atol=1f-6)

end

# Loss
function loss(H, X_prev_in, X_curr_in, Y_curr, Y_new)
    X_curr, X_new = H.forward(X_prev_in, X_curr_in)
    f = .5f0*norm(X_new - Y_new)^2 + .5*norm(X_curr - Y_curr)^2
    ΔX_curr = X_curr - Y_curr
    ΔX_new = X_new - Y_new
    ΔX_prev_, ΔX_curr_, X_prev_, X_curr_ = H.backward(ΔX_curr, ΔX_new, X_curr, X_new)

    # Check inverse is correct
    @test isapprox(norm(X_prev_ - X_prev_in)/norm(X_prev_in), 0f0; atol=1f-6)
    @test isapprox(norm(X_curr_ - X_curr_in)/norm(X_curr_in), 0f0; atol=1f-6)

    return f, ΔX_prev_, ΔX_curr_, H.W.grad, H.b.grad
end

function test_grad_X(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    H1 = HyperbolicLayer(n_in, k, s, p; action=action, α=.2f0, n_hidden=n_hidden)
    # Initial guess and residual
    X00 = randn(Float32, nx, ny, n_in, batchsize)
    X01 = randn(Float32, nx, ny, n_in, batchsize)
    dX0 = randn(Float32, nx, ny, n_in, batchsize)
    dX1 = randn(Float32, nx, ny, n_in, batchsize)

    # Observed data
    Y_curr, Y_new = H1.forward(X00, X01)

    # Gradient test for X
    maxiter = 10
    print("\nGradient test hyperbolic layer  input with action=$(action)\n")
    f0, ΔX_prev, ΔX_curr = loss(H1, X00, X01, Y_curr, Y_new)[1:3]
    h = .1f0
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)
    for j=1:maxiter
        f = loss(H1, X00 + h*dX0, X01 + h*dX1, Y_curr, Y_new)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX0, ΔX_prev) - h*dot(dX1, ΔX_curr))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

function test_grad_par(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    # Gradient test for W and b
    H0 = HyperbolicLayer(n_in, k, s, p; action=-1, α=.2f0, n_hidden=n_hidden)
    Hini = deepcopy(H0)

    dW = randn(Float32, size(H0.W.data))
    db = randn(Float32, size(H0.b.data))
    maxiter = 10

    # Data
    X0 = randn(Float32, nx, ny, n_in, batchsize)
    X1 = randn(Float32, nx, ny, n_in, batchsize)

    # Observed data
    Y_curr, Y_new = H0.forward(X0, X1)

    print("\nGradient test hyperbolic layer weights with action=$(action)\n")
    f0, ΔX_prev, ΔX_curr, ΔW, Δb = loss(H0, X0, X1, Y_curr, Y_new)

    h = .1f0
    err3 = zeros(Float32, maxiter)
    err4 = zeros(Float32, maxiter)

    for j=1:maxiter
        H0.W.data = Hini.W.data + h*dW
        H0.b.data = Hini.b.data + h*db
        f = loss(H0, X0, X1, Y_curr, Y_new)[1]
        err3[j] = abs(f - f0)
        err4[j] = abs(f - f0 - h*dot(dW, ΔW) - h*dot(db, Δb))
        print(err3[j], "; ", err4[j], "\n")
        h = h/2f0
    end
    @test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

end

for action in [1, -1, 0]
    test_inv(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    test_grad_X(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
    test_grad_par(nx, ny, n_in, batchsize, k, s, p, n_hidden, action)
end

###################################################################################################
# Jacobian-related tests

for action in [1, -1, 0]

    # Gradient test
    # Initialization
    HL = HyperbolicLayer(n_in, k, s, p; action=action, α=.2f0, n_hidden=8)
    θ = deepcopy(get_params(HL))
    HL0 = HyperbolicLayer(n_in, k, s, p; action=action, α=.2f0, n_hidden=8)
    θ0 = deepcopy(get_params(HL0))
    X1 = randn(Float32, nx, ny, n_in, batchsize)
    X2 = randn(Float32, nx, ny, n_in, batchsize)

    # Perturbation (normalized)
    dθ = θ-θ0; dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)
    dX1 = randn(Float32, nx, ny, n_in, batchsize); dX1 *= norm(X1)/norm(dX1)
    dX2 = randn(Float32, nx, ny, n_in, batchsize); dX2 *= norm(X2)/norm(dX2)

    # Jacobian eval
    dY1, dY2, Y1, Y2 = HL.jacobian(dX1, dX2, dθ, X1, X2)

    # Test
    print("\nJacobian test\n")
    h = 0.1f0
    maxiter = 5
    err5 = zeros(Float32, maxiter)
    err6 = zeros(Float32, maxiter)
    for j=1:maxiter
        set_params!(HL, θ+h*dθ)
        Y1_, Y2_ = HL.forward(X1+h*dX1, X2+h*dX2)
        err5[j] = sqrt(norm(Y1_ - Y1)^2f0+norm(Y2_ - Y2)^2f0)
        err6[j] = sqrt(norm(Y1_ - Y1 - h*dY1)^2f0+norm(Y2_ - Y2 - h*dY2)^2f0)
        print(err5[j], "; ", err6[j], "\n")
        h = h/2f0
    end

    @test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

    # Adjoint test
    set_params!(HL, θ)
    dY1, dY2, Y1, Y2 = HL.jacobian(dX1, dX2, dθ, X1, X2)
    dY1_ = randn(Float32, size(dY1)); dY2_ = randn(Float32, size(dY2));
    dX1_, dX2_, dθ_, _, _ = HL.adjointJacobian(dY1_, dY2_, Y1, Y2)
    a = dot(dY1, dY1_)+dot(dY2, dY2_)
    b = dot(dX1, dX1_)+dot(dX2, dX2_)+dot(dθ, dθ_)
    @test isapprox(a, b; rtol=1f-3)

end


text:
# adapted from LibIsolationForest
# see: https://github.com/msimms/LibIsolationForest/blob/master/julia/IsolationForest.jl
module IsolationForest

Feature = AbstractVector{T} where T <: AbstractFloat
Dataset = AbstractVector{<:Feature}

# Tree node, used internally.
mutable struct Node
    featureName::Int
    splitValue::AbstractFloat
    left::Union{Node, Nothing}
    right::Union{Node, Nothing}
end

# Isolation Forest.
mutable struct Forest
    numTrees::Int
    subSamplingSize::Int
    numFeatures::Int
    trees::AbstractArray{Union{Node, Nothing}}
    function Forest(numTrees, subSamplingSize, numFeatures, featureValues)
        forest = new(numTrees, subSamplingSize, numFeatures, [])
        for i = 1:forest.numTrees
            featureValues = deepcopy(featureValues)
            tree = create_tree(forest, featureValues, 0)
            if tree !== nothing
                push!(forest.trees, tree)
            end
        end
        return forest
    end
end

# Creates and returns a single tree. As this is a recursive function, depth indicates the current depth of the recursion.
function create_tree(forest::Forest, feature_values::Dataset, depth::Int)
    # Sanity check
    if forest.numFeatures <= 1
        return nothing
    end

    # If we've exceeded the maximum desired depth, then stop.
    if (forest.subSamplingSize > 0) && (depth >= forest.subSamplingSize)
        return nothing
    end

    # Randomly select a feature.
    randomly_selected_feature = rand(1:forest.numFeatures)

    # Randomly select a split value.
    feature_value_set = feature_values[randomly_selected_feature]
    feature_value_set_len = length(feature_value_set)

    if feature_value_set_len <= 1
        return nothing
    end
    split_value_index = rand(1:feature_value_set_len)
    split_value = feature_value_set[split_value_index]

    # Create a tree node to hold the split value.
    tree = Node(randomly_selected_feature, split_value, nothing, nothing)

    # Create two versions of the feature value set that we just used,
    # one for the left side of the tree and one for the right.
    temp_feature_values = feature_values

    # Create the left subtree.
    left_features = feature_value_set[1:split_value_index]

    temp_feature_values[randomly_selected_feature] = left_features
    tree.left = IsolationForest.create_tree(forest, temp_feature_values, depth + 1)

    # Create the right subtree.
    if split_value_index + 1 < feature_value_set_len
        right_features = feature_value_set[split_value_index + 1:feature_value_set_len]
        temp_feature_values[randomly_selected_feature] = right_features
        tree.right = IsolationForest.create_tree(forest, temp_feature_values, depth + 1)
    end

    return tree
end

# Scores the sample against the specified tree.
function score_sample_against_tree(tree::Node, features::Feature)
    depth = 0.0
    current_node = tree

    while current_node !== nothing
        found_feature = false

        # Find the next feature in the sample.
        for (current_feature_name, current_feature_value) in enumerate(features)

            # If the current node has the feature in question.
            if current_feature_name == current_node.featureName
                if current_feature_value < current_node.splitValue
                    current_node = current_node.left
                else
                    current_node = current_node.right
                end

                depth = depth + 1.0
                found_feature = true
                break
            end
        end

        # If the tree contained a feature not in the sample then take
        # both sides of the tree and average the scores together.
        if found_feature == false
            left_depth = depth + score_sample_against_tree(sample, current_node.left)
            right_depth = depth + score_sample_against_tree(sample, current_node.right)
            return (left_depth + right_depth) / 2.0
        end
    end

    return depth
end

# Scores the sample against the entire forest of trees. Result is the average path length.
function score_sample_against_forest(forest::Forest, features::Feature)
    num_trees = 0
    avg_path_len = 0.0

    for tree in forest.trees
        path_len = score_sample_against_tree(tree, features)
        if path_len > 0
            avg_path_len = avg_path_len + path_len
            num_trees = num_trees + 1
        end
    end

    if num_trees > 0
        avg_path_len = avg_path_len / num_trees
    end

    return avg_path_len
end

# Scores the sample against the entire forest of trees. Result is normalized so that values
# close to 1 indicate anomalies and values close to zero indicate normal values.
function H(i)
    return log(i) + 0.5772156649
end
function C(n)
    return 2 * H(n - 1) - (2 * (n - 1) / n)
end
function score_sample_against_forest_normalized(forest::Forest, features::Feature)

    # Compute the average path length for all valid trees.
    num_trees = 0
    avg_path_len = 0.0

    for tree in forest.trees
        path_len = score_sample_against_tree(tree, features)
        if path_len > 0
            avg_path_len = avg_path_len + path_len
            num_trees = num_trees + 1
        end
    end

    if num_trees > 0
        avg_path_len = avg_path_len / num_trees
    end

    # Normalize, per the original paper.
    score = 0.0
    if num_trees > 1.0
        exponent = -1.0 * (avg_path_len / C(num_trees))
        score = 2 ^ exponent
    end

    return score
end

end


text:
module DataLoader

    using CFTime, Dates
    using NCDatasets
    using DataStructures
    using Formatting

    mutable struct DataLoaderInfo
        data_dict :: Dict
        state_varnames  :: Array
        forcing_varnames :: Array
        joint_varnames   :: Array
        layers :: Union{Colon, UnitRange}
        year_rng  :: Array{Integer}

        function DataLoaderInfo(;
            hist_dir :: String,
            casename :: String,
            state_varnames :: Array,
            forcing_varnames :: Array,
            year_rng :: Array{Int64, 1},
            layers :: Union{Colon, UnitRange} = Colon(),
        )

            joint_varnames = Array{String}(undef,0)
            append!(joint_varnames, state_varnames, forcing_varnames)

            println(joint_varnames)

            data_dict = makeDataDict(
                hist_dir,
                casename,
                joint_varnames,
                year_rng;
                verbose = true,
            )

            return new(
                data_dict,
                state_varnames,
                forcing_varnames,
                joint_varnames,
                layers,
                year_rng,
            )

        end
    end

    function makeDataDict(
        hist_dir :: String,
        casename :: String,
        varnames :: Array{String, 1},
        year_rng :: Array{Int64, 1};
        verbose :: Bool = false,
    )

        expect_filenames = DataStructures.OrderedDict()

        println("I will also check a month prior to the begin year because we are using 12/31 of previous year as initial profile.")

        for y=(year_rng[1]-1):(year_rng[2]+1), m=1:12

            if (y == year_rng[1] - 1 && m != 12) || (y == year_rng[2]+1 && m != 1)
                continue
            end

            key = format("{:04d}-{:02d}", y, m)

            filename = format("$(casename).pop.h.daily.{:04d}-{:02d}-01.nc", y, m)
            full_path = joinpath(hist_dir, filename)
            verbose && println("Checking $(key)")

            if ! isfile(full_path)
                throw(ErrorException("Error: File $(filename) does not exist in folder $(hist_dir)."))
            end

            ds = Dataset(full_path, "r")
            dom = daysinmonth(DateTimeNoLeap, y, m)
            if ds.dim["time"] != dom
                throw(ErrorException("Error: File $(filename) does not have $(dom) records."))
            end

            for varname in varnames
                if ! haskey(ds, varname)
                    throw(ErrorException("Error: File $(filename) does not have variable `$(varname)`."))
                end
            end
            close(ds)

            expect_filenames[key] = Dict(
                :fp => full_path,
                :fn => filename,
                :y => y,
                :m => m,
            )

        end

        return expect_filenames

    end

    function getYMD(t::AbstractCFDateTime)
        return Dates.year(t), Dates.month(t), Dates.day(t)
    end

    function loadDataByDate(
        dli :: DataLoaderInfo,
        t   :: AbstractCFDateTime;
    )

        y, m, d = getYMD(t)

        key = format("{:04d}-{:02d}", y, m)

        full_path = dli.data_dict[key][:fp]

        ds = Dataset(full_path, "r")

        data = Dict()

        for varname in dli.joint_varnames

            v = ds[varname]
            dim = length(size(v))

            if dim == 4
                rng = (:, :, dli.layers, d)
            elseif dim == 3
                rng = (:, :, d)
            else
                throw(ErrorException("Unknown dimension: $(dim) of variable $(varname) in file $(full_path)"))
            end

            data[varname] = nomissing(ds[varname][rng...], NaN)

            if dim == 4
                data[varname] = permutedims(data[varname], [3, 1, 2])
            end
        end
        close(ds)

        return data

    end

    function loadInitDataAndForcing(
        dli :: DataLoaderInfo,
        t   :: AbstractCFDateTime,
    )

        t_curr = t
        data = loadDataByDate(dli, t_curr)

        t_prev = t_curr - Second(86400)
        data_prev = loadDataByDate(dli, t_prev)

        for varname in dli.state_varnames
            var      = data[varname]
            var_prev = data_prev[varname]

            @. var = (var + var_prev) / 2.0
        end

        return data

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module Factorions
using Memoize, ProgressMeter, Distributed

export is_factorion, fact_upper_bound, find_factorions, string_base
export find_in_range, find_in_parallel

@memoize function fast_fact(n::Int)
    return factorial(n)
end

"""
`sum_fact_digs(n)` returns the sum of the factorials of the
base-10 digits of `n`.
"""
function sum_fact_digs(n::Int, b::Int = 10)
    digs = digits(n, base = b)
    sum(map(fast_fact, digs))
end

"""
`is_factorion(n,b=10)` determines if `n` is a base-`b` *factorion*;
that is, if `n` is equal to the sum of the factorials of its base-`b` digits.
"""
function is_factorion(n::Int, b::Int = 10)::Bool
    if n <= 0
        return false
    end
    return n == sum_fact_digs(n, b)
end

"""
`fact_upper_bound(b=10)` returns an upper bound on the number of digits in
a base-`b` factorion.
"""
function fact_upper_bound(b::Int = 10)
    n::Int = 1
    fmax = fast_fact(b - 1)
    while b^(n - 1) < fmax * n
        n += 1
    end
    return n
end

function find_factorions(b::Int = 10)
    stop = b^fact_upper_bound(b)
    println("Trying $stop possibilities (up to $(fact_upper_bound(b)) digits)")
    list = Int[]
    PM = Progress(stop)
    for n = 1:stop
        if is_factorion(n, b)
            push!(list, n)
        end
        next!(PM)
    end
    for x in list
        println("$(string_base(x,b)) = $x")
    end

    return Set(list)

end

function find_in_range(rng, b::Int = 10)
    list = Int[]
    for x in rng
        if is_factorion(x, b)
            push!(list, x)
        end
    end
    return Set(list)
end

function find_in_parallel(b::Int = 10)
    stop = b^fact_upper_bound(b)
    println("Trying $stop possibilities (up to $(fact_upper_bound(b)) digits)")
    np = Threads.nthreads()
    #println(addprocs(np))
    results = Vector{Any}(undef, np)

    for j = 1:np
        rng = j:np:stop
        @info "Spawning process #$j"
        results[j] = @spawnat j + 1 find_in_range(rng, b)
    end
    sets = fetch.(results)
    X = union(sets...)
    for n in sort(collect(X))
        println("$n \t $(string_base(n,b))")
    end
    return union(sets...)
end

"""
`string_base(n,b)` renders the nonnegative integer `n` in
base `b` (we require `1 < b < 37`).
"""
function string_base(n::Int, b::Int = 10)
    symbs = "0123456789abcdefghijklmnopqrstuvwxyz"
    digs = digits(n, base = b)
    chars = [symbs[k+1] for k in digs]
    return join(reverse(chars)) * "_$b"
end

end # module


text:
function MultipleShooter(arg0::List, arg1::List, arg2::List, arg3::jdouble, arg4::jdouble)
    return MultipleShooter((List, List, List, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4)
end


text:
export rotz, simulation_neel

function rotz(x)
  #returns the rotation matrix that rotates the vector x onto the z axis.

  x = x/norm(x);
  n = cross(x,[0;0;1]);
  if norm(n)>0
    n = n/norm(n);
  else
    R = Matrix{Float64}(I, 3, 3)
    return R
  end
  a = acos(x[3]);
  n1 = n[1];
  n2 = n[2];
  n3 = n[3];
  R = [n1^2*(1-cos(a))+cos(a)  n1*n2*(1-cos(a))-n3*sin(a)  n1*n3*(1-cos(a))+n2*sin(a);
    n2*n1*(1-cos(a))+n3*sin(a)  n2^2*(1-cos(a))+cos(a)  n2*n3*(1-cos(a))-n1*sin(a);
    n3*n1*(1-cos(a))-n2*sin(a)  n3*n2*(1-cos(a))+n1*sin(a)  n3^2*(1-cos(a))+cos(a)];

  return R
end

struct NeelParams
  B::Function
  m_offset::SparseMatrixCSC{Complex{Float64},Int64}
  m_b3::SparseMatrixCSC{Complex{Float64},Int64}
  m_bp::SparseMatrixCSC{Complex{Float64},Int64}
  m_bm::SparseMatrixCSC{Complex{Float64},Int64}
  tau_N::Float64
  ytmp::Vector{ComplexF64}
  idx_offset::Vector{Int64}
  idx_b3::Vector{Int64}
  idx_bp::Vector{Int64}
  idx_bm::Vector{Int64}
end

function neel_odesys(y_out, y, p, t)
#@time begin
  B_ = p.B(t);
  B_1 = B_[1];
  B_2 = B_[2];
  B_3 = B_[3];

  #M = p.m_offset + B_3 .* p.m_b3 + (B_1 + 1im*B_2) .* p.m_bp + (B_1 - 1im*B_2) .* p.m_bm;
  #dydt = (1/(2*p.tau_N)) * M*y;

  y_out .= 0.0
  mul!(p.ytmp, p.m_bm, y)
  factor1 = (1.0/(2.0*p.tau_N)) *(B_1 - 1im*B_2)
  y_out .+= factor1 .* p.ytmp
  mul!(p.ytmp, p.m_bp, y)
  factor2 = (1.0/(2.0*p.tau_N)) *(B_1 + 1im*B_2)
  y_out .+= factor2 .* p.ytmp
  mul!(p.ytmp, p.m_b3, y)
  factor3 = (1.0/(2.0*p.tau_N)) * B_3
  y_out .+= factor3 .* p.ytmp
  mul!(p.ytmp, p.m_offset, y)
  factor4 = (1.0/(2.0*p.tau_N))
  y_out .+= factor4 .* p.ytmp

  #y_out .= dydt

#end
  return
end

function neel_odesys_jac(y_out, y, p, t)
#@time begin
  B_ = p.B(t);
  B_1 = B_[1];
  B_2 = B_[2];
  B_3 = B_[3];

  #=
  M = p.m_offset + B_3 .* p.m_b3 + (B_1 + 1im*B_2) .* p.m_bp + (B_1 - 1im*B_2).*p.m_bm;
  dydt = (1/(2*p.tau_N))*M;
  y_out .= dydt
  =#

  y_out.nzval .= 0
  factor1 = (1.0/(2.0*p.tau_N)) *(B_1 - 1im*B_2)
  y_out.nzval[p.idx_bm] .+= factor1 .* p.m_bm.nzval
  factor2 = (1.0/(2.0*p.tau_N)) *(B_1 + 1im*B_2)
  y_out.nzval[p.idx_bp] .+= factor2 .* p.m_bp.nzval
  factor3 = (1.0/(2.0*p.tau_N)) * B_3
  y_out.nzval[p.idx_b3] .+= factor3 .* p.m_b3.nzval
  factor4 = (1.0/(2.0*p.tau_N))
  y_out.nzval[p.idx_offset] .+= factor4 .* p.m_offset.nzval

  #end
  return
end

function simulation_neel(Bb, pr1, pr2, tau_N, alpha, t_vec, N)
  counter = 0;
  nz = 0;
  for r=0:N
    for q=-r:r
        counter = counter+1;

        nz = nz+1;
        if q!=-r
            nz = nz+1;
            if r!=0 && q!=r && q!=-r
                nz = nz+1;
            end
        end
        if q!=r
            nz = nz+1;
        end
        if r<(N-1)
            nz = nz+1;
        end
        if r>1 && q>(-r+1)
            nz = nz+1;
        end
        if q>(-r+1)
            nz = nz+1;
        end
        if q<(r-1)
            nz = nz+1;
        end
        if r<N
            nz = nz+1;
            nz = nz+1;
            nz = nz+1;
        end

    end
  end

  V = zeros(ComplexF64,nz);
  I = zeros(Int64,nz);
  J = zeros(Int64,nz);
  counter = 0;
  ind = 1;
  for r=0:N
    for q=-r:r
        counter = counter+1;
        I[ind] = counter;
        J[ind] = counter;
        V[ind] = -r*(r+1) + pr2* 2*(r^2+r-3*q^2)/((2*r+3)*(2*r-1));
        ind = ind+1;
        if q!=-r
            if r!=0 && q!=r
                I[ind] = counter-2*r;
                J[ind] = counter;
                V[ind] = 2*1im/alpha * pr2 * q*(r-q)/(2*r-1);
                ind = ind+1;
            end
        end
        if r<(N-1)
            I[ind] = counter+2*(r+1)+2*(r+2);
            J[ind] = counter;
            V[ind] = -pr2 * 2*r*(r+q+1)*(r+q+2)/((2*r+3)*(2*r+5));
            ind = ind+1;
        end
        if r>1 && q>(-r+1)
            I[ind] = counter-2*r-2*(r-1);
            J[ind] = counter;
            V[ind] = pr2 * 2*(r+1)*(r-q)*(r-q-1)/((2*r-3)*(2*r-1));
            ind = ind+1;
        end
        if r<N
            I[ind] = counter+2*(r+1);
            J[ind] = counter;
            V[ind] = 2*1im/alpha * pr2 * q*(r+q+1)/(2*r+3);
            ind = ind+1;
        end
    end
  end

  K = findfirst(isequal(0.0), I) - 1

  m_offset = sparse(J[1:K], I[1:K], V[1:K], (N+1)^2, (N+1)^2);
  dropzeros!(m_offset)

  V = zeros(ComplexF64,nz);
  I = zeros(Int64,nz);
  J = zeros(Int64,nz);
  counter = 0;
  ind = 1;
  for r=0:N
    for q=-r:r
        counter = counter+1;
        I[ind] = counter;
        J[ind] = counter;
        V[ind] = 1im/(2*alpha) * pr1 *2*q;
        ind = ind+1;
        if q!=-r
            if r!=0 && q!=r
                I[ind] = counter-2*r;
                J[ind] = counter;
                V[ind] = pr1*(r+1)*(r-q)/(2*r-1);
                ind = ind+1;
            end
        end
        if r<N
            I[ind] = counter+2*(r+1);
            J[ind] = counter;
            V[ind] = -pr1*r*(r+q+1)/(2*r+3);
            ind = ind+1;
        end
    end
  end

  K = findfirst(isequal(0.0), I) - 1
  m_b3 = sparse(J[1:K], I[1:K], V[1:K], (N+1)^2, (N+1)^2);
  dropzeros!(m_b3)

  V = zeros(ComplexF64,nz);
  I = zeros(Int64,nz);
  J = zeros(Int64,nz);
  counter = 0;
  ind = 1;
  for r=0:N
    for q=-r:r
        counter = counter+1;
        if q!=r
            I[ind] = counter+1;
            J[ind] = counter;
            V[ind] = 1im/(2*alpha) * pr1 *(r-q)*(r+q+1);
            ind = ind+1;
        end
        if q<(r-1)
            I[ind] = counter-2*r+1;
            J[ind] = counter;
            V[ind] = pr1 * (r+1)*(r-q)*(r-q-1)/(4*r-2);
            ind = ind+1;
        end
        if r<N
            I[ind] = counter+2*(r+1)+1;
            J[ind] = counter;
            V[ind] = pr1 * r*(r+q+1)*(r+q+2)/(4*r+6);
            ind = ind+1;
        end
    end
  end

  K = findfirst(isequal(0.0), I) - 1
  m_bp = sparse(J[1:K], I[1:K], V[1:K], (N+1)^2, (N+1)^2);
  dropzeros!(m_bp)

  V = zeros(ComplexF64,nz);
  I = zeros(Int64,nz);
  J = zeros(Int64,nz);
  counter = 0;
  ind = 1;
  for r=0:N
    for q=-r:r
        counter = counter+1;
        if q!=-r
            I[ind] = counter-1;
            J[ind] = counter;
            V[ind] = 1im/(2*alpha)*pr1;
            ind = ind+1;
        end
        if q>(-r+1)
            I[ind] = counter-2*r-1;
            J[ind] = counter;
            V[ind] = -pr1 * (r+1)/(4*r-2);
            ind = ind+1;

        end
        if r<N
            I[ind] = counter+2*(r+1)-1;
            J[ind] = counter;
            V[ind] = -pr1 * r/(4*r+6);
            ind = ind+1;
        end
    end
  end

  K = findfirst(isequal(0.0), I) - 1
  m_bm = sparse(J[1:K], I[1:K], V[1:K], (N+1)^2, (N+1)^2);
  dropzeros!(m_bm)

  # initial value
  y0 = zeros(ComplexF64, (N+1)^2);
  y0[1] = 1/(4*pi);
  ytmp = zeros(ComplexF64, (N+1)^2);

  # calculate the indices occuring in the jacobian
  tmp_off = deepcopy(m_offset)
  tmp_off.nzval .= 1
  tmp_b3 = deepcopy(m_b3)
  tmp_b3.nzval .= 1
  tmp_bp = deepcopy(m_bp)
  tmp_bp.nzval .= 1
  tmp_bm = deepcopy(m_bm)
  tmp_bm.nzval .= 1

  Mzero = tmp_off .+ tmp_b3 .+ tmp_bp .+ tmp_bm

  idx_offset = getIdxInM(Mzero, tmp_off)
  idx_b3 = getIdxInM(Mzero, tmp_b3)
  idx_bp = getIdxInM(Mzero, tmp_bp)
  idx_bm = getIdxInM(Mzero, tmp_bm)

  p = NeelParams(Bb, m_offset, m_b3, m_bp, m_bm, tau_N, ytmp, idx_offset, idx_b3, idx_bp, idx_bm)

  ff = ODEFunction(neel_odesys, jac = neel_odesys_jac, jac_prototype = Mzero)
  prob = ODEProblem(ff, y0, (t_vec[1],t_vec[end]), p)

  #@time sol = solve(prob, QNDF(), reltol=1e-3,abstol=1e-3)
  #@time sol = solve(prob, CVODE_BDF(), reltol=1e-3,abstol=1e-3)
  #@time sol = solve(prob, Rosenbrock23(autodiff=false), reltol=1e-3)
  #@time sol = solve(prob, Rodas5(autodiff=false), dt=1e-3, reltol=1e-3)
  @time sol = solve(prob, TRBDF2(autodiff=false), dt=1e-3, reltol=1e-3)
  #@time sol = solve(prob, Rodas5(autodiff=false), reltol=1e-3)

  @show sol.destats

  y = zeros(ComplexF64,length(t_vec), (N+1)^2)

  for ti=1:length(t_vec)
    y[ti, :] = sol(t_vec[ti])
  end
  return t_vec, y
end

function getIdxInM(M, V)
  M.nzval .= -1.0
  V.nzval .= 2.0
  M .+= V
  return findall(a-> real(a) > 0, M.nzval)
end


text:
@reexport module HARPSN

# Import types from EchelleTools
#import ..ManifestFormatEntry
import ..EchelleObservation, ..EchelleObservationSet, ..OrdersType
import ..AbstractInstrument #, ..HARPSNType
import ..instrument_symbol_list

using Reexport
@reexport using ..Blaze
#using ..Blaze.BlazeMap

struct HARPSNType <: AbstractInstrument end
push!(instrument_symbol_list,:HARPSN)

include("constants.jl")

# Set paths # TODO: Move to config file
base_path_default = "SET ME IN customizations.jl"
manifest_subdir_default = "."
data_subdir_default = "e2ds"
include("customizations.jl")
base_path = base_path_default

include("manifests.jl")
get_harpsn_base_path = Manifests.get_base_path
set_harpsn_base_path! = Manifests.set_base_path!
export get_harpsn_base_path, set_harpsn_base_path!
read_harpsn_manifest_file = Manifests.read
read_harpsn_manifest_files = Manifests.read
get_harpsn_data_filenames = Manifests.get_data_filenames
get_harpsn_data_filenames_and_bsrv = Manifests.get_data_filenames_and_bsrv
export read_harpsn_manifest_file, read_harpsn_manifest_files, get_harpsn_data_filenames
export get_harpsn_data_filenames_and_bsrv

include("readdata.jl")
read_harpsn_data = read_data
read_harpsn_blaze_list = read_blaze_list
export read_harpsn_data
export read_harpsn_blaze_list

end # module HARPSN


text:
export update_data!

"recursively applies new_data to data, overwriting information"
function update_data!(data::Dict{String,Any}, new_data::Dict{String,Any})
    if haskey(data, "per_unit") && haskey(new_data, "per_unit")
        if data["per_unit"] != new_data["per_unit"]
            error("update_data requires datasets in the same units, try make_per_unit and make_mixed_units")
        end
    else
        @warn "running update_data with data that does not include per_unit field, units may be incorrect"
    end
    _update_data!(data, new_data)
end

"recursive call of _update_data"
function _update_data!(data::Dict{String,Any}, new_data::Dict{String,Any})
    for (key, new_v) in new_data
        if haskey(data, key)
            v = data[key]
            if isa(v, Dict) && isa(new_v, Dict)
                _update_data!(v, new_v)
            else
                data[key] = new_v
            end
        else
            data[key] = new_v
        end
    end
end

"checks if a given network data is a multinetwork"
ismultinetwork(data::Dict{String,Any}) = (haskey(data, "multinetwork") && data["multinetwork"] == true)

"Transforms a single network into a multinetwork with several deepcopies of the original network"
function im_replicate(sn_data::Dict{String,Any}, count::Int; global_keys::Set{String} = Set{String}())
    @assert count > 0
    if ismultinetwork(sn_data)
        error("replicate can only be used on single networks")
    end

    if length(global_keys) <= 0
        @warn "deprecation warning, calls to replicate should explicitly specify a set of global_keys"
        # old default
        for (k,v) in sn_data
             if !(typeof(v) <: Dict)
                @warn "adding global key $(k)"
                push!(global_keys, k)
             end
        end
    end

    name = get(sn_data, "name", "anonymous")

    mn_data = Dict{String,Any}(
        "nw" => Dict{String,Any}()
    )

    mn_data["multinetwork"] = true

    sn_data_tmp = deepcopy(sn_data)
    for k in global_keys
        if haskey(sn_data_tmp, k)
            mn_data[k] = sn_data_tmp[k]
        end

        # note this is robust to cases where k is not present in sn_data_tmp
        delete!(sn_data_tmp, k)
    end

    mn_data["name"] = "$(count) replicates of $(name)"

    for n in 1:count
        mn_data["nw"]["$n"] = deepcopy(sn_data_tmp)
    end

    return mn_data
end

"builds a table of component data"
function component_table(data::Dict{String,Any}, component::String, fields::Vector{String})
    if ismultinetwork(data)
        return Dict((i, _component_table(nw_data, component, fields)) for (i,nw_data) in data["nw"])
    else
        return _component_table(data, component, fields)
    end
end
component_table(data::Dict{String,Any}, component::String, field::String) = component_table(data, component, [field])

function _component_table(data::Dict{String,Any}, component::String, fields::Vector{String})
    comps = data[component]
    if !_iscomponentdict(comps)
        @error "$(component) does not appear to refer to a component list"
    end

    items = []
    sorted_comps = sort(collect(comps); by=x->parse(Int, x[1]))
    for (i,comp) in sorted_comps
        push!(items, parse(Int, i))
    end
    for key in fields
        for (i,comp) in sorted_comps
            if haskey(comp, key)
                push!(items, comp[key])
            else
                push!(items, NaN)
            end
        end
    end

    return reshape(items, length(comps), length(fields)+1)
end

"prints the text summary for a data dictionary to stdout"
function print_summary(obj::Dict{String,Any}; kwargs...)
    summary(stdout, obj; kwargs...)
end

"prints the text summary for a data dictionary to IO"
function summary(io::IO, data::Dict{String,Any};
    float_precision::Int = 3,
    component_types_order = Dict(),
    component_parameter_order = Dict(),
    max_parameter_value = 999.0,
    component_status_parameters = Set(["status"])
    )

    if ismultinetwork(data)
        error("summary does not yet support multinetwork data")
    end

    component_types = []
    other_types = []

    println(io, _bold("Metadata"))
    for (k,v) in sort(collect(data); by=x->x[1])
        if typeof(v) <: Dict && _iscomponentdict(v)
            push!(component_types, k)
            continue
        end

        println(io, "  $(k): $(_value2string(v, float_precision))")
    end

    if length(component_types) > 0
        println(io, "")
        println(io, _bold("Table Counts"))
    end
    for k in sort(component_types, by=x->get(component_types_order, x, max_parameter_value))
        println(io, "  $(k): $(length(data[k]))")
    end

    for comp_type in sort(component_types, by=x->get(component_types_order, x, max_parameter_value))
        if length(data[comp_type]) <= 0
            continue
        end
        println(io, "")
        println(io, "")
        println(io, _bold("Table: $(comp_type)"))

        components = data[comp_type]

        display_components = Dict()
        active_components = Set()
        for (i, component) in components
            disp_comp = copy(component)

            status_found = false
            for (k, v) in disp_comp
                if k in component_status_parameters
                    status_found = true
                    if !(v == 0)
                        push!(active_components, i)
                    end
                end

                disp_comp[k] = _value2string(v, float_precision)
            end
            if !status_found
                push!(active_components, i)
            end

            display_components[i] = disp_comp
        end

        comp_key_sizes = Dict{String, Int}()
        default_values = Dict{String, Any}()
        for (i, component) in display_components
            # a special case for "index", for example when reading solution data
            if haskey(comp_key_sizes, "index")
                comp_key_sizes["index"] = max(comp_key_sizes["index"], length(i))
            else
                comp_key_sizes["index"] = length(i)
            end

            for (k, v) in component
                if haskey(comp_key_sizes, k)
                    comp_key_sizes[k] = max(comp_key_sizes[k], length(v))
                else
                    comp_key_sizes[k] = length(v)
                end

                if haskey(default_values, k)
                    if default_values[k] != v
                        default_values[k] = nothing
                    end
                else
                    default_values[k] = v
                end
            end
        end

        # when there is only one component nothing is default
        if length(display_components) == 1
            default_values = Dict{String, Any}()
        else
            default_values = Dict{String, Any}([x for x in default_values if !isa(x.second, Nothing)])
        end

        #display(default_values)

        # account for header width
        for (k, v) in comp_key_sizes
            comp_key_sizes[k] = max(length(k), v)
        end

        comp_id_pad = comp_key_sizes["index"] # not clear why this is offset so much
        delete!(comp_key_sizes, "index")
        comp_keys_ordered = sort([k for k in keys(comp_key_sizes) if !(haskey(default_values, k))], by=x->(get(component_parameter_order, x, max_parameter_value), x))

        header = join([lpad(k, comp_key_sizes[k]) for k in comp_keys_ordered], ", ")

        pad = " "^(comp_id_pad+2)
        println(io, "  $(pad)$(header)")
        for k in sort([k for k in keys(display_components)]; by=x->parse(Int, x))
            comp = display_components[k]
            items = []
            for ck in comp_keys_ordered
                if haskey(comp, ck)
                    push!(items, lpad("$(comp[ck])", comp_key_sizes[ck]))
                else
                    push!(items, lpad("--", comp_key_sizes[ck]))
                end
            end
            line = "  $(lpad(k, comp_id_pad)): $(join(items, ", "))"
            if k in active_components
                println(io, line)
            else
                println(io, _grey(line))
            end
        end

        if length(default_values) > 0
            println(io, "")
            println(io, "  default values:")
            for k in sort([k for k in keys(default_values)], by=x->(get(component_parameter_order, x, max_parameter_value), x))
                println(io, "    $(k): $(default_values[k])")
            end
        end
    end

end

"Attempts to determine if the given data is a component dictionary"
function _iscomponentdict(data::Dict)
    return all( typeof(comp) <: Dict for (i, comp) in data )
end

"Makes a string bold in the terminal"
function _bold(s::String)
    return "\033[1m$(s)\033[0m"
end

"""
Makes a string grey in the terminal, does not seem to work well on Windows terminals
more info can be found at https://en.wikipedia.org/wiki/ANSI_escape_code
"""
function _grey(s::String)
    return "\033[38;5;239m$(s)\033[0m"
end

"converts any value to a string, summarizes arrays and dicts"
function _value2string(v, float_precision::Int)
    if typeof(v) <: AbstractFloat
        return _float2string(v, float_precision)
    end
    if typeof(v) <: Array
        return "[($(length(v)))]"
    end
    if typeof(v) <: Dict
        return "{($(length(v)))}"
    end

    return "$(v)"
end

"""
converts a float value into a string of fixed precision

sprintf would do the job but this work around is needed because
sprintf cannot take format strings during runtime
"""
function _float2string(v::AbstractFloat, float_precision::Int)
    #str = "$(round(v; digits=float_precision))"
    str = "$(round(v; digits=float_precision))"
    lhs = length(split(str, '.')[1])
    return rpad(str, lhs + 1 + float_precision, "0")
end

"tests if two dicts are equal, up to floating point precision"
function compare_dict(d1, d2)
    for (k1,v1) in d1
        if !haskey(d2, k1)
            #println(k1)
            return false
        end
        v2 = d2[k1]

        if isa(v1, Number)
            if !compare_numbers(v1, v2)
                return false
            end
        elseif isa(v1, Array)
            if length(v1) != length(v2)
                return false
            end
            for i in 1:length(v1)
                if isa(v1[i], Number)
                    if !compare_numbers(v1[i], v2[i])
                        return false
                    end
                else
                    if v1 != v2
                        #println(v1, " ", v2)
                        return false
                    end
                end
            end
        elseif isa(v1, Dict)
            if !compare_dict(v1, v2)
                #println(v1, " ", v2)
                return false
            end
        else
            #println("2")
            if !isapprox(v1, v2)
                #println(v1, " ", v2)
                return false
            end
        end
    end
    return true
end

function Base.isapprox(a::Any, b::Any; kwargs...)
    return a == b
end

"tests if two numbers are equal, up to floating point precision"
function compare_numbers(v1, v2)
    if isnan(v1)
        #println("1.1")
        if !isnan(v2)
            #println(v1, " ", v2)
            return false
        end
    else
        #println("1.2")
        if !isapprox(v1, v2)
            #println(v1, " ", v2)
            return false
        end
    end
    return true
end


text:
push!(LOAD_PATH,"../src/")
using IntegratedOptics
using Documenter

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

makedocs(;
    modules=[IntegratedOptics],
    authors="Mohamed Kamal AbdElrahman",
    repo="https://github.com/MKAbdElrahman/IntegratedOptics.jl/blob/{commit}{path}#{line}",
    sitename="IntegratedOptics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MKAbdElrahman.github.io/IntegratedOptics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;devbranch = "main",
    repo="github.com/MKAbdElrahman/IntegratedOptics.jl",
)


text:
##############################
#
# obj-Files
#
##############################

function load(io::Stream{format"OBJ"}; facetype=GLTriangleFace,
        pointtype=Point3f, normaltype=Vec3f, uvtype=Vec2f)

    points, v_normals, uv, faces = pointtype[], normaltype[], uvtype[], facetype[]
    f_uv_n_faces = (faces, GLTriangleFace[], GLTriangleFace[])
    last_command = ""
    attrib_type  = nothing
    for full_line in eachline(stream(io))
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(full_line))
        !isascii(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !all(iscntrl, line) #ignore comments
            lines = split(line)
            command = popfirst!(lines) #first is the command, rest the data

            if "v" == command # mesh always has vertices
                push!(points, Point{3, Float32}(parse.(Float32, lines))) # should automatically convert to the right type in vertices(mesh)
            elseif "vn" == command
                push!(v_normals, Vec3f(parse.(Float32, lines)))
            elseif "vt" == command
                if length(lines) == 2
                    push!(uv, Vec2f(parse.(Float32, lines)))
                elseif length(lines) == 3
                    push!(uv, Vec3f(parse.(Float32, lines)))
                else
                    error("Unknown UVW coordinate: $lines")
                end
            elseif "f" == command # mesh always has faces
                if any(x-> occursin("//", x), lines)
                    fs = process_face_normal(lines)
                elseif any(x-> occursin("/", x), lines)
                    fs = process_face_uv_or_normal(lines)
                else
                    append!(faces, triangulated_faces(facetype, lines))
                    continue
                end
                for i = 1:length(first(fs))
                    append!(f_uv_n_faces[i], triangulated_faces(facetype, getindex.(fs, i)))
                end
            else
                #TODO
            end
        end
    end

    point_attributes = Dict{Symbol, Any}()
    non_empty_faces = filtertuple(!isempty, f_uv_n_faces)

    if length(non_empty_faces) > 1
        N = length(points)
        void = tuple((_typemax(eltype(facetype)) for _ in 1:length(non_empty_faces))...)
        vertices = fill(void, N)

        if !isempty(v_normals)
            point_attributes[:normals] = Vector{normaltype}(undef, N)
        end
        if !isempty(uv)
            point_attributes[:uv] = Vector{uvtype}(undef, N)
        end

        for (k, fs) in enumerate(zip(non_empty_faces...))
            f = collect(first(fs)) # position indices
            for i in eachindex(non_empty_faces)
                l = 2
                vertex = getindex.(fs, i) # one of each indices (pos/uv/normal)

                if vertices[vertex[1]] == void
                    # Replace void
                    vertices[vertex[1]] = vertex
                    f[i] = vertex[1]
                    if !isempty(uv)
                        point_attributes[:uv][vertex[1]] = uv[vertex[l]]
                        l += 1
                    end
                    if !isempty(v_normals)
                        point_attributes[:normals][vertex[1]] = v_normals[vertex[l]]
                    end
                elseif vertices[vertex[1]] == vertex
                    # vertex is correct, nothing to replace
                    f[i] = vertex[1]
                else
                    @views j = findfirst(==(vertex), vertices[N+1:end])
                    if j === nothing
                        # vertex is unique, add it as a new one and adjust
                        # points, uv, normals
                        push!(vertices, vertex)
                        f[i] = length(vertices)
                        push!(points, points[vertex[1]])
                        if !isempty(uv)
                            push!(point_attributes[:uv], uv[vertex[l]])
                            l += 1
                        end
                        if !isempty(v_normals)
                            push!(point_attributes[:normals], v_normals[vertex[l]])
                        end
                    else
                        # vertex has already been added, adjust face
                        # (points, uv, normals correct because they've been pushed)
                        f[i] = j + N
                    end
                end
            end
            # remap indices
            faces[k] = facetype(f)
        end
    else # we have vertex indexing - no need to remap
        if !isempty(v_normals)
            point_attributes[:normals] = v_normals
        end
        if !isempty(uv)
            point_attributes[:uv] = uv
        end
    end

    return Mesh(meta(points; point_attributes...), faces)
end

# of form "faces v1 v2 v3 ....""
process_face(lines::Vector) = (lines,) # just put it in the same format as the others
# of form "faces v1//vn1 v2//vn2 v3//vn3 ..."
process_face_normal(lines::Vector) = split.(lines, "//")
# of form "faces v1/vt1 v2/vt2 v3/vt3 ..." or of form "faces v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...."
process_face_uv_or_normal(lines::Vector) = split.(lines, ('/',))

function triangulated_faces(::Type{Tf}, vertex_indices::Vector) where {Tf}
    poly_face = NgonFace{length(vertex_indices), UInt32}(parse.(UInt32, vertex_indices))
    return convert_simplex(Tf, poly_face)
end

function _typemax(::Type{OffsetInteger{O, T}}) where {O, T}
    typemax(T)
end

function save(f::Stream{format"OBJ"}, mesh::AbstractMesh)
    io = stream(f)
    for p in decompose(Point3f, mesh)
        println(io, "v ", p[1], " ", p[2], " ", p[3])
    end

    if hasproperty(mesh, :uv)
        for uv in mesh.uv
            println(io, "vt ", uv[1], " ", uv[2])
        end
    end

    if hasproperty(mesh, :normals)
        for n in decompose_normals(mesh)
            println(io, "vn ", n[1], " ", n[2], " ", n[3])
        end
    end

    for f in decompose(GLTriangleFace, mesh)
        println(io, "f ", convert(Int, f[1]), " ", convert(Int, f[2]), " ", convert(Int, f[3]))
    end
end


text:
include("imports.jl")

@testset ExtendedTestSet "catafirst" begin
    x = XNode(:body, [XLeaf(1), XLeaf(2)])

    x_ = cata(x, SelectLeaf()) do leaf
        return XLeaf(leaf[] + 1)
    end
    @test children(x_)[1][] == 2
    @test children(x_)[2][] == 3

    x__ = catafirst(x, SelectLeaf()) do leaf
        return XLeaf(leaf[] + 1)
    end
    @test children(x__)[1][] == 2
    @test children(x__)[2][] == 2
end

@testset ExtendedTestSet "replace" begin
    x = XNode(:body, [XLeaf(1), XLeaf(2)])
    node = XNode(:body)
    @test Pollen.replace(x, node, SelectNode()) == node

    x_ = Pollen.replacefirst(x, node, SelectLeaf())
    @test tag(x_) == :body
    @test tag(children(x_)[1]) == :body
end

@testset ExtendedTestSet "insert" begin
    x = XNode(:body, [XLeaf(1), XLeaf(2)])
    @testset ExtendedTestSet "NthChild" begin
        x_ = insert(x, XLeaf(0), NthChild(1, SelectNode()))
        @test children(x_) == XLeaf.(0:2)
        @test insert(x, XLeaf(0), NthChild(1, SelectNode())) == insertfirst(x, XLeaf(0), NthChild(1, SelectNode()))
    end

    @testset ExtendedTestSet "Before" begin
        x_ = insert(x, XLeaf(0), Before(SelectLeaf()))
        @test children(x_) == XLeaf.(0:2)
    end

    @testset ExtendedTestSet "Before" begin
        x_ = insert(x, XLeaf(0), After(SelectLeaf()))
        @test children(x_) == XLeaf.([1, 0, 2])
    end
end

@testset ExtendedTestSet "gettext" begin
    x = XNode(:body, XLeaf.(["Hello", " ", "World"]))
    @test Pollen.gettext(x) == "Hello World"
end


text:
# diff from CosSpace -> SinSpace

function cosspacediff{T<:Number}(v::Vector{T})
    if length(v)==1
        w = zeros(T,1)
    else
        w = Array{T}(length(v)-1)
        for k=1:length(v)-1
            @inbounds w[k] = -k*v[k+1]
        end
    end

    w
end

# diff from SinSpace -> CosSpace

function sinspacediff{T<:Number}(v::Vector{T})
    w = Array{T}(length(v)+1)
    w[1] = zero(T)
    for k=1:length(v)
        @inbounds w[k+1] = k*v[k]
    end

    w
end

# diff from Fourier -> Fourier

function fourierdiff{T<:Number}(v::Vector{T})
    n = 2(length(v)÷2)+1
    w = Array{T}(n)
    w[1] = zero(T)
    n > 1 && (w[n-1] = zero(T))
    for k=1:n÷2-1
        @inbounds w[2k] = -k*v[2k+1]
        @inbounds w[2k+1] = k*v[2k]
    end
    n > 1 && (w[n] = (n÷2)*v[n-1]; n == length(v) && (w[n-1] = -(n÷2)*v[n]))

    w
end

# diff from Taylor -> Taylor

function taylor_diff{T<:Number}(v::Vector{T})
    w = Array{T}(length(v))
    for k=1:length(v)
        @inbounds w[k] = (k-1)*v[k]
    end

    w
end

# diff from Hardy{false} -> Hardy{false}

function hardyfalse_diff{T<:Number}(v::Vector{T})
    w = Array{T}(length(v))
    for k=1:length(v)
        @inbounds w[k] = -k*v[k]
    end

    w
end

# diff from Laurent -> Laurent

function laurentdiff{T<:Number}(v::Vector{T})
    w = Array{T}(length(v))
    w[1] = zero(T)
    for k=1:length(v)÷2
        @inbounds w[2k] = -k*v[2k]
        @inbounds w[2k+1] = k*v[2k+1]
    end

    w
end


text:
import LinearAlgebra.logdet

struct LogDetAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}

    function LogDetAtom(x::AbstractExpr)
        children = (x,)
        return new(:logdet, hash(children), children, (1, 1))
    end
end

function sign(x::LogDetAtom)
    return NoSign()
end

function monotonicity(x::LogDetAtom)
    return (NoMonotonicity(),)
end

function curvature(x::LogDetAtom)
    return ConcaveVexity()
end

function evaluate(x::LogDetAtom)
    return log(det(evaluate(x.children[1])))
end

function conic_form!(x::LogDetAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        A = x.children[1]
        D = Variable(size(A)) # diagonal matrix
        U = Variable(size(A)) # upper triangular matrix

        # objective given by the sum of the log of diagonal matrix D
        objective = conic_form!(sum(log(diag(D))), unique_conic_forms)

        # enforce D and U to be diagonal and upper triangular
        for i in 1:A.size[1]
            for j in 1:A.size[2]
                if i != j
                    conic_form!(D[i, j] == 0, unique_conic_forms)
                end
                if i > j
                    conic_form!(U[i, j] == 0, unique_conic_forms)
                end
            end
        end

        # diagonals of D and U are the same
        conic_form!(diag(D) == diag(U), unique_conic_forms)

        # A and [D U; U' A] need to be positive semidefinite
        conic_form!(isposdef(A), unique_conic_forms)
        conic_form!(isposdef([D U; U' A]), unique_conic_forms)

        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

logdet(x::AbstractExpr) = LogDetAtom(x)


text:
using HorizonSideRobots
include("Functions.jl")
include("Structs.jl")

function cross_border(robot::Robot)
    cross_robot = PutmarkersRobot(robot)
    sides = [Nord, West, Sud, Ost]
    for side in sides
        m = movements!(cross_robot, side)
        rside = reverse_side(side)
        movements!(cross_robot, rside, m)
    end
    putmarker!(robot)
end

max_stars_repo_name: Mananhina/mirea-progs

max_issues_repo_name: Mananhina/mirea-progs

max_forks_repo_name: Mananhina/mirea-progs


text:
using DelayDiffEq
using Test

@testset "in-place" begin
  # define functions (Hutchinson's equation)
  function f(du, u, h, p, t)
    du[1] = u[1] * (1  - h(p, t - 1)[1])
    nothing
  end

  njacs = Ref(0)
  function jac(J, u, h, p, t)
    njacs[] += 1
    J[1, 1] = 1 - h(p, t - 1)[1]
    nothing
  end

  nWfacts = Ref(0)
  function Wfact(W, u, h, p, dtgamma, t)
    nWfacts[] += 1
    W[1,1] = dtgamma * (1 - h(p, t - 1)[1]) - 1
    nothing
  end

  nWfact_ts = Ref(0)
  function Wfact_t(W, u, h, p, dtgamma, t)
    nWfact_ts[] += 1
    W[1,1] = 1 - h(p, t - 1)[1] - inv(dtgamma)
    nothing
  end

  h(p, t) = [0.0]

  # define problems
  prob = DDEProblem(DDEFunction{true}(f), [1.0], h, (0.0, 40.0); constant_lags = [1])
  prob_jac = remake(prob; f = DDEFunction{true}(f; jac = jac))
  prob_Wfact = remake(prob; f = DDEFunction{true}(f; Wfact = Wfact))
  prob_Wfact_t = remake(prob; f = DDEFunction{true}(f; Wfact_t = Wfact_t))

  # compute solutions
  for alg in (Rosenbrock23(), TRBDF2())
    sol = solve(prob, MethodOfSteps(alg))

    ## Jacobian
    njacs[] = 0
    sol_jac = solve(prob_jac, MethodOfSteps(alg))

    # check number of function evaluations
    @test !iszero(njacs[])
    @test njacs[] == sol_jac.destats.njacs
    if alg isa Rosenbrock23
      @test njacs[] == sol_jac.destats.nw
    else
      @test_broken njacs[] == sol_jac.destats.nw
    end

    # check resulting solution
    @test sol.t ≈ sol_jac.t
    @test sol.u ≈ sol_jac.u

    ## Wfact
    nWfacts[] = 0
    sol_Wfact = solve(prob_Wfact, MethodOfSteps(alg))

    # check number of function evaluations
    if alg isa Rosenbrock23
      @test !iszero(nWfacts[])
      @test nWfacts[] == njacs[]
      @test iszero(sol_Wfact.destats.njacs)
    else
      @test_broken !iszero(nWfacts[])
      @test_broken nWfacts[] == njacs[]
      @test_broken iszero(sol_Wfact.destats.njacs)
    end
    @test_broken nWfacts[] == sol_Wfact.destats.nw

    # check resulting solution
    @test sol.t ≈ sol_Wfact.t
    @test sol.u ≈ sol_Wfact.u

    ## Wfact_t
    nWfact_ts[] = 0
    sol_Wfact_t = solve(prob_Wfact_t, MethodOfSteps(alg))

    # check number of function evaluations
    if alg isa Rosenbrock23
      @test_broken !iszero(nWfact_ts[])
      @test_broken nWfact_ts[] == njacs[]
      @test_broken iszero(sol_Wfact_t.destats.njacs)
    else
      @test !iszero(nWfact_ts[])
      @test_broken nWfact_ts[] == njacs[]
      @test iszero(sol_Wfact_t.destats.njacs)
    end
    @test_broken nWfact_ts[] == sol_Wfact_t.destats.nw

    # check resulting solution
    if alg isa Rosenbrock23
      @test sol.t ≈ sol_Wfact_t.t
      @test sol.u ≈ sol_Wfact_t.u
    else
      @test_broken sol.t ≈ sol_Wfact_t.t
      @test_broken sol.u ≈ sol_Wfact_t.u
    end
  end
end

@testset "out-of-place" begin
  # define functions (Hutchinson's equation)
  f(u, h, p, t) = u[1] .* (1  .- h(p, t - 1))

  njacs = Ref(0)
  function jac(u, h, p, t)
    njacs[] += 1
    reshape(1 .- h(p, t - 1), 1, 1)
  end

  nWfacts = Ref(0)
  function Wfact(u, h, p, dtgamma, t)
    nWfacts[] += 1
    reshape(dtgamma .* (1 .- h(p, t - 1)) .- 1, 1, 1)
  end

  nWfact_ts = Ref(0)
  function Wfact_t(u, h, p, dtgamma, t)
    nWfact_ts[] += 1
    reshape((1 - inv(dtgamma)) .- h(p, t - 1), 1, 1)
  end

  h(p, t) = [0.0]

  # define problems
  prob = DDEProblem(DDEFunction{false}(f), [1.0], h, (0.0, 40.0); constant_lags = [1])
  prob_jac = remake(prob; f = DDEFunction{false}(f; jac = jac))
  prob_Wfact = remake(prob; f = DDEFunction{false}(f; Wfact = Wfact))
  prob_Wfact_t = remake(prob; f = DDEFunction{false}(f; Wfact_t = Wfact_t))

  # compute solutions
  for alg in (Rosenbrock23(), TRBDF2())
    sol = solve(prob, MethodOfSteps(alg))

    ## Jacobian
    njacs[] = 0
    sol_jac = solve(prob_jac, MethodOfSteps(alg))

    # check number of function evaluations
    @test !iszero(njacs[])
    @test_broken njacs[] == sol_jac.destats.njacs
    @test_broken njacs[] == sol_jac.destats.nw

    # check resulting solution
    @test sol.t ≈ sol_jac.t
    @test sol.u ≈ sol_jac.u

    ## Wfact
    nWfacts[] = 0
    sol_Wfact = solve(prob_Wfact, MethodOfSteps(alg))

    # check number of function evaluations
    @test_broken !iszero(nWfacts[])
    @test_broken nWfacts[] == njacs[]
    @test_broken iszero(sol_Wfact.destats.njacs)
    @test_broken nWfacts[] == sol_Wfact.destats.nw

    # check resulting solution
    @test sol.t ≈ sol_Wfact.t
    @test sol.u ≈ sol_Wfact.u

    ## Wfact_t
    nWfact_ts[] = 0
    sol_Wfact_t = solve(prob_Wfact_t, MethodOfSteps(alg))

    # check number of function evaluations
    @test_broken !iszero(nWfact_ts[])
    @test_broken nWfact_ts[] == njacs[]
    @test_broken iszero(sol_Wfact_ts.destats.njacs)
    @test_broken nWfact_ts[] == sol_Wfact_t.destats.nw

    # check resulting solution
    @test sol.t ≈ sol_Wfact_t.t
    @test sol.u ≈ sol_Wfact_t.u
  end
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using Pkg;
Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, ArviZ
import MCMCChains

makedocs(
    modules = [ArviZ],
    sitename = "ArviZ.jl",
    pages = [
        "Home" => "index.md",
        "Quickstart" => "quickstart.md",
        "API" => "api.md",
        "Reference" => "reference.md",
    ],
    checkdocs = :exports,
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "CI"),
        assets = ["assets/favicon.ico"],
        sidebar_sitename = false,
        canonical = "stable",
    ),
    linkcheck = true,
)

deploydocs(
    repo = "github.com/arviz-devs/ArviZ.jl.git",
    push_preview = true,
)


text:
@inline function calcpoint(blockIdx, blockDim, threadIdx, size)
    i = (blockIdx - 1) * blockDim + threadIdx
    return i, Float32(i)
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,1}, texture::CuDeviceTexture{<:Any,1})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    @inbounds dst[i] = texture[u]
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,2}, texture::CuDeviceTexture{<:Any,2})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    @inbounds dst[i,j] = texture[u,v]
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{<:Any,3}, texture::CuDeviceTexture{<:Any,3})
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    k, w = calcpoint(blockIdx().z, blockDim().z, threadIdx().z, size(dst)[3])
    @inbounds dst[i,j,k] = texture[u,v,w]
    return nothing
end

function fetch_all(texture)
    dims = size(texture)
    d_out = CuArray{eltype(texture)}(undef, dims...)
    @cuda threads = dims kernel_texture_warp_native(d_out, texture)
    d_out
end

@testset "CuTextureArray(::CuArray)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)
    d_a3D = CuArray(a3D)

    texarr1D = CuTextureArray(d_a1D)
    tex1D = CuTexture(texarr1D)
    @test fetch_all(tex1D) == d_a1D

    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D

    texarr3D = CuTextureArray(d_a3D)
    tex3D = CuTexture(texarr3D)
    @test fetch_all(tex3D) == d_a3D
end

@testset "CuTextureArray(::Array)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end

    texarr1D = CuTextureArray(a1D)
    copyto!(texarr1D, a1D)
    tex1D = CuTexture(texarr1D)
    @test Array(fetch_all(tex1D)) == a1D

    texarr2D = CuTextureArray(a2D)
    tex2D = CuTexture(texarr2D)
    @test Array(fetch_all(tex2D)) == a2D

    tex2D_dir = CuTexture(CuTextureArray(a2D))
    @test Array(fetch_all(tex2D_dir)) == a2D

    texarr3D = CuTextureArray(a3D)
    tex3D = CuTexture(texarr3D)
    @test Array(fetch_all(tex3D)) == a3D
end

@testset "CuTexture(::CuArray)" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)

    # Strangely, this is not working
    texwrap1D = CuTexture(d_a1D)
    @test_broken fetch_all(texwrap1D) == d_a1D

    # This works as long as d_a2D is well pitched
    texwrap2D = CuTexture(d_a2D)
    @test fetch_all(texwrap2D) == d_a2D
end

@testset "type support" begin
    for T in (Int32, UInt32, Int16, UInt16, Int8, UInt8, Float32, Float16)
        testheight, testwidth, testdepth = 32, 32, 4
        a2D = rand(T, testheight, testwidth)
        d_a2D = CuArray(a2D)

        # Using CuTextureArray
        tex_2D = CuTexture(d_a2D)
        @test fetch_all(tex_2D) == d_a2D

        # Wrapping CuArray
        # This works as long as d_a2D is well pitched
        texwrap_2D = CuTexture(d_a2D)
        @test fetch_all(texwrap_2D) == d_a2D
    end
end

@testset "multiple channels" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int32(i), Int32(j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D

    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int16(i), Int16(j), Int16(i + j), Int16(i - j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray(d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D
end


text:
"""
    violin(x, y; kwargs...)
Draw a violin plot.
# Arguments
- `x`: positions of the categories
- `y`: variables whose density is computed
# Keywords
- `orientation=:vertical`: orientation of the violins (`:vertical` or `:horizontal`)
- `width=0.8`: width of the violin
- `show_median=true`: show median as midline
- `side=:both`: specify `:left` or `:right` to only plot the violin on one side
- `datalimits`: specify values to trim the `violin`. Can be a `Tuple` or a `Function` (e.g. `datalimits=extrema`)
"""
@recipe(Violin, x, y) do scene
    Theme(;
        default_theme(scene, Poly)...,
        npoints = 200,
        boundary = automatic,
        bandwidth = automatic,
        side = :both,
        width = automatic,
        dodge = automatic,
        n_dodge = automatic,
        x_gap = 0.2,
        dodge_gap = 0.03,
        datalimits = (-Inf, Inf),
        max_density = automatic,
        show_median = false,
        mediancolor = theme(scene, :linecolor),
        medianlinewidth = theme(scene, :linewidth),
    )
end

conversion_trait(x::Type{<:Violin}) = SampleBased()

getuniquevalue(v, idxs) = v

function getuniquevalue(v::AbstractVector, idxs)
    u = view(v, idxs)
    f = first(u)
    msg = "Collection must have the same value across all indices"
    all(isequal(f), u) || throw(ArgumentError(msg))
    return f
end

function plot!(plot::Violin)
    x, y = plot[1], plot[2]
    args = @extract plot (width, side, color, show_median, npoints, boundary, bandwidth,
        datalimits, max_density, dodge, n_dodge, x_gap, dodge_gap)
    signals = lift(x, y, args...) do x, y, width, vside, color, show_median, n, bound, bw, limits, max_density, dodge, n_dodge, x_gap, dodge_gap
        x̂, violinwidth = xw_from_dodge(x, width, 1, x_gap, dodge, n_dodge, dodge_gap)

        # Allow `side` to be either scalar or vector
        sides = broadcast(x̂, vside) do _, s
            return s == :left ? - 1 : s == :right ? 1 : 0
        end

        sa = StructArray((x = x̂, side = sides))

        specs = map(StructArrays.finduniquesorted(sa)) do (key, idxs)
            v = view(y, idxs)
            k = KernelDensity.kde(v;
                npoints = n,
                (bound === automatic ? NamedTuple() : (boundary = bound,))...,
                (bw === automatic ? NamedTuple() : (bandwidth = bw,))...,
            )
            l1, l2 = limits isa Function ? limits(v) : limits
            i1, i2 = searchsortedfirst(k.x, l1), searchsortedlast(k.x, l2)
            kde = (x = view(k.x, i1:i2), density = view(k.density, i1:i2))
            c = getuniquevalue(color, idxs)
            return (x = key.x, side = key.side, color = to_color(c), kde = kde, median = median(v))
        end

        max = if max_density === automatic
            maximum(specs) do spec
                _, max = extrema_nan(spec.kde.density)
                return max
            end
        else
            max_density
        end

        vertices = Vector{Point2f}[]
        lines = Pair{Point2f, Point2f}[]
        colors = RGBA{Float32}[]

        for spec in specs
            scale = 0.5*violinwidth/max
            xl = reverse(spec.x .- spec.kde.density .* scale)
            xr = spec.x .+ spec.kde.density .* scale
            yl = reverse(spec.kde.x)
            yr = spec.kde.x

            x_coord, y_coord = if spec.side == -1 # left violin
                [spec.x; xl; spec.x], [yl[1]; yl; yl[end]]
            elseif spec.side == 1 # right violin
                [spec.x; xr; spec.x], [yr[1]; yr; yr[end]]
            else
                [spec.x; xr; spec.x; xl], [yr[1]; yr; yl[1]; yl]
            end
            verts = Point2f.(x_coord, y_coord)
            push!(vertices, verts)

            if show_median
                # interpolate median bounds between corresponding points
                xm = spec.median
                ip = findfirst(>(xm), spec.kde.x)
                ym₋, ym₊ = spec.kde.density[ip-1], spec.kde.density[ip]
                xm₋, xm₊ = spec.kde.x[ip-1], spec.kde.x[ip]
                ym = (xm * (ym₊ - ym₋) + xm₊ * ym₋ - xm₋ * ym₊) / (xm₊ - xm₋)
                median_left = Point2f(spec.side == 1 ? spec.x : spec.x - ym * scale, xm)
                median_right = Point2f(spec.side == -1 ? spec.x : spec.x + ym * scale, xm)
                push!(lines, median_left => median_right)
            end

            push!(colors, spec.color)
        end

        return (vertices = vertices, lines = lines, colors = colors)
    end

    poly!(
        plot,
        lift(s -> s.vertices, signals),
        color = lift(s -> s.colors, signals),
        strokecolor = plot[:strokecolor],
        strokewidth = plot[:strokewidth],
    )
    linesegments!(
        plot,
        lift(s -> s.lines, signals),
        color = plot[:mediancolor],
        linewidth = plot[:medianlinewidth],
        visible = plot[:show_median],
        inspectable = plot[:inspectable]
    )
end


text:
"""
GA0 Genetic Algorithm with 0 difficulty
Authors: Carmine Spagnuolo and Giuseppe D'Ambrosio

"""

using eqjl

import Pkg
Pkg.add("JSON")
import JSON
Pkg.add("Evolutionary")
import Evolutionary
Pkg.add("Random")
import Random
Pkg.add("Distributions")
import Distributions

function map(pops)
    eqjl.OUT_put(join(pops, ","))
    result = eqjl.IN_get()
    parse(Float64, result)
end

function custom_mutate(sigma::Float64)
    function mutation(recombinant::T) where {T<:Array}
        d = Distributions.Normal(0, sigma)
        println("Mutate: ", recombinant)
        recombinant .+= rand(d)
        println("to: ", recombinant)
        return recombinant
    end
    return mutation
end

#=
:param num_iter: number of generations
:param num_pop: size of population
:param seed: random seed
:param csv_file_name: csv file name (e.g., "params_for_deap.csv")
=#

function main()
    # Reading settings from file
    eqjl.OUT_put("Settings")
    settings_filename = eqjl.IN_get()
    settings = JSON.parsefile(settings_filename)

    Random.seed!(settings["seed"])
    sigma = settings["sigma"]
    num_iter = settings["num_iter"]
    num_pop = settings["num_pop"]
    mate_pb = settings["mate_pb"]
    mutate_pb = settings["mutate_pb"]
    crossover_pb = settings["crossover_pb"]

    # Creating population
    individuals = ()->Float64[rand(0.0:0.01:5.0), rand(0.0:0.01:5.0)]

    # Setting the GA algorithm
    algorithm = Evolutionary.GA(
        selection = Evolutionary.tournament(floor(Int, num_pop/2)),
        mutation =  custom_mutate(sigma),
        crossover = Evolutionary.uniformbin(mate_pb),
        mutationRate = mutate_pb,
        crossoverRate = 0.5,
        populationSize = num_pop,
        ε = 0.02
    )

    # Set "show_trace=false, store_trace=false" true to got more details
    opts = Evolutionary.Options(iterations=num_iter, show_every=1, show_trace=false, store_trace=false)

    res = Evolutionary.optimize(map, individuals, algorithm, opts)

    eqjl.OUT_put("FINAL")
    println("###FINAL sended")
    eqjl.OUT_put(res)
    # println(res)
end
@async main()

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
using QuantumPEPS
using Yao
using Flux

nx=4
ny=4
depth=1
nv=5
nbatch=1024
maxiter=200
J2=0.5
lr=0.1
periodic=false

model = J1J2(nx, ny; J2=J2, periodic=periodic)
config = QMPSConfig(nv, nx*ny-nv+1, depth)
optimizer = Flux.Optimise.ADAM(lr)
qpeps, history = train(config, model; maxiter=maxiter, nbatch=nbatch, optimizer=optimizer, use_cuda=false)