text:
struct IntSetDomain <: AbstractDomain
    values          ::Array{Int}
    indexes         ::Array{Int}
    min             ::Int
    max             ::Int
    requiring_index ::StateObject{Int}
    excluding_index ::StateObject{Int}
    trailer         ::Trailer
end

function IntSetDomain(trailer::Trailer, min::Int, max::Int)
    n = max - min + 1
    requiring_index = SeaPearl.StateObject{Int}(1, trailer)
    excluding_index = SeaPearl.StateObject{Int}(n, trailer)

    values = collect(1:n)
    indexes = collect(1:n)

    return IntSetDomain(values, indexes, min, max, requiring_index, excluding_index, trailer)
end

"""
    reset_domain!(dom::IntSetDomain)

Used in `reset_model!`.
"""
function reset_domain!(dom::IntSetDomain)
    setValue!(dom.requiring_index, 1)
    setValue!(dom.excluding_index, dom.max - dom.min + 1)
    sort!(dom.values)
    sort!(dom.indexes)
    dom
end

function exclude!(dom::IntSetDomain, v::Int)
    @assert dom.min <= v <= dom.max

    if !is_possible(dom, v)
        return nothing
    end
    if is_required(dom, v)
        throw(ErrorException("You can't exclude $(v) as it is required in $(dom)"))
    end
    exchangePositions!(dom, v, dom.values[dom.excluding_index.value] + dom.min - 1)
    setValue!(dom.excluding_index, dom.excluding_index.value - 1)
end

function require!(dom::IntSetDomain, v::Int)
    @assert dom.min <= v <= dom.max

    if is_required(dom, v)
        return nothing
    end
    if !is_possible(dom, v)
        throw(ErrorException("You can't require $(v) as it is excluded from $(dom)"))
    end
    exchangePositions!(dom, v, dom.values[dom.requiring_index.value] + dom.min - 1)
    setValue!(dom.requiring_index, dom.requiring_index.value + 1)
end

"""
    function exclude_all!(dom::IntSetDomain)

Exclude all possible not yet required values from the set domain.
"""
function exclude_all!(dom::IntSetDomain)
    excluded = dom.values[dom.requiring_index.value:dom.excluding_index.value] .+ (dom.min -1)
    setValue!(dom.excluding_index, dom.requiring_index.value - 1)
    return excluded
end

"""
    function require_all!(dom::IntSetDomain)

Require all possible not yet excluded values from the set domain.
"""
function require_all!(dom::IntSetDomain)
    required = dom.values[dom.requiring_index.value:dom.excluding_index.value] .+ (dom.min -1)
    setValue!(dom.requiring_index, dom.excluding_index.value + 1)
    return required
end

"""
    function is_possible(dom::IntSetDomain, v::Int)

Return true if `v` is a possible value for `dom`.
"""
function is_possible(dom::IntSetDomain, v::Int)
    if dom.min > v || v > dom.max
        return false
    end
    return dom.indexes[v - dom.min + 1] <= dom.excluding_index.value
end

"""
    function is_possible(dom::IntSetDomain, v::Int)

Return true if `v` is a required value for `dom`.
"""
function is_required(dom::IntSetDomain, v::Int)
    if dom.min > v || v > dom.max
        return false
    end
    return dom.indexes[v - dom.min + 1] < dom.requiring_index.value
end

"""
    function required_values(dom::IntSetDomain)::Set{Int}

Return a set containing all the required values of `dom`.
"""
function required_values(dom::IntSetDomain)
    return Set{Int}(dom.values[1:(dom.requiring_index.value - 1)] .+ (dom.min - 1))
end

"""
    function possible_not_required_values(dom::IntSetDomain)::Set{Int}

Return a set containing the possible values, but not the required ones, for `dom`.
"""
function possible_not_required_values(dom::IntSetDomain)
    return Set{Int}(dom.values[dom.requiring_index.value:dom.excluding_index.value] .+ (dom.min - 1))
end

"""
    exchangePositions!(dom::IntDomain, v1::Int, v2::Int)

Intended for internal use only, exchange the position of `v1` and `v2` in the array of the domain.
"""
function exchangePositions!(dom::IntSetDomain, v1::Int, v2::Int)
    @assert dom.min <= v1 <= dom.max
    @assert dom.min <= v2 <= dom.max

    v1, v2 = v1 - dom.min + 1, v2 - dom.min + 1

    i1, i2 = dom.indexes[v1], dom.indexes[v2]

    dom.values[i1] = v2
    dom.values[i2] = v1
    dom.indexes[v1] = i2
    dom.indexes[v2] = i1

    return dom
end

function Base.show(io::IO, dom::IntSetDomain)
    print(io, "required_values = ", required_values(dom), "; possibles_values = ", possible_not_required_values(dom))
end

function Base.show(io::IO, ::MIME"text/plain", dom::IntSetDomain)
    println(io, typeof(dom), ", min = ", dom.min, ", max = ", dom.max)
    println(io, "   required_values = ", required_values(dom))
    print(io, "   possible_values = ", possible_not_required_values(dom))
end


text:
## path to TRPL folder --> set it before running the code!
path_to_trpl = "D:\\Projects\\TRPL" # some path to the TRPL folder (just an example here)

## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## loading packages
using DataFrames, CSV, XLSX, LinearAlgebra,
      Distributions, Measurements, ForwardDiff, NLopt

## directory
cd(joinpath(path_to_trpl,"Data\\Real"))

## output results? (true or false)
out = false

## name of measurement
name_measurement = "sample1"

## load measurement data
data_measurement = Matrix(CSV.read(name_measurement*".csv",DataFrame,skipto=2))

## function for data processing (data structre is explained in README.md)
function data_processing(data)
    t = data[4:end,1]
    idx_end = sum(t .< 200)
    num_traj = size(data,2)-1
    dims = 9 + num_traj
    P = round.(data[1,2:end],sigdigits=3)
    ND = 10 .^ data[2,2:end]
    background = data[3,2:end]
    IPL = data[4:end,2:num_traj+1]
    tspan = (t[1], t[idx_end])
    N = length(t)
    idx_IPL = zeros(Bool,size(IPL))
    for i = 1:num_traj
	    idx_IPL[:,i] = IPL[:,i] .> 0
    end
    C_scaling = maximum(IPL)
    return t[1:idx_end], tspan, num_traj, dims, IPL[1:idx_end,:], idx_IPL[1:idx_end,:], P, ND, background, N, C_scaling
end

## process measurement data
t, tspan, num_traj, dims, IPL, idx_IPL, P, ND, background, N, C_scaling = data_processing(data_measurement)
IPL = IPL[:,1] # we use the transient at the lowest laser power for the mono-exponential fit
background = background[1]

## mono-exponential model
model(p) = @. p[1] * exp(-t/p[2]) + background

## mle loss function
function loss(par)
    sim = model(par)
    idx_used = (sim .> 0) .& (IPL .> 0)
    l = sum(sim[idx_used] .- IPL[idx_used] .* log.(sim[idx_used]))
    return l
end

## definition nlopt loss function
function loss_nlopt(par::Vector,grad::Vector)
    if length(grad) > 0
	    grad[:] = ForwardDiff.gradient(loss,par)
    end
    println("l = $(loss(par))")
    return loss(par)
end

## LBFGS (local optmization) via NLopt
opt = Opt(:LD_LBFGS,2)
opt.min_objective = loss_nlopt
opt.ftol_abs = 1e-14
opt.maxtime = 600
(optf,p_opt,ret) = optimize(opt,[1000, 10])

## calculate standard errors
H_opt = ForwardDiff.hessian(loss,p_opt)
SE_opt = sqrt.(diag(inv(H_opt)))

## print results
println("\nτ_mono = $(round(p_opt[2],sigdigits=3)) ± $(round(SE_opt[2],sigdigits=3)) ns")

## confidence level (1-α)
alpha = 0.05 # corresponds to a 95% confidence level
q = quantile(Normal(),1-alpha/2)

## get fit parameters with standard erros
cd(string(path_to_trpl,"\\Results\\Real"))
files = filter(x->endswith(x, ".xlsx"), readdir())
for i = 1:length(files)
    if occursin(Regex("$name_measurement"),files[i])
        global idx = i
    end
end
xf = XLSX.readxlsx(files[idx])
sh = xf["Results"]
p_fit = Float64.(sh["C3:K3"])
CIHW_fit = Float64.(sh["C4:K4"])
p_fit = p_fit .± (CIHW_fit ./ q)

## compair to effective lifetime
tau_eff = round(sum([1/p_fit[i] for i = 1:4]) .^ -1,sigdigits=3)
println("\nτ_eff  = $tau_eff ns")

## write output
if out
    file_name = "mono_exp_fit_$(name_measurement).txt"
    if isfile(file_name)
        rm(file_name)
    end
    open(file_name,"w") do io
        println(io,"\nτ_mono = $(round(p_opt[2],sigdigits=3)) ± $(round(SE_opt[2],sigdigits=3)) ns")
        println(io,"\nτ_eff  = $tau_eff ns")
    end
end

ext: jl

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
using Q.K_KeyTable
@testset "table tests" begin
  @test_throws ArgumentError K_KeyTable(kj(0))
  @test_throws ArgumentError K_KeyTable(xD(ktn(0,0),ktn(0,0)))
  let x = K_Table(a=Int[]), p = kpointer(x)
    @test_throws ArgumentError K_KeyTable(xD(r1(p),ktn(0,0)))
    @test (kt = K_KeyTable(xD(r1(p), r1(p))); xt(kpointer(kt)) == XD)
  end
  @test (x = K_KeyTable(1; a=Int[], b=Int[]); size(x) == (0, 2))
  @test begin
    t = K_KeyTable(1; a=[0, 1, 2], b=[0., 10., 20.])
    show_to_string(t) == """
    3×2 Q.K_KeyTable
    │ Row │ a │ b    │
    ├─────┼───┼──────┤
    │ 1   │ 0 │ 0.0  │
    │ 2   │ 1 │ 10.0 │
    │ 3   │ 2 │ 20.0 │"""
  end
  @test (x = K_KeyTable(1;a=Int[],b=Int[]); x[:a] == x[:b] == Int[])
end


text:
using Random
using GroupsCore

import AbstractAlgebra: AbstractPermutationGroup, AbstractPerm, degree

# to be moved to AA:
AbstractAlgebra.degree(G::AbstractAlgebra.Generic.SymmetricGroup) = G.n

"""
    WreathProduct(G::Group, P::AbstractPermutationGroup) <: Group
Return the wreath product of a group `G` by permutation group `P`, usually
written as `G ≀ P`.

As set `G ≀ P` is the same as `Gᵈ × P` and the group can be understood as a
semi-direct product of `P` acting on `d`-fold cartesian product of `G` by
permuting coordinates. To be more precise, the multiplication inside wreath
product is defined as
>  `(n, σ) * (m, τ) = (n*(m^σ), σ*τ)`
where `m^σ` denotes the action (from the right) of the permutation `σ` on
`d`-tuples of elements from `G`.
"""
struct WreathProduct{DP<:DirectPower,PGr<:AbstractPermutationGroup} <:
       GroupsCore.Group
    N::DP
    P::PGr

    function WreathProduct(G::Group, P::AbstractPermutationGroup)
        N = DirectPower{degree(P)}(G)
        return new{typeof(N),typeof(P)}(N, P)
    end
end

struct WreathProductElement{
    DPEl<:DirectPowerElement,
    PEl<:AbstractPerm,
    Wr<:WreathProduct,
} <: GroupsCore.GroupElement
    n::DPEl
    p::PEl
    parent::Wr

    function WreathProductElement(
        n::DirectPowerElement,
        p::AbstractPerm,
        W::WreathProduct,
    )
        new{typeof(n),typeof(p),typeof(W)}(n, p, W)
    end
end

Base.one(W::WreathProduct) = WreathProductElement(one(W.N), one(W.P), W)

Base.eltype(::Type{<:WreathProduct{DP,PGr}}) where {DP,PGr} =
    WreathProductElement{eltype(DP),eltype(PGr),WreathProduct{DP,PGr}}

function Base.iterate(G::WreathProduct)
    itr = Iterators.product(G.N, G.P)
    res = iterate(itr)
    @assert res !== nothing
    elt = WreathProductElement(first(res)..., G)
    return elt, (iterator = itr, state = last(res))
end

function Base.iterate(G::WreathProduct, state)
    itr, st = state.iterator, state.state
    res = iterate(itr, st)
    res === nothing && return nothing
    elt = WreathProductElement(first(res)..., G)
    return elt, (iterator = itr, state = last(res))
end

function Base.IteratorSize(::Type{<:WreathProduct{DP,PGr}}) where {DP,PGr}
    dpI = Base.IteratorSize(DP)
    pgI = Base.IteratorSize(PGr)

    if dpI isa Base.IsInfinite || pgI isa Base.IsInfinite
        return Base.IsInfinite()
    elseif dpI isa Base.SizeUnknown || pgI isa Base.SizeUnknown
        return Base.SizeUnknown()
    else
        return Base.HasShape{2}()
    end
end

Base.size(G::WreathProduct) = (length(G.N), length(G.P))

GroupsCore.order(::Type{I}, G::WreathProduct) where {I<:Integer} =
    convert(I, order(I, G.N) * order(I, G.P))

function GroupsCore.gens(G::WreathProduct)
    N_gens = [WreathProductElement(n, one(G.P), G) for n in gens(G.N)]
    P_gens = [WreathProductElement(one(G.N), p, G) for p in gens(G.P)]
    return [N_gens; P_gens]
end

Base.isfinite(G::WreathProduct) = isfinite(G.N) && isfinite(G.P)

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{<:WreathProduct},
)

    G = rs[]
    return WreathProductElement(rand(rng, G.N), rand(rng, G.P), G)
end

GroupsCore.parent(g::WreathProductElement) = g.parent
GroupsCore.parent_type(::Type{WreathProductElement{S,T,Wr}}) where {S,T,Wr} = Wr

Base.:(==)(g::WreathProductElement, h::WreathProductElement) =
    parent(g) === parent(h) && g.n == h.n && g.p == h.p

Base.hash(g::WreathProductElement, h::UInt) =
    hash(g.n, hash(g.p, hash(g.parent, h)))

function Base.deepcopy_internal(g::WreathProductElement, stackdict::IdDict)
    return WreathProductElement(
        Base.deepcopy_internal(g.n, stackdict),
        Base.deepcopy_internal(g.p, stackdict),
        parent(g),
    )
end

_act(p::AbstractAlgebra.AbstractPerm, n::DirectPowerElement) =
    DirectPowerElement(n.elts[p.d], parent(n))

function Base.inv(g::WreathProductElement)
    pinv = inv(g.p)
    return WreathProductElement(_act(pinv, inv(g.n)), pinv, parent(g))
end

function Base.:(*)(g::WreathProductElement, h::WreathProductElement)
    @assert parent(g) === parent(h)
    return WreathProductElement(g.n * _act(g.p, h.n), g.p * h.p, parent(g))
end

Base.isone(g::WreathProductElement) = isone(g.n) && isone(g.p)

Base.show(io::IO, G::WreathProduct) =
    print(io, "Wreath Product of $(G.N.group) by $(G.P)")
Base.show(io::IO, g::WreathProductElement) = print(io, "( $(g.n)≀$(g.p) )")


text:
using StatsBase
using DataFrames
using CSV

#trainingSizeInput = parse(Int64, ARGS[1])
trainingSize = 12

#path = "/zhome/9f/d/88706/SpecialeCode/Thesis/ML/Lasso_Test"
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/ML"
cd("$(homedir())/Documents/GitHub/Thesis/Data")
path = "$(homedir())/Documents/GitHub/Thesis/Data"
@everywhere include("ParallelModelGeneration.jl")
include("SupportFunction.jl")
include("DataLoad.jl")
println("Leeeeroooy Jenkins")

possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

modelMatrix = zeros(industriesTotal, possibilities)
testModel = [0 1 0 0 0]
for i=2:industriesTotal
    modelMatrix[i, :] = testModel
end
##START OF A METHOD

path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexDataDiff/"

XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

riskAversions = linspace(0, 2.4, 10)
riskAversions = linspace(0, 4, 16)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

nRows = size(XArrays[1])[1]
w1N = repeat([0.1], outer = 10) #1/N weights
return1NMatrix = zeros(nRows-trainingSize)
returnSAAMatrix = zeros(nRows-trainingSize)

weightsSAA = zeros(nRows-trainingSize, 10)
forecastRows = zeros(nRows-trainingSize, 10)

startPoint = 241 #194608
endPoint = 1080 #201607

g=1
fileName = "Results"
gammaRisk = riskAversions[g] #riskAversion in MV optimization
total = (endPoint-trainingSize)

t=1
println("time $t / $total, gammaRisk $g / 10 ")
trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)
expectedReturns = zeros(industriesTotal)
for i=1:10
    expectedReturns[i] = mean(trainingYArrays[i][:])
end

Sigma =  cov(trainingXArrays[1][:,1:10])

a = w1N'*Sigma*w1N*4
b = w1N'*expectedReturns
a/b
#A=U^(T)U where U is upper triangular with real positive diagonal entries
F = lufact(Sigma)

U = F[:U]  #Cholesky factorization of Sigma

#getting the actual Y values for each industry
valY = zeros(10)
for i = 1:10
    valY[i] = validationY[i][1]
end
indexes = 10
M = JuMP.Model(solver = GurobiSolver(OutputFlag = 0))
@variables M begin
        w[1:indexes]
        u[1:indexes]
        z
        y
end
gammaRisk=2.4
@objective(M,Min, gammaRisk*y - expectedReturns'*w)
@constraint(M, 0 .<= w)
@constraint(M, sum(w[i] for i=1:indexes) == 1)
@constraint(M, norm([2*U'*w;y-1]) <= y+1)
solve(M)
obj24 = getobjectivevalue(M)
wStar0 = getvalue(w)
wStar24 = getvalue(w)

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
facts("A simple square") do
    m = Model()
    @variable(m, 0 ≤ x[1:2] ≤ 1)

    v, r = get_extrema(m)

    @fact length(v) --> 4
    @fact length(r) --> 0
    @fact [0,0] in v --> true
    @fact [1,0] in v --> true
    @fact [0,1] in v --> true
    @fact [1,1] in v --> true
end


text:
# Convert the buffer forward rates from  1/(Molar sec) to 1/(molecule sec)
# Arguments:  f = Float64, the forward rate in 1/(Molar sec)
#             v = the volume in metres cubed
# Returns: Float64, the forward rate in per molecule per second

function perMolarToPerMolecule(f, v)
  N_av = 6.02214085774e+23
  v_litres = 1000 * v
  return f/(N_av * v_litres) # factor of 1000 to convert from litres to m^3
end


text:
using Turing, Random
using Turing: ParticleContainer, getweights, getweight,
    effectiveSampleSize, Trace, current_trace, VarName, VarInfo,
    Sampler, consume, produce, fork, getlogp
using Turing.Core: logZ, reset_logweights!, increase_logweight!,
    resample_propagate!, reweight!
using Test

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "container.jl" begin
    @turing_testset "copy particle container" begin
        pc = ParticleContainer(Trace[])
        newpc = copy(pc)

        @test newpc.logWs == pc.logWs
        @test typeof(pc) === typeof(newpc)
    end

    @turing_testset "particle container" begin
        # Create a resumable function that always yields `logp`.
        function fpc(logp)
            f = let logp = logp
                () -> begin
                    while true
                        produce(logp)
                    end
                end
            end
            return f
        end

        # Dummy sampler that is not actually used.
        sampler = Sampler(PG(5), empty_model())

        # Create particle container.
        logps = [0.0, -1.0, -2.0]
        particles = [Trace(fpc(logp), empty_model(), sampler, VarInfo()) for logp in logps]
        pc = ParticleContainer(particles)

        # Initial state.
        @test all(iszero(getlogp(particle.vi)) for particle in pc.vals)
        @test pc.logWs == zeros(3)
        @test getweights(pc) == fill(1/3, 3)
        @test all(getweight(pc, i) == 1/3 for i in 1:3)
        @test logZ(pc) ≈ log(3)
        @test effectiveSampleSize(pc) == 3

        # Reweight particles.
        reweight!(pc)
        @test all(iszero(getlogp(particle.vi)) for particle in pc.vals)
        @test pc.logWs == logps
        @test getweights(pc) ≈ exp.(logps) ./ sum(exp, logps)
        @test all(getweight(pc, i) ≈ exp(logps[i]) / sum(exp, logps) for i in 1:3)
        @test logZ(pc) == log(sum(exp, logps))

        # Reweight particles.
        reweight!(pc)
        @test all(iszero(getlogp(particle.vi)) for particle in pc.vals)
        @test pc.logWs == 2 .* logps
        @test getweights(pc) == exp.(2 .* logps) ./ sum(exp, 2 .* logps)
        @test all(getweight(pc, i) ≈ exp(2 * logps[i]) / sum(exp, 2 .* logps) for i in 1:3)
        @test logZ(pc) == log(sum(exp, 2 .* logps))

        # Resample and propagate particles.
        resample_propagate!(pc)
        @test all(iszero(getlogp(particle.vi)) for particle in pc.vals)
        @test pc.logWs == zeros(3)
        @test getweights(pc) == fill(1/3, 3)
        @test all(getweight(pc, i) == 1/3 for i in 1:3)
        @test logZ(pc) ≈ log(3)
        @test effectiveSampleSize(pc) == 3

        # Reweight particles.
        reweight!(pc)
        @test all(iszero(getlogp(particle.vi)) for particle in pc.vals)
        @test pc.logWs ⊆ logps
        @test getweights(pc) == exp.(pc.logWs) ./ sum(exp, pc.logWs)
        @test all(getweight(pc, i) ≈ exp(pc.logWs[i]) / sum(exp, pc.logWs) for i in 1:3)
        @test logZ(pc) == log(sum(exp, pc.logWs))

        # Increase unnormalized logarithmic weights.
        logws = copy(pc.logWs)
        increase_logweight!(pc, 2, 1.41)
        @test pc.logWs == logws + [0, 1.41, 0]

        # Reset unnormalized logarithmic weights.
        logws = pc.logWs
        reset_logweights!(pc)
        @test pc.logWs === logws
        @test all(iszero, pc.logWs)
    end

    @turing_testset "trace" begin
        n = Ref(0)

        alg = PG(5)
        spl = Sampler(alg, empty_model())
        dist = Normal(0, 1)
        function f2()
            t = TArray(Int, 1);
            t[1] = 0;
            while true
                ct = current_trace()
                vn = @varname x[n]
                Turing.assume(Random.GLOBAL_RNG, spl, dist, vn, ct.vi)
                n[] += 1
                produce(t[1])
                vn = @varname x[n]
                Turing.assume(Random.GLOBAL_RNG, spl, dist, vn, ct.vi)
                n[] += 1
                t[1] = 1 + t[1]
            end
        end

        # Test task copy version of trace
        tr = Trace(f2, empty_model(), spl, VarInfo())

        consume(tr); consume(tr)
        a = fork(tr);
        consume(a); consume(a)

        @test consume(tr) == 2
        @test consume(a) == 4
    end
end


text:
# returns Sl(n,x) using the naive series expansion
function sl_series(n::Integer, x::Float64)
    kmax = ceil(Int64, eps(Float64)^(-inv(n)))

    if iseven(n)
        co = cos(x)
        co2 = one(x) # cos((n-2)*x)
        co1 = co     # cos((n-1)*x)
        sum = co
        for k in 2:kmax
            con = 2*co*co1 - co2 # cos(n*x)
            co2 = co1
            co1 = con
            sum += con/Float64(k)^n
        end
        sum
    else
        si, co = sincos(x)
        si2 = zero(x) # sin((n-2)*x)
        si1 = si      # sin((n-1)*x)
        sum = si
        for k in 2:kmax
            si = 2*co*si1 - si2 # sin(n*x)
            si2 = si1
            si1 = si
            sum += si/Float64(k)^n
        end
        sum
    end
end

"""
    sl(n::Integer, x::Float64)::Float64

Returns the value of the Glaisher-Clausen function
``\\operatorname{Sl}_n(x)`` for integers ``n > 0`` and a real angle
``x`` of type `Float64`.  This function is defined as

```math
\\begin{aligned}
\\operatorname{Sl}_n(x) &= \\Re[\\operatorname{Li}_n(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\cos(kx)}{k^n}, \\qquad \\text{for}~n~\\text{even}, \\\\
\\operatorname{Sl}_n(x) &= \\Im[\\operatorname{Li}_n(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\sin(kx)}{k^n}, \\qquad \\text{for}~n~\\text{odd}.
\\end{aligned}
```

Note: We set ``\\operatorname{Sl}_1(0) = 0`` for consistency with the
series expansion.

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> sl(10, 1.0)
0.5398785706335891
```
"""
function sl(n::Integer, x::Float64)::Float64
    n < 1 && throw(DomainError(n, "sl(n,x) undefined for n < 1"))

    (x, sgn) = range_reduce(n + 1, x)

    n ==  1 && x == zero(Float64) && return zero(Float64)
    n ==  1 && return sgn*(pi/2 - 1/2*x)
    n ==  2 && return sgn*(pi^2/6 + (-1/2*pi + 1/4*x)*x)
    n ==  3 && return sgn*(x*(pi^2/6 + (-1/4*pi + 1/12*x)*x))

    x2 = x*x

    n ==  4 && return sgn*(pi^4/90 + x2*(-1/12*pi^2 + (pi/12 - 1/48*x)*x))
    n ==  5 && return sgn*(x*(pi^4/90 + x2*(-1/36*pi^2 + (pi/48 - 1/240*x)*x)))
    n ==  6 && return sgn*(pi^6/945 + x2*(-1/180*pi^4 + x2*(pi^2/144 + (-1/240*pi + 1/1440*x)*x)))
    n ==  7 && return sgn*(x*(pi^6/945 + x2*(-1/540*pi^4 + x2*(pi^2/720 + (-1/1440*pi + 1/10080*x)*x))))
    n ==  8 && return sgn*(pi^8/9450 + x2*(-1/1890*pi^6 + x2*(pi^4/2160 + x2*(-1/4320*pi^2 + (pi/10080 - 1/80640*x)*x))))
    n ==  9 && return sgn*(x*(pi^8/9450 + x2*(-1/5670*pi^6 + x2*(pi^4/10800 + x2*(-1/30240*pi^2 + (pi/80640 - 1/725760*x)*x)))))
    n == 10 && return sgn*(pi^10/93555 + x2*(-1/18900*pi^8 + x2*(pi^6/22680 + x2*(-1/64800*pi^4 + x2*(pi^2/241920 + (-1/725760*pi + 1/7257600*x)*x)))))
    n == 11 && return sgn*(x*(pi^10/93555 + x2*(-1/56700*pi^8 + x2*(pi^6/113400 + x2*(-1/453600*pi^4 + x2*(pi^2/2177280 + (-1/7257600*pi + 1/79833600*x)*x))))))
    n == 12 && return sgn*((691*pi^12)/638512875 + x2*(-1/187110*pi^10 + x2*(pi^8/226800 + x2*(-1/680400*pi^6 + x2*(pi^4/3628800 + x2*(-1/21772800*pi^2 + (pi/79833600 - 1/958003200*x)*x))))))
    n == 13 && return sgn*(x*((691*pi^12)/638512875 + x2*(-1/561330*pi^10 + x2*(pi^8/1134000 + x2*(-1/4762800*pi^6 + x2*(pi^4/32659200 + x2*(-1/239500800*pi^2 + (pi/958003200 - 1/12454041600*x)*x)))))))
    n == 14 && return sgn*((2*pi^14)/18243225 + x2*((-691*pi^12)/1277025750 + x2*(pi^10/2245320 + x2*(-1/6804000*pi^8 + x2*(pi^6/38102400 + x2*(-1/326592000*pi^4 + x2*(pi^2/2874009600 + (-1/12454041600*pi + 1/174356582400*x)*x)))))))
    n == 15 && return sgn*(x*((2*pi^14)/18243225 + x2*((-691*pi^12)/3831077250 + x2*(pi^10/11226600 + x2*(-1/47628000*pi^8 + x2*(pi^6/342921600 + x2*(-1/3592512000*pi^4 + x2*(pi^2/37362124800 + (-1/174356582400*pi + 1/2615348736000*x)*x))))))))
    n == 16 && return sgn*((3617*pi^16)/325641566250 + x2*(-1/18243225*pi^14 + x2*((691*pi^12)/15324309000 + x2*(-1/67359600*pi^10 + x2*(pi^8/381024000 + x2*(-1/3429216000*pi^6 + x2*(pi^4/43110144000 + x2*(-1/523069747200*pi^2 + (pi/2615348736000 - 1/41845579776000*x)*x))))))))
    n == 17 && return sgn*(x*((3617*pi^16)/325641566250 + x2*(-1/54729675*pi^14 + x2*((691*pi^12)/76621545000 + x2*(-1/471517200*pi^10 + x2*(pi^8/3429216000 + x2*(-1/37721376000*pi^6 + x2*(pi^4/560431872000 + x2*(-1/7846046208000*pi^2 + (pi/41845579776000 - 1/711374856192000*x)*x)))))))))
    n == 18 && return sgn*((43867*pi^18)/38979295480125 + x2*((-3617*pi^16)/651283132500 + x2*(pi^14/218918700 + x2*((-691*pi^12)/459729270000 + x2*(pi^10/3772137600 + x2*(-1/34292160000*pi^8 + x2*(pi^6/452656512000 + x2*(-1/7846046208000*pi^4 + x2*(pi^2/125536739328000 + (-1/711374856192000*pi + 1/12804747411456000*x)*x)))))))))
    n == 19 && return sgn*(x*((43867*pi^18)/38979295480125 + x2*((-3617*pi^16)/1953849397500 + x2*(pi^14/1094593500 + x2*((-691*pi^12)/3218104890000 + x2*(pi^10/33949238400 + x2*(-1/377213760000*pi^8 + x2*(pi^6/5884534656000 + x2*(-1/117690693120000*pi^4 + x2*(pi^2/2134124568576000 + (-1/12804747411456000*pi + 1/243290200817664000*x)*x))))))))))
    n == 20 && return sgn*((174611*pi^20)/1531329465290625 + x2*((-43867*pi^18)/77958590960250 + x2*((3617*pi^16)/7815397590000 + x2*(-1/6567561000*pi^14 + x2*((691*pi^12)/25744839120000 + x2*(-1/339492384000*pi^10 + x2*(pi^8/4526565120000 + x2*(-1/82383485184000*pi^6 + x2*(pi^4/1883051089920000 + x2*(-1/38414242234368000*pi^2 + (pi/243290200817664000 - 1/4865804016353280000*x)*x))))))))))
    n == 21 && return sgn*(x*((174611*pi^20)/1531329465290625 + x2*((-43867*pi^18)/233875772880750 + x2*((3617*pi^16)/39076987950000 + x2*(-1/45972927000*pi^14 + x2*((691*pi^12)/231703552080000 + x2*(-1/3734416224000*pi^10 + x2*(pi^8/58845346560000 + x2*(-1/1235752277760000*pi^6 + x2*(pi^4/32011868528640000 + x2*(-1/729870602452992000*pi^2 + (pi/4865804016353280000 - 1/102181884343418880000*x)*x)))))))))))
    n == 22 && return sgn*((155366*pi^22)/13447856940643125 + x2*((-174611*pi^20)/3062658930581250 + x2*((43867*pi^18)/935503091523000 + x2*((-3617*pi^16)/234461927700000 + x2*(pi^14/367783416000 + x2*((-691*pi^12)/2317035520800000 + x2*(pi^10/44812994688000 + x2*(-1/823834851840000*pi^8 + x2*(pi^6/19772036444160000 + x2*(-1/576213633515520000*pi^4 + x2*(pi^2/14597412049059840000 + (-1/102181884343418880000*pi + 1/2248001455555215360000*x)*x)))))))))))
    n == 23 && return sgn*(x*((155366*pi^22)/13447856940643125 + x2*((-174611*pi^20)/9187976791743750 + x2*((43867*pi^18)/4677515457615000 + x2*((-3617*pi^16)/1641233493900000 + x2*(pi^14/3310050744000 + x2*((-691*pi^12)/25487390728800000 + x2*(pi^10/582568930944000 + x2*(-1/12357522777600000*pi^8 + x2*(pi^6/336124619550720000 + x2*(-1/10948059036794880000*pi^4 + x2*(pi^2/306545653030256640000 + (-1/2248001455555215360000*pi + 1/51704033477769953280000*x)*x))))))))))))
    n == 24 && return sgn*((236364091*pi^24)/201919571963756521875 + x2*((-77683*pi^22)/13447856940643125 + x2*((174611*pi^20)/36751907166975000 + x2*((-43867*pi^18)/28065092745690000 + x2*((3617*pi^16)/13129867951200000 + x2*(-1/33100507440000*pi^14 + x2*((691*pi^12)/305848688745600000 + x2*(-1/8155965033216000*pi^10 + x2*(pi^8/197720364441600000 + x2*(-1/6050243151912960000*pi^6 + x2*(pi^4/218961180735897600000 + x2*(-1/6744004366665646080000*pi^2 + (pi/51704033477769953280000 - 1/1240896803466478878720000*x)*x))))))))))))
    n == 25 && return sgn*(x*((236364091*pi^24)/201919571963756521875 + x2*((-77683*pi^22)/40343570821929375 + x2*((174611*pi^20)/183759535834875000 + x2*((-43867*pi^18)/196455649219830000 + x2*((3617*pi^16)/118168811560800000 + x2*(-1/364105581840000*pi^14 + x2*((691*pi^12)/3976032953692800000 + x2*(-1/122339475498240000*pi^10 + x2*(pi^8/3361246195507200000 + x2*(-1/114954619886346240000*pi^6 + x2*(pi^4/4598184795453849600000 + x2*(-1/155112100433309859840000*pi^2 + (pi/1240896803466478878720000 - 1/31022420086661971968000000*x)*x)))))))))))))

    sgn*sl_series(n, x)
end


text:
# Rcpp derivatives (RcppSundials)
function create_derivatives_rcpp(language, sorted_model, name)
  code = ""
  for level in 1:length(sorted_model.SortedEquations)
    for (lhs, rhs) in sorted_model.SortedEquations[level]
        if level == 1
          if(ismatch(r"params|states|forcs", string(rhs.Expr)))
           number = match(r"(?<=\[)[\d]+(?=\])",string(rhs.Expr))
           new_expr = replace(string(rhs.Expr), "[$(number.match)]", "[$(int(number.match) - 1)]")
           code *= "const double $(new_expr);\n"
         else
           code *= "const double $(rhs.Expr);\n"
         end
        else
          mod_expr = copy(rhs.Expr)
          mod_expr = sub_product(mod_expr)
          mod_expr = substitute_power(mod_expr)
          mod_expr = sub_minmax(mod_expr)
          code *= "const double $lhs" * " = " * replace(replace(string(mod_expr), ".*", "*"), ":", "") * ";\n"
        end
    end
  end
  # Return line
  return_line = create_return_line_Rcpp(sorted_model.NamesDerivatives,sorted_model.NamesObserved)
  # Return the output
  code = replace(code, "1.0 *", "")
  up_boiler_plate =
"""
  array<vector<double>, 2> $(name)(const double& t, const vector<double>& states,
            const vector<double>& params, const vector<double>& forcs) { \n
"""
  low_boiler_plate =
"""
}
"""
  return paste("\n", up_boiler_plate, code, return_line,low_boiler_plate)
end

# Rcpp jacobian (RcppSundials)
function create_jacobian_rcpp(language, compressed_model, jacobian_matrix, names_derivatives, name)
  # Substitute * in AST to avoid pretty printing of scalar product
  for i = 1:size(jacobian_matrix, 1), j = 1:size(jacobian_matrix, 2)
      if isa(jacobian_matrix[i,j], Expr)
        jacobian_matrix[i,j] = sub_product(jacobian_matrix[i,j])
      end
  end
  code = ""
  for (lhs, rhs) in compressed_model.SortedEquations[1]
      if (ismatch(r"params|states|forcs", string(rhs.Expr)))
       number = match(r"(?<=\[)[\d]+(?=\])",string(rhs.Expr))
       new_expr = replace(string(rhs.Expr), "[$(number.match)]", "[$(int(number.match) - 1)]")
       code *= "const double $(new_expr);\n"
     else
       code *= "const double $(rhs.Expr);\n"
     end
  end
  for i = 1:size(jacobian_matrix, 1), j = 1:size(jacobian_matrix, 2)
    if isa(jacobian_matrix[i,j], Number) && eval(jacobian_matrix[i,j]) == 0
      continue
    else
      rhs = substitute_power(jacobian_matrix[i,j])
      rhs = sub_minmax(rhs)
      rhs = replace(replace(string(rhs), ".*", "*"), ":", "")
      code *= "output.at($(i-1),$(j - 1)) = $rhs;\n"
    end
  end

  up_boiler_plate =
"""
arma::mat $(name)_jacobian(const double& t, const vector<double>& states,
          const vector<double>& params, const vector<double>& forcs) {
  arma::mat output = arma::eye(states.size(), states.size());\n
"""
  low_boiler_plate =
"""
    return output;
  }
"""

  return paste("\n", up_boiler_plate, code,low_boiler_plate)

end

# Write the file with all the model functions
function write_code_rcpp!(dynamic_type, model_function, name, file)

  up_boiler_plate =
"""
#include <RcppArmadillo.h>
#include <array>
#include <vector>
#include <math.h>
using namespace std;

double ifelse(bool condition, const double& result1, const double& result2) {
  if(condition) {
    return result1;
  } else {
    return result2;
  }
}

inline double heaviside(const double& arg) {
  return arg <= 0.0 ? 0.0 : 1.0;
}

inline double dirac(const double& arg) {
  return abs(arg) <= numeric_limits<double>::epsilon()  ? 0 : numeric_limits<double>::infinity();
}

inline double Min(const double& arg1) {
  return arg1;
}

inline double Max(const double& arg1) {
  return arg1;
}

inline double Min(const double& arg1, const double& arg2) {
  return arg1 <= arg2 ? arg1 : arg2;
}

inline double Max(const double& arg1, const double& arg2) {
  return arg1 >= arg2 ? arg1 : arg2;
}

extern "C" {
"""
  low_boiler_plate =
"""
};
"""

  f = open("$(file).cpp","w")
  println(f, up_boiler_plate)
  println(f, model_function)
  println(f, low_boiler_plate)
  close(f)

  nothing

end

############################ AUXILLIARY FUNCTIONS ##############################

# Because Julia will use "pretty printing" for scalar product, we need to fool the parser
# by temporarily substituying by another binary operator and then doing string replacement
function sub_product(ex::Expr)
  for i in 1:length(ex.args)
  if ex.args[i] == :*
    ex.args[i] = :.*
  elseif isa(ex.args[i], Expr)
     ex.args[i] = sub_product(ex.args[i])
  end
end
return ex
end

sub_product(ex::Symbol) = ex

# Substitute the annoying x^y for pow(x,y)
function substitute_power(ex::Expr)
  new_ex = deepcopy(ex)
  for i in 1:length(new_ex.args)
    if isa(new_ex.args[i], Expr)
      new_ex.args[i] = substitute_power(new_ex.args[i])
    elseif new_ex.args[i] == :^
      new_ex.args[i] = :pow
    end
  end
  return new_ex
end

substitute_power(ex::Any) = ex

function sub_minmax(ex::Expr)
  new_ex = deepcopy(ex)
  for i in 1:length(new_ex.args)
    if new_ex.args[i] == :min || new_ex.args[i] == :max
        new_ex.args[i] == :min ?  (new_ex.args[i] = :Min) : (new_ex.args[i] = :Max)
    elseif isa(new_ex.args[i], Expr)
        new_ex.args[i] = sub_minmax(new_ex.args[i])
  end
end
return new_ex
end

sub_minmax(ex::Symbol) = ex

# Create a return line with an array containing the derivatives and observed functions (for the STL version of the model)
function create_return_line_Rcpp(states, observed)
    return_line = "vector<double> derivatives{"
    for i in states
      if i != states[end]
        return_line = return_line  * i * ","
      else
        return_line = return_line  * i
      end
    end
    return_line = return_line * "};\n vector<double> observed{"
    for i in observed
      if i != observed[end]
        return_line = return_line * i * ", "
      else
        return_line = return_line  * i
      end
    end
    return_line = return_line * "};\n array<vector<double>,2> output{derivatives, observed};\n return output;"
end


text:
### Compute the gradients using a gradient function and matrices Js ###

base_kernel(k::Kernel) = eval(nameof(typeof(k)))

function compute_hyperparameter_gradient(k::Kernel,gradient_function::Function,J::Vector{<:AbstractMatrix})
    return map(gradient_function,J)
end

function compute_hyperparameter_gradient(k::Kernel,gradient_function::Function,Jmm::Vector{<:AbstractMatrix},Jnm::Vector{<:AbstractMatrix},Jnn::Vector{<:AbstractVector},∇E_μ::AbstractVector{T},∇E_Σ::AbstractVector{T},i::Inference,opt::AbstractOptimizer) where {T<:Real}
    return map(gradient_function,Jmm,Jnm,Jnn,[∇E_μ],[∇E_Σ],i,[opt])
end

function apply_gradients_lengthscale!(opt::Optimizer,k::Kernel{T,<:ScaleTransform{<:Base.RefValue{<:Real}}},g::AbstractVector) where {T}
    ρ = first(KernelFunctions.params(k))
    logρ = log(ρ) + update(opt,first(g)*ρ)
    KernelFunctions.set_params!(k,exp(logρ))
end

function apply_gradients_lengthscale!(opt::Optimizer,k::Kernel{T,<:ScaleTransform{<:AbstractVector}},g::AbstractVector) where {T}
    ρ = first(KernelFunctions.params(k))
    logρ = log.(ρ) .+ update(opt,g.*ρ)
    KernelFunctions.set_params!(k,exp.(logρ))
end

function apply_gradients_variance!(gp::Abstract_GP,g::Real)
    logσ = log(gp.σ_k)
    logσ += update(gp.opt_σ,g*gp.σ_k)
    gp.σ_k = exp(logσ)
end

function apply_gradients_mean_prior!(μ::PriorMean,g::AbstractVector,X::AbstractMatrix)
    update!(μ,g,X)
end

function kernelderivative(kernel::Kernel{T,<:ScaleTransform{<:Base.RefValue}},X) where {T}
    p = collect(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kernelmatrix(base_kernel(kernel)(x...),X,obsdim=1),p),size(X,1),size(X,1),1)
    return [J[:,:,1]]
end

function kernelderivative(kernel::Kernel{T,<:ScaleTransform{<:AbstractVector}},X) where {T}
    p = first(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kernelmatrix(base_kernel(kernel)(x),X,obsdim=1),p),size(X,1),size(X,1),length(p))
    return [J[:,:,i] for i in 1:length(p)]
end

function kernelderivative(kernel::Kernel{T,<:ScaleTransform{<:Base.RefValue}},X,Y) where {T}
    p = collect(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kernelmatrix(base_kernel(kernel)(x...),X,Y,obsdim=1),p),size(X,1),size(Y,1),1)
    return [J[:,:,1]]
end

function kernelderivative(kernel::Kernel{T,<:ScaleTransform{<:AbstractVector}},X,Y) where {T}
    p = first(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kernelmatrix(base_kernel(kernel)(x),X,Y,obsdim=1),p),size(X,1),size(Y,1),length(p))
    return [J[:,:,i] for i in 1:length(p)]
end

function kerneldiagderivative(kernel::Kernel{T,<:ScaleTransform{<:Base.RefValue}},X) where {T}
    p = collect(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kerneldiagmatrix(base_kernel(kernel)(x...),X,obsdim=1),p),size(X,1),length(p))
    return [J[:,1]]
end

function kerneldiagderivative(kernel::Kernel{T,<:ScaleTransform{<:AbstractVector}},X) where {T}
    p = first(KernelFunctions.params(kernel))
    J = reshape(ForwardDiff.jacobian(x->kerneldiagmatrix(base_kernel(kernel)(x),X,obsdim=1),p),size(X,1),length(p))
    return [J[:,i] for i in 1:length(p)]
end

function indpoint_derivative(kernel::Kernel,Z::InducingPoints)
    reshape(ForwardDiff.jacobian(x->kernelmatrix(kernel,x,obsdim=1),Z.Z),size(Z,1),size(Z,1),size(Z,1),size(Z,2))
end

function indpoint_derivative(kernel::Kernel,X,Z::InducingPoints)
    reshape(ForwardDiff.jacobian(x->kernelmatrix(kernel,X,x,obsdim=1),Z.Z),size(X,1),size(Z,1),size(Z,1),size(Z,2))
end


text:
"""
Gives a reward of +1 if sp is the target lane
"""
struct TargetLaneReward <: AbstractMLRewardModel
    target_lane::Int
end
function POMDPs.reward(mdp::MLMDP{MLState, MLAction, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState) where D<:AbstractMLDynamicsModel
    return nothing === (s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

function POMDPs.reward(mdp::MLPOMDP{MLState, MLAction, MLPhysicalState, D, TargetLaneReward},
          s::MLState,
          ::MLAction,
          sp::MLState) where D<:AbstractMLDynamicsModel
    return nothing === (s.terminal) && sp.cars[1].y == mdp.rmodel.target_lane
end

"""
Reward of +1 on transition INTO target lane, -lambda on unsafe transitions, -lane_change_cost if sp between lanes
"""
@with_kw struct SuccessReward <: AbstractMLRewardModel
    lambda::Float64                 = 1.0  # always positive
    target_lane::Int                = 4
    brake_penalty_thresh::Float64   = 4.0  # always positive
    speed_thresh::Float64           = 15.0 # always positive
    lane_change_cost::Float64       = 0.0 # always positive
end

function POMDPs.reward(p::NoCrashProblem{SuccessReward}, s::MLState, ::MLAction, sp::MLState)
    if sp.cars[1].y == p.rmodel.target_lane && s.cars[1].y != p.rmodel.target_lane
        r = 1.0
    else
        r = 0.0
    end
    min_speed = minimum(c.vel for c in sp.cars)
    nb_brakes = detect_braking(p, s, sp)
    if nb_brakes > 0 || min_speed < p.rmodel.speed_thresh
        r -= p.rmodel.lambda
    end
    if !isinteger(sp.cars[1].y)
        r -= p.rmodel.lane_change_cost
    end
    return r
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -25.780717722472463 <= q <= 30.482165922896684)

                     add_NL_constraint(m, :(gelu(0.6681236842141707 + 0.4190368949234773*gelu(-0.3928529255603399 + 0.13822049866887154*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.7912756220587216*gelu(-0.5063192893501953 + 0.4673927080857103*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.030565575094480213*gelu(-0.34808094533003686 + 0.7360164287507445*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) + gelu(-0.037313120688143186 + -0.8220585173358868*gelu(-0.3928529255603399 + 0.13822049866887154*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.9374338955858046*gelu(-0.5063192893501953 + 0.4673927080857103*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + -0.1479923276913464*gelu(-0.34808094533003686 + 0.7360164287507445*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) + gelu(0.5136304015740518 + -0.41595218402127365*gelu(-0.3928529255603399 + 0.13822049866887154*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.4562953764921889*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.9204155811731316*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + -0.01955723268679055*gelu(-0.5063192893501953 + 0.4673927080857103*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.2643988862479607*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.8391655243989531*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5]))) + 0.13292433203840748*gelu(-0.34808094533003686 + 0.7360164287507445*gelu(-0.9290908758950072 + -0.9885689195552865*$(x[1]) + 0.968971534656105*$(x[2]) + 0.0899903458574629*$(x[3]) + 0.3621960275577596*$(x[4]) + 0.6428275846795262*$(x[5])) + -0.8890315695148892*gelu(-0.6068769875537994 + 0.23936570066126128*$(x[1]) + -0.8813597146688212*$(x[2]) + 0.7513604373502818*$(x[3]) + -0.8801053294405703*$(x[4]) + -0.09321218509736573*$(x[5])) + -0.7245715776404307*gelu(-0.4931215033246312 + -0.7619598530221103*$(x[1]) + -0.2076077765621065*$(x[2]) + -0.31174640065939396*$(x[3]) + -0.2011575506947998*$(x[4]) + -0.9004282460591928*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export upsample_nearest, ∇upsample_nearest,
    upsample_bilinear, ∇upsample_bilinear,
    pixel_shuffle

"""
    upsample_nearest(x, scale::NTuple{S,Int})
    upsample_nearest(x; size::NTuple{S,Int})

Upsamples the array `x` by integer multiples along the first `S` dimensions.
Subsequent dimensions of `x` are not altered.

Either the `scale` factors or the final output `size` can be specified.

See also [`upsample_bilinear`](@ref), for two dimensions of an `N=4` array.

# Example
```jldoctest
julia> upsample_nearest([1 2 3; 4 5 6], (2, 3))
4×9 Array{$Int,2}:
 1  1  1  2  2  2  3  3  3
 1  1  1  2  2  2  3  3  3
 4  4  4  5  5  5  6  6  6
 4  4  4  5  5  5  6  6  6

julia> ans == upsample_nearest([1 2 3; 4 5 6]; size=(4, 9))  # equivalent
true

julia> upsample_nearest([1 2 3; 4 5 6], (2,))
4×3 Array{$Int,1}:
 1  2  3
 1  2  3
 4  5  6
 4  5  6

julia> ans == upsample_nearest([1 2 3; 4 5 6], size=(4,))
true
```
"""
function upsample_nearest(x::AbstractArray; size::NTuple{S,Int}) where S
    xsize = Base.size(x)[1:S]
    all(size .% xsize .== 0) || throw(ArgumentError("expected output size divisible by input size"))
    return upsample_nearest(x, size .÷ xsize)
end

function upsample_nearest(x::AbstractArray{T,N}, scales::NTuple{S, <:Integer}) where {T,N,S}
    S in 1:N || throw(ArgumentError("can't upsample ndims(x)=$N with scale=$scales"))
    outsize = ntuple(d -> d<=S ? scales[d] * size(x,d) : size(x,d), N)
    out = similar(x, T, outsize)
    writesize = ntuple(N+S) do d
        d > 2S && return size(x, d-S)
        isodd(d) ? scales[cld(d,2)] : size(x, cld(d,2))
    end
    readsize = ntuple(N+S) do d
        d > 2S && return size(x, d-S)
        isodd(d) ? 1 : size(x, cld(d,2))
    end
    reshape(out, writesize) .= reshape(x, readsize)
    out
end

function ∇upsample_nearest(x::AbstractArray{T,N}, scales::NTuple{S, <:Integer}) where {T,N,S}
    outsize = ntuple(N) do d
        d > S && return size(x,d)
        rem(size(x,d), scales[d]) == 0 || throw(ArgumentError("expected input array evenly divisible by scale=$scales, got size(x)=$(size(x))"))
        div(size(x,d), scales[d])
    end
    tempsize = ntuple(N+S) do d
        d > 2S && return size(x, d-S)
        s = scales[cld(d,2)]
        isodd(d) ? s : div(size(x, cld(d,2)),s)
    end
    mid = sum(reshape(x, tempsize), dims=ntuple(d -> 2d-1, S))
    reshape(mid, outsize)
end

function rrule(::typeof(upsample_nearest), x::AbstractArray, s::Tuple)
    Ω = upsample_nearest(x, s)
    upsample_nearest_pullback(Δ) = (NO_FIELDS, ∇upsample_nearest(Δ, s), DoesNotExist())
    return Ω, upsample_nearest_pullback
end

# utility function
@inline function compute_source_index_and_lambda(
    ratio, # 0 < ratio < 1
    output_index,
    input_size,
    output_size
)
    real_input_index = ratio*output_index
    input_index0 = floor(Int, real_input_index) # typecast to int was here in C++
    offset = (input_index0 < input_size - 1) ? 1 : 0
    input_index1 = input_index0 + offset
    lambda1 = real_input_index - input_index0
    lambda0 = 1 - lambda1
    return input_index0, input_index1, lambda0, lambda1
end

"""
    upsample_bilinear(x::AbstractArray{T,4}, scale::NTuple{2,Real})
    upsample_bilinear(x::AbstractArray{T,4}; size::NTuple{2,Integer})

Upsamples the first 2 dimensions of the array `x` by the upsample factors stored in `scale`,
using bilinear interpolation. As an alternative to using `scale`, the resulting image `size`
can be directly specified with a keyword argument.

The size of the output is equal to
`(scale[1]*S1, scale[2]*S2, S3, S4)`, where `S1, S2, S3, S4 = size(x)`.

# Examples

```jldoctest
julia> x = reshape(Float32[1 2 3; 4 5 6], (2,3,1,1))
2×3×1×1 Array{Float32,4}:
[:, :, 1, 1] =
 1.0  2.0  3.0
 4.0  5.0  6.0

julia> upsample_bilinear(x, (2, 3))
4×9×1×1 Array{Float32,4}:
[:, :, 1, 1] =
 1.0  1.25  1.5  1.75  2.0  2.25  2.5  2.75  3.0
 2.0  2.25  2.5  2.75  3.0  3.25  3.5  3.75  4.0
 3.0  3.25  3.5  3.75  4.0  4.25  4.5  4.75  5.0
 4.0  4.25  4.5  4.75  5.0  5.25  5.5  5.75  6.0

julia> ans == upsample_bilinear(x; size=(4, 9))  # specify ouput size instead
true

julia> upsample_bilinear(x, (2.5, 3.5))  # non-integer scaling factors are allowed
5×10×1×1 Array{Float32,4}:
[:, :, 1, 1] =
 1.0   1.22222  1.44444  1.66667  1.88889  2.11111  2.33333  2.55556  2.77778  3.0
 1.75  1.97222  2.19444  2.41667  2.63889  2.86111  3.08333  3.30556  3.52778  3.75
 2.5   2.72222  2.94444  3.16667  3.38889  3.61111  3.83333  4.05556  4.27778  4.5
 3.25  3.47222  3.69444  3.91667  4.13889  4.36111  4.58333  4.80556  5.02778  5.25
 4.0   4.22222  4.44444  4.66667  4.88889  5.11111  5.33333  5.55556  5.77778  6.0
```
"""
function upsample_bilinear(x::AbstractArray{<:Any,4}, scale::NTuple{2,Real})
    outsize = ntuple(i -> floor(Int, scale[i] * Base.size(x, i)), 2)
    return upsample_bilinear(x; size=outsize)
end

upsample_bilinear(x, scale::Real) = upsample_bilinear(x, (scale,scale))

function upsample_bilinear(x::AbstractArray{T,4}; size::NTuple{2,Integer}) where T
    w,h,c,n = Base.size(x)
    if (w,h) == size
        return x
    end
    y = similar(x, T, size..., c, n)
    return upsample_bilinear_whcn!(y, x)
end

function upsample_bilinear(x::AbstractArray{T,4}; size::NTuple{2,Integer}) where T<:Integer
    y = float.(x)
    res = upsample_bilinear(y; size=size)
    return round.(T, res)
end

# this is the core function which works on arrays of arbitrary size
# the implementation is a translation of https://github.com/pytorch/pytorch/blob/master/aten/src/ATen/native/cpu/UpSampleMoreKernel.cpp
# which implements open-cv style linear interpolation / upsampling
# for simplicity, corners are aligned and all logic for other behaviour has been stripped
# - whcn because there is also a cwhn implementation
# - the function is parallelized using @threads
# - RGB types could be supported via reinterpreting
# - integer types need to be converted to Float and back
# - rationals work, but are slow
function upsample_bilinear_whcn!(output::AbstractArray{T,4}, input::AbstractArray{T,4}) where T
    size(input)[3:4] == size(output)[3:4] || error("Number of input and output channels and batches must match. Got input $(size(input)) and output $(size(output))")
    in_w, in_h, channels, batches = size(input)
    # treat batch and channel dimension as one for better parallelization granularity
    channels *= batches
    out_w, out_h, _, _ = size(output)
    output_slice_size = out_h * out_w

    # T() and // so that we can handle rationals (super slow)
    width_scale  = T((in_w - 1) // (out_w - 1))
    height_scale = T((in_h - 1) // (out_h - 1))

    @inline idx(c, h, w) = c * in_h * in_w + h * in_w + w + 1

    @inbounds Threads.@threads for c in 0:channels-1
        for oh in 0:out_h-1
            ih0, ih1, h0lambda, h1lambda = compute_source_index_and_lambda(height_scale, oh, in_h, out_h)
            for ow in 0:out_w-1
                iw0, iw1, w0lambda, w1lambda = compute_source_index_and_lambda(width_scale, ow, in_w, out_w)
                output_offset = c * output_slice_size + oh * out_w + ow + 1
                output[output_offset] =
                    (h0lambda * w0lambda * input[idx(c, ih0, iw0)] + # h0 * w0 * i00
                     h0lambda * w1lambda * input[idx(c, ih0, iw1)] + # h0 * w1 * i01
                     h1lambda * w0lambda * input[idx(c, ih1, iw0)] + # h1 * w0 * i10
                     h1lambda * w1lambda * input[idx(c, ih1, iw1)])  # h1 * w1 * i11
            end
        end
    end
    return output
end

"""
    ∇upsample_bilinear(Δ::AbstractArray{T,4}; size::NTuple{2,Integer}) where T

# Arguments
- `Δ`: Incoming gradient array, backpropagated from downstream layers
- `size`: Lateral (W,H) size of the image upsampled in the first place

# Outputs
- `dx`: Downsampled version of `Δ`
"""
function ∇upsample_bilinear(Δ::AbstractArray{T,4}; size::NTuple{2,Integer}) where T
    w, h, c, n = Base.size(Δ)
    out_w, out_h = size
    if (w,h) == (out_w, out_h)
        return Δ
    end
    dx = zero(similar(Δ, T, out_w, out_h, c, n))
    return ∇upsample_bilinear_whcn!(dx, Δ)
end

function ∇upsample_bilinear_whcn!(dx::AbstractArray{T,4}, Δ::AbstractArray{T,4}) where T
    size(dx)[3:4] == size(Δ)[3:4] || error("Number of input and output channels and batches must match. Got input $(size(input)) and output $(size(output))")
    in_w, in_h, channels, batches = size(dx)

    # treat batch and channel dimension as one for better parallelization granularity
    channels *= batches
    out_w, out_h, _, _ = size(Δ)
    output_slice_size = out_h * out_w

    width_scale  = T((in_w - 1) // (out_w - 1))
    height_scale = T((in_h - 1) // (out_h - 1))

    @inline idx(c, h, w) = c * in_h * in_w + h * in_w + w + 1

    @inbounds Threads.@threads for c in 0:channels-1
        for oh in 0:out_h-1
            ih0, ih1, h0lambda, h1lambda = compute_source_index_and_lambda(height_scale, oh, in_h, out_h)
            for ow in 0:out_w-1
                iw0, iw1, w0lambda, w1lambda = compute_source_index_and_lambda(width_scale, ow, in_w, out_w)
                output_offset = c * output_slice_size + oh * out_w + ow + 1
                Δ_value = Δ[output_offset]
                dx[idx(c, ih0, iw0)] += h0lambda * w0lambda * Δ_value # i00
                dx[idx(c, ih0, iw1)] += h0lambda * w1lambda * Δ_value # i01
                dx[idx(c, ih1, iw0)] += h1lambda * w0lambda * Δ_value # i10
                dx[idx(c, ih1, iw1)] += h1lambda * w1lambda * Δ_value # i11
            end
        end
    end
    return dx
end

function rrule(::typeof(upsample_bilinear), x; size)
    Ω = upsample_bilinear(x; size=size)
    function upsample_bilinear_pullback(Δ)
        (NO_FIELDS, ∇upsample_bilinear(Δ; size=(Base.size(x,1),Base.size(x,2))))
    end
    return Ω, upsample_bilinear_pullback
end

"""
    pixel_shuffle(x, r::Integer)

Pixel shuffling operation, upscaling by a factor `r`.

For 4-arrays representing `N` images, the operation converts input `size(x) == (W, H, r^2*C, N)`
to output of size `(r*W, r*H, C, N)`. For `D`-dimensional data, it expects `ndims(x) == D+2`
with channel and batch dimensions, and divides the number of channels by `r^D`.

Used in super-resolution networks to upsample towards high resolution features.
Reference: Shi et. al., "Real-Time Single Image and Video Super-Resolution ...", CVPR 2016, https://arxiv.org/abs/1609.05158

# Examples

```jldoctest
julia> x = [10i + j + channel/10 for i in 1:2, j in 1:3, channel in 1:4, batch in 1:1]
2×3×4×1 Array{Float64,4}:
[:, :, 1, 1] =
 11.1  12.1  13.1
 21.1  22.1  23.1

[:, :, 2, 1] =
 11.2  12.2  13.2
 21.2  22.2  23.2

[:, :, 3, 1] =
 11.3  12.3  13.3
 21.3  22.3  23.3

[:, :, 4, 1] =
 11.4  12.4  13.4
 21.4  22.4  23.4

julia> pixel_shuffle(x, 2)  # 4 channels used up as 2x upscaling of image dimensions
4×6×1×1 Array{Float64,4}:
[:, :, 1, 1] =
 11.1  11.3  12.1  12.3  13.1  13.3
 11.2  11.4  12.2  12.4  13.2  13.4
 21.1  21.3  22.1  22.3  23.1  23.3
 21.2  21.4  22.2  22.4  23.2  23.4

julia> y = [i + channel/10 for i in 1:3, channel in 1:6, batch in 1:1]
3×6×1 Array{Float64, 3}:
[:, :, 1] =
 1.1  1.2  1.3  1.4  1.5  1.6
 2.1  2.2  2.3  2.4  2.5  2.6
 3.1  3.2  3.3  3.4  3.5  3.6

julia> pixel_shuffle(y, 2)  # 1D image, with 6 channels reduced to 3
6×3×1 Array{Float64,3}:
[:, :, 1] =
 1.1  1.3  1.5
 1.2  1.4  1.6
 2.1  2.3  2.5
 2.2  2.4  2.6
 3.1  3.3  3.5
 3.2  3.4  3.6
```
"""
function pixel_shuffle(x::AbstractArray, r::Integer)
    ndims(x) > 2 || throw(ArgumentError("expected x with at least 3 dimensions"))
    d = ndims(x) - 2
    sizein = size(x)[1:d]
    cin, n = size(x, d+1), size(x, d+2)
    cin % r^d == 0 || throw(ArgumentError("expected channel dimension to be divisible by r^d = $(
        r^d), where d=$d is the number of spatial dimensions. Given r=$r, input size(x) = $(size(x))"))
    cout = cin ÷ r^d
    x = reshape(x, sizein..., ntuple(i->r, d)..., cout, n)
    perm = hcat(d+1:2d, 1:d) |> transpose |> vec  # = [d+1, 1, d+2, 2, ..., 2d, d]
    x = permutedims(x, (perm..., 2d+1, 2d+2))
    return reshape(x, map(s -> s*r, sizein)..., cout, n)
end


text:
@io struct MachOSymtabCmd{H <: MachOHandle} <: MachOLoadCmd{H}
    symoff::UInt32
    nsyms::UInt32
    stroff::UInt32
    strsize::UInt32
end

show(io::IO, lc::MachOSymtabCmd) = write(io, "SymtabCmd")

# Accessors for symbols info
symtab_symbols_offset(cmd::MachOSymtabCmd) = cmd.symoff
symtab_num_symbols(cmd::MachOSymtabCmd) = cmd.nsyms

# Accessors for strtab info
symtab_strtab_offset(cmd::MachOSymtabCmd) = cmd.stroff
symtab_strtab_length(cmd::MachOSymtabCmd) = cmd.strsize

@derefmethod symtab_symbols_offset(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}
@derefmethod symtab_num_symbols(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}
@derefmethod symtab_strtab_offset(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}
@derefmethod symtab_strtab_length(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}) where {H <: MachOHandle}

function strtab_lookup(cmd::MachOLoadCmdRef{H,MachOSymtabCmd{H}}, index) where {H <: MachOHandle}
    oh = handle(cmd)
    seek(oh, symtab_strtab_offset(cmd) + index)
    return strip(readuntil(oh, '\0'), '\0')
end

# We don't actually use this yet, but it doesn't hurt to have it here
@io struct MachODySymtabCmd{H <: MachOHandle} <: MachOLoadCmd{H}
    ilocalsym::UInt32
    nlocalsym::UInt32
    iextdefsym::UInt32
    nextdefsym::UInt32
    iundefsym::UInt32
    nundefsym::UInt32
    tocoff::UInt32
    ntoc::UInt32
    modtaboff::UInt32
    nmodtab::UInt32
    extrefsymoff::UInt32
    nextrefsyms::UInt32
    indirectsymoff::UInt32
    nindirectsyms::UInt32
    extreloff::UInt32
    nextrel::UInt32
    locreloff::UInt32
    nlocrel::UInt32
end


text:
mutable struct LinearizedStep{T}
	z::AbstractVector{T}
	θ::AbstractVector{T}
	κ::T
	r::AbstractVector{T}
	rz::AbstractMatrix{T}
	rθ::AbstractMatrix{T}
end

function LinearizedStep(s::Simulation, z::AbstractVector{T}, θ::AbstractVector{T}, κ::T) where T
	model = s.model
	env = s.env

	nz = num_var(model, env)
	nθ = num_data(model)

	z0 = SizedVector{nz,T}(z)
	θ0 = SizedVector{nθ,T}(θ)
	κ0 = κ
	r0 = zeros(SizedVector{nz,T})
	rz0 = zeros(nz, nz)
	rθ0 = zeros(nz, nθ)

	s.res.r!(r0, z0, θ0, κ0)
	s.res.rz!(rz0, z0, θ0)
	s.res.rθ!(rθ0, z0, θ0)

	return LinearizedStep(z0, θ0, κ0, r0, rz0, rθ0)
end

"""
	Create LinearizedStep.
"""
function LinearizedStep(s::Simulation)
	model = s.model
	env = s.env

	nz = num_var(model, env)
	nθ = num_data(model)
	z0 = zeros(SizedVector{nz})
	θ0 = zeros(SizedVector{nθ})
	κ0 = 0.0

	return LinearizedStep(s, z0, θ0, κ0)
end

function update!(lin::LinearizedStep, s::Simulation, z, θ)
	lin.z .= z
	lin.θ .= θ
	s.res.r!(lin.r, z, θ, lin.κ)
	s.res.rz!(lin.rz, z, θ)
	s.res.rθ!(lin.rθ, z, θ)
	return nothing
end


text:
function lookup_function(mod::LLVM.Module, func_name::String)
    LLVM.Function(LLVM.API.LLVMGetNamedFunction(mod, func_name))
end

lookup_function(ee::ExecutionEngine, func_name::String) = functions(ee)[func_name]

link(lib::AbstractString) = LLVM.load_library_permantly(lib)

link_crt(ee::ExecutionEngine) = LLVM.API.LLVMRunStaticConstructors(ee.ref)

function get_buffer(x::String, name="", copy=true)
    data = unsafe_wrap(Vector{UInt8}, x)
    return MemoryBuffer(data, name, copy)
end


text:
Base.@kwdef struct Class
    domain::String
    system::String
    value::String
end

function DiagnosisClass(; system, value)
    domain = "diagnosis"
    class = Class(; domain, system, value)
    return class
end


text:
"""
    Module EntropyStableEuler

Includes general math tools
"""

module EntropyStableEuler

#const γ=1.4
#export γ
export logmean
export u_vfun, v_ufun, betafun, Sfun,sfun,pressure_fun
export euler_fluxes, wavespeed
export vortex, primitive_to_conservative

include("./logmean.jl")
include("./euler_fluxes.jl")
include("./euler_variables.jl")

# 2D isentropic vortex solution for testing. assumes domain around [0,20]x[-5,5]
function vortex(x,y,t,γ=1.4)

    x0 = 5
    y0 = 0
    beta = 5
    r2 = @. (x-x0-t)^2 + (y-y0)^2

    u = @. 1 - beta*exp(1-r2)*(y-y0)/(2*pi)
    v = @. beta*exp(1-r2)*(x-x0-t)/(2*pi)
    rho = @. 1 - (1/(8*γ*pi^2))*(γ-1)/2*(beta*exp(1-r2))^2
    rho = @. rho^(1/(γ-1))
    p = @. rho^γ

    return (rho, u, v, p)
end

end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
struct binned1dDensity
    grid::Array{Float64,1}
    cumarr::Array{Float64}
    density
end

function getbinned1dDensity(g, lim, Nbins)
    grid = [σi for σi=LinRange(lim[1], lim[2], Nbins)]
    #
    weights = [g(v...) for v in (grid[2:end] .+ grid[1:end-1]) ./ 2]
    weights ./= sum(weights)
    cumarr = [0;cumsum(vcat(weights...), dims=1)]
    return binned1dDensity(grid,cumarr,g)
end

import Base:rand
function rand(bD::binned1dDensity)
    # get
    binind = findfirst(bD.cumarr .> rand())-1
    σl, σr = bD.grid[binind], bD.grid[binind+1]
    σ = σl+rand()*(σr-σl)
    bD.density(σ) == 0.0 && return rand(bD)
    return σ
end

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

struct binned2dDensity
    grid::Matrix{Array{Float64,1}}
    cumarr::Array{Float64}
    density
end

arg1_lims(binned2dDensity) = (binned2dDensity.grid[1,1][1],binned2dDensity.grid[end,end][1])
arg2_lims(binned2dDensity) = (binned2dDensity.grid[1,1][2],binned2dDensity.grid[end,end][2])

function getbinned2dDensity(g, arg1_lims, arg2_lims,  N1, N2)
    grid = hcat([[[x1,x2]
        for x2=range(arg2_lims[1], arg2_lims[2], length=N1)]
        for x1=range(arg1_lims[1], arg1_lims[2], length=N2)]...)
    #
    weights = [g(v...) for v in (grid[2:end,2:end] .+ grid[1:end-1,1:end-1]) ./ 2]
    # normalize by the cell size
    for i=1:size(grid,1)-1, j=1:size(grid,2)-1
        weights[i,j] *= (grid[i+1,j][2] - grid[i,j][2])*(grid[i,j+1][1] - grid[i,j][1])
    end
    weights ./= sum(weights)
    cumarr = [0;cumsum(vcat(weights...), dims=1)]
    return binned2dDensity(grid,cumarr, g)
end

import Base:rand
function rand(bD::binned2dDensity)
    # get
    binind = findfirst(bD.cumarr .> rand())-2
    # pars back
    Nrows, Ncols = size(bD.grid)
    indCol = div(binind, Nrows-1) + 1
    indRow = mod(binind, Nrows-1) + 1
    indRow+1 > Nrows && error("indRow+1 > Nrows: binind = $binind")
    indCol+1 > Ncols && error("indCol+1 > Ncols: binind = $binind")
    s1, σ1 = bD.grid[indRow,indCol]
    s2, σ2 = bD.grid[indRow+1,indCol+1]
#         s1 != s2 && error("Something is wrong!")
    s = s1+rand()*(s2-s1)
    σ = σ1+rand()*(σ2-σ1)
    bD.density(s,σ) == 0.0 && return rand(bD)
    return [s, σ]
end

function gridded_density_function(x, bD::binned2dDensity)
    (s,σ) = x
    (bD.grid[1,1][1] > s || s > bD.grid[1,end][1]) && return 0.0
    (bD.grid[1,1][2] > σ || σ > bD.grid[end,1][2]) && return 0.0
    #
    mp_s = map(x->x[1], bD.grid[1,:])
    i0 = findfirst(mp_s .> s)
    s0 = (mp_s[i0]+mp_s[i0-1])/2
    #
    mp_σ = map(x->x[2], bD.grid[:,1])
    j0 = findfirst(mp_σ .> σ)
    σ0 = (mp_σ[j0]+mp_σ[j0-1])/2
    #
    return bD.density(s0,σ0)!=0.0 ?  bD.density(s0,σ0) : 0.0;
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
is_flu_outlier() = nothing

 # Flu specific
const lineages = ["h3n2" , "h1n1pdm", "vic", "yam"]
# From genebank file of strain A/Beijing/32/1992
const CDS = Dict(("h3n2","na") => [(4,1410)],
	("h3n2","ha") => [(1,48),(49,1035),(1036,1698)],
	("h1n1pdm", "na") => [(9,1418)],
	("h1n1pdm", "ha") => [(21,71), (72,1052), (1053,1718)])

const gene_positions = Dict(("h3n2","na") => Dict("NA" => (4,1410)),
	("h3n2","ha") => Dict("SigPep" => (1,48), "HA1" => (49,1035), "HA2" => (1036,1698)),
	("h1n1pdm", "na") => Dict("NA" => (9,1418)),
	("h1n1pdm", "ha") => Dict("SigPep"=>(21,71), "HA1"=>(72,1052), "HA2"=>(1053, 1718)))
# Offsets for going from local to global coordinates on proteins (mainly HA).
# Somehow, nextstrain misses one AA between SigPep and HA1, and another betw HA1 and HA2. Maybe stop codons that I translated? Anyway, that's why the 17 and 344 = 327 + 17 for h1n1pdm for example
const gene_offsets = Dict(("h1n1pdm","ha") => Dict("SigPep" => 0, "HA1" => 17, "HA2" => 344),
						("h1n1pdm","na") => Dict("NA" => 0),
						("h3n2","ha") => Dict("SigPep" => 0, "HA1" => 16, "HA2" => 345),
						("h3n2","na") => Dict("NA" => 0))

const substitution_rate = Dict(("h3n2","ha")=>4e-3, ("h3n2","na")=>3e-3) # Number of substitutions per site per year for NUCLEOTIDES
# Note: it's quite similar for AAs
global coalescence_time = Year(6)
global lbi_integration_time = Day(100)
# Settings from the seasonal-flu nextstrain module
const segments = ["ha","na"]
const nextstrain_regions = ["north_america", "south_america", "europe", "china", "oceania","southeast_asia", "japan_korea", "south_asia", "west_asia", "africa"]
#
const flu_usual_header_fields = ["strain", "virus", "", "date", "region", "country", "", "", "", "segment"]
const augur_all_header_fields = ["strain", "virus", "isolate_id", "date", "region", "country", "division", "location", "passage", "authors", "age", "gender"]
flu_usual_filters(flulineage) = [!is_flu_outlier(flulineage), BioTools.hasdate, s->BioTools.gapfilter(s,threshold=0.05)]
#
function get_standard_coordinates(i::Int64; lineage="h3n2", gene="ha")
	idx = sort([(k,x[2]/3) for (k,x) in Flu.gene_positions["h3n2", "ha"]], by=x->x[2])
	offset=0
	for (k,(ref,ix)) in enumerate(idx)
		if i <= ix
			return (ref, Int64(i-offset))
		end
		offset = ix
	end
	error("AA position $i too large for gene ($lineage,$gene)")
end

outliers = Dict{Union{String,Missing},Array{String}}(missing=>String[])
for l in lineages
	p = dirname(pathof(FluPredictibility)) * "/config/outliers_$(l).txt"
	if isfile(p)
		outliers[l] = vec(readdlm(p, String))
	else
		outliers[l] = String[]
	end
	p = dirname(pathof(FluPredictibility)) * "/config/my_outliers_$(l).txt"
	if isfile(p)
		append!(outliers[l], vec(readdlm(p,String)))
	end
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -16.424759279035918 <= q <= 20.23452089871623)

                     add_NL_constraint(m, :(softplus(0.8420310557117703 + 0.34198141401583015*softplus(-0.5389217244927917 + 0.5798025520034438*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.1887519896617813*softplus(-0.23921485861490321 + 0.5030319708107833*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.2882490310085508*softplus(-0.6545847926105979 + 0.9774628615021363*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) + softplus(-0.35188345522271236 + -0.11415509244934263*softplus(-0.5389217244927917 + 0.5798025520034438*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.8088728193024375*softplus(-0.23921485861490321 + 0.5030319708107833*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + 0.8271426964720727*softplus(-0.6545847926105979 + 0.9774628615021363*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) + softplus(0.4798469380896484 + 0.8134725364327324*softplus(-0.5389217244927917 + 0.5798025520034438*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + 0.4428632765655749*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + 0.9033224385695577*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + -0.12585836968848696*softplus(-0.23921485861490321 + 0.5030319708107833*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.2081481501163327*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.2939216839670169*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4]))) + 0.613011911779838*softplus(-0.6545847926105979 + 0.9774628615021363*softplus(0.8931287257144485 + 0.7044247236000336*$(x[1]) + 0.9746565459142618*$(x[2]) + -0.37297510015038515*$(x[3]) + 0.2762326932574686*$(x[4])) + -0.21180887804098658*softplus(-0.08155088834170021 + 0.2657824238897022*$(x[1]) + -0.7808635499515542*$(x[2]) + -0.8071446449659851*$(x[3]) + 0.29051627947702796*$(x[4])) + -0.7013494059895882*softplus(0.5572058350137663 + -0.028961919981738138*$(x[1]) + -0.957730339736961*$(x[2]) + -0.8718568123342201*$(x[3]) + -0.17639541567112582*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function drawdot(pos)
    gsave()
    sethue("red")
    circle(pos, 5, :fill)
    grestore()
end

function text_rotation_tests(fname)
    Drawing(800, 800, fname)
    origin()
    background("white")
    fontsize(10)

    for (n, θ) in enumerate(range(0, 2π, length=11))
        pt = polar(100, θ)
        @layer begin
            sethue("red")
            circle(pt, 5, :fill)
        end
        text("LEFT BOTTOM",  pt, halign=:left, valign=:bottom, angle=θ)
        text("RIGHT BOTTOM", pt, halign=:right, valign=:bottom, angle=θ)
        text("LEFT TOP",     pt, halign=:left, valign=:top, angle=θ)
        text("RIGHT TOP",    pt, halign=:right, valign=:top, angle=θ)
    end

    for (n, θ) in enumerate(range(0, 2π, length=15))
        pt = polar(300, θ)
        @layer begin
            sethue("red")
            circle(pt, 5, :fill)
        end
        text("CENTER BOTTOM",  pt, halign=:center, valign=:bottom, angle=θ)
        text("CENTER TOP",     pt, halign=:center, valign=:top, angle=θ)
    end

    @test finish() == true
    println("...finished test: output in $(fname)")
end

text_rotation_tests("text-rotation-tests.svg")


text:
# This file defines midlevel api wrappers. We include name normalization for methods that are
# applicable to different hdf5 api-layers. We still try to adhere close proximity to the underlying
# method name in the hdf5-library.

"""
    HDF5.set_extent_dims(dset::HDF5.Dataset, new_dims::Dims)

Change the current dimensions of a dataset to `new_dims`, limited by
`max_dims = get_extent_dims(dset)[2]`. Reduction is possible and leads to loss of truncated data.
"""
function set_extent_dims(dset::Dataset, size::Dims)
    checkvalid(dset)
    API.h5d_set_extent(dset, API.hsize_t[reverse(size)...])
end

"""
    HDF5.set_extent_dims(dspace::HDF5.Dataspace, new_dims::Dims, max_dims::Union{Dims,Nothing} = nothing)

Change the dimensions of a dataspace `dspace` to `new_dims`, optionally with the maximum possible
dimensions `max_dims` different from the active size `new_dims`. If not given, `max_dims` is set equal
to `new_dims`.
"""
function set_extent_dims(dspace::Dataspace, size::Dims, max_dims::Union{Dims,Nothing} = nothing)
    checkvalid(dspace)
    rank = length(size)
    current_size = API.hsize_t[reverse(size)...]
    maximum_size = isnothing(max_dims) ? C_NULL : [reverse(max_dims .% API.hsize_t)...]
    API.h5s_set_extent_simple(dspace, rank, current_size, maximum_size)
    return nothing
end

"""
    HDF5.get_extent_dims(obj::Union{HDF5.Dataspace, HDF5.Dataset, HDF5.Attribute}) -> dims, maxdims

Get the array dimensions from a dataspace, dataset, or attribute and return a tuple of `dims` and `maxdims`.
"""
function get_extent_dims(obj::Union{Dataspace,Dataset,Attribute})
    dspace = obj isa Dataspace ? checkvalid(obj) : dataspace(obj)
    h5_dims, h5_maxdims = API.h5s_get_simple_extent_dims(dspace)
    # reverse dimensions since hdf5 uses C-style order
    N = length(h5_dims)
    dims = ntuple(i -> @inbounds(Int(h5_dims[N-i+1])), N)
    maxdims = ntuple(i -> @inbounds(h5_maxdims[N-i+1]) % Int, N) # allows max_dims to be specified as -1 without triggering an overflow
    obj isa Dataspace || close(dspace)
    return dims, maxdims
end

"""
    HDF5.get_chunk_offset(dataset_id, index)

Get 0-based offset of chunk from 0-based `index`. The offsets are returned in Julia's column-major order rather than hdf5 row-major order.
For a 1-based API, see `HDF5.ChunkStorage`.
"""
function get_chunk_offset(dataset_id, index)
    extent = size(dataset_id)
    chunk = get_chunk(dataset_id)
    chunk_indices = CartesianIndices(ntuple(i -> 0:extent[i]÷chunk[i]-1, length(extent)))
    offset = API.hsize_t.(chunk_indices[index + 1].I .* chunk)
    return offset
end

"""
    HDF5.get_chunk_index(dataset_id, offset)

Get 0-based index of chunk from 0-based `offset` returned in Julia's column-major order.
For a 1-based API, see `HDF5.ChunkStorage`.
"""
function get_chunk_index(dataset_id, offset)
    extent = size(dataset_id)
    chunk = get_chunk(dataset_id)
    chunk_indices = LinearIndices(ntuple(i->0:extent[i]÷chunk[i]-1, length(extent)))
    chunk_indices[(offset .÷ chunk .+ 1)...] - 1
end

"""
    HDF5.get_num_chunks_per_dim(dataset_id)

Get the number of chunks in each dimension in Julia's column-major order.
"""
function get_num_chunks_per_dim(dataset_id)
    extent = size(dataset_id)
    chunk = get_chunk(dataset_id)
    return extent .÷ chunk
end

"""
    HDF5.get_num_chunks(dataset_id)

Returns the number of chunks in a dataset. Equivalent to `API.h5d_get_num_chunks(dataset_id, HDF5.H5S_ALL)`.
"""
function get_num_chunks(dataset_id)
    @static if v"1.10.5" ≤ API._libhdf5_build_ver
        API.h5d_get_num_chunks(dataset_id)
    else
        prod(get_num_chunks_per_dim(dataset_id))
    end
end

"""
    HDF5.get_chunk_length(dataset_id)

Retrieves the chunk size in bytes. Equivalent to `API.h5d_get_chunk_info(dataset_id, index)[:size]`.
"""
function get_chunk_length(dataset_id)
    type = API.h5d_get_type(dataset_id)
    chunk = get_chunk(dataset_id)
    return Int(API.h5t_get_size(type) * prod(chunk))
end

function vlen_get_buf_size(dataset_id)
    type = API.h5d_get_type(dataset_id)
    space = API.h5d_get_space(dataset_id)
    API.h5d_vlen_get_buf_size(dataset_id, type, space)
end

"""
    HDF5.read_chunk(dataset_id, offset, [buf]; dxpl_id = HDF5.API.H5P_DEFAULT, filters = Ref{UInt32}())

Helper method to read chunks via 0-based offsets in a `Tuple`.

Argument `buf` is optional and defaults to a `Vector{UInt8}` of length determined by `HDF5.get_chunk_length`.
Argument `dxpl_id` can be supplied a keyword and defaults to `HDF5.API.H5P_DEFAULT`.
Argument `filters` can be retrieved by supplying a `Ref{UInt32}` value via a keyword argument.

This method returns `Vector{UInt8}`.
"""
function read_chunk(dataset_id, offset,
        buf::Vector{UInt8} = Vector{UInt8}(undef, get_chunk_length(dataset_id));
        dxpl_id = API.H5P_DEFAULT,
        filters = Ref{UInt32}()
    )
    API.h5d_read_chunk(dataset_id, dxpl_id, offset, filters, buf)
    return buf
end

"""
    HDF5.read_chunk(dataset_id, index::Integer, [buf]; dxpl_id = HDF5.API.H5P_DEFAULT, filters = Ref{UInt32}())

Helper method to read chunks via 0-based integer `index`.

Argument `buf` is optional and defaults to a `Vector{UInt8}` of length determined by `HDF5.h5d_get_chunk_info`.
Argument `dxpl_id` can be supplied a keyword and defaults to `HDF5.API.H5P_DEFAULT`.
Argument `filters` can be retrieved by supplying a `Ref{UInt32}` value via a keyword argument.

This method returns `Vector{UInt8}`.
"""
function read_chunk(dataset_id, index::Integer,
        buf::Vector{UInt8} = Vector{UInt8}(undef, get_chunk_length(dataset_id));
        dxpl_id = API.H5P_DEFAULT,
        filters = Ref{UInt32}()
    )
    offset = [reverse(get_chunk_offset(dataset_id, index))...]
    read_chunk(dataset_id, offset, buf; dxpl_id = dxpl_id, filters = filters)
end

"""
    HDF5.write_chunk(dataset_id, offset, buf::Vector{UInt8}; dxpl_id = HDF5.API.H5P_DEFAULT, filter_mask = 0)

Helper method to write chunks via 0-based offsets `offset` as a `Tuple`.
"""
function write_chunk(dataset_id, offset, buf::Vector{UInt8}; dxpl_id = API.H5P_DEFAULT, filter_mask=0)
    API.h5d_write_chunk(dataset_id, dxpl_id, filter_mask, offset, length(buf), buf)
end

"""
    HDF5.write_chunk(dataset_id, index::Integer, buf::Vector{UInt8}; dxpl_id = API.H5P_DEFAULT, filter_mask = 0)

Helper method to write chunks via 0-based integer `index`.
"""
function write_chunk(dataset_id, index::Integer, buf::Vector{UInt8}; dxpl_id = API.H5P_DEFAULT, filter_mask = 0)
    offset = [reverse(get_chunk_offset(dataset_id, index))...]
    write_chunk(dataset_id, offset, buf; dxpl_id = dxpl_id, filter_mask = filter_mask)
end


text:
using JML
using Test
using  PrettyPrint

@testset "snooping" begin
    @test 2 == JML.second([1, 2, 3])
    @test "2" == JML.second(["1", "2", "3"])
    @test "3" == JML.LConst("3")._1
    @test 3 == JML.LConst(3)._1
    @test 1.0 == JML.LConst(1.0)._1
    @test nothing == JML.LConst(nothing)._1
    err1 = JML.SimpleMessage("aaaaa")
    err2 = JML.ModulePathNotFound("a.pml", "b")
    err3 = JML.NameUnresolved("var")
    err4 = JML.Positioned(1, 2, err3)
    err5 = JML.ComposeExc(err1, JML.ComposeExc(err2, err4))
    JML.print_exc(stdout, err5)
end

@testset "JML.jl" begin
    ENV["RUPYPATH"] = "./path"
    println(codegen(raw"""
       module S where
       import A
       def y = 1 `fst` 2 `fst` 3
       """, "a", :py))
    # Write your own tests here.
end

max_forks_repo_name: thautwarm/Rupy


text:
abstract type Spectrum end
abstract type Multispectral <: Spectrum end # Multispectral types
abstract type CMatch <: Multispectral end #Color Matching Functions (x_bar, y_bar, z_bar)
abstract type ConeFund <: Multispectral end # Cone fundamentals (l, m, s)
abstract type COpp <: Multispectral end # Opponent color channels
# abstract type MSpec <: Multispectral end # Generalized multispectral type, maybe implemented later

"""
    LSpec, Luminance Spectrum

Fields:

- `λ::Vector{Real}`, vector containing wavelengths
- `l::Vector{Real}`, vector containing luminance values
"""
struct LSpec <: Spectrum
    λ::Vector{Real}  #wavelength vector
    l::Vector{Real}  #reflectance vector
    function LSpec(λ, l)
        new(λ, l)
    end
end

"""
    RSpec, Reflectance Spectrum

Fields:

- `λ::Vector{Real}`, vector containing wavelengths
- `r::Vector{Real}, vector containing reflectanceS values
"""
struct RSpec <: Spectrum
    λ::Vector{Real}  #wavelength vector
    r::Vector{Real}  #reflectance vector
    function RSpec(λ, r)
        new(λ, r)
    end
end

"""
    TSpec, Transmittance spectrum

Fields:

- `λ::Vector{Real}`, vector containing wavelengths
- `t::Vector{Real}`, vector containing transmittance values
- `x::Real`, unit thickness value
"""
struct TSpec <: Spectrum
    λ::Vector{Real}  #start wavelength
    t::Vector{Real}  #transmittance vector at unit thickness
    x::Real          #thickness in terms of unit thickness
    function TSpec(λ, t, x)
        new(λ, t, x)
    end
end

#=  Multispectral types
    These include color matching functions and cone fundamentals defined by the CIE
=#

"""
    CIE31

Type of CIE 1931 2° observer color matching function.
"""
struct CIE31 <: CMatch
    λ::Vector{Real} #wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE31(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    LMS31

Type of CIE 1931 2° observer cone fundamentals.
"""
struct LMS31 <: ConeFund
    λ::Vector{Real}
    l::Vector{Real}
    m::Vector{Real}
    s::Vector{Real}
    function LMS31(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    CIE31_J

Type of CIE 1931 2° observer color matching function with Judd corrections.
"""
struct CIE31_J <: CMatch
    λ::Vector{Real} #wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE31_J(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    CIE31_JV

Type of CIE 1931 2° observer color matching function with corrections by Judd and Vos.
"""
struct CIE31_JV <: CMatch
    λ::Vector{Real} #wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE31_JV(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    CIE31_J

Type of CIE 1964 10° observer color matching function.
"""
struct CIE64 <: CMatch
    λ::Vector{Real} #wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE64(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    LMS64

Type of CIE 1964 10° observer cone fundamentals.
"""
struct LMS64 <: ConeFund
    λ::Vector{Real}
    l::Vector{Real}
    m::Vector{Real}
    s::Vector{Real}
    function LMS64(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    CIE12_2

Type of CIE 2012 2° observer color matching function. Current CIE 2° standard observer.
"""
struct CIE12_2 <: CMatch
    λ::Vector{Real} # wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE12_2(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    CIE12_10

Type of CIE 2012 10° observer color matching function. Current CIE 10° standard observer.
"""
struct CIE12_10 <: CMatch
    λ::Vector{Real} # wavelength vector
    x::Vector{Real} #x_bar cone matching function
    y::Vector{Real} #y_bar cone matching function
    z::Vector{Real} #z_bar cone matching function
    function CIE12_10(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    LMS06_2

Type of CIE 2006 2° observer cone fundamentals.
"""
struct LMS06_2 <: ConeFund
    λ::Vector{Real} #wavelength vector
    l::Vector{Real} #l cone response function
    m::Vector{Real} #m cone response function
    s::Vector{Real} #s cone response function
    function LMS06_2(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    LMS06_10

Type of CIE 2006 10° observer cone fundamentals.
"""
struct LMS06_10 <: ConeFund
    λ::Vector{Real} #wavelength vector
    l::Vector{Real} #l cone response function
    m::Vector{Real} #m cone response function
    s::Vector{Real} #s cone response function
    function LMS06_10(λ, l, m, s)
        new(λ, l, m, s)
    end
end

"""
    VYBRG

Type of Ingling-Tsou color opponent space functions with `Vλ` (lightness), `yb` (yellow-blue), `rg` (red-green) color channels.
"""
struct VYBRG <: COpp
    λ::Vector{Real}
    Vλ::Vector{Real}
    yb::Vector{Real}
    rg::Vector{Real}
    function VYBRG(λ, Vλ, yb, rg)
        new(λ, Vλ, yb, rg)
    end
end


text:
const WORLDCLIM_URL = "https://biogeo.ucdavis.edu/data/worldclim/v2.1"

rasterpath(::Type{WorldClim{T}}) where T =
    joinpath(rasterpath(), "WorldClim", string(nameof(T)))
rasterpath(T::Type{<:WorldClim}, layer) = joinpath(rasterpath(T), string(layer))
rasterpath(T::Type{<:WorldClim}, layer, x) =
    joinpath(rasterpath(T, layer), rastername(T, layer, x))

_zipfile_to_read(raster_name, zf) = first(filter(f -> f.name == raster_name, zf.files))


text:
using StaticArrays
using JuLIP, JuLIP.Potentials
using JuLIP.Potentials: site_energy, site_energy_d
using JuLIP: AbstractCalculator

function evec(i)
   if i == 1
      return SVector(1.0,0.0,0.0)
   elseif i == 2
      return SVector(0.0,1.0,0.0)
   elseif i == 3
      return SVector(0.0,0.0,1.0)
   end
   @error("`evec`: input must be 1,2 or 3")
end

function eval_dataset!(D::Dict, calc::AbstractCalculator; key = "train")
    at = Atoms(D["at"])
    for dat in D["data"]   # TODO: parallelise this loop!
       # get some info about the local lsq system
       dt = dat["datatype"]
       # evaluate this data point
       val = eval_dat(Val(Symbol(dt)), dat, calc, at)
       dat[key] = val
    end
    return D
end

"""
`eval_dat(::Val{SYM}, dat, calc, at) -> Vector{Float64}`

The following `SYM`s  are registered so far:
- `:E` : total energy
- `:F` : forces
- `:V` : virial
- `:Es` : site energy
- `:dEs` : site energy gradient
- `:d2Esh` : site energy hessian, FD approximation
- `:FC` : force-constants, i.e. a representation of the hessian
"""
function eval_dat end

eval_dat(::Val{:Es}, dat, calc, at) =
        site_energy(calc, at, 1)

eval_dat(::Val{:dEs}, dat, calc, at) =
        site_energy_d(calc, at, 1)

eval_dat(::Val{:d2Esh}, dat, calc, at) =
        _d2Esh(calc, at, 1, dat["l"], dat["i"], dat["h"])

eval_dat(::Val{:E}, dat, calc, at) =
        energy(calc, at)

eval_dat(::Val{:F}, dat, calc, at) =
        collect(mat(forces(calc, at))[:])

eval_dat(::Val{:V}, dat, calc, at) =
        virial(calc, at)[SVector(1,2,3,5,6,9)]

eval_dat(::Val{:FC}, dat, calc, at) =
        _force_constants(calc, at, dat["i0"], dat["h"], dat["i"])

# -----------------------------------------------------------------------------
#   SPECIAL EVALUATION ROUTINES FOR TB/QM MODELS
# -----------------------------------------------------------------------------

function eval_dataset_tb!(D::Dict, calc::AbstractCalculator; key = "train")
    at = Atoms(D["at"])
    DTs = [ dat["datatype"]  for dat in D["data"] ]
    for dt in unique(DTs)
        Idt = findall(DTs .== dt)
        eval_dataset_tb!(Val(Symbol(dt)), D.data[Idt], calc, at; key=key)
    end
    return D
end

function eval_dataset_tb!(valdt::Union{Val{:Es}, Val{:dEs}}, data, calc, at; key=key)
    @assert length(data) == 1
    data[1][key] = eval_dat(valdt, data[1], calc, at)
    return nothing
end

function eval_dataset_tb!(::Val{:d2Esh}, data, calc, at; key=key)
    h = data[1]["h"]
    # TODO: test all h are the same!
    for i = 1:3, sig in [1, -1]
        # perturb X

        # solve the eval problem

        # compute dEs on all neighbours we care about
        # as determined by data sketch
        for dat in data
            if dat["i"] == i
                l = dat["l"]
                @assert dat["h"] == h
                # compute dEs

                # write dEs into the data point
                calc[key] += dEs * sig
            end
        end
    end
    for dat in data
        calc[key] /= (2*h)
    end
end

"""
l0 : site that we are perturbing
 l : site at which the site energy is evaluated
 i : direction of the perturbation (E1, E2, E3)
 h : finite-difference step
"""
function _d2Esh(calc, at::Atoms,
                l0::Integer, l::Integer, i::Integer, h::Real)
   X = positions(at)
   X[l0] += h * evec(i)
   set_positions!(at, X)
   dVp = site_energy_d(calc, at, l)
   X[l0] -= 2*h * evec(i)
   set_positions!(at, X)
   dVm = site_energy_d(calc, at, l)
   X[l0] += h * evec(i)
   set_positions!(at, X)
   return (dVp - dVm) / (2*h)
end

function _force_constants(calc, at::Atoms, i0, h, i)
    at[i0] += h * evec(i)
    Fp = forces(calc, at)
    at[i0] -= 2 * h * evec(i)
    Fm = forces(calc, at)
    at[i0] +=  h * evec(i0)
    return  mat( (Fp - Fm) / (2*h) )[:] |> collect
end


text:
push!(LOAD_PATH, homedir()*"/Projects/QuIPS/src")

# Quantum Virtual Machine Library

module QVM

export QCircuit, run!, step!, reset!, operate!

using Operators

const C32 = Complex{Float32}

mutable struct QCircuit
    wfn::Vector{C32}
    ops::Vector{Operator}
    out::BitVector
    pos::Int
    N::Int

    function QCircuit(N::Int)
        wfn = zeros(C32, 2^N); wfn[1] = 1
        out = BitVector(undef, N)
        new(wfn, [], out, 0, N)
    end

    function QCircuit(quip::Vector, N::Int)
        wfn = zeros(C32, 2^N); wfn[1] = 1
        ops = compile(quip)
        out = BitVector(undef, N)
        new(wfn, ops, out, 0, N)
    end
end

function compile(quip::Vector, verbose=true)
    println("\ncompiling quip to circuit\n")
    Vs = []
    for tag in quip
        if tag[1] == :MEASURE || tag[1] == :M
            push!(Vs, Measurement(tag[2]))
        elseif String(tag[1])[1] == 'C'
            push!(Vs, ControlGate(tag))
        else
            push!(Vs, Gate(tag))
        end
        if verbose
            V = Vs[end]
            print("  compiling "); show(V)
        end
    end
    println("\ncompiled\n")
    Vs
end

function run!(QC::QCircuit, verbose=true)
    if verbose
        println("\nbeginning run\n")
    end
    for V in QC.ops
        QC.pos += 1
        evolve!(QC, V)
        if verbose
            if typeof(V) == Measurement
                print("  "); show(V)
                println("   -> ", QC.out[V.k])
            else
                print("  "); show(V)
            end
        end
    end
    if verbose
        println("\nfinished run\n")
    end
end

function step!(QC::QCircuit)
    if QC.pos < length(QC.ops)
        QC.pos += 1
        V = QC.ops[QC.pos]
        evolve!(QC, V)
    end
end

function reset!(QC::QCircuit)
    QC.wfn = zero.(QC.wfn); QC.wfn[1] = 1
    QC.out = BitVector(undef, QC.N)
    QC.pos = 0
end

function evolve!(QC::QCircuit, V::Operator)
    if typeof(V) == Measurement
        ψ, outcome = V(QC.wfn, QC.N)
        QC.out[V.k] = outcome
        QC.wfn .= ψ
    else
        ψ = V(QC.wfn, QC.N)
        QC.wfn .= ψ
    end
end

function operate!(QC::QCircuit, V::Operator)
    insert!(QC.ops, QC.pos + 1, V)
    step!(QC)
end

operate!(QC::QCircuit, Vs::Vector) = for V in Vs operate!(QC, V) end

end

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
using SQLite

db = SQLite.DB("pinglog-sqlite3.db")
df = SQLite.query(db,"SELECT date, dest, ping FROM pinglog;")

df[:date] = map(x->Dates.unix2datetime(parse(Float32, x)), df[:date])
df[:ping] = map(x->parse(Float32, x), df[:ping])
df = df[df[:date] .> DateTime("2017-12-08"), :]
pings_to_google = df[df[:dest] .== "google.com", [:date, :ping]]
pings_to_router = df[df[:dest] .== "192.168.0.1", [:date, :ping]]

using PlotlyJS
using Plots
plotlyjs()

p1 = Plots.plot(
    pings_to_google[:date],
    pings_to_google[:ping],
    title="Ping logs from rpi3",
    label="google.com",
    color="red")
p2 = Plots.plot(
    pings_to_router[:date],
    pings_to_router[:ping],
    xlabel="Time [DateTime]",
    label="router",
    color="blue")
Plots.plot(p1, p2, layout=(2,1),
    link=:x,
    ylabel="Ping [ms]")

max_stars_repo_name: cmey/pinglog

max_issues_repo_name: cmey/pinglog

max_forks_repo_name: cmey/pinglog


text:
module exclusionprocess

export ExclusionProcess, update!

import ..epstate: EPState
import ..boundarycondition: BoundaryCondition

abstract type ExclusionProcess end

"""
Generic update function. Needs methods added for each process and boundary condition.
"""
update!(::Type{S}, ::Type{T}, state::EPState) where {S<:ExclusionProcess, T<:BoundaryCondition}  =
throw(ArgumentError("Update not defined for exclusion process $S with boundary condition $T."))

abstract type Update end

end

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
struct MeshPointNM{U,D,C,N,T}
    patch::Simplex{U,D,C,N,T}
    bary::SVector{D,T}
    cart::SVector{U,T}
end

Base.length(m::MeshPointNM) = length(m.cart)
Base.getindex(p::MeshPointNM, i::Int) = p.cart[i]

"""
    cartesian(neighborhood) -> point

Return the cartesian coordinates of the point where `neighborhood` is located.
"""
function cartesian(mp::MeshPointNM)
    mp.cart
end

"""
    parametric(neighborhood) -> point

Return the parameters where `neighborhood` is located.
"""
function parametric(mp::MeshPointNM)
    mp.bary
end

chart(nbd::MeshPointNM) = mp.patch

"Return the barycentric coordinates of `mp`"
barycentric(mp::MeshPointNM) = SVector(mp.bary[1], mp.bary[2], 1-mp.bary[1]-mp.bary[2])

"""
A number defines a neighborhood in euclidian space
"""
jacobian(x::Number) = one(x)

jacobian(mp::MeshPointNM) = volume(mp.patch) * factorial(dimension(mp.patch))

"""
    tangents(neighborhod, i) -> tangent_i

Return the i-th tangent vector at the neighborhood.
"""
function tangents(mp::MeshPointNM, i)
    mp.patch.tangents[i]
end

"""
    normal(neighborhood)

Return the normal at a neighborhood on a surface.
"""
function normal(mp::MeshPointNM)
    mp.patch.normals[1]
end

"""
    neighborhood(chart, params)

Create a neighborhood from a chart and a set of parameter values.
"""
function neighborhood(p::Simplex, bary)
  D = dimension(p)
  T = coordtype(p)
  P = SVector{D,T}
  cart = barytocart(p, bary)
  MeshPointNM(p, P(bary), cart)
end

"""
    center(simplex) -> neighborhood

Create the neighborhood at the center of the simplex, i.e. the point corresponding
to parameter `(1/(D+1), 1/(D+1), ...)` where `D` is the simplex dimension.
"""
@generated function center(p::Simplex{U,D,C,N,T}) where {U,D,C,N,T}
    uv = ones(T,D)/(D+1)
    :(neighborhood(p, $uv))
end


text:
"""
    function spearman_corr(;
        players::Array{Player},
        ranking_scheme1::String,
        ranking_scheme2::String
    )::Float64

Calculates the Spearman Correlation between two ranks.

Arguments
=========
- `players`:
    the players whose rankings to consider.
- `ranking_scheme1`:
    ranking scheme 1 to use.
- `ranking_scheme2`:
    ranking scheme 2 to use.
"""

function spearman_corr(;
    players::Array{Player},
    ranking_scheme1::String,
    ranking_scheme2::String
)::Float64
    rankings1 = map(player -> player.rankings[ranking_scheme1], players)
    rankings2 = map(player -> player.rankings[ranking_scheme2], players)
    return corspearman(rankings1, rankings2)
end

"""
    function kendtau_corr(;
        players::Array{Player},
        ranking_scheme1::String,
        ranking_scheme2::String
    )::Float64

Calculates the Kendall Tau Correlation between two ranks.

Arguments
=========
- `players`:
    the players whose rankings to consider.
- `ranking_scheme1`:
    ranking scheme 1 to use.
- `ranking_scheme2`:
    ranking scheme 2 to use.
"""

function kendtau_corr(;
    players::Array{Player},
    ranking_scheme1::String,
    ranking_scheme2::String
)::Float64
    rankings1 = map(player -> player.rankings[ranking_scheme1], players)
    rankings2 = map(player -> player.rankings[ranking_scheme2], players)
    return corkendall(rankings1, rankings2)
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking