text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    This function computes the sunlight, penumbra, and umbra times of the
#    satellite for one orbit every day in a year.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Longo, C. R. O., Rickman, S. L (1995). Method for the Calculation of
#       Spacecraft Umbra and Penumbra Shadow Terminator Points. NASA Technical
#       Paper 3547.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export eclipse_time_summary

"""
    function eclipse_time_summary(JD₀::Number, a::Number, e::Number, i::Number, RAAN::Number, w::Number, Δd::Integer, relative::Bool = false, Δt₀::AbstractFloat = -1.0)

Compute the eclipse time of an orbit with semi-major axis `a` [m], eccentricity
`e`, inclination `i` [rad], initial right ascension of the ascending node `RAAN`
[rad], and initial argument of perigee `w` [rad]. The orbit epoch, which is also
the day in which the analysis will begin, is `JD₀` [Julian Day]. The analysis
will be performed for each day during `Δd` days.

This function will compute the eclipse time of one orbit per day.

If the argument `relative` is `true`, then the computed times will be relative
to the nodal period [%]. Otherwise, they will be computed in seconds. By
default, `relative = false`.

The argument `Δt₀` can be used to select the time step in which the orbit will
be propagated. Notice that this algorithm performs a numerical search to find
the beginning of each section (sunlight, penumbra, and umbra) with millisecond
precision. Thus, selecting a high number for `Δt₀` will make the analysis
faster, but the accuracy is lost if a region time span is smalled than `Δt₀`. If
this parameter is omitted or if it is negative, then the time step will be
selected automatically to match a mean anomaly step of 5°.

All the analysis is performed using a J2 orbit propagator.

# Returns

The following table:

        day | Sunlight Time | Penumbra Time | Umbra Time
       -----+---------------+---------------+------------

"""
function eclipse_time_summary(JD₀::Number, a::Number, e::Number, i::Number,
                              RAAN::Number, w::Number, Δd::Integer,
                              relative::Bool = false, Δt₀::AbstractFloat = -1.0)

    # Get the period of the orbit.
    T = period(a, e, i, :J2)

    # If the user did not specify the integration step, then compute based on
    # a mean anomaly step.
    if Δt₀ <= 0
        # Step of the orbit propagation (mean anomaly) [rad].
        Δm = 0.5*π/180

        # Time step related to the desired mean anomaly step [s].
        n   = angvel(a, e, i, :J2)
        Δt₀ = Δm/n
    end

    # Vector of the days in which the beta angle will be computed.
    days = collect(0:1:Δd-1)

    # Preallocate the output variables.
    p_time = zeros(Δd)  # Penumbra time.
    u_time = zeros(Δd)  # Umbra time.
    s_time = zeros(Δd)  # Sunlight time.

    # Configure the orbit propagator.
    orbp = init_orbit_propagator(Val{:J2}, JD₀, a, e, i, RAAN, w, 0)

    # Lambda functions
    # ==========================================================================

    # Function to compute the lightning condition given an instant of the day
    # `t`, the day from orbit epoch `d`, and the Sun vector `s_i`.
    f(t,d,s_i)::Int = begin
        ~, r_i, ~ = propagate!(orbp, 86400d + t)
        return satellite_lighting_condition(r_i, s_i)
    end

    # Return `true` if the lightning condition given an instant of the day `t`,
    # the day from orbit epoch `d`, and the Sun vector `s_i` is equal `state`.
    fb(t,d,s_i,state)::Bool = f(t,d,s_i) == state

    # Accumulate the time step `Δt` according to the state `state`.
    accum(Δt, state, ind, s, p, u)::Nothing = begin
        @inbounds if state == SAT_LIGHTING_SUNLIGHT
            s[ind] += Δt
        elseif state == SAT_LIGHTING_PENUMBRA
            p[ind] += Δt
        elseif state == SAT_LIGHTING_UMBRA
            u[ind] += Δt
        end

        return nothing
    end

    # Loop
    # ==========================================================================

    @inbounds for d in days
        # Get the sun position represented in the Inertial coordinate frame.
        s_i = sun_position_i(JD₀+d)

        # Initial state.
        state = f(0,d,s_i)

        # Compute the eclipse time during one orbit.
        Δt = Δt₀
        k  = Δt

        while true
            new_state = f(k,d,s_i)

            # Check if the state has changed.
            if new_state != state
                # Refine to find the edge.
                k₀ = k - Δt
                k₁ = k
                kc = find_crossing(fb, k₀, k₁, true, false, d, s_i, state)

                # Times to be added in the previous and current states.
                Δ₀ = kc - k₀
                Δ₁ = k₁ - kc

                accum(Δ₀, state,     d+1, s_time, p_time, u_time)
                accum(Δ₁, new_state, d+1, s_time, p_time, u_time)

            # If not, just add the time step to the current state.
            else
                accum(Δt, state, d+1, s_time, p_time, u_time)
            end

            state = new_state

            abs(k - T) < 1e-3 && break

            # Make sure that the last interval will have the exact size so that
            # the end of analysis is the end of the orbit.
            (k + Δt > T) && (Δt = T - k)

            k += Δt
        end
    end

    if (!relative)
        [days s_time p_time u_time]
    else
        total_time = s_time + p_time + u_time
        [days s_time./total_time p_time./total_time u_time./total_time]
    end
end


text:
module Calculus
    import Compat
    import Base.ctranspose
    export check_derivative,
           check_gradient,
           check_hessian,
           check_second_derivative,
           deparse,
           derivative,
           differentiate,
           gradient,
           hessian,
           integrate,
           jacobian,
           second_derivative

    # TODO: Debate type system more carefully
    # abstract BundledFunction
    # abstract ScalarFunction
    # abstract VectorFunction
    # abstract ForwardDifference
    # abstract CentralDifference
    # abstract ComplexDifference
    # abstract GradientEstimator
    # abstract HessianEstimator

    # typealias NonDifferentiableFunction Function
    # type DifferentiableFunction
    #   f::Function
    #   g::Function
    # end
    # type TwiceDifferentiableFunction
    #   f::Function
    #   g::Function
    #   h::Function
    # end
    # type NonDifferentiableBundledFunction <: BundledFunction
    #   f::Function
    #   fstorage::Any
    # end
    # type DifferentiableBundledFunction <: BundledFunction
    #   f::Function
    #   g::Function
    #   fstorage::Any
    #   gstorage::Any
    # end
    # type TwiceDifferentiableBundledFunction <: BundledFunction
    #   f::Function
    #   g::Function
    #   h::Function
    #   fstorage::Any
    #   gstorage::Any
    #   hstorage::Any
    # end

    include("finite_difference.jl")
    include("derivative.jl")
    include("check_derivative.jl")
    include("integrate.jl")
    include("symbolic.jl")
    include("differentiate.jl")
    include("deparse.jl")
end

ext: jl


text:
"""
Compute the upper bounds on `A_i ⋅ x under common linear constraints`
"""
function compute_bounds(bp::BilevelProblem; optimizer, silent=true, iterated=true)
    m = Model(optimizer)
    silent && MOI.set(m, MOI.Silent(), true)
    (_, x, v, λ, σ, s, _) = high_point_relaxation(m, bp)
    optimize!(m)
    st = termination_status(m)
    @assert st == MOI.OPTIMAL
    function bound_vector()
        bounds_val = Float64[]
        for i in 1:bp.ml
            @objective(m, Max, bp.A[i,:] ⋅ x)
            optimize!(m)
            st = termination_status(m)
            if st == MOI.INFEASIBLE || st == MOI.INFEASIBLE_OR_UNBOUNDED
                @warn "Status $st"
                return (st, [NaN])
            elseif st == MOI.DUAL_INFEASIBLE
                push!(bounds_val, Inf)
            else
                st == MOI.OPTIMAL || error("Status $st")
                push!(bounds_val, objective_value(m))
            end
        end
        return (MOI.OPTIMAL, bounds_val)
    end
    (st, Abounds) = bound_vector()
    if st != MOI.OPTIMAL
        return (st, Abounds)
    end
    niter = 1
    while iterated && all(isfinite, Abounds)
        @info niter
        Aprev = copy(Abounds)
        add_bound_inequalities(m, bp, x, v, λ, Aprev)
        (st, Abounds) = bound_vector()
        if st != MOI.OPTIMAL
            return (st, Abounds)
        end
        if all(tup -> tup[1] ≈ tup[2], zip(Abounds, Aprev))
            break
        end
        niter += 1
    end
    return (MOI.OPTIMAL, Abounds)
end

function add_bound_inequalities(m::JuMP.Model, bp::BilevelProblem, x, v, λ, bounds::AbstractVector{<:Real})
    @constraint(m, λ ⋅ bp.b + v ⋅ bp.d ≤ bounds ⋅ λ)
    return nothing
end

function add_bound_inequalities(m::JuMP.Model, bp::BilevelProblem, x, v, λ; optimizer, silent=true, verbose=false)
    (st, Aplus) = compute_bounds(bp; optimizer=optimizer, silent=silent)
    if st != MOI.OPTIMAL
        return false
    end
    if all(isfinite, Aplus)
        add_bound_inequalities(m, bp, x, v, λ, Aplus)
        return true
    end
    verbose && @warn("Not all bounds are finite: # inf: $(count(isfinite, Aplus))")
    return false
end


text:
#=
   IMPORTANT

   Julia v0.4.0 does not provide elementary functions that always round to nearest.
   For example:  exp( 0.05489354701373524 ) evaluates 1.0564281487582248
                                            actually  1.0564281487582246

   Consequently, the bounds provided below must be expanded to allow for this.

   This file should be used with functions that do not always evaluate to nearest.
=#

#=
   functions y=f(x) such that f(x+dx) > f(x) for dx>0 and f(x) defined for -Inf..x..Inf
=#

for (fn) in (:exp, :expm1, :log, :log1p,
             :sin, :cos, :tan, :csc, :sec, :cot,
             :asin, :acos, :atan, :acsc, :asec, :acot,
             :sinh, :cosh, :tanh, :csch, :sech, :coth,
             :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
             :erf, :erfinv)
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            lo = ($fn)(lo)
            hi = ($fn)(hi)
            if lo > hi
               lo,hi = hi,lo
            end
            lo = prevFloat(lo)
            hi = nextFloat(hi)

            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end

#= in arith
for (fn, domainMin) in ((:sqrt, 0.0), )
    @eval begin
        function ($fn){S<:Sculpt,Q<:Qualia,C<:Clay}(x::Flex{S,Q,C})
            if x.lo < ($domainMin)
                throw(ErrorException("DomainError: $($fn) expected x.lo>=($($domainMin)), got $(x)."))
            end
            loIsOpen, hiIsOpen = boundries(S)
            # round to nearest, make sure lo,hi are min(this,that), max(this,that) respectively
            lo, hi = value(x)
            lo = ($fn)(lo)
            hi = ($fn)(hi)
            if loIsOpen
                with_rounding(C, RoundDown) do
                    lo = min(lo, ($fn)(x.lo))
                end
            end
            if hiIsOpen
                with_rounding(C, RoundUp) do
                    hi = max(hi, ($fn)(x.hi))
                end
            end
            Flex{CLCL,Q,C}(lo,hi)
        end
    end
end
=#


text:
using RigidBodyDynamics

function forwarddiff_compatible(::typeof(mass_matrix!),
        statecache::StateCache, resultcache::DynamicsResultCache, normalize = Ref(false))
    function (out::AbstractVector, q::AbstractVector{T}) where T
        state = statecache[T]
        result = resultcache[T]
        set_configuration!(state, q)
        if normalize[]
            normalize_configuration!(state)
        end
        mass_matrix!(result, state)
        copyto!(out, result.massmatrix)
    end
end

function forwarddiff_compatible(::typeof(dynamics_bias!),
        statecache::StateCache, resultcache::DynamicsResultCache, v::AbstractVector, normalize = Ref(false))
    function (out::AbstractVector, q::AbstractVector{T}) where T
        state = statecache[T]
        result = resultcache[T]
        set_configuration!(state, q)
        set_velocity!(state, v)
        if normalize[]
            normalize_configuration!(state)
        end
        dynamics_bias!(result, state)
        copyto!(out, result.dynamicsbias)
    end
end

function forwarddiff_compatible(::typeof(inverse_dynamics!),
        statecache::StateCache, resultcache::DynamicsResultCache, v::AbstractVector, v̇::AbstractVector, normalize = Ref(false))
    function (out::AbstractVector, q::AbstractVector{T}) where T
        state = statecache[T]
        result = resultcache[T]
        set_configuration!(state, q)
        set_velocity!(state, v)
        if normalize[]
            normalize_configuration!(state)
        end
        inverse_dynamics!(result.dynamicsbias, result.jointwrenches, result.accelerations, state, v̇, result.totalwrenches)
        copyto!(out, result.dynamicsbias)
    end
end

function forwarddiff_compatible(::typeof(dynamics!),
    statecache::StateCache, resultcache::DynamicsResultCache, v::AbstractVector, τ::AbstractVector, normalize = Ref(false))
    function (out::AbstractVector, q::AbstractVector{T}) where T
        state = statecache[T]
        result = resultcache[T]
        set_configuration!(state, q)
        set_velocity!(state, v)
        if normalize[]
            normalize_configuration!(state)
        end
        dynamics!(result, state, τ)
        copyto!(out, result.v̇)
    end
end


text:
function createhighordermesh(mesh,app)

display("run createnodes...");
#[mesh.dgnodes,mesh.elemtype,perm] = createnodes(mesh.p,mesh.t,app.porder);
mesh.dgnodes, elemtype, perm = createnodes(mesh.p,mesh.t,app.porder);

display("run facenumbering...");
mesh.f, mesh.tprd = facenumbering(mesh.p,mesh.t,elemtype,mesh.boundaryexpr,mesh.periodicexpr);

# Project nodes on the curved boundaries
if length(mesh.curvedboundaryexpr)>0 && app.porder>1 && maximum(mesh.curvedboundary[:])>0
    display("Project dgnodes onto the curved boundaries...");
    fd = mesh.curvedboundaryexpr;
    nd = size(mesh.p,1);
    nfe,ne = size(mesh.f);
    if nd == 2
        for i = 1:ne
            for j = 1:nfe
                if mesh.f[j,i]!=0 # boundary element
                    k = abs(mesh.f[j,i]); # get boundary index
                    if mesh.curvedboundary[k]==1 # if this boundary is curved
                        p = mesh.dgnodes[perm[:,j],:,i];
                        deps = 1e-8*maximum(maximum(p,dim=1)-minimum(p,dim=1));
                        d = fd[k](p);
                        dgradx = (fd[k]([p[:,1]+deps p[:,2]])-d)/deps;
                        dgrady = (fd[k]([p[:,1] p[:,2]+deps])-d)/deps;
                        dgrad2 = dgradx.*dgradx + dgrady.*dgrady;
                        dgrad2[dgrad2==0] = 1;
                        p = p-[d.*dgradx./dgrad2 d.*dgrady./dgrad2];
                        mesh.dgnodes[perm[:,j],:,i] = p;
                    end
                end
            end
        end
    elseif nd==3
        for i = 1:ne
            for j = 1:nfe
                if mesh.f[j,i]!=0 # boundary element
                    k = abs(mesh.f[j,i]); # get boundary index
                    if mesh.curvedboundary[k]==1 # if this boundary is curved
                        p = mesh.dgnodes[perm[:,j],:,i];
                        deps = 1e-8*maximum(maximum(p,dim=1)-minimum(p,dim=1));
                        d = fd[k](p);
                        dgradx = (fd[k]([p[:,1]+deps,p[:,2],p[:,3]])-d)/deps;
                        dgrady = (fd[k]([p[:,1],p[:,2]+deps,p[:,3]])-d)/deps;
                        dgradz = (fd[k]([p[:,1],p[:,2],p[:,3]+deps])-d)/deps;
                        dgrad2 = dgradx.*dgradx+dgrady.*dgrady+dgradz.*dgradz;
                        dgrad2[dgrad2==0] = 1;
                        p = p-[d.*dgradx./dgrad2 d.*dgrady./dgrad2 d.*dgradz./dgrad2];
                        mesh.dgnodes[perm[:,j],:,i] = p;
                    end
                end
            end
        end
    end
end

return mesh

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
function unsafepointers()
    intspace = [42]
    address = pointer_from_objref(intspace)
    println("The address of intspace is $address")
    anotherint = unsafe_pointer_to_objref(address)
    println("intspace is $(intspace[1]), memory at $address, reference value $(anotherint[1])")
    intspace[1] = 123456
    println("Now, intspace is $(intspace[1]), memory at $address, reference value $(anotherint[1])")
    anotherint[1] = 7890
    println("Now, intspace is $(intspace[1]), memory at $(pointer_from_objref(anotherint)), reference value $(anotherint[1])")
end

unsafepointers()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# import Pkg; Pkg.add("Documenter")
using Documenter, NNHelferlein

makedocs(modules = [NNHelferlein],
         clean = false,
         assets = ["assets/favicon.ico"],
         sitename = "NNHelferlein.jl",
         authors = "Andreas Dominik",
         pages = [
                  "Introduction" => "index.md",
                  "Overview" => "overview.md",
                  "Examples" => "examples.md",
                  "API Reference" => "api.md",
                  "License" => "license.md"
                  ],
         #          # Use clean URLs, unless built as a "local" build
          html_prettyurls = false, #!("local" in ARGS),
          html_canonical = "https://andreasdominik.github.io/NNHelferlein.jl/stable/"
         )

deploydocs(
    # root   = "<current-directory>",
    target = "build",
    repo   = "github.com/andreasdominik/NNHelferlein.jl.git",
    branch = "gh-pages",
    # deps   = nothing | <Function>,
    # make   = nothing | <Function>,
    devbranch = "main",
    devurl = "dev",
    # versions = ["stable" => "v^", "v#.#", devurl => "dev"],
    push_preview    = false,
    # repo_previews   = repo,
    # branch_previews = branch
)


text:
using AtomicLevels
using UnicodeFun
using Test

@testset "Terms" begin
    @testset "Construction" begin
        @test T"1S" == Term(0, 0, p"even")
        @test T"1S" == Term(0, 0, 1)
        @test T"1Se" == Term(0, 0, p"even")
        @test T"1So" == Term(0, 0, p"odd")
        @test T"1So" == Term(0, 0, -1)
        @test T"2So" == Term(0, 1//2, p"odd")
        @test T"4P" == Term(1, 3//2, p"even")
        @test T"3D" == Term(2, 1, p"even")
        @test T"3Do" == Term(2, 1, p"odd")
        @test T"1[54]" == Term(54, 0, p"even")
        @test T"1[3/2]" == Term(3//2, 0, p"even")
        @test T"2[3/2]o" == Term(3//2, 1//2, p"odd")

        @test_throws ArgumentError parse(Term, "1[4/3]")
        @test_throws ArgumentError parse(Term, "1[43/]")
        @test_throws ArgumentError parse(Term, "1[/43/]")
        @test_throws ArgumentError parse(Term, "1[/43]")
        @test_throws ArgumentError parse(Term, "P")
        @test_throws ArgumentError parse(Term, "asdf")
    end

    @testset "Properties" begin
        @test multiplicity(T"1S") == 1
        @test multiplicity(T"1So") == 1
        @test multiplicity(T"2S") == 2
        @test multiplicity(T"2So") == 2
        @test multiplicity(T"3S") == 3
        @test multiplicity(T"3So") == 3

        @test weight(T"1S") == 1
        @test weight(T"1P") == 3
        @test weight(T"2S") == 2
        @test weight(T"2P") == 6

        @test T"1S" > T"1So"
        @test T"1So" < T"1S"
        @test T"1S" < T"2S"
        @test T"1S" < T"1P"
        @test T"1S" < T"3S"
        @test T"1P" < T"3S"
    end

    @testset "Pretty printing" begin
        map([T"1S" => "¹S",
             T"2So" => "²Sᵒ",
             T"4[3/2]" => "⁴[3/2]"]) do (t,s)
            @test "$(t)" == s
        end
    end

    @testset "Orbital terms" begin
        function test_single_orbital_terms(orb::Orbital, occ::Int, ts::Vector{Term})
            cts = sort(terms(orb, occ))
            ts = sort(ts)
            ccts = copy(cts)
            for t in cts
                if t in ts
                    @test count(e -> e == t, ccts) == count(e -> e == t, ts)
                    ccts = filter!(e -> e != t, ccts)
                    ts = filter!(e -> e != t, ts)
                end
            end
            if length(ts) != 0
                println("fail, terms: ", join(string.(cts), ", "))
                println("missing: ", join(string.(ts), ", "))
                println("should not be there: ", join(string.(ccts), ", "))
                println("==========================================================")
            end
            @test length(ts) == 0
        end

        function test_single_orbital_terms(orb::Orbital, occs::Tuple{Int,Int}, ts::Vector{Term})
            map(occs) do occ
                test_single_orbital_terms(orb, occ, ts)
            end
        end

        function get_orbital(o::AbstractString)
            n = something(findfirst(isequal(o[1]), AtomicLevels.spectroscopic), 0)
            ℓ = n - 1
            orb = Orbital(n, ℓ)
            g = 2(2ℓ + 1)
            occ = length(o)>1 ? parse(Int, o[2:end]) : 1
            if g != occ && g/2 != occ
                occ = (occ, g-occ)
            end
            orb, occ
        end

        function test_orbital_terms(o_ts::Pair{<:ST,<:ST}) where {ST<:AbstractString}
            orb,occ = get_orbital(o_ts[1])
            m = match(r"([0-9]+)([A-Z])", o_ts[2])
            L = something(findfirst(isequal(lowercase(m[2])[1]), AtomicLevels.spectroscopic), 0)-1
            S = (parse(Int, m[1])-1)//2
            test_single_orbital_terms(orb, occ, [Term(L, S, parity(orb)^occ[1])])
        end

        function test_orbital_terms(o_ts::Pair{<:ST,<:Vector{ST}}) where {ST<:AbstractString}
            orb,occ = get_orbital(o_ts[1])

            p = parity(orb)^occ[1]
            ts = o_ts[2]

            p1 = r"([0-9]+)\(((?:[A-Z][0-9]*)+)\)"
            p2 = r"([0-9]+)([A-Z])"
            toL = s -> something(findfirst(isequal(lowercase(s)[1]), AtomicLevels.spectroscopic), 0) - 1
            ts = map(ts) do t
                m = match(p1, t)
                if m != nothing
                    S = (parse(Int, m[1]) - 1)//2
                    map(eachmatch(r"([A-Z])([0-9]*)", m[2])) do mmm
                        mm = mmm.match
                        [Term(toL(mm[1]),S,p)
                         for j in 1:(length(mm)>1 ? parse(Int, mm[2:end]) : 1)]
                    end
                else
                    m = match(p2, t)
                    S = (parse(Int, m[1]) - 1)//2
                    Term(toL(m[2]), S, p)
                end
            end
            test_single_orbital_terms(orb, occ, vcat(vcat(ts...)...))
        end

        # Table taken from Cowan, p. 110
        # Numbers following term symbols indicate the amount of times
        # different terms with the same (L,S) occur.
        test_orbital_terms("s" => "2S")
        for o in ["s2", "p6", "d10", "f14"]
            test_orbital_terms(o => "1S")
        end
        test_orbital_terms("p" => "2P")
        test_orbital_terms("p2" => ["1(SD)", "3P"])
        test_orbital_terms("p3" => ["2(PD)", "4S"])
        test_orbital_terms("d" => "2D")
        test_orbital_terms("d2" => ["1(SDG)", "3(PF)"])
        test_orbital_terms("d3" => ["2(PD2FGH)", "4(PF)"])
        test_orbital_terms("d4" => ["1(S2D2FG2I)", "3(P2DF2GH)", "5D"])
        test_orbital_terms("d5" => ["2(SPD3F2G2HI)", "4(PDFG)", "6S"])
        test_orbital_terms("f" => "2F")
        test_orbital_terms("f2" => ["1(SDGI)", "3(PFH)"])
        test_orbital_terms("f3" => ["2(PD2F2G2H2IKL)", "4(SDFGI)"])
        test_orbital_terms("f4" => ["1(S2D4FG4H2I3KL2N)", "3(P3D2F4G3H4I2K2LM)", "5(SDFGI)"])
        test_orbital_terms("f5" => ["2(P4D5F7G6H7I5K5L3M2NO)", "4(SP2D3F4G4H3I3K2LM)", "6(PFH)"])
        test_orbital_terms("f6" => ["1(S4PD6F4G8H4I7K3L4M2N2Q)", "3(P6D5F9G7H9I6K6L3M3NO)", "5(SPD3F2G3H2I2KL)", "7F"])
        test_orbital_terms("f7" => ["2(S2P5D7F10G10H9I9K7L5M4N2OQ)", "4(S2P2D6F5G7H5I5K3L3MN)", "6(PDFGHI)", "8S"])
        # # Data below are from Xu2006
        # test_orbital_terms("g9" => [2(S8 P19 D35 F40 G52 H54 I56 K53 L53 M44 N40 O32 Q26 R19 T15 U9 V7 W4 X2 YZ) 4(S6 P16 D24 F34 G38 H40 I42 K39 L35 M32 N26 O20 Q16 R11 T7 U5 V3 WX) 6(S3P3D9F8G12H10I12K9L9M6N6O3Q3RT) 8(PDFGHIKL) 10(S)])
        # test_orbital_terms("h11" => [2(S36 P107 D173 F233 G283 H325 I353 K370 L376 M371 N357 O335 Q307 R275 T241 U207 V173 W142 X114 Y88 Z68 2150 2236 2325 2417 2511 267 274 282 29 30) 4 (S37 P89 D157 F199 G253 H277 I309 K313 L323 M308 N300 O271 Q251 R216 T190 U155 V131 W101 X81 Y59 Z45 2130 2222 2313 249 255 263 27 28) 6 (S12 P35 D55 F76 G90 H101 I109 K111 L109 M105 N97 O87 Q77 R65 T53 U43 V33 W24 X18 Y12 Z8 215 223 2324) 8(S4 P4 D12 F11 G17 H15 I19 K16 L18 M14 N14 O10 Q10 R6 T6 U3 V3WX) 10(PDFGHIKLMN) 12S])

        @test_throws ArgumentError terms(o"2p", 7)
    end

    @testset "Count terms" begin
        @test count_terms(o"1s", 1, T"2S") == 1
        @test count_terms(o"1s", 2, T"1S") == 1
        @test count_terms(o"3d", 1, T"2D") == 1
        @test count_terms(o"3d", 3, T"2D") == 2
        @test count_terms(o"3d", 5, T"2D") == 3
    end

    @testset "Intermediate terms" begin
        @test_throws ArgumentError IntermediateTerm(T"2S", 2)

        @test string(IntermediateTerm(T"2D", 3)) == "₃²D"

        @testset "Seniority" begin
            # Table taken from p. 25 of
            #
            # - Froese Fischer, C., Brage, T., & Jönsson, P. (1997). Computational
            #   Atomic Structure : An MCHF Approach. Bristol, UK Philadelphia, Penn:
            #   Institute of Physics Publ.
            subshell_terms = [
                o"1s" => (1 => "2S1",
                         2 => "1S0"),
                o"2p" => (1 => "2P1",
                         2 => ("1S0", "1D2", "3P2"),
                         3 => ("2P1", "2D3", "4S3")),
                o"3d" => (1 => "2D1",
                         2 => ("1S0", "1D2", "1G2", "3P2", "3F2"),
                         3 => ("2D1", "2P3", "2D3", "2F3", "2G3", "2H3", "4P3", "4F3"),
                         4 => ("1S0", "1D2", "1G2", "3P2", "3F2", "1S4", "1D4", "1F4", "1G4", "1I4", "3P4", "3D4", "3F4", "3G4", "3H4", "5D4"),
                         5 => ("2D1", "2P3", "2D3", "2F3", "2G3", "2H3", "4P3", "4F3", "2S5", "2D5", "2F5", "2G5", "2I5", "4D5", "4G5", "6S5")),
                o"4f" => (1 => "2F1",
                         2 => ("1S0", "1D2", "1G2", "1I2", "3P2", "3F2", "3H2"))
            ]

            foreach(subshell_terms) do (orb, data)
                foreach(data) do (occ, expected_its)
                    p = parity(orb)^occ
                    expected_its = map(expected_its isa String ? (expected_its,) : expected_its) do ei
                        t_str = "$(ei[1:2])$(isodd(p) ? "o" : "")"
                        IntermediateTerm(parse(Term, t_str), parse(Int, ei[end]))
                    end
                    its = intermediate_terms(orb, occ)
                    @test its == [expected_its...]
                end
            end
        end
    end
end


text:
"""
    Frechet(α,θ)

The *Fréchet distribution* with shape `α` and scale `θ` has probability density function

```math
f(x; \\alpha, \\theta) = \\frac{\\alpha}{\\theta} \\left( \\frac{x}{\\theta} \\right)^{-\\alpha-1}
e^{-(x/\\theta)^{-\\alpha}}, \\quad x > 0
```

```julia
Frechet()        # Fréchet distribution with unit shape and unit scale, i.e. Frechet(1, 1)
Frechet(a)       # Fréchet distribution with shape a and unit scale, i.e. Frechet(a, 1)
Frechet(a, b)    # Fréchet distribution with shape a and scale b

params(d)        # Get the parameters, i.e. (a, b)
shape(d)         # Get the shape parameter, i.e. a
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Fréchet_distribution on Wikipedia](http://en.wikipedia.org/wiki/Fréchet_distribution)

"""
struct Frechet{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    θ::T
    Frechet{T}(α::T, θ::T) where {T<:Real} = new{T}(α, θ)
end

function Frechet(α::T, θ::T; check_args=true) where {T <: Real}
    check_args && @check_args(Frechet, α > zero(α) && θ > zero(θ))
    return Frechet{T}(α, θ)
end

Frechet(α::Real, θ::Real) = Frechet(promote(α, θ)...)
Frechet(α::Integer, θ::Integer) = Frechet(float(α), float(θ))
Frechet(α::T) where {T <: Real} = Frechet(α, one(T))
Frechet() = Frechet(1.0, 1.0, check_args=false)

@distr_support Frechet 0.0 Inf

#### Conversions
function convert(::Type{Frechet{T}}, α::S, θ::S) where {T <: Real, S <: Real}
    Frechet(T(α), T(θ))
end
function convert(::Type{Frechet{T}}, d::Frechet{S}) where {T <: Real, S <: Real}
    Frechet(T(d.α), T(d.θ), check_args=false)
end

#### Parameters

shape(d::Frechet) = d.α
scale(d::Frechet) = d.θ
params(d::Frechet) = (d.α, d.θ)
partype(::Frechet{T}) where {T} = T

#### Statistics

function mean(d::Frechet{T}) where {T}
    α = d.α
    return α > 1 ? d.θ * gamma(1 - 1 / α) : T(Inf)
end

median(d::Frechet) = d.θ * logtwo^(-1 / d.α)

mode(d::Frechet) = (iα = -1/d.α; d.θ * (1 - iα) ^ iα)

function var(d::Frechet{T}) where {T<:Real}
    if d.α > 2
        iα = 1 / d.α
        return d.θ^2 * (gamma(1 - 2 * iα) - gamma(1 - iα)^2)
    else
        return T(Inf)
    end
end

function skewness(d::Frechet{T}) where T<:Real
    if d.α > 3
        iα = 1 / d.α
        g1 = gamma(1 - iα)
        g2 = gamma(1 - 2 * iα)
        g3 = gamma(1 - 3 * iα)
        return (g3 - 3g2 * g1 + 2 * g1^3) / ((g2 - g1^2)^1.5)
    else
        return T(Inf)
    end
end

function kurtosis(d::Frechet{T}) where T<:Real
    if d.α > 3
        iα = 1 / d.α
        g1 = gamma(1 - iα)
        g2 = gamma(1 - 2iα)
        g3 = gamma(1 - 3iα)
        g4 = gamma(1 - 4iα)
        return (g4 - 4g3 * g1 + 3 * g2^2) / ((g2 - g1^2)^2) - 6
    else
        return T(Inf)
    end
end

function entropy(d::Frechet)
    1 + MathConstants.γ / d.α + MathConstants.γ + log(d.θ / d.α)
end

#### Evaluation

function logpdf(d::Frechet{T}, x::Real) where T<:Real
    (α, θ) = params(d)
    if x > 0
        z = θ / x
        return log(α / θ) + (1 + α) * log(z) - z^α
    else
        return -T(Inf)
    end
end

cdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? exp(-((d.θ / x) ^ d.α)) : zero(T)
ccdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? -expm1(-((d.θ / x) ^ d.α)) : one(T)
logcdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? -(d.θ / x) ^ d.α : -T(Inf)
logccdf(d::Frechet{T}, x::Real) where {T<:Real} = x > 0 ? log1mexp(-((d.θ / x) ^ d.α)) : zero(T)

quantile(d::Frechet, p::Real) = d.θ * (-log(p)) ^ (-1 / d.α)
cquantile(d::Frechet, p::Real) = d.θ * (-log1p(-p)) ^ (-1 / d.α)
invlogcdf(d::Frechet, lp::Real) = d.θ * (-lp)^(-1 / d.α)
invlogccdf(d::Frechet, lp::Real) = d.θ * (-log1mexp(lp))^(-1 / d.α)

function gradlogpdf(d::Frechet{T}, x::Real) where T<:Real
    (α, θ) = params(d)
    insupport(Frechet, x) ? -(α + 1) / x + α * (θ^α) * x^(-α-1) : zero(T)
end
function heslogpdf(d::Frechet{T}, x::Real) where T<:Real
    (α, θ) = params(d)
    insupport(Frechet, x) ? (α + 1) / x^2 - (α + α^2 ) * θ^α * x^(-α-2) : zero(T)
end

## Sampling

rand(rng::AbstractRNG, d::Frechet) = d.θ * randexp(rng) ^ (-1 / d.α)


text:
using StaticArrays
using ViewerGL
GL = ViewerGL

const rgb  = [77 202 137;
       130 70 88;
       198 180 71;
       133 185 98;
       175 115 69;
       227 184 128;
       144 121 86;
       203 107 113;
       231 157 134;
       146 230 136;
       79 153 125;
       187 154 38]

const RGBMaya = [[rgb[k,:]' 255]./ 255 for k=1:12]
const MayaColors = GL.Point4d.(RGBMaya)

"""
Palette of twelve predefined colors
"""
const WHITE   = Point4d([1.0, 1.0, 1.0, 1.0])
const RED     = Point4d([1.0, 0.0, 0.0, 1.0])
const GREEN   = Point4d([0.0, 1.0, 0.0, 1.0])
const BLUE    = Point4d([0.0, 0.0, 1.0, 1.0])
const CYAN    = Point4d([0.0, 1.0, 1.0, 1.0])
const MAGENTA = Point4d([1.0, 0.0, 1.0, 1.0])
const YELLOW  = Point4d([1.0, 1.0, 0.0, 1.0])
const ORANGE  = Point4d([1.0, 0.65, 1.0, 1.0])
const PURPLE  = Point4d([0.5, 0.0, 0.5, 1.0])
const BROWN   = Point4d([0.65, 0.16, 0.16, 1.0])
const GRAY    = Point4d([0.5, 0.5, 0.5, 1.0])
const BLACK   = Point4d([0.0, 0.0, 0.0, 1.0])

const palette = Dict(1=>WHITE, 2=>RED, 3=>GREEN, 4=>BLUE,
        5=>CYAN, 6=>MAGENTA, 7=>YELLOW, 8=>ORANGE,
        9=>PURPLE, 10=>BROWN, 11=>GRAY, 12=>BLACK )

#COLORS = palette
const COLORS = MayaColors

"""
	mutable struct GLColor

"""
mutable struct GLColor

    choice::Int   # to select between per_body, per_cell, and per_vertex

    color_per_body::Vector{Point4d}
    color_per_cell::Vector{Point4d}
    color_per_vertex::Vector{Point4d}

    colors::GL.GLVertexBuffer

    # constructor
    GLColor() = GLColor(1,COLORS[1],Point4d[],Point4d[],Point4d[],Point4d[])

    # constructor
    GLColor(choice) = GLColor(choice,COLORS[1],Point4d[],Point4d[],Point4d[],Point4d[])

    # constructor
    function GLColor(choice,color_per_body,color_per_cell,color_per_vertex,colors)
        ret = new(choice,color_per_body,color_per_cell,color_per_vertex,colors)
        return ret
    end

end


text:
module ShammaModel
using DataFrames
using Unitful
using Unitful: ms, s, Hz, kHz
using ProgressMeter
using Requires
using PlotAxes

export ms, s, Hz, kHz

next!(x::Progress) = ProgressMeter.next!(x)
next!(x::Nothing) = nothing

struct Filter
  A::Vector{Float64}
  B::Vector{Float64}
end
struct CochFilters
  filters::Vector{Filter}
  norm::Float64
end
const cochlear = Ref{CochFilters}()

include("AxisMeta.jl")
include("audiospect.jl")
include("cortical.jl")
include("plot_axes.jl")

const localunits = Unitful.basefactors
const localpromotion = Unitful.promotion
function __init__()
  cochlear[] = jldopen(joinpath(@__DIR__,"..","data","cochba.jld2"),"r") do file
    filters = map(keys(file["filters"])) do ch
      A = file["filters/$ch/A"]
      B = file["filters/$ch/B"]
      Filter(A,B)
    end
    CochFilters(filters,file["norm"])
  end

  merge!(Unitful.basefactors, localunits)
  merge!(Unitful.promotion, localpromotion)
  Unitful.register(ShammaModel)
end

end

max_stars_repo_name: haberdashPI/ShammaModel

max_issues_repo_name: haberdashPI/ShammaModel

max_forks_repo_name: haberdashPI/ShammaModel


text:
# Load data to use for tests
path = dirname(@__FILE__)
fred = CSV.read("$path/../reference/fred_160812.csv", DataFrame)

# Specify vintage and dates
custom_settings = Dict{Symbol, Setting}(
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :n_anticipated_shocks     => Setting(:n_anticipated_shocks, 6),
    :population_forecast      => Setting(:population_forecast, false))
m = AnSchorfheide(testing = true, custom_settings = custom_settings)
m <= Setting(:saveroot, "$path/../reference/")
m <= Setting(:rate_expectations_source, :ois)

pop_forecast = DataFrame(CNP16OV = cumsum(fred[!,:CNP16OV][end] .* (ones(8) .* 1.03)))
pop_forecast[!,:date] = DSGE.get_quarter_ends(DSGE.next_quarter(fred[!,:date][end]), fred[!,:date][end] + Year(2))

mb_means = load("$path/../reference/means_bands_out.jld2")["exp_full_means"] # this is AnSchorfheide
m_dates = DSGE.quarter_range(DSGE.quartertodate("1960-Q1"), DSGE.quartertodate("2015-Q3"))
start_date = m_dates[1]
histobs_df = DataFrame(hcat(m_dates, mb_means[:histobs]'),
                       vcat(:date, collect(keys(m.observables))))
histpseudo_df = DataFrame(hcat(m_dates, mb_means[:histpseudo]'),
                          vcat(:date, collect(keys(m.pseudo_observables))))

revobs1 = Matrix{Float64}(reverse_transform(m, mb_means[:histobs], start_date,
                                            collect(keys(m.observables)), :obs, verbose = :none)[!,2:end])
revobs2 = Matrix{Float64}(reverse_transform(m, histobs_df, :obs, verbose = :none)[:,2:end])
revpseudo1 = Matrix{Float64}(reverse_transform(m, mb_means[:histpseudo], start_date,
                                               collect(keys(m.pseudo_observables)), :pseudo,
                                               verbose = :none)[:,2:end])
revpseudo2 = Matrix{Float64}(reverse_transform(m, histpseudo_df, :pseudo, verbose = :none)[!,2:end])
q4revobs1 = Matrix{Float64}(reverse_transform(m, mb_means[:histobs], start_date,
                                            collect(keys(m.observables)), :obs, verbose = :none,
                                              fourquarter = true)[4:end,2:end])
q4revobs2 = Matrix{Float64}(reverse_transform(m, histobs_df, :obs, verbose = :none,
                                              fourquarter = true)[4:end,2:end])
q4revpseudo1 = Matrix{Float64}(reverse_transform(m, mb_means[:histpseudo], start_date,
                                                 collect(keys(m.pseudo_observables)), :pseudo,
                                                 verbose = :none,
                                                 fourquarter = true)[4:end,2:end])
q4revpseudo2 = Matrix{Float64}(reverse_transform(m, histpseudo_df, :pseudo, verbose = :none,
                                                 fourquarter = true)[4:end,2:end])

exp_revobs, exp_revpseudo, exp_q4revobs, exp_q4revpseudo =
    JLD2.jldopen("$path/../reference/reverse_transform_out.jld2", "r") do file
    read(file, "exp_revobs"), read(file, "exp_revpseudo"),
    read(file, "exp_q4revobs"), read(file, "exp_q4revpseudo")
end

@testset "Check reverse transform returns expected output" begin
    @test @test_matrix_approx_eq revobs1 exp_revobs
    @test @test_matrix_approx_eq revobs2 exp_revobs
    @test @test_matrix_approx_eq revpseudo1 exp_revpseudo
    @test @test_matrix_approx_eq revpseudo2 exp_revpseudo
    @test @test_matrix_approx_eq q4revobs1 exp_q4revobs
    @test @test_matrix_approx_eq q4revobs2 exp_q4revobs
    @test @test_matrix_approx_eq q4revpseudo1 exp_q4revpseudo
    @test @test_matrix_approx_eq q4revpseudo2 exp_q4revpseudo
end

nothing


text:
#----------------------------------------------------------------------------------------------#
#                               MOO.jl -- Moodle engine backend                                #
#----------------------------------------------------------------------------------------------#

#----------------------------------------------------------------------------------------------#
#                                       Engine Metadata                                        #
#----------------------------------------------------------------------------------------------#

function meta(::Type{MOO})
    (
        name = "Moodle",
        symb = "MOO",
    )
end


text:
using Walker2014Sampling
using Random
using Test

##### Example 1

k = 3;
x = 5;
p(x) = 1.0
s = Walker2014Sampler(k);
rng = MersenneTwister(0);
rand(rng, p, s, x)

@test sum(s.r) ≈ 1.0
@test s.q[1] ≈ (s.p[1] + s.p[2] + s.p[3])^(-1)
@test s.q[2] ≈ (s.p[2] + s.p[3] + s.p[4])^(-1)
@test s.q[3] ≈ (s.p[3] + s.p[4] + s.p[5])^(-1)

##### Example 2

k = 3;
x = 2;
p(x) = 1.0
s = Walker2014Sampler(k);
rng = MersenneTwister(0);
rand(rng, p, s, x)

@test sum(s.r) ≈ 1.0
@test s.q[1] ≈ (s.p[1] + s.p[2] + s.p[3])^(-1)
@test s.q[2] ≈ (s.p[2] + s.p[3] + s.p[4])^(-1)
@test s.q[3] ≈ (s.p[3] + s.p[4] + s.p[5])^(-1)

##### Example 3

k = 3;
x = 5;
p(x) = Float64(x >= 2)
s = Walker2014Sampler(k);
rng = MersenneTwister(0);
rand(rng, p, s, x)

@test sum(s.r) ≈ 1.0
@test s.q[1] ≈ (s.p[1] + s.p[2] + s.p[3])^(-1)
@test s.q[2] ≈ (s.p[2] + s.p[3] + s.p[4])^(-1)
@test s.q[3] ≈ (s.p[3] + s.p[4] + s.p[5])^(-1)


text:
@testset "Testing files with unusually large aspects" begin

     @testset "sample15 - More than 127 points" begin
         @test_nowarn readc3d(joinpath(datadir, "sample15", "FP1.C3D"))
         @test_nowarn readc3d(joinpath(datadir, "sample15", "FP2.C3D"))
     end

     @testset "sample17 - More than 128 channels" begin
         @test_nowarn readc3d(joinpath(datadir, "sample17",  "128analogchannels.c3d"))
     end

     @testset "sample19 - 34672 frames analog data" begin
         @test_nowarn readc3d(joinpath(datadir, "sample19",  "sample19.c3d"))
     end

     if isdir(joinpath(datadir, "sample31"))
         @testset "sample31 - more than 65535 frames" begin
             @test_nowarn readc3d(joinpath(datadir, "sample31", "large01.c3d"))
             @test_nowarn readc3d(joinpath(datadir, "sample31", "large02.c3d"))
         end
     end

     @testset "sample33 - DATA_START greater than 127" begin
         @test_nowarn readc3d(joinpath(datadir, "sample33",  "bigparlove.c3d"))
     end

end


text:
# step one, generate optimizer info
function rehearse_contraction(tn::CuTensorNetwork, max_workspace_size::Integer, optimizer_conf::OptimizerConfig=OptimizerConfig())
    info = CuTensorNetworkContractionOptimizerInfo(tn.desc)
    # set optimizer options
    config = CuTensorNetworkContractionOptimizerConfig(optimizer_conf)
    # perform optimization
    cutensornetContractionOptimize(handle(), tn.desc, config, max_workspace_size, info)
    plan = CuTensorNetworkContractionPlan(tn.desc, info, max_workspace_size)
    return info, plan
end

function num_slices(info::CuTensorNetworkContractionOptimizerInfo)
    n = Ref{Int64}()
    cutensornetContractionOptimizerInfoGetAttribute(handle(), info, CUTENSORNET_CONTRACTION_OPTIMIZER_INFO_NUM_SLICES, n, sizeof(Int64))
    return n[]
end

function num_sliced_modes(info::CuTensorNetworkContractionOptimizerInfo)
    n = Ref{Int32}()
    cutensornetContractionOptimizerInfoGetAttribute(handle(), info, CUTENSORNET_CONTRACTION_OPTIMIZER_INFO_NUM_SLICED_MODES, n, sizeof(Int32))
    return n[]
end

function flop_count(info::CuTensorNetworkContractionOptimizerInfo)
    n = Ref{Float64}()
    cutensornetContractionOptimizerInfoGetAttribute(handle(), info, CUTENSORNET_CONTRACTION_OPTIMIZER_INFO_FLOP_COUNT, n, sizeof(Float64))
    return n[]
end

# step 2, contract
function perform_contraction!(tn::CuTensorNetwork, info, plan, workspace_size::Integer, ::NoAutoTune; prefs::AutotunePreferences=AutotunePreferences(), stream::CuStream=stream())
    with_workspace(workspace_size) do workspace
        for slice_ix in 0:num_slices(info)-1
            cutensornetContraction(handle(), plan, pointer.(tn.input_arrs),
                                   pointer(tn.output_arr), workspace, workspace_size,
                                   slice_ix, stream)
        end
    end
    return tn
end

function perform_contraction!(tn::CuTensorNetwork, info, plan, workspace_size::Integer, ::AutoTune; prefs::AutotunePreferences=AutotunePreferences(), stream::CuStream=stream())
    ctn_prefs = CuTensorNetworkAutotunePreference(prefs)
    with_workspace(workspace_size) do workspace
        cutensornetContractionAutotune(handle(), plan, pointer.(tn.input_arrs),
                                       pointer(tn.output_arr), workspace, workspace_size,
                                       ctn_prefs, stream)
        for slice_ix in 0:num_slices(info)-1
            cutensornetContraction(handle(), plan, pointer.(tn.input_arrs),
                                   pointer(tn.output_arr), workspace, workspace_size,
                                   slice_ix, stream)
        end
    end
    return tn
end


text:
"""
The first attempt Nathan and I made at IO
"""
using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using DelimitedFiles

"""
Plot the performance statistics of greedy_multi_path for 1 random end-user pair
over a range of graph sizes and for different numbers of max_path (1,2,3) (ie. the
maximum number of paths that can be purified by a given end-user pair.)
"""
file_path = "data/maxpaths.txt"
generate_new_data = true
if generate_new_data

    # Params
    num_trials = 50::Int64
    min_size = 2::Int64
    max_size = 20::Int64

    @assert min_size < max_size

    perf_data1 = []
    perf_data2 = []
    perf_data3 = []
    perf_data4 = []

    size_list = collect(min_size:1:max_size)
    for (j, data_array) in enumerate([perf_data1, perf_data2, perf_data3, perf_data4])
        println("Collecting for max_paths: $j")
        for i in size_list
            println("Collecting for gridsize: $i")
            # Generate ixi graph:
            net = GridNetwork(i, i)

            # Collect performance statistics
            performance, dummy, dummy, dummy = net_performance(net, num_trials, 1, max_paths=j)
            push!(data_array, performance)
        end
    end

    # Get values for x axis
    x = collect(min_size:1:max_size)

    # Extract data from performance data
    loss_arr1 = collect(map(x->x["loss"], perf_data1))
    z_arr1 = collect(map(x->x["Z"], perf_data1))
    loss_arr2 = collect(map(x->x["loss"], perf_data2))
    z_arr2 = collect(map(x->x["Z"], perf_data2))
    loss_arr3 = collect(map(x->x["loss"], perf_data3))
    z_arr3 = collect(map(x->x["Z"], perf_data3))
    loss_arr4 = collect(map(x->x["loss"], perf_data4))
    z_arr4 = collect(map(x->x["Z"], perf_data4))

    # Save data to txt
    save_variables = ["min_size" "max_size" "num_trials" "num_pairs" "perf_data1" "perf_data2" "perf_data3" "perf_data4"]
    open(file_path, "w") do io
        writedlm(io, save_variables)
        writedlm(io, [min_size,
                      max_size,
                      num_trials, # num_pairs
                      perf_data1,
                      perf_data2,
                      perf_data3,
                      perf_data4])
    end
# else
#     read_data = DelimitedFiles.readdlm(file_path)
#     variable_names = read_data[1]
#     min_size = read_data[2]
#     max_size = read_data[3]
#     num_trials = read_data[4]
#     num_pairs = read_data[5]
#     perf_data1 = read_data[6]
#     perf_data2 = read_data[7]
#     perf_data3 = read_data[8]
#     perf_data4 = read_data[9]

end

# Plot
plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), color = :red, legend=:left)
plot!(x, z_arr1, linewidth=2, label=L"$F_1$", linestyle=:dash, color =:red)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", color =:blue)
plot!(x, z_arr2, linewidth=2, label=L"$F_2$", linestyle=:dash, color =:blue)
plot!(x, loss_arr3, linewidth=2, label=L"$\eta_3$", color =:green)
plot!(x, z_arr3, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:green)
plot!(x, loss_arr4, linewidth=2, label=L"$\eta_3$", color =:purple)
plot!(x, z_arr4, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:purple)
xaxis!(L"$\textrm{Grid Size}$")

savefig("plots/cost_maxpaths.png")
savefig("plots/cost_maxpaths.pdf")

# Plot
plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), yaxis=:log, color = :red, legend=:left)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", yaxis=:log, color =:blue)
plot!(x[3:end], loss_arr3[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:green)
plot!(x[3:end], loss_arr4[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:purple)
xaxis!(L"$\textrm{Grid Size}$")

savefig("plots/cost_maxpathslog.png")
savefig("plots/cost_maxpathslog.pdf")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
using Test: @test, @testset, @inferred
using Statistics: mean, std, var, cor
using Random: randn!, bitrand
using LogExpFunctions: softmax
using RestrictedBoltzmannMachines: BinaryRBM, energy, free_energy, metropolis, metropolis!
import RestrictedBoltzmannMachines as RBMs

@testset "metropolis β=$β" for β in [0.5, 1.0, 2.0]
    N = 5
    M = 2
    rbm = BinaryRBM(randn(N), randn(M), randn(N, M) / √N)
    T = 10000
    B = 100
    v = bitrand(N, B, T)
    for t in 2:T
        v[:,:,t] .= metropolis(rbm, v[:,:,t - 1]; β)
    end
    counts = Dict{BitVector, Int}()
    for t in 1000:T, n in 1:B
        counts[v[:,n,t]] = get(counts, v[:,n,t], 0) + 1
    end
    freqs = Dict(v => c / sum(values(counts)) for (v,c) in counts)
    𝒱 = [BitVector(digits(Bool, x; base=2, pad=N)) for x in 0:(2^N - 1)]
    @test cor([get(freqs, v, 0.0) for v in 𝒱], softmax(-β * free_energy.(Ref(rbm), 𝒱))) > 0.99
end

@testset "metropolis! β=$β" for β in [0.5, 1.0, 2.0]
    N = 5
    M = 2
    rbm = BinaryRBM(randn(N), randn(M), randn(N, M) / √N)
    T = 10000
    B = 100
    v = bitrand(N, B, T)
    metropolis!(v, rbm; β)
    counts = Dict{BitVector, Int}()
    for t in 1000:T, n in 1:B
        counts[v[:,n,t]] = get(counts, v[:,n,t], 0) + 1
    end
    freqs = Dict(v => c / sum(values(counts)) for (v,c) in counts)
    𝒱 = [BitVector(digits(Bool, x; base=2, pad=N)) for x in 0:(2^N - 1)]
    @test cor([get(freqs, v, 0.0) for v in 𝒱], softmax(-β * free_energy.(Ref(rbm), 𝒱))) > 0.99
end


text:
@testset "ff f4 reduce sort invariant" begin
    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:lex)

    @test Groebner.groebner([x, y]) == Groebner.groebner([y, x]) == [y, x]

    fs = [x^2 + y, x*y]
    @test Groebner.groebner(fs) == [y^2, x*y, x^2 + y]

end

@testset "ff f4 yesreduce" begin

    root = Groebner.change_ordering(Groebner.rootn(3, ground=GF(2^31 - 1)), :degrevlex)
    x1, x2, x3 = gens(parent(first(root)))
    gb = Groebner.groebner(root, reduced=true)
    @test gb == [x1 + x2 + x3, x2^2 + x2*x3 + x3^2, x3^3 + 2147483646]

    root = Groebner.change_ordering(Groebner.rootn(6, ground=GF(2^31 - 1)), :degrevlex)
    x1, x2, x3, x4, x5, x6 = gens(parent(first(root)))
    gb = Groebner.groebner(root, reduced=true)
    @test gb == [
        x1 + x2 + x3 + x4 + x5 + x6,
        x2^2 + x2*x3 + x3^2 + x2*x4 + x3*x4 + x4^2 +
        x2*x5 + x3*x5 + x4*x5 + x5^2 + x2*x6 + x3*x6 + x4*x6 +
        x5*x6 + x6^2, x3^3 + x3^2*x4 + x3*x4^2 + x4^3 +
        x3^2*x5 + x3*x4*x5 + x4^2*x5 + x3*x5^2 + x4*x5^2 +
        x5^3 + x3^2*x6 + x3*x4*x6 + x4^2*x6 + x3*x5*x6 +
        x4*x5*x6 + x5^2*x6 + x3*x6^2 + x4*x6^2 + x5*x6^2 +
        x6^3, x4^4 + x4^3*x5 + x4^2*x5^2 + x4*x5^3 + x5^4 +
        x4^3*x6 + x4^2*x5*x6 + x4*x5^2*x6 + x5^3*x6 + x4^2*x6^2 +
        x4*x5*x6^2 + x5^2*x6^2 + x4*x6^3 + x5*x6^3 + x6^4, x5^5 +
        x5^4*x6 + x5^3*x6^2 + x5^2*x6^3 + x5*x6^4 + x6^5,
        x6^6 + 2147483646]

    ku = Groebner.change_ordering(Groebner.ku10(ground=GF(2^31-1)), :degrevlex)
    x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = gens(parent(first(ku)))
    gb = Groebner.groebner(ku, reduced=true)

    @test gb == [
        x9 + 1272065637*x10 + 875418006,
        x8 + 1529540685*x10 + 617942964,
        x7 + 1539832471*x10 + 607651173,
        x6 + 1314302432*x10 + 833181218,
        x5 + 1453635454*x10 + 693848197,
        x4 + 673118236*x10 + 1474365406,
        x3 + 269783061*x10 + 1877700587,
        x2 + 1042807874*x10 + 1104675778,
        x1 + 389079675*x10 + 1758403970,
        x10^2 + 1222705397*x10 + 924778249]

end


text:
# TODO: add waveguide sites for matched conditions
"""
	module Simulations

Package data used to define simulation and generate the relevant differential operators.
Exports are `Simulation`, `smooth!`, `update_dielectric!`, `update_pump`
"""
module Simulations

export Simulation
export smooth!
export update_dielectric!
export update_pump!

files = (
	"1D/Simulations1D.jl",
	"2D/Symmetric/Simulations2D.jl",
	# "3D/Simulations3D.jl"
	)

using ..Curlcurls
using ..Dispersions
using ..Domains
using ..Laplacians
using ..Lattices
using ..Points
using ..SelfEnergies
using ..Shapes
using ..VectorFields
# using LinearAlgebra
using SparseArrays
using Statistics
using RecipesBase

import ..AbstractDomain
import ..Symmetric, ..Unsymmetric

"""
	Simulation(ω₀, domains...; [k₂₀, k₃₀, k₁₀]) -> simulation

Build a simulation object which can be passed to the various solvers to find
resonance frequencies, etc.

There are three types of `domains`: `LatticeDomain`, `NondispersiveDomain`,
`DispersiveDomain`, which can be provided in any order, though at least one
`LatticeDomain` must be provided.
Where domains of the same type overlap, the one that appears first takes priority.
"""
struct Simulation{N,CLASS,T,TLDOM,TNDOM,TDDOM,TSE}
	lattice_domains::TLDOM
	nondispersive_domains::TNDOM
	dispersive_domains::TDDOM

	x::Vector{Point{N,Cartesian}}
	lattice_domain_indices::Vector{Int}
	nondispersive_domain_indices::Vector{Int}
	dispersive_domain_indices::Vector{Int}
	ε::Vector{ComplexF64}
	F::Vector{Float64}
	χ::Vector{AbstractDispersion}
	Fs::Vector{Vector{Float64}}

	laplacian::Laplacian{N,CLASS}
	# curlcurl::Curlcurl{N,CLASS}
	self_energy::TSE
	α::NTuple{N,Vector{ComplexF64}}
	σ::NTuple{N,Vector{ComplexF64}}
	x_half::NTuple{N,Vector{Point{N,Cartesian}}}
	σ_half::NTuple{N,Vector{ComplexF64}}
	ω₀::Float64
	k₁₀::Float64
	k₂₀::Float64
	k₃₀::Float64
	smoothed::Base.RefValue{Bool}
end

# for arbitrary domain order
function Simulation(
			ω₀::Real,
			domains::Vararg{AbstractDomain{N}};
			kwargs...
			) where N

	lattice_domains = _get_lattice_domains(domains...)
	nondispersive_domains = _get_nondispersive_domains(domains...)
	dispersive_domains = _get_dispersive_domains(domains...)

	# check
	length(lattice_domains)>0 || throw(ArgumentError("must provide at least one `LatticeDomain`, but 0 provided instead"))
	return Simulation(ω₀, lattice_domains, nondispersive_domains, dispersive_domains; kwargs...)
end

# functions to extract lattice_domains, etc, from arbitrary list of domains
fnames = (:_get_lattice_domains, :_get_nondispersive_domains, :_get_dispersive_domains)
types = (LatticeDomain, NondispersiveDomain, DispersiveDomain)
for i ∈ eachindex(types)
	@eval begin
		$(fnames[i])(ldom::$(types[i]), args...) = (ldom,$(fnames[i])(args...)...)
		$(fnames[i])(arg, args...) = $(fnames[i])(args...)
		$(fnames[i])(ldom1::$(types[i]), ldom2::$(types[i]), args...) = (ldom1,ldom2,$(fnames[i])(args...)...)
		$(fnames[i])(arg, ldom::$(types[i]), args...) = (ldom,$(fnames[i])(args...)...)
		$(fnames[i])(ldom::$(types[i]), arg, args...) = (ldom,$(fnames[i])(args...)...)
		$(fnames[i])(arg1, arg2, args...) = $(fnames[i])(args...)
		$(fnames[i])() = ()
	end
end

Base.length(sim::Simulation) = length(sim.x)
Base.ndims(sim::Simulation{N}) where N = N

# Base.conj(sim::Simulation) = Simulation(sim.ω₀,map(conj,sim.domains)...)

import ..NUM_SUBPIXELS

"""
	smooth!(sim,[num_sub_pixel])

Use sub-pixel sampling along domain walls to smooth the dieletric and pump profiles.

`num_sub_pixel` defaults to `NUM_SUB_PIXEL` in module `Defaults`, which is currently set to 5.
"""
function smooth!(sim::Simulation, num_sub_pixel::Integer = NUM_SUBPIXELS)
    smooth_dielectric!(sim, num_sub_pixel)
	smooth_pump!(sim, num_sub_pixel)
	sim.smoothed[] = true
    return sim
end

"""
	update!(sim)

update simulation `sim` by recomputing the dielectric and pump functions
"""
function VectorFields.update!(sim::Simulation)
	update_dielectric!(sim)
	update_pump!(sim)
	return nothing
end

"""
	update_dielectric!(sim,[ε])

update simulation `sim` by recomputing the dielectric function or by giving it an explicit array `ε`
"""
function update_dielectric!(sim::Simulation)
	simulation_dielectric!(sim.ε, sim, sim.x)
	sim.smoothed[] = false
	return nothing
end
function update_dielectric!(sim::Simulation,ε::AbstractVecOrMat)
	copyto!(sim.ε,ε)
	sim.smoothed[] = false
	return nothing
end

"""
	update_pump!(sim,[F])

update simulation `sim` by recomputing the pump function or by giving it an explicit array `F`
"""
function update_pump!(sim::Simulation)
	simulation_pump!(sim.F, sim, sim.x)
	sim.smoothed[] = false
	return nothing
end
function update_pump!(sim::Simulation,F::AbstractVecOrMat)
	copyto!(sim.F, F)
	sim.smoothed[] = false
	return nothing
end

# the dielectric as a function of position (x,y)
# function simulation_dielectric(sim::Simulation, x::Vector)
	# ε = Vector{ComplexF64}(undef,length(x))
	# simulation_dielectric!(ε, sim , x)
	# return ε
# end
function simulation_dielectric!(ε::Vector, sim::Simulation, x::Vector, backgrounds=map(d->d.ε, sim.lattice_domains), dielectrics=map(d->d.dielectric, sim.nondispersive_domains))
    for i ∈ eachindex(x) ε[i] = simulation_dielectric(sim, x[i], backgrounds, dielectrics) end
    return nothing
end
@inline function simulation_dielectric(sim::Simulation, x::Point, backgrounds=map(d->d.ε, sim.lattice_domains), dielectrics=map(d->d.dielectric, sim.nondispersive_domains))
    lattice_index = which_domain(sim.lattice_domains, x)
	nondispersive_index = which_domain(sim.nondispersive_domains, x)
	if nondispersive_index==0
		return backgrounds[lattice_index]
	else
	    return dielectrics[nondispersive_index](x)
	end
end

# the dielectric as a function of position (x,y)
# function simulation_pump(sim::Simulation, x::Vector)
# 	F = Vector{Float64}(undef,length(x))
# 	simulation_pump!(F, sim , x)
# 	return F
# end
function simulation_pump!(F::Vector, sim::Simulation, x::Vector, pumps=map(d->d.pump, sim.dispersive_domains))
    for i ∈ eachindex(x) F[i] = simulation_pump(sim, x[i], pumps) end
    return nothing
end
@inline function simulation_pump(sim::Simulation, x::Point, pumps=map(d->d.pump, sim.dispersive_domains))
	dispersive_index = which_domain(sim.dispersive_domains, x)
	if dispersive_index==0
		return 0.0
	else
	    return pumps[dispersive_index](x)
	end
end

################################################################################
# utilities

# identify which domain is associated with a given spatial position (x)
function which_domains(domains, x::Vector)
    indices = Vector{Int}(undef, length(x))
	shapes = map(z->z.shape, domains)
    which_domains!(indices, domains, shapes, x)
    return indices
end
function which_domains!(indices::Vector{Int}, domains, shapes, x)
    @inbounds for i ∈ eachindex(x) indices[i] = which_domain(domains, shapes, x[i]) end
    return nothing
end
# above is for arrays of points, below is for single points
function which_domain(domains,x::Point)
	shapes = map(z->z.shape,domains)
    return which_domain(domains, shapes, x)
end
@inline function which_domain(domains, shapes, x::Point)
	index = findfirst(z->z(x),shapes)
    isnothing(index) ? index=0 : nothing
    iszero(index) ? nothing : (isvoid(domains[index]) ? index=0 : nothing)
    return index
end

isvoid(domain::AbstractDomain) = domain.type==:Void

smooth_dielectric!() = nothing
smooth_pump!() = nothing

foreach(include,files)

# take individually defined domains and concatenate all relevant data
function append_domains!(domain_index,nnm,nnp,name,x,indices,ε,F,interior,bulk,surface,corner,domains)
	count = 1
	for d ∈ domains
		append!(domain_index,fill(count,length(d.x)))
		for i ∈ eachindex(d.nnm)
			append!(nnm[i],d.nnm[i].+length(nnm[i]))
			append!(nnp[i],d.nnp[i].+length(nnp[i]))
		end
		append!(name,fill(d.name,length(d.x)))
		append!(x,d.x)
		append!(indices,d.indices)
		append!(ε,d.ε)
		append!(F,d.F)
		append!(interior,d.interior)
		append!(bulk,d.bulk)
		append!(surface,d.surface)
		append!(corner,d.corner)
		count += 1
	end
	return nothing
end

# identify which sites are to be removed by determining which sites belong to the topmost domain.
get_sites_to_remove(domains,dom,x) = which_domains(domains,x).!==dom

# remove the pre-computed nearest neighbors, which makes the bulk nn calculation very fast
function trim_nn!(nnm,nnp,removed,bulk)
	cr = cumsum(removed)
	for j ∈ eachindex(nnm)
		for i ∈ eachindex(nnm[j])
			if bulk[i]
				nnm[j][i] -= cr[nnm[j][i]]
				nnp[j][i] -= cr[nnp[j][i]]
			end
		end
	end
end

# remove the sites determined by get_sites_to_remove
remove_sites!(sites::Tuple,removed) = map(z->remove_sites!(z,removed),sites)
remove_sites!(site::NTuple{N,Array},removed) where N = map(z->deleteat!(z,findall(removed)),site)
remove_sites!(site,removed) = deleteat!(site,findall(removed))

################################################################################

VectorField{M}(sim::Simulation{N}, args...) where {N,M} = VectorField{N,M}(sim.x, args...)
VectorField{N,M}(sim::Simulation{N}, args...) where {N,M} = VectorField{N,M}(sim.x, args...)

################################################################################
#Pretty Printing

import ..PRINTED_COLOR_GOOD
import ..PRINTED_COLOR_WARN
import ..PRINTED_COLOR_DARK

function Base.show(io::IO,sim::Simulation{N,CLASS}) where {N,CLASS}
	suffix = length(sim.domains)>1 ? "s" : ""
	if sim.smoothed[]
		printstyled(io,"Smoothed ",color=PRINTED_COLOR_GOOD)
		print(io,"$(N)D $CLASS ")
		printstyled(io," Simulation",color=PRINTED_COLOR_DARK)
		println(io, " with ", length(sim.domains)," domain",suffix)
	else
		printstyled(io,"Unsmoothed ",color=PRINTED_COLOR_WARN)
		print(io,N, "D ", CLASS)
		printstyled(io," Simulation",color=PRINTED_COLOR_DARK)
		println(io, " with ", length(sim.domains)," domain",suffix)
	end
    print(io,"\tn sites: ")
	printstyled(io,length(sim),"\n",color=:light_cyan)
	println(io,"\t\t====================")
	for d ∈ eachindex(sim.domains)
		printstyled(io,"\t\tDomain ",color=PRINTED_COLOR_DARK)
		print(io, d," (",sim.domains[d].type,"): ")
		print(io,sim.domains[d].name)
		d < length(sim.domains) ? println(io) : nothing
	end
end

@recipe f(v::VectorField, sim::Simulation) = sim, v
@recipe f(e::VectorField,sim::Simulation,by::Function) = sim,e,by
@recipe f(sim::Simulation,by::Function,e::VectorField) = sim,e,by
@recipe f(e::VectorField,by::Function,sim::Simulation) = sim,e,by
@recipe f(by::Function,sim::Simulation,e::VectorField) = sim,e,by
@recipe f(by::Function,e::VectorField,sim::Simulation) = sim,e,by

end # module


text:
struct Data{T <: AbstractVector{<:Real}}
  yC::T
  yT::T
  yC_finite::T
  yT_finite::T
  NC::Int
  NT::Int
  NC_finite::Int
  NT_finite::Int
  ZC::Int
  ZT::Int
end

function Data(yC::AbstractVector{<:Real}, yT::AbstractVector{<:Real})
  yC_finite = yC[isfinite.(yC)]
  yT_finite = yT[isfinite.(yT)]
  NC = length(yC)
  NT = length(yT)
  NC_finite = length(yC_finite)
  NT_finite = length(yT_finite)
  ZC = NC - NC_finite
  ZT = NT - NT_finite
  Data(yC, yT, yC_finite, yT_finite, NC, NT, NC_finite, NT_finite, ZC, ZT)
end

const samplenames = ('C', 'T')

ref_yfinite(data::Data, i::Char) = (i == 'C') ? data.yC_finite : data.yT_finite
ref_N(data::Data, i::Char) = (i == 'C') ? data.NC : data.NT
ref_Nfinite(data::Data, i::Char) = (i == 'C') ? data.NC_finite : data.NT_finite
ref_Z(data::Data, i::Char) = (i == 'C') ? data.ZC : data.ZT

function subsample_data(data::Data, nc::Int, nt::Int)
  if nc > 0 && nt >0
    return Data(sample(data.yC, nc, replace=false),
                sample(data.yT, nt, replace=false))
  else
    return data
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
struct Matern{D,L} <: UnknownAcvTimeSeriesModel{D,Float64}
    σ::SHermitianCompact{D,Float64,L}
    ν::SHermitianCompact{D,Float64,L}
    a::SHermitianCompact{D,Float64,L}
    σ²::SHermitianCompact{D,Float64,L}
    a²::SHermitianCompact{D,Float64,L}
    νplushalf::SHermitianCompact{D,Float64,L}
    sdfconst::SHermitianCompact{D,Float64,L}
    ∂ν_part::SHermitianCompact{D,Float64,L}
    ∂a_part1::SHermitianCompact{D,Float64,L}
    ∂a_part2::SHermitianCompact{D,Float64,L}
    ∂ν²_part::SHermitianCompact{D,Float64,L}
    ∂ν∂a_part::SHermitianCompact{D,Float64,L}
    ∂a²_part::SHermitianCompact{D,Float64,L}
    function Matern{D,L}(θ) where {D,L}
        @boundscheck checkparameterlength(θ, MaternAcv{D,L})
        L == triangularnumber(D) || error("Matern{D,L} should satisfy L == D*(D+1)÷2")
        all(x->x>0, θ) || throw(ArgumentError("all parameters of Matern should be > 0."))
        σ = @views SHermitianCompact(SVector{L,Float64}(θ[1:L]))
        all(i==j ? true : σ[i,j] < 1 for i in 1:size(σ,1) for j in 1:i) || throw(ArgumentError("ρ parameters must be < 1."))
        ν = @views SHermitianCompact(SVector{L,Float64}(θ[L+1:2L]))
        a = @views SHermitianCompact(SVector{L,Float64}(θ[2L+1:end]))
        σ² = SHermitianCompact((σ.lowertriangle).^2)
        a² = SHermitianCompact((a.lowertriangle).^2)
        νplushalf = SHermitianCompact((ν.lowertriangle).+0.5)
        variance_part = SHermitianCompact(SMatrix{D,D,Float64}(i==j ? σ[i,i]^2 : σ[i,i]*σ[j,j]*σ[i,j] for i in 1:D, j in 1:D))
        sdfconst = SHermitianCompact(matern_sdf_normalising.(ν.lowertriangle,a.lowertriangle).*variance_part.lowertriangle)
        ∂ν_part  = SHermitianCompact(digamma.(νplushalf.lowertriangle).-digamma.(ν.lowertriangle).+2.0.*log.(a.lowertriangle))
        ∂a_part1 = SHermitianCompact(2.0.*ν.lowertriangle./a.lowertriangle)
        ∂a_part2 = SHermitianCompact(a.lowertriangle.*(2.0.*ν.lowertriangle.+1))
        ∂ν²_part  = SHermitianCompact(trigamma.(νplushalf.lowertriangle).-trigamma.(ν.lowertriangle))
        ∂ν∂a_part = SHermitianCompact(1.0./a.lowertriangle)
        ∂a²_part = SHermitianCompact(ν.lowertriangle./(a².lowertriangle))
        new{D,L}(σ,ν,a,σ²,a²,νplushalf,sdfconst,∂ν_part,∂a_part1,∂a_part2,∂ν²_part,∂ν∂a_part,∂a²_part)
    end
end
matern_sdf_normalising(ν,a) = gamma(ν+0.5)*a^(2ν) / (gamma(ν)*sqrt(π))

const Matern2D = Matern{2,3}
const Matern3D = Matern{3,6}
const Matern4D = Matern{4,10}

npars(::Type{Matern{D,L}}) where {D,L} = 3triangularnumber(D)
minbins(::Type{Matern{D,L}}) where {D,L} = 4096
nalias(model::Matern) = 3
lowerbounds(m::Type{Matern{D,L}}) where {D,L} = upperbounds(m)
upperbounds(::Type{Matern{D,L}}) where {D,L} =
    error("Bounds not specified for the multivariate matern model here as parameter space is complicated in
    this formulation. This should be treated with greater care than fit currently does.")

function parameternames(::Type{Matern{D,L}}) where {D,L}
    σ = reduce(vcat,[i==j ? "σ_$i" : "ρ_$i$j" for j in i:D] for i in 1:D)
    ν = reduce(vcat,[i==j ? "ν_$i" : "ν_$i$j" for j in i:D] for i in 1:D)
    a = reduce(vcat,[i==j ? "a_$i" : "a_$i$j" for j in i:D] for i in 1:D)
    return Tuple(Symbol.([σ;ν;a]))
end

function add_sdf!(out, model::Matern{D,L}, ω) where {D,L}
    ω² = ω^2
    count = 1
    for i ∈ 1:D, j ∈ i:D
        out[count] += model.sdfconst[i,j]/((model.a²[i,j]+ω²)^(model.νplushalf[i,j]))
        count += 1
    end
    return nothing
end

function grad_add_sdf!(out, model::Matern{D,L}, ω) where {D,L}
    ω² = ω^2
    count = 1
    for i ∈ 1:D, j ∈ i:D
        asq_plus_omsq = model.a²[i,j]+ω²
        sdf = model.sdfconst[i,j]/((asq_plus_omsq)^(model.νplushalf[i,j]))
        # σ/ρ
        if i == j
            out[count,count] += 2sdf/model.σ[i,i]
        else
            out[count,count] += sdf/model.σ[i,j]
            out[count,indexLT(i,i,D)] += sdf/model.σ[i,i] # derivative with respect to σᵢ
            out[count,indexLT(j,j,D)] += sdf/model.σ[j,j] # derivative with respect to σⱼ
        end
        # ν ## digamma(model.νplushalf[i,j])-digamma(model.ν[i,j])+2log(model.a[i,j]) # precomputed as no dependency on ω
        out[count,L+count] += sdf * (model.∂ν_part[i,j]-log(asq_plus_omsq))
        # a ## 2model.ν[i,j]/model.a[i,j] - model.a[i,j]*(2model.ν[i,j]+1)/... # precomputed as no dependency on ω
        out[count,2L+count] += sdf * (model.∂a_part1[i,j] - model.∂a_part2[i,j]/asq_plus_omsq)
        count += 1
    end
    return nothing
end

function hess_add_sdf!(out, model::Matern{D,L}, ω) where {D,L}
    ω² = ω^2
    count = 1
    for i ∈ 1:D, j ∈ i:D
        asq_plus_omsq = model.a²[i,j]+ω²
        sdf = model.sdfconst[i,j]/((asq_plus_omsq)^(model.νplushalf[i,j]))
        sdf∂ν = sdf * (model.∂ν_part[i,j]-log(asq_plus_omsq))
        sdf∂a = sdf * (model.∂a_part1[i,j] - (model.∂a_part2[i,j])/asq_plus_omsq)
        # σ/ρ
        if i == j
            out[count,indexLT(count,count,   3L)]   += 2sdf/model.σ²[i,i]  # ∂σᵢ²
            out[count,indexLT(count,count+L, 3L)]   += 2sdf∂ν/model.σ[i,i] # ∂σᵢ∂νᵢⱼ
            out[count,indexLT(count,count+2L,3L)]   += 2sdf∂a/model.σ[i,i] # ∂σᵢ∂aᵢⱼ
        else
            iind = indexLT(i,i,D)
            jind = indexLT(j,j,D)
            out[count,indexLT(iind, count,   3L)] += sdf/(model.σ[i,i]*model.σ[i,j]) # ∂σᵢ∂ρᵢⱼ
            out[count,indexLT(iind, jind,    3L)] += sdf/(model.σ[i,i]*model.σ[j,j]) # ∂σᵢ∂σⱼ
            out[count,indexLT(jind, count,   3L)] += sdf/(model.σ[i,j]*model.σ[j,j]) # ∂σⱼ∂ρᵢⱼ
            out[count,indexLT(iind, count+L, 3L)] += sdf∂ν/model.σ[i,i] # ∂σᵢ∂νᵢⱼ
            out[count,indexLT(iind, count+2L,3L)] += sdf∂a/model.σ[i,i] # ∂σᵢ∂aᵢⱼ
            out[count,indexLT(count,count+L, 3L)] += sdf∂ν/model.σ[i,j] # ∂ρᵢⱼ∂νᵢⱼ
            out[count,indexLT(count,count+2L,3L)] += sdf∂a/model.σ[i,j] # ∂ρᵢⱼ∂aᵢⱼ
            out[count,indexLT(jind, count+L, 3L)] += sdf∂ν/model.σ[j,j] # ∂σⱼ∂νᵢⱼ
            out[count,indexLT(jind, count+2L,3L)] += sdf∂a/model.σ[j,j] # ∂σⱼ∂aᵢⱼ
        end
        # ∂νᵢⱼ²
        out[count,indexLT(count+L, count+L, 3L)] += sdf∂ν^2/sdf + sdf*model.∂ν²_part[i,j]
        # ∂νᵢⱼ∂aᵢⱼ
        out[count,indexLT(count+L, count+2L,3L)] += sdf∂ν*sdf∂a/sdf + 2sdf*(model.∂ν∂a_part[i,j]-model.a[i,j]/asq_plus_omsq)
        # ∂aᵢⱼ²
        out[count,indexLT(count+2L,count+2L,3L)] += sdf∂a^2/sdf + 2sdf*(model.νplushalf[i,j]*(2model.a²[i,j]/(asq_plus_omsq^2) - 1/asq_plus_omsq) - model.∂a²_part[i,j])
        count += 1
    end
    return nothing
end

## Univariate

struct Matern1D <: UnknownAcvTimeSeriesModel{1,Float64}
    σ::Float64
    ν::Float64
    a::Float64
    σ²::Float64
    a²::Float64
    νplushalf::Float64
    sdfconst::Float64
    ∂ν_part::Float64
    ∂a_part1::Float64
    ∂a_part2::Float64
    ∂ν²_part::Float64
    ∂ν∂a_part::Float64
    ∂a²_part::Float64
    function Matern1D(σ,ν,a)
        σ > 0 || throw(ArgumentError("Matern1D process requires 0 < σ."))
        ν > 0 || throw(ArgumentError("Matern1D process requires 0 < ν."))
        a > 0 || throw(ArgumentError("Matern1D process requires 0 < a."))
        sdfconst = matern_sdf_normalising(ν,a)*σ^2
        ∂ν_part  = digamma(ν+0.5)-digamma(ν)+2.0*log(a)
        ∂a_part1 = 2.0ν/a
        ∂a_part2 = a*(2.0*ν+1)
        ∂ν²_part  = trigamma(ν+0.5)-trigamma(ν)
        ∂ν∂a_part = 1.0/a
        ∂a²_part = ν/(a^2)
        new(σ, ν, a, σ^2, a^2, ν+0.5, sdfconst, ∂ν_part, ∂a_part1, ∂a_part2, ∂ν²_part, ∂ν∂a_part, ∂a²_part)
    end
    function Matern1D(x::AbstractVector{Float64})
        length(x) == npars(Matern1D) || throw(ArgumentError("Matern1D process has $(npars(Matern1D)) parameters, but $(length(x)) were provided."))
        @inbounds Matern1D(x[1], x[2], x[3])
    end
end

npars(::Type{Matern1D}) = 3
minbins(::Type{Matern1D}) = 4096
nalias(model::Matern1D) = 3
lowerbounds(::Type{Matern1D}) = [0,0,0]
upperbounds(::Type{Matern1D}) = [Inf,Inf,Inf]

sdf(model::Matern1D, ω) = model.sdfconst/((model.a²+ω^2)^(model.νplushalf))

function grad_add_sdf!(out, model::Matern1D, ω)
    ω² = ω^2
    asq_plus_omsq = model.a²+ω²
    sdf = model.sdfconst/((asq_plus_omsq)^(model.νplushalf))
    # σ
    out[1] += 2sdf/model.σ
    # ν
    out[2] += sdf * (model.∂ν_part-log(asq_plus_omsq))
    # a
    out[3] += sdf * (model.∂a_part1 - model.∂a_part2/asq_plus_omsq)

    return nothing
end

function hess_add_sdf!(out, model::Matern1D, ω)
    ω² = ω^2
    asq_plus_omsq = model.a²+ω²
    sdf = model.sdfconst/((asq_plus_omsq)^(model.νplushalf))
    sdf∂ν = sdf * (model.∂ν_part-log(asq_plus_omsq))
    sdf∂a = sdf * (model.∂a_part1 - (model.∂a_part2)/asq_plus_omsq)

    # ∂σ²
    out[1] += 2sdf/model.σ²
    # ∂σ∂ν
    out[2] += 2sdf∂ν/model.σ
    # ∂σ∂a
    out[3] += 2sdf∂a/model.σ
    # ∂ν²
    out[4] += sdf∂ν^2/sdf + sdf*model.∂ν²_part
    # ∂ν∂a
    out[5] += sdf∂ν*sdf∂a/sdf + 2sdf*(model.∂ν∂a_part-model.a/asq_plus_omsq)
    # ∂a²
    out[6] += sdf∂a^2/sdf + 2sdf*(model.νplushalf*(2model.a²/(asq_plus_omsq^2) - 1/asq_plus_omsq) - model.∂a²_part)

    return nothing
end


text:
using Compat
import LibGit2
import JSON
const OS_ARCH_CoolProp = (Sys.WORD_SIZE == 64) ? "64bit" : (Sys.iswindows() ? "32bit__cdecl" : "32bit");
const destpathbase = abspath(@__FILE__, "..", "lib");

# const branchname = LibGit2.branch(LibGit2.GitRepo(abspath(@__FILE__, "..", "..")));
const branchname = "master"
@info "On $branchname";
latestVersion_CoolProp = "";
_download(s, d) = isfile(d) ? throw("file exists ...") : begin
    @info "Downloading $s => $d ."
    download(s, d);
    @info "Done.";
end
if isdir(destpathbase)
    @warn "$destpathbase exists make sure to remove old library files.";
else
    mkdir(destpathbase);
end
try
    latestVersion_CoolProp = JSON.parse(read(download("https://sourceforge.net/projects/coolprop/best_release.json"), String))["release"]["filename"][11:15];
catch err
    latestVersion_CoolProp = "6.1.0";
    @warn "unable to download may be a windows machine firewall.. , set latestVersion_CoolProp = $latestVersion_CoolProp";
end
latestVersion_CoolProp = JSON.parse(read(download("https://sourceforge.net/projects/coolprop/best_release.json"), String))["release"]["filename"][11:15];
coolpropurlbase = "http://netix.dl.sourceforge.net/project/coolprop/CoolProp/$latestVersion_CoolProp/";

(branchname == "nightly") && (global coolpropurlbase = "http://www.coolprop.dreamhosters.com/binaries/");
(branchname == "master") && (global coolpropurlbase = "http://netix.dl.sourceforge.net/project/coolprop/CoolProp/$latestVersion_CoolProp/");
try
    println("CoolProp latestVersion = $latestVersion_CoolProp ...")
    # _download(coolpropurlbase * "Julia/CoolProp.jl", joinpath(destpathbase,"CoolProp.jl"));
    @info "I'm Getting CoolProp Binaries...";
    @static if  Sys.iswindows()
        urlbase = coolpropurlbase * "shared_library/Windows/$OS_ARCH_CoolProp/";
        _download(joinpath(urlbase,"CoolProp.dll"), joinpath(destpathbase,"CoolProp.dll"));
        _download(joinpath(urlbase,"CoolProp.lib"), joinpath(destpathbase,"CoolProp.lib"));
        _download(joinpath(urlbase,"exports.txt"), joinpath(destpathbase,"exports.txt"));
    end
    @static if  Sys.islinux()
        (branchname == "nightly") && (global coolpropurlbase = "http://netix.dl.sourceforge.net/project/coolprop/CoolProp/nightly/");
        urlbase = coolpropurlbase * "shared_library/Linux/$OS_ARCH_CoolProp/libCoolProp.so.$latestVersion_CoolProp";
        _download(urlbase, joinpath(destpathbase,"CoolProp.so"));
    end
    @static if  Sys.isapple()
        urlbase = coolpropurlbase * "shared_library/Darwin/$OS_ARCH_CoolProp/libCoolProp.dylib";
        _download(urlbase, joinpath(destpathbase,"CoolProp.dylib"));
    end
    (branchname == "nightly") && begin
        @info "Building help tables..."
        touch(joinpath(destpathbase, "fluids.table"))
        touch(joinpath(destpathbase, "parameters.table"))
        import CoolProp: buildfluids , buildparameters
        buildfluids();
        buildparameters();
    end
catch err
    println("Build error: $err")
end


text:
@testset "not equal" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, y[1:3], CS.Integers([-3, 1, 2, 3]))
    @constraint(m, sum(y) + 1 != 5)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]

    # doesn't check the length
    # 1+2+1 + constant (1) == 5
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 1])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 2])

    constr_indices = constraint.indices
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        -3,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        1,
    )
    CS.fix!(com, com.search_space[constr_indices[1]], 2)
    CS.fix!(com, com.search_space[constr_indices[2]], 1)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        1,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1
    # now setting it to 1 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        1,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, y[1:3], CS.Integers([-3, 1, 2, 3]))
    @constraint(m, sum(y) + 1 != 5)
    optimize!(m)

    com = CS.get_inner_model(m)
    constraint = get_constraints_by_type(com, CS.LinearConstraint)[1]

    CS.fix!(com, com.search_space[constr_indices[1]], 2)
    CS.fix!(com, com.search_space[constr_indices[2]], 1)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(com.search_space[constr_indices[3]])) == [-3, 2, 3]
end

@testset "not equal is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) != 7)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 5; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) != 7)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 5; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end


text:
# Este é um comentário em Julia, como em Python
print("Olá mundo!")

max_stars_repo_name: thecrypters/julia-demo

max_issues_repo_name: thecrypters/julia-demo

max_forks_repo_name: thecrypters/julia-demo


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

# Uniform random generation

# This file contains the creation of Sampler objects and the associated generation of
# random values from them. More specifically, given the specification S of a set
# of values to pick from (e.g. 1:10, or "a string"), we define
#
# 1) Sampler(rng, S, ::Repetition) -> sampler
# 2) rand(rng, sampler) -> random value
#
# Note that the 1) is automated when the sampler is not intended to carry information,
# i.e. the default fall-backs SamplerType and SamplerTrivial are used.

## from types: rand(::Type, [dims...])

### random floats

Sampler(::Type{RNG}, ::Type{T}, n::Repetition) where {RNG<:AbstractRNG,T<:AbstractFloat} =
    Sampler(RNG, CloseOpen01(T), n)

# generic random generation function which can be used by RNG implementors
# it is not defined as a fallback rand method as this could create ambiguities

rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float16}}) =
    Float16(reinterpret(Float32,
                        (rand(r, UInt10(UInt32)) << 13)  | 0x3f800000) - 1)

rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01{Float32}}) =
    reinterpret(Float32, rand(r, UInt23()) | 0x3f800000) - 1

rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen12_64}) =
    reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt52()))

rand(r::AbstractRNG, ::SamplerTrivial{CloseOpen01_64}) = rand(r, CloseOpen12()) - 1.0

#### BigFloat

const bits_in_Limb = sizeof(Limb) << 3
const Limb_high_bit = one(Limb) << (bits_in_Limb-1)

struct SamplerBigFloat{I<:FloatInterval{BigFloat}} <: Sampler{BigFloat}
    prec::Int
    nlimbs::Int
    limbs::Vector{Limb}
    shift::UInt

    function SamplerBigFloat{I}(prec::Int) where I<:FloatInterval{BigFloat}
        nlimbs = (prec-1) ÷ bits_in_Limb + 1
        limbs = Vector{Limb}(undef, nlimbs)
        shift = nlimbs * bits_in_Limb - prec
        new(prec, nlimbs, limbs, shift)
    end
end

Sampler(::Type{<:AbstractRNG}, I::FloatInterval{BigFloat}, ::Repetition) =
    SamplerBigFloat{typeof(I)}(precision(BigFloat))

function _rand(rng::AbstractRNG, sp::SamplerBigFloat)
    z = BigFloat()
    limbs = sp.limbs
    rand!(rng, limbs)
    @inbounds begin
        limbs[1] <<= sp.shift
        randbool = iszero(limbs[end] & Limb_high_bit)
        limbs[end] |= Limb_high_bit
    end
    z.sign = 1
    GC.@preserve limbs unsafe_copyto!(z.d, pointer(limbs), sp.nlimbs)
    (z, randbool)
end

function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::CloseOpen12{BigFloat})
    z = _rand(rng, sp)[1]
    z.exp = 1
    z
end

function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::CloseOpen01{BigFloat})
    z, randbool = _rand(rng, sp)
    z.exp = 0
    randbool &&
        ccall((:mpfr_sub_d, :libmpfr), Int32,
              (Ref{BigFloat}, Ref{BigFloat}, Cdouble, Int32),
              z, z, 0.5, Base.MPFR.ROUNDING_MODE[])
    z
end

# alternative, with 1 bit less of precision
# TODO: make an API for requesting full or not-full precision
function _rand(rng::AbstractRNG, sp::SamplerBigFloat, ::CloseOpen01{BigFloat}, ::Nothing)
    z = _rand(rng, sp, CloseOpen12(BigFloat))
    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Culong, Int32),
          z, z, 1, Base.MPFR.ROUNDING_MODE[])
    z
end

rand(rng::AbstractRNG, sp::SamplerBigFloat{T}) where {T<:FloatInterval{BigFloat}} =
    _rand(rng, sp, T())

### random integers

#### UniformBits

rand(r::AbstractRNG, ::SamplerTrivial{UInt10Raw{UInt16}}) = rand(r, UInt16)
rand(r::AbstractRNG, ::SamplerTrivial{UInt23Raw{UInt32}}) = rand(r, UInt32)

rand(r::AbstractRNG, ::SamplerTrivial{UInt52Raw{UInt64}}) =
    _rand52(r, rng_native_52(r))

_rand52(r::AbstractRNG, ::Type{Float64}) = reinterpret(UInt64, rand(r, CloseOpen12()))
_rand52(r::AbstractRNG, ::Type{UInt64})  = rand(r, UInt64)

rand(r::AbstractRNG, ::SamplerTrivial{UInt104Raw{UInt128}}) =
    rand(r, UInt52Raw(UInt128)) << 52 ⊻ rand(r, UInt52Raw(UInt128))

rand(r::AbstractRNG, ::SamplerTrivial{UInt10{UInt16}})   = rand(r, UInt10Raw())  & 0x03ff
rand(r::AbstractRNG, ::SamplerTrivial{UInt23{UInt32}})   = rand(r, UInt23Raw())  & 0x007fffff
rand(r::AbstractRNG, ::SamplerTrivial{UInt52{UInt64}})   = rand(r, UInt52Raw())  & 0x000fffffffffffff
rand(r::AbstractRNG, ::SamplerTrivial{UInt104{UInt128}}) = rand(r, UInt104Raw()) & 0x000000ffffffffffffffffffffffffff

rand(r::AbstractRNG, sp::SamplerTrivial{<:UniformBits{T}}) where {T} =
        rand(r, uint_default(sp[])) % T

#### BitInteger

# rand_generic methods are intended to help RNG implementors with common operations
# we don't call them simply `rand` as this can easily contribute to create
# amibuities with user-side methods (forcing the user to resort to @eval)

rand_generic(r::AbstractRNG, T::Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32}) =
    rand(r, UInt52Raw()) % T[]

rand_generic(r::AbstractRNG, ::Type{UInt64}) =
    rand(r, UInt52Raw()) << 32 ⊻ rand(r, UInt52Raw())

rand_generic(r::AbstractRNG, ::Type{UInt128}) = _rand128(r, rng_native_52(r))

_rand128(r::AbstractRNG, ::Type{UInt64}) =
    ((rand(r, UInt64) % UInt128) << 64) ⊻ rand(r, UInt64)

function _rand128(r::AbstractRNG, ::Type{Float64})
    xor(rand(r, UInt52Raw(UInt128))  << 96,
        rand(r, UInt52Raw(UInt128))  << 48,
        rand(r, UInt52Raw(UInt128)))
end

rand_generic(r::AbstractRNG, ::Type{Int128}) = rand(r, UInt128) % Int128
rand_generic(r::AbstractRNG, ::Type{Int64})  = rand(r, UInt64) % Int64

### random complex numbers

rand(r::AbstractRNG, ::SamplerType{Complex{T}}) where {T<:Real} =
    complex(rand(r, T), rand(r, T))

### random characters

# returns a random valid Unicode scalar value (i.e. 0 - 0xd7ff, 0xe000 - # 0x10ffff)
function rand(r::AbstractRNG, ::SamplerType{T}) where {T<:AbstractChar}
    c = rand(r, 0x00000000:0x0010f7ff)
    (c < 0xd800) ? T(c) : T(c+0x800)
end

## Generate random integer within a range

### BitInteger

# there are two implemented samplers for unit ranges, which assume that Float64 (i.e.
# 52 random bits) is the native type for the RNG:
# 1) "Fast", which is the most efficient when the underlying RNG produces rand(Float64)
#     "fast enough". The tradeoff is faster creation of the sampler, but more
#     consumption of entropy bits
# 2) "Default" which tries to use as few entropy bits as possible, at the cost of a
#    a bigger upfront price associated with the creation of the sampler

#### helper functions

uint_sup(::Type{<:Base.BitInteger32}) = UInt32
uint_sup(::Type{<:Union{Int64,UInt64}}) = UInt64
uint_sup(::Type{<:Union{Int128,UInt128}}) = UInt128

#### Fast

struct SamplerRangeFast{U<:BitUnsigned,T<:BitInteger} <: Sampler{T}
    a::T      # first element of the range
    bw::UInt  # bit width
    m::U      # range length - 1
    mask::U   # mask generated values before threshold rejection
end

SamplerRangeFast(r::AbstractUnitRange{T}) where T<:BitInteger =
    SamplerRangeFast(r, uint_sup(T))

function SamplerRangeFast(r::AbstractUnitRange{T}, ::Type{U}) where {T,U}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = (last(r) - first(r)) % unsigned(T) % U # % unsigned(T) to not propagate sign bit
    bw = (sizeof(U) << 3 - leading_zeros(m)) % UInt # bit-width
    mask = ((1 % U) << bw) - (1 % U)
    SamplerRangeFast{U,T}(first(r), bw, m, mask)
end

function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt32,T}) where T
    a, bw, m, mask = sp.a, sp.bw, sp.m, sp.mask
    # below, we don't use UInt32, to get reproducible values, whether Int is Int64 or Int32
    x = rand(rng, LessThan(m, Masked(mask, UInt52Raw(UInt32))))
    (x + a % UInt32) % T
end

function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt64,T}) where T
    a, bw, m, mask = sp.a, sp.bw, sp.m, sp.mask
    x = bw <= 52 ? rand(rng, LessThan(m, Masked(mask, UInt52Raw()))) :
                   rand(rng, LessThan(m, Masked(mask, uniform(UInt64))))
    (x + a % UInt64) % T
end

function rand(rng::AbstractRNG, sp::SamplerRangeFast{UInt128,T}) where T
    a, bw, m, mask = sp.a, sp.bw, sp.m, sp.mask
    x = bw <= 52  ?
        rand(rng, LessThan(m % UInt64, Masked(mask % UInt64, UInt52Raw()))) % UInt128 :
    bw <= 104 ?
        rand(rng, LessThan(m, Masked(mask, UInt104Raw()))) :
        rand(rng, LessThan(m, Masked(mask, uniform(UInt128))))
    x % T + a
end

#### Default

# remainder function according to Knuth, where rem_knuth(a, 0) = a
rem_knuth(a::UInt, b::UInt) = a % (b + (b == 0)) + a * (b == 0)
rem_knuth(a::T, b::T) where {T<:Unsigned} = b != 0 ? a % b : a

# maximum multiple of k <= sup decremented by one,
# that is 0xFFFF...FFFF if k = (typemax(T) - typemin(T)) + 1 and sup == typemax(T) - 1
# with intentional underflow
# see http://stackoverflow.com/questions/29182036/integer-arithmetic-add-1-to-uint-max-and-divide-by-n-without-overflow

# sup == 0 means typemax(T) + 1
maxmultiple(k::T, sup::T=zero(T)) where {T<:Unsigned} =
    (div(sup - k, k + (k == 0))*k + k - one(k))::T

# similar but sup must not be equal to typemax(T)
unsafe_maxmultiple(k::T, sup::T) where {T<:Unsigned} =
    div(sup, k + (k == 0))*k - one(k)

struct SamplerRangeInt{T<:Integer,U<:Unsigned} <: Sampler{T}
    a::T      # first element of the range
    bw::Int   # bit width
    k::U      # range length or zero for full range
    u::U      # rejection threshold
end

SamplerRangeInt(r::AbstractUnitRange{T}) where T<:BitInteger =
    SamplerRangeInt(r, uint_sup(T))

function SamplerRangeInt(r::AbstractUnitRange{T}, ::Type{U}) where {T,U}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    a = first(r)
    m = (last(r) - first(r)) % unsigned(T) % U
    k = m + one(U)
    bw = (sizeof(U) << 3 - leading_zeros(m)) % Int
    mult = if U === UInt32
        maxmultiple(k)
    elseif U === UInt64
        bw <= 52 ? unsafe_maxmultiple(k, one(UInt64) << 52) :
                   maxmultiple(k)
    else # U === UInt128
        bw <= 52  ? unsafe_maxmultiple(k, one(UInt128) << 52) :
        bw <= 104 ? unsafe_maxmultiple(k, one(UInt128) << 104) :
                    maxmultiple(k)
    end

    SamplerRangeInt{T,U}(a, bw, k, mult) # overflow ok
end

Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{T},
        ::Repetition) where {T<:BitInteger} = SamplerRangeInt(r)

rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt32}) where {T<:BitInteger} =
    (unsigned(sp.a) + rem_knuth(rand(rng, LessThan(sp.u, UInt52Raw(UInt32))), sp.k)) % T

# this function uses 52 bit entropy for small ranges of length <= 2^52
function rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt64}) where T<:BitInteger
    x = sp.bw <= 52 ? rand(rng, LessThan(sp.u, UInt52())) :
                      rand(rng, LessThan(sp.u, uniform(UInt64)))
    return ((sp.a % UInt64) + rem_knuth(x, sp.k)) % T
end

function rand(rng::AbstractRNG, sp::SamplerRangeInt{T,UInt128}) where T<:BitInteger
    x = sp.bw <= 52  ? rand(rng, LessThan(sp.u, UInt52(UInt128))) :
        sp.bw <= 104 ? rand(rng, LessThan(sp.u, UInt104(UInt128))) :
                       rand(rng, LessThan(sp.u, uniform(UInt128)))
    return ((sp.a % UInt128) + rem_knuth(x, sp.k)) % T
end

### BigInt

struct SamplerBigInt <: Sampler{BigInt}
    a::BigInt         # first
    m::BigInt         # range length - 1
    nlimbs::Int       # number of limbs in generated BigInt's (z ∈ [0, m])
    nlimbsmax::Int    # max number of limbs for z+a
    mask::Limb        # applied to the highest limb
end

function Sampler(::Type{<:AbstractRNG}, r::AbstractUnitRange{BigInt}, ::Repetition)
    m = last(r) - first(r)
    m < 0 && throw(ArgumentError("range must be non-empty"))
    nd = ndigits(m, base=2)
    nlimbs, highbits = divrem(nd, 8*sizeof(Limb))
    highbits > 0 && (nlimbs += 1)
    mask = highbits == 0 ? ~zero(Limb) : one(Limb)<<highbits - one(Limb)
    nlimbsmax = max(nlimbs, abs(last(r).size), abs(first(r).size))
    return SamplerBigInt(first(r), m, nlimbs, nlimbsmax, mask)
end

function rand(rng::AbstractRNG, sp::SamplerBigInt)
    x = MPZ.realloc2(sp.nlimbsmax*8*sizeof(Limb))
    GC.@preserve x begin
        limbs = UnsafeView(x.d, sp.nlimbs)
        while true
            rand!(rng, limbs)
            limbs[end] &= sp.mask
            MPZ.mpn_cmp(x, sp.m, sp.nlimbs) <= 0 && break
        end
        # adjust x.size (normally done by mpz_limbs_finish, in GMP version >= 6)
        x.size = sp.nlimbs
        while x.size > 0
            limbs[x.size] != 0 && break
            x.size -= 1
        end
    end
    MPZ.add!(x, sp.a)
end

## random values from AbstractArray

Sampler(::Type{RNG}, r::AbstractArray, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerSimple(r, Sampler(RNG, firstindex(r):lastindex(r), n))

rand(rng::AbstractRNG, sp::SamplerSimple{<:AbstractArray,<:Sampler}) =
    @inbounds return sp[][rand(rng, sp.data)]

## random values from Dict

function Sampler(::Type{RNG}, t::Dict, ::Repetition) where RNG<:AbstractRNG
    isempty(t) && throw(ArgumentError("collection must be non-empty"))
    # we use Val(Inf) below as rand is called repeatedly internally
    # even for generating only one random value from t
    SamplerSimple(t, Sampler(RNG, LinearIndices(t.slots), Val(Inf)))
end

function rand(rng::AbstractRNG, sp::SamplerSimple{<:Dict,<:Sampler})
    while true
        i = rand(rng, sp.data)
        Base.isslotfilled(sp[], i) && @inbounds return (sp[].keys[i] => sp[].vals[i])
    end
end

## random values from Set

Sampler(::Type{RNG}, t::Set{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
    SamplerTag{Set{T}}(Sampler(RNG, t.dict, n))

rand(rng::AbstractRNG, sp::SamplerTag{<:Set,<:Sampler}) = rand(rng, sp.data).first

## random values from BitSet

function Sampler(RNG::Type{<:AbstractRNG}, t::BitSet, n::Repetition)
    isempty(t) && throw(ArgumentError("collection must be non-empty"))
    SamplerSimple(t, Sampler(RNG, minimum(t):maximum(t), Val(Inf)))
end

function rand(rng::AbstractRNG, sp::SamplerSimple{BitSet,<:Sampler})
    while true
        n = rand(rng, sp.data)
        n in sp[] && return n
    end
end

## random values from AbstractDict/AbstractSet

# we defer to _Sampler to avoid ambiguities with a call like Sampler(rng, Set(1), Val(1))
Sampler(RNG::Type{<:AbstractRNG}, t::Union{AbstractDict,AbstractSet}, n::Repetition) =
    _Sampler(RNG, t, n)

# avoid linear complexity for repeated calls
_Sampler(RNG::Type{<:AbstractRNG}, t::Union{AbstractDict,AbstractSet}, n::Val{Inf}) =
    Sampler(RNG, collect(t), n)

# when generating only one element, avoid the call to collect
_Sampler(::Type{<:AbstractRNG}, t::Union{AbstractDict,AbstractSet}, ::Val{1}) =
    SamplerTrivial(t)

function nth(iter, n::Integer)::eltype(iter)
    for (i, x) in enumerate(iter)
        i == n && return x
    end
end

rand(rng::AbstractRNG, sp::SamplerTrivial{<:Union{AbstractDict,AbstractSet}}) =
    nth(sp[], rand(rng, 1:length(sp[])))

## random characters from a string

# we use collect(str), which is most of the time more efficient than specialized methods
# (except maybe for very small arrays)
Sampler(RNG::Type{<:AbstractRNG}, str::AbstractString, n::Val{Inf}) = Sampler(RNG, collect(str), n)

# when generating only one char from a string, the specialized method below
# is usually more efficient
Sampler(RNG::Type{<:AbstractRNG}, str::AbstractString, ::Val{1}) =
    SamplerSimple(str, Sampler(RNG, 1:_lastindex(str), Val(Inf)))

isvalid_unsafe(s::String, i) = !Base.is_valid_continuation(GC.@preserve s unsafe_load(pointer(s), i))
isvalid_unsafe(s::AbstractString, i) = isvalid(s, i)
_lastindex(s::String) = sizeof(s)
_lastindex(s::AbstractString) = lastindex(s)

function rand(rng::AbstractRNG, sp::SamplerSimple{<:AbstractString,<:Sampler})::Char
    str = sp[]
    while true
        pos = rand(rng, sp.data)
        isvalid_unsafe(str, pos) && return str[pos]
    end
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


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

# Float32/Float64 based on C implementations from FDLIBM (http://www.netlib.org/fdlibm/)
# and FreeBSD:
#
## ====================================================
## Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
##
## Developed at SunPro, a Sun Microsystems, Inc. business.
## Permission to use, copy, modify, and distribute this
## software is freely granted, provided that this notice
## is preserved.
## ====================================================
## Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
## Debugged and optimized by Bruce D. Evans.

"""
    cbrt(x::Real)

Return the cube root of `x`, i.e. ``x^{1/3}``. Negative values are accepted
(returning the negative real root when ``x < 0``).

The prefix operator `∛` is equivalent to `cbrt`.

# Examples
```jldoctest
julia> cbrt(big(27))
3.0

julia> cbrt(big(-27))
-3.0
```
"""
cbrt(x::Real) = cbrt(float(x))
cbrt(x::AbstractFloat) = x < 0 ? -(-x)^(1//3) : x^(1//3)

"""
    _approx_cbrt(x)

Approximate `cbrt` to 5 bits precision

    cbrt(2^e * (1+m)) ≈ 2^(e÷3) * (1 + (e%3+m)÷3)

where:
 - `e` is integral and >= 0
 - `m` is real and in [0, 1),
 - `÷` is integer division
 - `%` is integer remainder

The RHS is always >= the LHS and has a maximum relative error of about 1 in 16.
Adding a bias of -0.03306235651 to the `(e%3+m)÷3` term reduces the error to about 1 in
32.

With the IEEE floating point representation, for finite positive normal values, ordinary
integer division of the value in bits magically gives almost exactly the RHS of the above
provided we first subtract the exponent bias and later add it back.  We do the
subtraction virtually to keep e >= 0 so that ordinary integer division rounds towards
minus infinity; this is also efficient. All operations can be done in 32-bit.

These implementations assume that NaNs, infinities and zeros have already been filtered.
"""
@inline function _approx_cbrt(x::T) where {T<:Union{Float32,Float64}}
    # floor(UInt32, adj * exp2(k)) should be evaluated to 2 constants.
    adj = exponent_bias(T)*2/3 - 0.03306235651
    k = significand_bits(T) - (8*sizeof(T) - 32)

    u = highword(x) & 0x7fff_ffff
    if u >= Base.Math.highword(floatmin(T))
        v = div(u, UInt32(3)) + floor(UInt32, adj * exp2(k))
    else
        # subnormal
        x *= maxintfloat(T)
        adj -= exponent(maxintfloat(T))/3
        u = highword(x) & 0x7fff_ffff
        v = div(u, UInt32(3)) + floor(UInt32, adj * exp2(k))
    end
    return copysign(fromhighword(T, v), x)
end

@inline function _improve_cbrt(x::Float32, t::Float32)
    # Newton iterations solving
    #   t^2 - x/t == 0
    # with update
    #   t <- t*(t^3 + 2*x)/(2*t^3 + x)

    # Use double precision so that its terms can be arranged for efficiency
    # without causing overflow or underflow.
    xx = Float64(x)
    tt = Float64(t)

    # 1st step: 16 bits accuracy
    tt3 = tt^3
    tt *= (2*xx + tt3)/(x + 2*tt3)

    # 2nd step: 47 bits accuracy
    tt3 = tt^3
    tt *= (2*xx + tt3)/(x + 2*tt3)

    return Float32(tt)
end

@inline function _improve_cbrt(x::Float64, t::Float64)
    # cbrt to 23 bits:
    #
    #    cbrt(x) = t * cbrt(x / t^3) ~= t * P(t^3 / x)
    #
    # where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
    # to within 2^-23.5 when |r - 1| < 1/10.  The rough approximation
    # has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
    # gives us bounds for r = t^3/x.

    r = (t*t)*(t/x)
    t *= (@horner(r, 1.87595182427177009643, -1.88497979543377169875, 1.621429720105354466140) +
          r^3 * @horner(r, -0.758397934778766047437, 0.145996192886612446982))

    # Round t away from zero to 23 bits (sloppily except for ensuring that
    # the result is larger in magnitude than cbrt(x) but not much more than
    # 2 23-bit ulps larger).  With rounding towards zero, the error bound
    # would be ~5/6 instead of ~4/6.  With a maximum error of 2 23-bit ulps
    # in the rounded t, the infinite-precision error in the Newton
    # approximation barely affects third digit in the final error
    # 0.667; the error in the rounded t can be up to about 3 23-bit ulps
    # before the final error is larger than 0.667 ulps.

    u = reinterpret(UInt64, t)
    u = (u + 0x8000_0000) & UInt64(0xffff_ffff_c000_0000)
    t = reinterpret(Float64, u)

    # one step Newton iteration solving
    #   t^3 - x == 0
    # with update
    #   t <- t + t * (x/t^2 - t) / (3*t)

    # to 53 bits with error < 0.667 ulps
    s = t*t             # t*t is exact
    r = x/s             # error <= 0.5 ulps; |r| < |t|
    w = t+t             # t+t is exact
    r = (r - t)/(w + r) # r-t is exact; w+r ~= 3*t
    t = muladd(t, r, t) # error <= 0.5 + 0.5/3 + epsilon
    return t
end

function cbrt(x::Union{Float32,Float64})
    if !isfinite(x) || iszero(x)
        return x
    end
    t = _approx_cbrt(x)
    return _improve_cbrt(x, t)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia