text:
module PWAHelper

using amplitudes_compass

using DelimitedFiles
using Juno

export precalculate_compass_basis, read_precalc_basis
# export precalculate_compass_basis_txt, read_precalc_basis_txt
export get_npars
export get_parameter_map, make_pblock_inds
export extnd, extnd!, shrnk, cohsq, cohts, cohts!
export contract_to_intensity, get_intesity
export get_parameter_ranges, normalize_pars!
export calculate_integrals_from_precalc_basis

function precalculate_compass_basis(basis,fin,fout)
    io = open(fin,"r")
    Nd = read(io,Int32)
    variable_mat = Array{Float64}(undef,6,Nd)
    read!(io,variable_mat)
    #variable_mat = readdlm(fin); Nd = size(variable_mat,1)
    function_mat = fill(0.0im, Nd, length(basis))
    @progress for e in 1:Nd
        for (i,b) in enumerate(basis)
            function_mat[e,i] = b(@view(variable_mat[2:end,e])..., variable_mat[1,e]);
        end
    end
    io = open(fout,"w")
    write(io,trunc(Int64, Nd))
    write(io,trunc(Int64, length(basis)+length(basis)))
    write(io,[real(function_mat) imag(function_mat)])
    close(io)
end

function read_precalc_basis(fname)
    io = open(fname,"r")
    r=read(io,Int64)
    c=read(io,Int64)
    ld = Array{Float64}(undef,r,c)
    read!(io,ld)
    Nh = div(size(ld,2),2)
    return ld[:,1:Nh]+1im*ld[:,(Nh+1):end]
end

function get_parameter_map(block_inds, Nw)
    temp = []; numb = []
    for (i,bl) in enumerate(block_inds)
        # push false for the first parameter, true for others
        push!(temp,false,[true for i=1:(length(bl)-1)]...);
        push!(numb,bl...)
    end
    Tmap = fill(0,2,sum(temp .+ 1))
    count=1
    for (i,b) in enumerate(temp)
        Tmap[1,count] = numb[i]
        count+=1
        if b
            Tmap[2,count] = numb[i];
            count+=1
        end
    end
    return Tmap
end

function get_npars(block_inds)
    sum(2*length(bl)-1 for bl in block_inds)
end

function make_pblock_inds(block_inds)
    Nps = [2*length(bl)-1 for bl in block_inds]
    counter = 1
    pbls = [let v = collect(counter:(counter+Np-1))
                counter += Np
                v
            end for Np in Nps]
    return pbls
end

extnd(Ψ, tmap)  = [((tmap[2,i]==0) ? Ψ[tmap[1,i]] : 1im*Ψ[tmap[2,i]]) for i in 1:size(tmap,2)]
function extnd!(X, Ψ, tmap)
    for i in 1:size(tmap,2)
        @inbounds X[i] = (tmap[2,i]==0) ? Ψ[tmap[1,i]] : 1im*Ψ[tmap[2,i]]
    end
end
function shrnk(p, tmap)
    Nw = tmap[2,end] # work around
    outv = fill(0.0im,Nw)
    for i in 1:size(tmap,2)
        (tmap[1,i] != 0) && (outv[tmap[1,i]] += p[i]);
        (tmap[2,i] != 0) && (outv[tmap[2,i]] += 1.0im*p[i]);
    end
    outv
end

function cohsq(X, pblocks)
    sum(abs2, (sum(@view(X[bl]))) for bl in pblocks)
end
function cohts(X, pblocks)
    Np = length(X)
    vblocks = [let v = fill(0.0,Np)
        v[bl] .= 1.0
        v
    end for bl in pblocks]
    sum(sum(X[bl])*pl for (pl,bl) in zip(vblocks, pblocks))
end
function cohts!(X, pblocks)
    # Np = length(X)
    # vblocks = [
    #     let v = fill(0.0,Np)
    #         v[bl] .= 1.0
    #         v
    #     end for bl in pblocks]
    # sum(sum(X[bl])*pl for (pl,bl) in zip(vblocks, pblocks))
    for bl in pblocks
        @inbounds X[bl] .= sum(@view(X[bl]))
    end
end

function contract_to_intensity(ΨΨstar, block_inds)
    Tmap = get_parameter_map(block_inds, size(ΨΨstar,1))
    Np = size(Tmap,2)
    pblocks = make_pblock_inds(block_inds)
    COH = fill(0.0, Np, Np)
    for bl in pblocks
        COH[bl,bl] .= 1.0
    end
    get(v,w) = (v==0 || w==0) ? 0.0im : ΨΨstar[v,w]
    [COH[i,j]*( # meaning is ([1]-i[2]) * ([1]+i[2])
            get(Tmap[1,i],Tmap[1,j]) +
            get(Tmap[2,i],Tmap[2,j]) -
        1im*get(Tmap[2,i],Tmap[1,j]) +
        1im*get(Tmap[1,i],Tmap[2,j])
            ) for j in 1:Np, i in 1:Np]
end

function get_intesity(pars, form, block_inds)
    pblocks = make_pblock_inds(block_inds)
    BM = real.(contract_to_intensity(form,block_inds))
    Np = length(pars)
    sum(pars[i]*BM[i,j]*pars[j] for i=1:Np, j=1:Np)
end

function get_parameter_ranges(form, block_inds)
    Np = get_npars(block_inds)
    [let p = fill(0.0, Np)
        p[i] = 1.0
        1/sqrt(get_intesity(p, form, block_inds))
     end for i in 1:Np]
end

function normalize_pars!(pars, form, block_inds)
    Intens = get_intesity(pars, form, block_inds)
    pars ./=  sqrt(Intens)
    pars
end

end


text:
function FieldHansenTesseralLinear(arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::Field)
    return FieldHansenTesseralLinear((jint, jint, jint, jint, jint, Field), arg0, arg1, arg2, arg3, arg4, arg5)
end

function compute_init_values(obj::FieldHansenTesseralLinear, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(obj, "computeInitValues", void, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_derivative(obj::FieldHansenTesseralLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getDerivative", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end

function get_value(obj::FieldHansenTesseralLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getValue", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineAssessPatchesResult <: SwaggerModel
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    assessmentActivityId::Any # spec type: Union{ Nothing, String } # spec name: assessmentActivityId
    rebootPending::Any # spec type: Union{ Nothing, Bool } # spec name: rebootPending
    criticalAndSecurityPatchCount::Any # spec type: Union{ Nothing, Int32 } # spec name: criticalAndSecurityPatchCount
    otherPatchCount::Any # spec type: Union{ Nothing, Int32 } # spec name: otherPatchCount
    startDateTime::Any # spec type: Union{ Nothing, DateTime } # spec name: startDateTime
    patches::Any # spec type: Union{ Nothing, Vector{VirtualMachineSoftwarePatchProperties} } # spec name: patches
    error::Any # spec type: Union{ Nothing, ApiError } # spec name: error

    function VirtualMachineAssessPatchesResult(;status=nothing, assessmentActivityId=nothing, rebootPending=nothing, criticalAndSecurityPatchCount=nothing, otherPatchCount=nothing, startDateTime=nothing, patches=nothing, error=nothing)
        o = new()
        validate_property(VirtualMachineAssessPatchesResult, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("assessmentActivityId"), assessmentActivityId)
        setfield!(o, Symbol("assessmentActivityId"), assessmentActivityId)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("rebootPending"), rebootPending)
        setfield!(o, Symbol("rebootPending"), rebootPending)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("criticalAndSecurityPatchCount"), criticalAndSecurityPatchCount)
        setfield!(o, Symbol("criticalAndSecurityPatchCount"), criticalAndSecurityPatchCount)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("otherPatchCount"), otherPatchCount)
        setfield!(o, Symbol("otherPatchCount"), otherPatchCount)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("startDateTime"), startDateTime)
        setfield!(o, Symbol("startDateTime"), startDateTime)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("patches"), patches)
        setfield!(o, Symbol("patches"), patches)
        validate_property(VirtualMachineAssessPatchesResult, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type VirtualMachineAssessPatchesResult

const _property_map_VirtualMachineAssessPatchesResult = Dict{Symbol,Symbol}(Symbol("status")=>Symbol("status"), Symbol("assessmentActivityId")=>Symbol("assessmentActivityId"), Symbol("rebootPending")=>Symbol("rebootPending"), Symbol("criticalAndSecurityPatchCount")=>Symbol("criticalAndSecurityPatchCount"), Symbol("otherPatchCount")=>Symbol("otherPatchCount"), Symbol("startDateTime")=>Symbol("startDateTime"), Symbol("patches")=>Symbol("patches"), Symbol("error")=>Symbol("error"))
const _property_types_VirtualMachineAssessPatchesResult = Dict{Symbol,String}(Symbol("status")=>"String", Symbol("assessmentActivityId")=>"String", Symbol("rebootPending")=>"Bool", Symbol("criticalAndSecurityPatchCount")=>"Int32", Symbol("otherPatchCount")=>"Int32", Symbol("startDateTime")=>"DateTime", Symbol("patches")=>"Vector{VirtualMachineSoftwarePatchProperties}", Symbol("error")=>"ApiError")
Base.propertynames(::Type{ VirtualMachineAssessPatchesResult }) = collect(keys(_property_map_VirtualMachineAssessPatchesResult))
Swagger.property_type(::Type{ VirtualMachineAssessPatchesResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineAssessPatchesResult[name]))}
Swagger.field_name(::Type{ VirtualMachineAssessPatchesResult }, property_name::Symbol) =  _property_map_VirtualMachineAssessPatchesResult[property_name]

const _allowed_VirtualMachineAssessPatchesResult_status = ["InProgress", "Failed", "Succeeded", "CompletedWithWarnings"]

function check_required(o::VirtualMachineAssessPatchesResult)
    true
end

function validate_property(::Type{ VirtualMachineAssessPatchesResult }, name::Symbol, val)
    if name === Symbol("status")
        Swagger.validate_param(name, "VirtualMachineAssessPatchesResult", :enum, val, _allowed_VirtualMachineAssessPatchesResult_status)
    end
end


text:
#@testset "BigMat.jl" begin
    # tests for MPT_SVD, MPT_SVD, contract
    # vector to MPS
    v1      = [1.0; 4.5; 2; 7; 3; 7; 3.3; 10; 6; 8; 1.1; 4.6; 1; 6; 2; 7; 12; 4.5; 1; 0; 0; 3; 6; 2.1];
    v2      = [1.0; 4.5; 2; 7; 3; 7; 3.3; 10; 6; 8; 1.1; 4.6; 1; 6; 2; 7; 12; 4.5; 1; 0; 0; 3; 6; 2.1];
    M1      = randn(64,64)
    M2      = Symmetric(M1*M1')

    tt1,err = MPT_SVD(v1,[2 3 4],0.07);
    tt2     = TT_ALS(v2,[2 3 4],[1,2,4,1]);
    ttm1    = TTm_ALS(M1,[4 4 4; 4 4 4],[1,16,16,1]);
    ttm2    = TTm_ALS(M2,[4 4 4; 4 4 4],[1,16,16,1])
    ttm2    = TTm_ALS(Symmetric(M1*M1'),[4 4 4; 4 4 4],[1,3,3,1])

    @test (norm(tt1) - norm(tt2)) < 1e-13
    @test norm(mps2vec(tt1) - v1)/norm(v1) < 1e-13
    @test norm(mpo2mat(ttm1) - M1)/norm(M1) < 1e-13
    @test norm(mpo2mat(ttm2) - mpo2mat(ttm2)') < 1e-10

    # matrix to MPO
    M1 = v1*v1';
    M2 = v2*v2';
    mpo1,err = MPT_SVD(M1,[2 3 4; 2 3 4],0.0);
    mpo2 = MPT_ALS(M2,[2 3 4; 2 3 4],[1,4,16,1]);
    @test (norm(mpo1) - norm(mpo2)) < 1e-14

    @run mpo2mat(mpo1)
    @test norm(mpo2mat(mpo1) - M1)/norm(M1) < 1e-14
    @test norm(mpo1)-norm(mpo2mat(mpo1)) < 1e-14

    tt3 = MPT_SVD(v1,[2 3 4],0.07);
    @test norm(mps2vec(tt3)-v1)/norm(v1) <= 0.07

    # TT-ALS without orthogonalization needs to be checked. Wihtout pinv, the error is very high
    # maxiter + residulas need to be implemented

    # test of rktp2tn
    Ũ = [rand(50,40),rand(50,40),rand(50,40)];
    ϵ = 1.0;
    Umpo,err2 = rkrp2tn(Ũ,ϵ);
    Uapprox = mpo2mat(Umpo);
    U12 = zeros(2500,40);
    @inbounds @simd for j = 1:40
        @views kron!(U12[:,j],Ũ[1][:,j],Ũ[2][:,j])
    end
    Utrue = zeros(2500*50,40);
    @inbounds @simd for j = 1:40
        @views kron!(Utrue[:,j],U12[:,j],Ũ[3][:,j])
    end
    @test norm(Utrue-Uapprox)/norm(Utrue) < ϵ

    # get matrix entry from mpo2mat
    Ktest = zeros(24,24);
    linear = LinearIndices((1:2, 1:3, 1:4, 1:2, 1:3, 1:4));
    for i1 = 1:2
        for i2 = 1:3
            for i3 = 1:4
                for j1 = 1:2
                    for j2 = 1:3
                        for j3 = 1:4
                            tmp = mpo2[1][:,i1,j1,:]*mpo2[2][:,i2,j2,:]*mpo2[3][:,i3,j3,:];
                            Ktest[linear[i1,i2,i3,j1,j2,j3]] = tmp[1];
                        end
                    end
                end
            end
        end
    end
    @test norm(Ktest-M2)/norm(M2) < 1e-10

    # TT-SVD for symmetrical matrix
    X, y, f, K  = gensynthdata(1024,1,0.01,1.0,1.0)
    Kttm,Lttm,Uttm,S,err = MPT_SVD(Symmetric(K),[4 4 4 4 4; 4 4 4 4 4],0.1)
    k = mps2vec(Kttm)
    Kt = reshape(k,(4,4,4,4,4,4,4,4,4,4))
    Kt = permutedims(Kt,[1 2 3 4 5 10 9 8 7 6])
    Kapprox = reshape(Kt,(1024, 1024))
    norm(K[:]-Kt[:])/norm(K[:])
    rank(Kttm)

    L  = (Kttm[1]*Kttm[2]*Kttm[3]*Kttm[4]*Kttm[5])[1,:,:]
    L2 = mpo2mat(Lttm)
    norm(L*L' - K)/norm(K)
    norm(L2*L2' - K)/norm(K)
    U = mpo2mat(Uttm)
    norm(U*Diagonal(S)*U'-K)/norm(K)

    # test Hierarchical Tucker function
    tensor = rand(5,5,5);
    ranks  = [2,2,2];
    @run leaves2roottrunc(tensor,ranks)

#end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
"""
A file with structs and functions for sampling reactions and updating reaction rates in spatial SSAs
"""

### spatial rx rates ###
struct RxRates{F,M}
    "rx_rates[i,j] is rate of reaction i at site j"
    rates::Matrix{F}

    "rx_rates_sum[j] is sum of reaction rates at site j"
    sum_rates::Vector{F}

    "AbstractMassActionJump"
    ma_jumps::M
end

"""
    RxRates(num_sites::Int, ma_jumps::M) where {M}

initializes RxRates with zero rates
"""
function RxRates(num_sites::Int, ma_jumps::M) where {M}
    numrxjumps = get_num_majumps(ma_jumps)
    rates = zeros(Float64, numrxjumps, num_sites)
    RxRates{Float64,M}(rates, vec(sum(rates, dims=1)), ma_jumps)
end

num_rxs(rx_rates::RxRates) = get_num_majumps(rx_rates.ma_jumps)

"""
    reset!(rx_rates::RxRates)

make all rates zero
"""
function reset!(rx_rates::RxRates)
    fill!(rx_rates.rates, zero(eltype(rx_rates.rates)))
    fill!(rx_rates.sum_rates, zero(eltype(rx_rates.sum_rates)))
    nothing
end

"""
    total_site_rx_rate(rx_rates::RxRates, site)

return total reaction rate at site
"""
function total_site_rx_rate(rx_rates::RxRates, site)
    @inbounds rx_rates.sum_rates[site]
end

"""
    update_rx_rates!(rx_rates, rxs, u, site)

update rates of all reactions in rxs at site
"""
function update_rx_rates!(rx_rates::RxRates{F,M}, rxs, u, site) where {F, M <: MassActionJump}
    ma_jumps = rx_rates.ma_jumps
    @inbounds for rx in rxs
        set_rx_rate_at_site!(rx_rates, site, rx, evalrxrate((@view u[:,site]), rx, ma_jumps))
    end
end

function update_rx_rates!(rx_rates::RxRates{F,M}, rxs, u, site) where {F, M <: SpatialMassActionJump}
    ma_jumps = rx_rates.ma_jumps
    @inbounds for rx in rxs
        set_rx_rate_at_site!(rx_rates, site, rx, evalrxrate(u, rx, ma_jumps, site))
    end
end

"""
    sample_rx_at_site(rx_rates::RxRates, site, rng)

sample a reaction at site, return reaction index
"""
function sample_rx_at_site(rx_rates::RxRates, site, rng)
    linear_search((@view rx_rates.rates[:,site]), rand(rng) * total_site_rx_rate(rx_rates, site))
end

# helper functions
function set_rx_rate_at_site!(rx_rates::RxRates, site, rx, rate)
    @inbounds old_rate = rx_rates.rates[rx, site]
    @inbounds rx_rates.rates[rx, site] = rate
    @inbounds rx_rates.sum_rates[site] += rate - old_rate
    old_rate
end

function Base.show(io::IO, ::MIME"text/plain", rx_rates::RxRates)
    num_rxs, num_sites = size(rx_rates.rates)
    println(io, "RxRates with $num_rxs reactions and $num_sites sites")
end


text:
function vc_integral(
            vc::LogisticSingle{FT},
            p_dos::FT,
            p_ups::FT
) where {FT<:AbstractFloat}
    @assert p_ups <= 0 && p_dos <= 0;

    # unpack data from VC
    @unpack a,b = vc;
    _t_dos = log(a * exp(b*p_dos) + 1);
    _t_ups = log(a * exp(b*p_ups) + 1);

    return (a+1) / (a*b) * (_t_ups - _t_dos)
end

function vc_integral(
            vc::LogisticSingle{FT},
            p_dos::FT,
            p_ups::FT,
            h::FT,
            E::FT,
            Kmax::FT
) where {FT<:AbstractFloat}
    @assert p_ups <= 0 && p_dos <= 0;

    # unpack data from VC
    @unpack a,b = vc;
    _krghe = Kmax * ρg_MPa(FT) * h * (a+1) / a + E;
    _lower = b * _krghe;
    _multi = Kmax * (a+1) / a * E;
    _upper_dos = log(a * _krghe * exp(b*p_dos) + E);
    _upper_ups = log(a * _krghe * exp(b*p_ups) + E);

    return _multi * (_upper_ups - _upper_dos) / _lower
end

function vc_integral(
            vc::PowerSingle{FT},
            p_dos::FT,
            p_ups::FT
) where {FT<:AbstractFloat}
    @assert p_ups <= 0 && p_dos <= 0;

    # unpack data from vc
    @unpack a,b = vc;
    _f_dos = p_dos * _₂F₁(1, 1/b, 1+1/b, -a*(-p_dos)^b);
    _f_ups = p_ups * _₂F₁(1, 1/b, 1+1/b, -a*(-p_ups)^b);

    return _f_ups - _f_dos
end

function vc_integral(
            vc::PowerSingle{FT},
            p_dos::FT,
            p_ups::FT,
            h::FT,
            E::FT,
            Kmax::FT
) where {FT<:AbstractFloat}
    @assert p_ups <= 0 && p_dos <= 0;

    # unpack data from VC
    @unpack a,b = vc;
    _krghe = Kmax * ρg_MPa(FT) * h + E;
    _multi = Kmax * E / _krghe;
    _f_dos = p_dos * _₂F₁(1, 1/b, 1+1/b, -E * a * (-p_dos)^b / _krghe);
    _f_ups = p_ups * _₂F₁(1, 1/b, 1+1/b, -E * a * (-p_ups)^b / _krghe);

    return _multi * (_f_ups - _f_dos)
end

function vc_integral(
            vc::WeibullSingle{FT},
            p_dos::FT,
            p_ups::FT
) where {FT<:AbstractFloat}
    @assert p_ups <= 0 && p_dos <= 0;

    # unpack data from VC
    @unpack b,c = vc;

    # compute the incomplete gamma function
    _Γ_dos = gamma(1/c, (-p_dos/b)^c);
    _Γ_ups = gamma(1/c, (-p_ups/b)^c);

    return b/c * (_Γ_ups - _Γ_dos)
end

max_stars_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land


text:
export
    wncard,
    wncomd,
    wncond,
    wncond!,
    wndifd,
    wnelmd,
    wnexpd,
    wnexpd!,
    wnextd,
    wnextd!,
    wnfetd,
    wnfild,
    wnfild!,
    wnfltd,
    wnfltd!,
    wnincd,
    wninsd,
    wninsd!,
    wnintd,
    wnreld,
    wnsumd,
    wnunid,
    wnvald,
    wnvald!

"""
    wncard(window)

Return the cardinality (number of intervals) of a double precision window.

### Arguments ###

- `window`: Input window

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wncard_c.html)
"""
function wncard(window)
    ccall((:wncard_c, libcspice), SpiceInt, (Ref{Cell{SpiceDouble}},), window.cell)
end

"""
    wncomd(window, left, right)

Determine the complement of a double precision window with respect to a specified interval.

### Arguments ###

- `window`: Input window
- `left`:  Left endpoint of the complement interval
- `right`:  Right endpoint of the complement interval

### Output ###

Returns the complement of `window` with respect to `[left,right]`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wncomd_c.html)
"""
function wncomd(window, left, right)
    result = SpiceDoubleCell(window.cell.size + 1)
    ccall((:wncomd_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
          left, right, window.cell, result.cell)
    handleerror()
    result
end

"""
    wncond!(window, left, right)

Contract each of the intervals of a double precision window.

### Arguments ###

- `window`: Window to be contracted
- `left`:  Amount added to each left endpoint
- `right`: Amount subtracted from each right endpoint

### Output ###

Returns the contracted window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wncond_c.html)
"""
function wncond!(window, left, right)
    ccall((:wncond_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, Ref{Cell{SpiceDouble}}),
          left, right, window.cell)
    handleerror()
    window
end

@deprecate wncond wncond!

"""
    wndifd(a, b)

Place the difference of two double precision windows into a third window.

### Arguments ###

- `a`: Input window
- `b`: Input window

### Output ###

Returns a window containing the difference of `a` and `b`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wndifd_c.html)
"""
function wndifd(a, b)
    c = SpiceDoubleCell(2 * (wncard(a) + wncard(b)))
    ccall((:wndifd_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
          a.cell, b.cell, c.cell)
    handleerror()
    c
end

"""
    wnelmd(window, point)

Determine whether a point is an element of a double precision window.

### Arguments ###

- `window`: Input window
- `point`: Input point

### Output ###

Returns `true` if `point` is an element of `window`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnelmd_c.html)
"""
function wnelmd(window, point)
    res = ccall((:wnelmd_c, libcspice), SpiceBoolean, (SpiceDouble, Ref{Cell{SpiceDouble}}),
                point, window.cell)
    handleerror()
    Bool(res)
end

"""
    wnexpd(window, left, right)

Expand each of the intervals of a double precision window.

### Arguments ###

- `left`: Amount subtracted from each left endpoint
- `right`: Amount added to each right endpoint

### Output ###

Returns the expanded window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnexpd_c.html)
"""
function wnexpd!(window, left, right)
    ccall((:wnexpd_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, Ref{Cell{SpiceDouble}}),
          left, right, window.cell)
    handleerror()
    window
end

@deprecate wnexpd wnexpd!

"""
    wnextd!(window, side)

Extract the left or right endpoints from a double precision window.

### Arguments ###

- `window`: Window to be extracted
- `side`: Extract left (`:L`) or right (`:R`) endpoints

### Output ###

Returns the extracted window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnextd_c.html)
"""
function wnextd!(window, side)
    ccall((:wnextd_c, libcspice), Cvoid, (SpiceChar, Ref{Cell{SpiceDouble}}),
          first(string(side)), window.cell)
    handleerror()
    window
end

@deprecate wnextd wnextd!

"""
    wnfetd(window, n)

Fetch a particular interval from a double precision window.

### Arguments ###

- `window`: Input window
- `n`: Index of interval to be fetched

### Output ###

Returns a tuple consisting of the left and right endpoints of the n-th interval
in the input window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnfetd_c.html)
"""
function wnfetd(window, n)
    left = Ref{SpiceDouble}()
    right = Ref{SpiceDouble}()
    ccall((:wnfetd_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, SpiceInt, Ref{SpiceDouble}, Ref{SpiceDouble}),
          window.cell, n - 1, left, right)
    handleerror()
    left[], right[]
end

"""
    wnfild!(window, small)

Fill small gaps between adjacent intervals of a double precision window.

### Arguments ###

- `window`: Window to be filled
- `small`: Limiting measure of small gaps

### Output ###

Returns the updated window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnfild_c.html)
"""
function wnfild!(window, small)
    ccall((:wnfild_c, libcspice), Cvoid, (SpiceDouble, Ref{Cell{SpiceDouble}}),
          small, window.cell)
    handleerror()
    window
end

"""
    wnfild!(window, small)

Filter (remove) small intervals from a double precision window.

### Arguments ###

- `window`: Window to be filtered
- `small`: Limiting measure of small intervals

### Output ###

Returns the updated window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnfltd_c.html)
"""
function wnfltd!(window, small)
    ccall((:wnfltd_c, libcspice), Cvoid, (SpiceDouble, Ref{Cell{SpiceDouble}}),
          small, window.cell)
    handleerror()
    window
end

@deprecate wnfltd wnfltd!

"""
    wnincd(window, left, right)

Determine whether an interval is included in a double precision window.

### Arguments ###

- `window`: Input window
- `left`: Left endpoint of the input interval
- `right`: Right endpoint of the input interval

### Output ###

Returns `true` when `(left, right)` is contained in `window`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnincd_c.html)
"""
function wnincd(window, left, right)
    res = ccall((:wnincd_c, libcspice), SpiceBoolean,
                (SpiceDouble, SpiceDouble, Ref{Cell{SpiceDouble}}),
                left, right, window.cell)
    handleerror()
    Bool(res)
end

"""
    wninsd!(window, left, right)

Insert an interval into a double precision window.

### Arguments ###

- `window`: Input window
- `left`: Left endpoint of the new interval
- `right`: Right endpoint of the new interval

### Output ###

Returns the updated windows.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wninsd_c.html)
"""
function wninsd!(window, left, right)
    ccall((:wninsd_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, Ref{Cell{SpiceDouble}}),
          left, right, window.cell)
    handleerror()
    window
end

@deprecate wninsd wninsd!

"""
    wnintd(a, b)

Place the intersection of two double precision windows into a third window.

### Arguments ###

- `a`: Input window
- `b`: Input window

### Output ###

Returns a window containing the intersection of `a` and `b`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnintd_c.html)
"""
function wnintd(a, b)
    c = SpiceDoubleCell(2 * (wncard(a) + wncard(b)))
    ccall((:wnintd_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
          a.cell, b.cell, c.cell)
    handleerror()
    c
end

"""
    wnreld(a, op, b)

Compare two double precision windows.

!!! note

    Consider using overloaded operators instead, i.e. `a == b`, `a ⊆ b`, and `a ⊊ b`.

### Arguments ###

- `a`: First window
- `op`: Comparison operator
- `b`: Second window

### Output ###

Returns the result of comparison `a (op) b`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnreld_c.html)
"""
function wnreld(a, op, b)
    if op == :(==)
        op = "="
    elseif op == :(!=)
        op = "<>"
    else
        op = string(op)
    end
    res = ccall((:wnreld_c, libcspice), SpiceBoolean,
                (Ref{Cell{SpiceDouble}}, Cstring, Ref{Cell{SpiceDouble}}),
                a.cell, op, b.cell)
    handleerror()
    Bool(res)
end

Base.:(==)(a::SpiceDoubleCell, b::SpiceDoubleCell) = wnreld(a, :(==), b)
Base.:(⊆)(a::SpiceDoubleCell, b::SpiceDoubleCell) = wnreld(a, :<=, b)
Base.:(⊊)(a::SpiceDoubleCell, b::SpiceDoubleCell) = wnreld(a, :<, b)

"""
    wnsumd(window)

Summarize the contents of a double precision window.

### Arguments ###

- `window`: Window to be summarized

### Output ###

Returns a tuple consisting of:

- `meas`: Total measure of intervals in window
- `avg`: Average measure
- `stddev`: Standard deviation
- `shortest`: Location of shortest interval
- `longest`: Location of longest interval

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnsumd_c.html)
"""
function wnsumd(window)
    meas = Ref{SpiceDouble}()
    avg = Ref{SpiceDouble}()
    stddev = Ref{SpiceDouble}()
    shortest = Ref{SpiceInt}()
    longest = Ref{SpiceInt}()
    ccall((:wnsumd_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, Ref{SpiceDouble}, Ref{SpiceDouble},
           Ref{SpiceDouble}, Ref{SpiceInt}, Ref{SpiceInt}),
          window.cell, meas, avg, stddev, shortest, longest)
    handleerror()
    meas[], avg[], stddev[], shortest[] + 1, longest[] + 1
end

"""
    wnunid(a, b)

Place the union of two double precision windows into a third window.

### Arguments ###

- `a`: Input window
- `b`: Input window

### Output ###

Returns a window containing the union of `a` and `b`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnunid_c.html)
"""
function wnunid(a, b)
    c = SpiceDoubleCell(2 * (wncard(a) + wncard(b)))
    ccall((:wnunid_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
          a.cell, b.cell, c.cell)
    handleerror()
    c
end

"""
    wnvald!(window)

Form a valid double precision window from the contents of a window array.

### Arguments ###

- `window`: A (possibly uninitialized) `SpiceDoubleCell` containing endpoints of
    (possibly unordered and non-disjoint) intervals.

### Output ###

Returns the validated window.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/wnvald_c.html)
"""
function wnvald!(window)
    ccall((:wnvald_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{Cell{SpiceDouble}}),
          window.cell.size, card(window), window.cell)
    handleerror()
    window
end

@deprecate wnvald wnvald!


text:
module MultiFidelityABC
include("simmethods.jl")
include("mfabc.jl")
include("produce_figs.jl")
end
include("Viral.jl")
include("Repressilator.jl")
using .MultiFidelityABC
using LaTeXStrings

### PRODUCE HUGE SIMULATIONS (only needed once)
#=
Best to begin Julia with parallel workers set up using:
% julia -p 10
for 10 parallel workers, for example.
=#
# MakeBenchmarkCloud(Repressilator.mf_prob, 5000000, "./repressilator/output/")
# MakeBenchmarkCloud(Viral.mf_prob, 100000, "./repressilator/output/")

#=
## SET UP

# Import the output from the repressilator simulations
bm = MakeBenchmarkCloud("./repressilator/output/")
sample_size = 10^4

epsilons = (50.0, 50.0)
F1(k) = (1.9<k[2]<2.1)
F2(k) = (1.5<k[2]<1.6)
F3(k) = k[2]
F = [F1, F2, F3]
budgets = [100.0,200.0,300.0]
=#

#=
## PRODUCE FIGURES

# Fig 1a
fig1a = view_distances(sim_set, epsilons)
# Fig 1b
fig1b = view_distances(sim_set, epsilons, 2, L"n")

# Fig 2a
fig2a = compare_efficiencies(bm, sim_sets, epsilons, output="theory")
# Fig 2b, Table 1
fig2b, t1, t1_latex = compare_efficiencies(bm, sim_sets, epsilons, output="plot")

# Table 2, 3 (and table 2 insert)
t2,t3,t2etas = observed_variances(bm, sim_sets, epsilons, F, budgets)

## SET UP
# Import the output from the viral dynamics simulations
allsims = BenchmarkCloud("./output/")
splitsims = Iterators.partition(allsims, 10^4)
bm = collect(first(splitsims))
simset = Iterators.drop(splitsims,1)

epsilons = (??)

## PRODUCE FIGURES
# Fig 3a
fig3a = view_distances(bm, epsilons)
# Fig 3b
fig3b = view_distances(bm, epsilons, 1, L"k_1")

# Fig 4
fig4a = efficiency_histogram(bm, simset, epsilons, method = "mf")
fig4b = efficiency_histogram(bm, simset, epsilons, method = "ed")
fig4c = efficiency_histogram(bm, simset, epsilons, method = "er")

=#

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text:
include("ImpedanceAnalyzer_Functions.jl")

function PlotStats_AgainstHeightOffPM(ZMagVec::Vector,ZPhaseVec::Vector,Height::AbstractFloat,newPlot,CVal)

    HeightVec = repeat([Height], Int(length(ZMagVec)))

    if newPlot
        Tmp = display(dotplot([HeightVec HeightVec],[ZMagVec ZPhaseVec], layout=(2,1),color=CVal,legend=[false false],xlabel = [" " "Pancake dist. from steel [mm]"],ylabel=["|Z|" "∠Z [deg]"]))
        # display(StatsPlots.scatter(ZMagVec,KnownMagVec))
    else
        display(dotplot!([HeightVec HeightVec],[ZMagVec ZPhaseVec], layout=(2,1),color=CVal,legend=[false false]))
        # display(StatsPlots.scatter!(ZMagVec,KnownMagVec))
    end

end

function PlotManyElements(PathArray,HeightList,ColorVec)

    for i in 1:length(KnownImpedanceVec)

        StatsData = UnpackStatsRead(PathArray[i])
        ZPhaseVec = StatsData.ZPhaseArray
        ZMagVec = StatsData.ZMagArray

        if i==1
            PlotStats_AgainstHeightOffPM(ZMagVec,ZPhaseVec,HeightList[i],true,ColorVec[i])
        else
            PlotStats_AgainstHeightOffPM(ZMagVec,ZPhaseVec,HeightList[i],false,ColorVec[i])
        end

    end

end

function PlotNTimetraces(PathArray,LabelList,ColorVec)
    for i in 1:length(PathArray)

        StatsData = UnpackStatsRead(PathArray[i])
        ZPhaseVec = StatsData.ZPhaseArray
        ZMagVec = StatsData.ZMagArray
        if i==1
            display(plot([ZMagVec ZPhaseVec], layout=(2,1),color=ColorVec[i],label = string(LabelList[i]) ,legend=[true false],xlabel = [" " "Acquisition number"],ylabel=["|Z|" "∠Z [deg]"]))
        else
            display(plot!([ZMagVec ZPhaseVec], layout=(2,1),color=ColorVec[i],label = string(LabelList[i]) ,legend=[true false],xlabel = [" " "Acquisition number"],ylabel=["|Z|" "∠Z [deg]"]))
        end
    end
end

function PlotNTimetraces_Norm(PathArray,LabelList,ColorVec)
    for i in 1:length(PathArray)

        StatsData = UnpackStatsRead(PathArray[i])
        ZPhaseVec = StatsData.ZPhaseArray
        ZMagVec = StatsData.ZMagArray
        if i==1
            display(plot([ZMagVec./ZMagVec[1] ZPhaseVec./ZPhaseVec[1]], layout=(2,1),color=ColorVec[i],label = string(LabelList[i]) ,legend=[true false],xlabel = [" " "Acquisition number"],ylabel=["|Z|" "∠Z [deg]"]))
        else
            display(plot!([ZMagVec./ZMagVec[1] ZPhaseVec./ZPhaseVec[1]], layout=(2,1),color=ColorVec[i],label = string(LabelList[i]) ,legend=[true false],xlabel = [" " "Acquisition number"],ylabel=["|Z|" "∠Z [deg]"]))
        end
    end
end

function PlotN_Mean_Stdeviation(PathArray,HeightList,msize=8)
    MeanLVec = zeros(length(PathArray))
    StdLVec = zeros(length(PathArray))
    for i in 1:length(PathArray)

        StatsData = UnpackStatsRead(PathArray[i])

        MeanLVec[i] = StatsData.ImpedanceData["Mean_Ls"].*1e6
        StdLVec[i] = StatsData.ImpedanceData["Std_Ls"].*1e6
    end
    display(Plots.scatter([HeightList HeightList],[MeanLVec log10.(StdLVec)], layout=(2,1),legend=[false false],xlabel = [" " "Pancake dist. from steel [mm]"],ylabel=["Mean(L) [μH]" "log₁₀(Stdev(L)) [μH]"],ms=msize))

end

PathArray = [
"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_13_15_32_47_PM_Meas_Pancake_165mm.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_13_15_27_33_PM_Meas_Pancake_133mm.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_13_15_15_56_PM_Meas_Pancake_38mm.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_13_15_21_55_PM_Meas_Pancake_88mm.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_13_15_08_52_PM_Meas_Pancake_0mm.tdms"]

HeightList = [165. , 133.0 , 38.0 , 88.0, 0.0]

ColorVec = ["blue" , "red" , "green" , "purple" ,"yellow" ]

PathArray = [
"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_09_49_47_PM_Meas_Saddle-0mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_09_53_21_PM_Meas_Saddle-50mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_09_56_42_PM_Meas_Saddle-100mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_00_32_PM_Meas_Saddle-150mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_03_56_PM_Meas_Saddle-200mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_09_00_PM_Meas_Saddle-250mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_13_03_PM_Meas_Saddle-300mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_19_50_PM_Meas_Saddle-350mm-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_10_28_09_PM_Meas_Saddle-400mm-10Hz.tdms"

]

HeightList = 0:50:400

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

## 10 Hz ##

PathArray = [
"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_16_09_19_Demagnetized-MagnetizedonTop-10Hz-long.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_14_46_17_DemagnetizedNdFeB-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_14_08_06_CoilOnly-10Hz.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_13_53_49_NdFeB-1Inch-PreHeating-10Hz.tdms"]

LabelList = [PathArray[i][68:end-5] for i in 1:length(PathArray)]

## 100 Hz ##
PathArray = [
"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_15_16_41_DemagnetizedNdFeB-100Hz-long.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_15_41_39_CoilOnly-100Hz-long.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_15_55_51_Demagnetized-MagnetizedonTop-100Hz-long.tdms"

"C:\\Users\\MPI\\Documents\\Pancake Impedance on PM\\2022_01_14_14_00_20_NdFeB-1Inch-PreHeating-100Hz.tdms"]

LabelList = [PathArray[i][68:end-5] for i in 1:length(PathArray)]


text:
function definelogicalset(logicset::LogicalCombo, command::String)::LogicalCombo
    counter!()
    m = match(r"^\s*(.+?)(?:\b|\s)(?:in|∈)(?:\b|\s)([a-zA-Z0-9,._ :\"']+)(?:\|\|(.*?)){0,1}$", command)
    left, right, condition = strip.((x -> x === nothing ? "" : x).(m.captures))
    vars   = split(left, ",")  .|> strip
    values = split(right, ",") .|> strip

    occursin(r"\[*unique(\s|\])*$", command) && return defineuniquelogicalset(logicset, command)

    (logicset.type == "UniquePermutation") && throw("Unique sets cannot have variables added to them!")

    for v in vars
        m = match.(r"(\s|\"|\')", v)
        (m != nothing) && throw("Variable names cannot have {$(m.captures[1])}.")
    end

    for v in values
        m = match.(r"(\"|\')", v)
        (m != nothing) && throw("Variable values cannot be defined with {$(m.captures[1])}.")
    end

    varnames = string.(vars)

    if length(values) == 1 && occursin(r"^[0-9]+:[0-9]+$", values[1])
        values = collect(range(values[1]))
    elseif all([occursin(r"^[0-9]+$", i) for i in values])
        values =  [integer(i) for i in values]
    else
        values = string.(values)
    end

    if condition == ""
        logicset  = expand(logicset, varnames, values)
    else # Loops through conditionals and sets restrictions on data as it is generated
        for i in 1:length(varnames)
          logicset  = expand(logicset, [varnames[i]], values)
          try
            templogicset = operatorspawn(condition, logicset, verbose=false)
            logicset = templogicset
          catch
          end
        end
    end

  logicset
end

definelogicalset(command::String) = definelogicalset(LogicalCombo(), command)

function defineuniquelogicalset(logicset::LogicalCombo, command::String)::LogicalCombo
  (nfeasible(logicset) > 0) && throw("Unique Permutation Sets must be generated from scratch")

  command = replace(command, r"\|*\s*unique\s*$" => "")

  m = match(r"^\s*(.+?)(?:\b|\s)(?:in|∈)(?:\b|\s|)([a-zA-Z0-9,._ :\"']*)(?:\|\|(.*?)){0,1}$", command)
  left, right, condition = strip.((x -> x === nothing ? "" : x).(m.captures))
  vars   = split(left , ",")  .|> strip
  values = split(right, ",")  .|> strip

  if length(values) == 1 && occursin(r"^[0-9]+:[0-9]+$", values[1])
      values = collect(range(values[1]))
  elseif all([occursin(r"^[0-9]+$", i) for i in values])
      values =  [integer(i) for i in values]
  else
      values = string.(values)
  end

  mykeys    = [Symbol(v) for v in vars]

  (right == "") && (mydomain = 1:length(mykeys))
  (length(values) > 1)  && (mydomain = values)

  (length(mydomain) != length(mykeys)) && throw("Variable Lengths Need to Equal Number of Variables")

  mylogical = fill(true, factorial(length(mydomain)))
  LogicalCombo(mykeys, mydomain, mylogical, permutationuniquelookup, "UniquePermutation")
end


text:
using POMDPModels
using Test

let
    problem = BabyPOMDP()
    solver = RandomSolver(rng=Random.MersenneTwister(1))
    policy = solve(solver, problem)
    sim = RolloutSimulator(max_steps=10, rng=Random.MersenneTwister(1))
    r1 = @inferred simulate(sim, problem, policy, updater(policy), initialstate_distribution(problem))

    sim = RolloutSimulator(max_steps=10, rng=Random.MersenneTwister(1))
    dummy = @inferred simulate(sim, problem, policy, updater(policy), nothing, true)

    problem = LegacyGridWorld()
    solver = RandomSolver(rng=Random.MersenneTwister(1))
    policy = solve(solver, problem)
    sim = RolloutSimulator(max_steps=10, rng=Random.MersenneTwister(1))
    r2 = @inferred simulate(sim, problem, policy, initialstate(problem, sim.rng))

    problem = LegacyGridWorld()
    solver = RandomSolver(rng=Random.MersenneTwister(1))
    policy = solve(solver, problem)
    sim = RolloutSimulator(Random.MersenneTwister(1), 10) # new constructor
    r2 = @inferred simulate(sim, problem, policy, initialstate(problem, sim.rng))

    @test isapprox(r1, -27.27829, atol=1e-3)
    @test isapprox(r2, 0.0, atol=1e-3)
end


text:
import ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid

abstract type AbstractDiscretizedDomain end

"""
    DiscretizedDomain
"""
struct DiscretizedDomain{𝒜, ℬ, 𝒞} <: AbstractDiscretizedDomain
    domain::𝒜
    resolution::ℬ
    numerical::𝒞
end

function DiscretizedDomain(
    domain::ProductDomain;
    elements = nothing,
    polynomial_order = nothing,
    overintegration_order = nothing,
    FT = Float64,
    mpicomm = MPI.COMM_WORLD,
    array = ClimateMachine.array_type(),
    topology = StackedBrickTopology,
    brick_builder = uniform_brick_builder,
)

    grid = DiscontinuousSpectralElementGrid(
        domain,
        elements = elements,
        polynomialorder = polynomial_order .+ overintegration_order,
        FT = FT,
        mpicomm = mpicomm,
        array = array,
        topology = topology,
        brick_builder = brick_builder,
    )
    return DiscretizedDomain(domain, (; elements, polynomial_order, overintegration_order), grid)
end

function DiscretizedDomain(
    domain::SphericalShell;
    elements = nothing,
    polynomial_order = nothing,
    overintegration_order = nothing,
    FT = Float64,
    mpicomm = MPI.COMM_WORLD,
    array = ClimateMachine.array_type(),
)
    new_polynomial_order = convention(polynomial_order, Val(2))
    new_polynomial_order = new_polynomial_order .+ convention(overintegration_order, Val(2))
    horizontal, vertical = new_polynomial_order

    grid = DiscontinuousSpectralElementGrid(
        domain,
        elements = elements,
        polynomialorder = (; vertical, horizontal),
        FT = FT,
        mpicomm = mpicomm,
        array = array,
    )
    return DiscretizedDomain(domain, (; elements, polynomial_order, overintegration_order), grid)
end

function DiscretizedDomain(
    domain::SingleStack;
    elements = nothing,
    polynomial_order = nothing,
    overintegration_order = nothing,
    FT = Float64,
    mpicomm = MPI.COMM_WORLD,
    array = ClimateMachine.array_type(),
)
    new_polynomial_order = convention(polynomial_order, Val(2))
    new_polynomial_order = new_polynomial_order .+ convention(overintegration_order, Val(2))
    horizontal, vertical = new_polynomial_order

    grid = DiscontinuousSpectralElementGrid(
        domain,
        elements = elements,
        polynomialorder = (; vertical, horizontal),
        FT = FT,
        mpicomm = mpicomm,
        array = array,
    )
    return DiscretizedDomain(domain, (; elements, polynomial_order, overintegration_order), grid)
end

"""
function DiscontinuousSpectralElementGrid(domain::ProductDomain; elements = nothing, polynomialorder = nothing)
# Description
Computes a DiscontinuousSpectralElementGrid as specified by a product domain
# Arguments
-`domain`: A product domain object
# Keyword Arguments
-`elements`: A tuple of integers ordered by (Nx, Ny, Nz) for number of elements
-`polynomialorder`: A tupe of integers ordered by (npx, npy, npz) for polynomial order
-`FT`: floattype, assumed Float64 unless otherwise specified
-`topology`: default = StackedBrickTopology
-`mpicomm`: default = MPI.COMM_WORLD
-`array`: default = ClimateMachine.array_type()
-`brickbuilder`: default = uniform_brick_builder,
  brickrange=uniform_brick_builder(domain, elements)
# Return
A DiscontinuousSpectralElementGrid object
"""
function DiscontinuousSpectralElementGrid(
    domain::ProductDomain;
    elements,
    polynomialorder,
    FT = Float64,
    mpicomm = MPI.COMM_WORLD,
    array = ClimateMachine.array_type(),
    topology = StackedBrickTopology,
    brick_builder = uniform_brick_builder,
)

    if elements === nothing
        error_message = "Please specify the number of elements as a tuple whose size is commensurate with the domain,"
        error_message *= " e.g., a 3 dimensional domain would need a specification like elements = (10,10,10)."
        error_message *= " or elements = (vertical = 8, horizontal = 5)"

        @error(error_message)
        return nothing
    end

    if polynomialorder === nothing
        error_message = "Please specify the polynomial order as a tuple whose size is commensurate with the domain,"
        error_message *= "e.g., a 3 dimensional domain would need a specification like polynomialorder = (3,3,3)."
        error_message *= " or polynomialorder = (vertical = 8, horizontal = 5)"

        @error(error_message)
        return nothing
    end

    dimension = ndims(domain)

    if (dimension < 2) || (dimension > 3)
        error_message = "SpectralElementGrid only works with dimensions 2 or 3. "
        error_message *= "The current dimension is " * string(ndims(domain))

        println("The domain is ", domain)
        @error(error_message)
        return nothing
    end

    elements = convention(elements, Val(dimension))
    if ndims(domain) != length(elements)
        @error("Incorrectly specified elements for the dimension of the domain")
        return nothing
    end

    polynomialorder = convention(polynomialorder, Val(dimension))
    if ndims(domain) != length(polynomialorder)
        @error("Incorrectly specified polynomialorders for the dimension of the domain")
        return nothing
    end

    brickrange = brick_builder(domain, elements, FT = FT)

    if dimension == 2
        boundary = ((1, 2), (3, 4))
    else
        boundary = ((1, 2), (3, 4), (5, 6))
    end

    periodicity = domain.periodicity
    connectivity = dimension == 2 ? :face : :full

    topl = topology(mpicomm, brickrange; periodicity = periodicity, boundary = boundary, connectivity = connectivity)

    grid =
        DiscontinuousSpectralElementGrid(topl, FloatType = FT, DeviceArray = array, polynomialorder = polynomialorder)

    return grid
end

function DiscontinuousSpectralElementGrid(
    domain::SphericalShell;
    elements,
    polynomialorder,
    mpicomm = MPI.COMM_WORLD,
    boundary = (1, 2),
    FT = Float64,
    array = Array,
)
    Rrange = grid1d(domain.radius, domain.radius + domain.height, nelem = elements.vertical)

    topl = StackedCubedSphereTopology(mpicomm, elements.horizontal, Rrange; boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = array,
        polynomialorder = (polynomialorder.horizontal, polynomialorder.vertical),
        meshwarp = equiangular_cubed_sphere_warp,
    )
    return grid
end

"""
function DiscontinuousSpectralElementGrid(
    domain::SingleStack;
"""
function DiscontinuousSpectralElementGrid(
    domain::SingleStack;
    elements,
    polynomialorder,
    mpicomm = MPI.COMM_WORLD,
    boundary = ((0, 0), (0, 0), (5, 6)),
    periodicity = (true, true, false),
    meshwarp = (x...) -> identity(x),
    FT = Float64,
    array = Array,
)
    xmin, xmax = zero(FT), domain.xmax
    ymin, ymax = zero(FT), domain.ymax

    zmin, zmax = (domain.zmin, domain.zmax)

    brickrange = (
        grid1d(xmin, xmax, nelem = elements.horizontal),
        grid1d(ymin, ymax, nelem = elements.horizontal),
        grid1d(zmin, zmax, nelem = elements.vertical),
    )

    topl = StackedBrickTopology(mpicomm, brickrange, periodicity = periodicity, boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = array,
        polynomialorder = (polynomialorder.horizontal, polynomialorder.vertical),
        meshwarp = meshwarp,
    )
    return grid
end

"""
    Brick builder
"""
function uniform_brick_builder(domain::ProductDomain, elements; FT = Float64)
    dimension = ndims(domain)

    tuple_ranges = []
    for i in 1:dimension
        push!(tuple_ranges, range(FT(domain[i].min); length = elements[i] + 1, stop = FT(domain[i].max)))
    end

    brickrange = Tuple(tuple_ranges)
    return brickrange
end

"""
    Conventions for polynomial order and overintegration order
"""
function convention(a::NamedTuple{(:vertical, :horizontal), T}, ::Val{3}) where {T}
    return (a.horizontal, a.horizontal, a.vertical)
end

function convention(a::Number, ::Val{3})
    return (a, a, a)
end

function convention(a::NamedTuple{(:vertical, :horizontal), T}, ::Val{2}) where {T}
    return (a.horizontal, a.vertical)
end

function convention(a::Number, ::Val{2})
    return (a, a)
end

function convention(a::Tuple, b)
    return a
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
#===============================================================================

    thermal_performance_diagnostic.jl

===============================================================================#

import DataStructures.OrderedDict
import Memento.debug

#===============================================================================
    preset values
===============================================================================#

const deltaTSetPoint    = 1
const highTemperature   = 95
const lowTemperature    = 50
const HIGH_TEMP_mDELTA  = highTemperature - deltaTSetPoint
const LOW_TEMP_pDELTA   = lowTemperature + deltaTSetPoint
const MIN_AVG_RAMP_RATE = 2         ## °C/s
const MAX_TOTAL_TIME    = 22.5e3    ## ms
const MAX_BLOCK_DELTA   = 2         ## °C
const MIN_HEATING_RATE  = 1         ## °C/s
const MAX_TIME_TO_HEAT  = 90e3      ## ms

#===============================================================================
    function definition >>
===============================================================================#

## called by dispatch()
function act(
    ::Type{Val{thermal_performance_diagnostic}},
    req             ::Associative;
    out_format      ::OutputFormat = pre_json_output
    ## remove MySql dependency
    #
    # db_conn ::MySQL.MySQLHandle,
    # exp_id ::Integer, # really used
    # calib_info ::Union{Integer,OrderedDict} ## not used for computation
)
    debug(logger, "at act(::Type{Val{thermal_performance_diagnostic}})")

    ## remove MySql dependency
    #
    # queryTemperatureData = "SELECT * FROM temperature_logs WHERE experiment_id = $exp_id ORDER BY elapsed_time"
    # temperatureData = MySQL.mysql_execute(db_conn, queryTemperatureData)[1]
    # num_dp = length(temperatureData[1]) ## dp = data points

    elapsed_times = req[ELAPSED_TIME_KEY]
    num_dp = length(elapsed_times)
    #
    ## add a new column (not row) that is the average of the two heat block zones
    hbzt_avg =
        map(range(1, num_dp)) do i
            mean(
                map([HEAT_BLOCK_ZONE_1_TEMP_KEY, HEAT_BLOCK_ZONE_2_TEMP_KEY]) do zone
                    req[zone][i]
                end) ## do zone
        end ## do i
    #
    ## calculate average ramp rates up and down of the heat block
    ## first, calculate the time the heat block reaches the high temperature
    ## this is also the time the ramp up ends and the ramp down starts
    (apprxRampUpEndTime, apprxRampDownStartTime) =
        extrema(elapsed_times[hbzt_avg .> HIGH_TEMP_mDELTA])
    ## second, calculate the time the ramp up starts and the ramp down ends
    hbzt_lower = hbzt_avg .< LOW_TEMP_pDELTA
    # const elapsed_times_low_temp  = elapsed_times[hbzt_lower]
    # const apprxRampDownEndTime, apprxRampUpStartTime = extrema(elapsed_times_low_temp)
    apprxRampUpStartTime =
        try maximum(elapsed_times[hbzt_lower .& (elapsed_times .< apprxRampUpEndTime)])
        catch
            -Inf_T
        end ## try maximum
    apprxRampDownEndTime =
        try minimum(elapsed_times[hbzt_lower .& (elapsed_times .> apprxRampDownStartTime)])
        catch
            Inf_T
        end ## try minimum
    #
    temp_range_adj = (HIGH_TEMP_mDELTA - LOW_TEMP_pDELTA) * 1000
    #
    ## calculate the average ramp rate up and down in °C per second
    Heating_TotalTime   = apprxRampUpEndTime   - apprxRampUpStartTime
    Heating_AvgRampRate = round(temp_range_adj / Heating_TotalTime, JSON_DIGITS)
    Cooling_TotalTime   = apprxRampDownEndTime - apprxRampDownStartTime
    Cooling_AvgRampRate = round(temp_range_adj / Cooling_TotalTime, JSON_DIGITS)
    ## calculate maximum temperature difference between heat block zones during ramp up and down
    Heating_MaxBlockDeltaT, Cooling_MaxBlockDeltaT =
        map((
            [apprxRampUpStartTime,   apprxRampUpEndTime],
            [apprxRampDownStartTime, apprxRampDownEndTime]
        )) do time_vec
            elapsed_time_idc = find(elapsed_times) do elapsed_time
                time_vec[1] < elapsed_time < time_vec[2]
            end ## do elapsed_time
            round(
                maximum(
                    abs.(
                        req[HEAT_BLOCK_ZONE_1_TEMP_KEY][elapsed_time_idc] .-
                            req[HEAT_BLOCK_ZONE_2_TEMP_KEY][elapsed_time_idc])),
                JSON_DIGITS)
        end ## do time_vec
    ## calculate the average ramp rate of the lid heater in degrees °C per second
    lidHeaterStartRampTime =
        minimum(elapsed_times[
            req[LID_TEMP_KEY] .> LOW_TEMP_pDELTA])
    lidHeaterStopRampTime =
        maximum(elapsed_times[
            req[LID_TEMP_KEY] .< HIGH_TEMP_mDELTA])
    Lid_TotalTime = lidHeaterStopRampTime - lidHeaterStartRampTime
    Lid_HeatingRate = round(temp_range_adj / Lid_TotalTime, JSON_DIGITS)
    #
    output =
        OrderedDict(
            :Heating => OrderedDict(
                :AvgRampRate    => (Heating_AvgRampRate,    Heating_AvgRampRate    .>= MIN_AVG_RAMP_RATE),
                :TotalTime      => (Heating_TotalTime,      Heating_TotalTime      .<= MAX_TOTAL_TIME),
                :MaxBlockDeltaT => (Heating_MaxBlockDeltaT, Heating_MaxBlockDeltaT .<= MAX_BLOCK_DELTA)
            ),
            :Cooling => OrderedDict(
                :AvgRampRate    => (Cooling_AvgRampRate,    Cooling_AvgRampRate    .>= MIN_AVG_RAMP_RATE),
                :TotalTime      => (Cooling_TotalTime,      Cooling_TotalTime      .<= MAX_TOTAL_TIME),
                :MaxBlockDeltaT => (Cooling_MaxBlockDeltaT, Cooling_MaxBlockDeltaT .<= MAX_BLOCK_DELTA)
            ),
            :Lid => OrderedDict(
                :HeatingRate    => (Lid_HeatingRate, Lid_HeatingRate .>= MIN_HEATING_RATE),
                :TotalTime      => (Lid_TotalTime,   Lid_TotalTime   .<= MAX_TIME_TO_HEAT)
            ),
            :valid => true)
    return output |> out(out_format)
end ## act(::Type{Val{thermal_performance_diagnostic}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using DocumenterGitHubActionsExample: foo
using Test

@testset "DocumenterGitHubActionsExample.jl" begin
    # Write your tests here.
    @test foo(1) == 2
    @test foo(3) == 10
end

max_stars_repo_name: mortenpi/DocumenterGitHubActionsExample

max_issues_repo_name: mortenpi/DocumenterGitHubActionsExample

max_forks_repo_name: mortenpi/DocumenterGitHubActionsExample


text:
module Scheduler
# Write your package code here.
export Qualification, Aircrew, NonAircrew, Event, FlightEvent, GroundEvent, Position

include("Aircrew.jl")
include("Event.jl")

end


text:
module Quadratures

import GaussQuadrature
import StaticArrays: SVector, SMatrix, MMatrix
import LinearAlgebra: Diagonal

export QuadratureStyle,
    GLL, GL, polynomial_degree, degrees_of_freedom, quadrature_points

"""
   QuadratureStyle

Quadrature style supertype. See sub-types:
 - [`GLL`](@ref)
 - [`GL`](@ref)
 - [`Uniform`](@ref)
"""
abstract type QuadratureStyle end

"""
    polynomial_degree(QuadratureStyle) -> Int

Returns the polynomial degree of the `QuadratureStyle` concrete type
"""
function polynomial_degree end

"""
    degrees_of_freedom(QuadratureStyle) -> Int

Returns the degrees_of_freedom of the `QuadratureStyle` concrete type
"""
function degrees_of_freedom end

"""
    points, weights = quadrature_points(::Type{FT}, quadrature_style)

The points and weights of the quadrature rule in floating point type `FT`.
"""
function quadrature_points end

"""
    GLL{Nq}()

Gauss-Legendre-Lobatto quadrature using `Nq` quadrature points.
"""
struct GLL{Nq} <: QuadratureStyle end

@inline polynomial_degree(::GLL{Nq}) where {Nq} = Int(Nq - 1)
@inline degrees_of_freedom(::GLL{Nq}) where {Nq} = Int(Nq)

@generated function quadrature_points(::Type{FT}, ::GLL{Nq}) where {FT, Nq}
    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.both)
    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
end

"""
    GL{Nq}()

Gauss-Legendre quadrature using `Nq` quadrature points.
"""
struct GL{Nq} <: QuadratureStyle end

@inline polynomial_degree(::GL{Nq}) where {Nq} = Int(Nq - 1)
@inline degrees_of_freedom(::GL{Nq}) where {Nq} = Int(Nq)

@generated function quadrature_points(::Type{FT}, ::GL{Nq}) where {FT, Nq}
    points, weights = GaussQuadrature.legendre(FT, Nq, GaussQuadrature.neither)
    :($(SVector{Nq}(points)), $(SVector{Nq}(weights)))
end

"""
    Uniform{Nq}()

Uniformly-spaced quadrature.
"""
struct Uniform{Nq} <: QuadratureStyle end

@inline polynomial_degree(::Uniform{Nq}) where {Nq} = Int(Nq - 1)
@inline degrees_of_freedom(::Uniform{Nq}) where {Nq} = Int(Nq)

@generated function quadrature_points(::Type{FT}, ::Uniform{Nq}) where {FT, Nq}
    points = SVector{Nq}(range(-1 + 1 / Nq, step = 2 / Nq, length = Nq))
    weights = SVector{Nq}(ntuple(i -> 2 / Nq, Nq))
    :($points, $weights)
end

"""
    ClosedUniform{Nq}()

Uniformly-spaced quadrature including boundary.
"""
struct ClosedUniform{Nq} <: QuadratureStyle end

@inline polynomial_degree(::ClosedUniform{Nq}) where {Nq} = Int(Nq - 1)
@inline degrees_of_freedom(::ClosedUniform{Nq}) where {Nq} = Int(Nq)

@generated function quadrature_points(
    ::Type{FT},
    ::ClosedUniform{Nq},
) where {FT, Nq}
    points = SVector{Nq}(range(FT(-1), FT(1), length = Nq))
    weights = SVector{Nq}(
        1 / (Nq - 1),
        ntuple(i -> 2 / (Nq - 1), Nq - 2)...,
        1 / (Nq - 1),
    )
    :($points, $weights)
end

"""
    barycentric_weights(x::SVector{Nq}) where {Nq}

The barycentric weights associated with the array of point locations `x`:

```math
w_j = \\frac{1}{\\prod_{k \\ne j} (x_i - x_j)}
```

See [Berrut2004](@cite), equation 3.2.
"""
function barycentric_weights(r::SVector{Nq, T}) where {Nq, T}
    SVector{Nq}(ntuple(Nq) do i
        w = one(T)
        for j in 1:Nq
            if j != i
                w *= (r[j] - r[i])
            end
        end
        inv(w)
    end)
end
@generated function barycentric_weights(
    ::Type{FT},
    quadstyle::QuadratureStyle,
) where {FT}
    barycentric_weights(quadrature_points(FT, quadstyle())[1])
end

"""
    interpolation_matrix(x::SVector, r::SVector{Nq})

The matrix which interpolates the Lagrange polynomial of degree `Nq-1` through
the points `r`, to points `x`. The matrix coefficients are computed using the
Barycentric formula of [Berrut2004](@cite), section 4:
```math
I_{ij} = \\begin{cases}
1 & \\text{if } x_i = r_j, \\\\
0 & \\text{if } x_i = r_k \\text{ for } k \\ne j, \\\\
\\frac{\\displaystyle \\frac{w_j}{x_i - r_j}}{\\displaystyle \\sum_k \\frac{w_k}{x_i - r_k}} & \\text{otherwise,}
\\end{cases}
```
where ``w_j`` are the barycentric weights, see [`barycentric_weights`](@ref).
"""
function interpolation_matrix(
    points_to::SVector{Nto},
    points_from::SVector{Nfrom},
) where {Nto, Nfrom}
    T = eltype(points_to)
    bw = barycentric_weights(points_from)
    M = zeros(MMatrix{Nto, Nfrom, T, Nto * Nfrom})
    for i in 1:Nto
        x_to = points_to[i]
        skip_row = false
        for j in 1:Nfrom
            if x_to == points_from[j]
                # assign to one to avoid singularity condition
                M[i, j] = one(T)
                # skip over the equal boundry condition
                skip_row = true
            end
            skip_row && break
        end
        skip_row && continue
        w = bw ./ (x_to .- points_from)
        M[i, :] .= w ./ sum(w)
    end
    return SMatrix(M)
end

@generated function interpolation_matrix(
    ::Type{FT},
    quadto::QuadratureStyle,
    quadfrom::QuadratureStyle,
) where {FT}
    interpolation_matrix(
        quadrature_points(FT, quadto())[1],
        quadrature_points(FT, quadfrom())[1],
    )
end

"""
    V = orthonormal_poly(points, quad)

`V_{ij}` contains the `j-1`th Legendre polynomial evaluated at `points[i]`.
i.e. it is the mapping from the modal to the nodal representation.
"""
function orthonormal_poly(
    points::SVector{Np, FT},
    quad::GLL{Nq},
) where {FT, Np, Nq}
    N = Nq - 1
    a, b = GaussQuadrature.legendre_coefs(FT, N)
    if N == 0
        return SMatrix{Np, 1}(ntuple(x -> b[1], Np))
    end
    return SMatrix{Np, Nq}(GaussQuadrature.orthonormal_poly(points, a, b))
end

function spectral_filter_matrix(
    quad::GLL{Nq},
    Σ::SVector{Nq, FT},
) where {Nq, FT}
    points, _ = quadrature_points(FT, quad)
    V = orthonormal_poly(points, quad)
    return V * Diagonal(Σ) / V
end

function cutoff_filter_matrix(
    ::Type{FT},
    quad::GLL{Nq},
    Nc::Integer,
) where {FT, Nq}
    Σ = SVector(ntuple(i -> i <= Nc ? FT(1) : FT(0), Nq))
    return spectral_filter_matrix(quad, Σ)
end

"""
    differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}

The spectral differentiation matrix for the Lagrange polynomial of degree `Nq-1`
interpolating at points `r`.

The matrix coefficients are computed using the [Berrut2004](@cite), section 9.3:
```math
D_{ij} = \\begin{cases}
    \\displaystyle
    \\frac{w_j}{w_i (x_i - x_j)} &\\text{ if } i \\ne j \\\\
    -\\sum_{k \\ne j} D_{kj} &\\text{ if } i = j
\\end{cases}
```
where ``w_j`` are the barycentric weights, see [`barycentric_weights`](@ref).
"""
function differentiation_matrix(r::SVector{Nq, T}) where {Nq, T}
    wb = barycentric_weights(r)
    SMatrix{Nq, Nq, T, Nq * Nq}(
        begin
            if i == j
                D = zero(T)
                for l in 1:Nq
                    if l != i
                        D += one(T) / (r[i] - r[l])
                    end
                end
                D
            else
                (wb[i] / wb[j]) / (r[j] - r[i])
            end
        end for j in 1:Nq, i in 1:Nq
    )
end

"""
    differentiation_matrix(FT, quadstyle::QuadratureStyle)

The spectral differentiation matrix at the quadrature points of `quadstyle`,
using floating point types `FT`.
"""
@generated function differentiation_matrix(
    ::Type{FT},
    quadstyle::QuadratureStyle,
) where {FT}
    differentiation_matrix(quadrature_points(FT, quadstyle())[1])
end

end # module


text:
mutable struct RunEval
    rep_nr::Int64
    density::Float64
    outdegree_sd::Float64
    outdegree_mean::Float64
    indegree_sd::Float64
    indegree_mean::Float64
    supernode_centrality::Float64
    clust_coeff::Float64
    community_count::Float64
    opinionsd::Float64
    opchange_delta_mean::Float64

    function RunEval(simulation::Simulation)

        Random.seed!(0)
        label_prop = label_propagation(simulation.final_state[1])[1]
        n_communities = maximum(label_prop)
        # if n_communities == 1
        #     community_opinion_mean_sds = 0
        # else
        #     community_opinion_mean_sds = std([mean([agent.opinion for agent in simulation.final_state[2] if agent.id in findall(x->x==j, label_prop)]) for j in 1:maximum(label_prop)])
        # end
        new(
            simulation.repnr,
            density(simulation.final_state[1]),
            std(outdegree(simulation.final_state[1])),
            mean(outdegree(simulation.final_state[1])),
            std(indegree(simulation.final_state[1])),
            mean(indegree(simulation.final_state[1])),
            closeness_centrality(simulation.final_state[1])[findmax(outdegree(simulation.final_state[1]))[2]],
            global_clustering_coefficient(simulation.final_state[1]),
            n_communities,
            std([agent.opinion for agent in simulation.final_state[2]]),
            mean([abs(simulation.init_state[2][i].opinion - simulation.final_state[2][i].opinion) for i in 1:simulation.config.network.agent_count])
        )

    end

    function RunEval(
        rep_nr::Int64,
        density::Float64,
        outdegree_sd::Float64,
        outdegree_mean::Float64,
        indegree_sd::Float64,
        indegree_mean::Float64,
        supernode_centrality::Float64,
        clust_coeff::Float64,
        community_count::Float64,
        opinionsd::Float64,
        opchange_delta_mean::Float64
        )

        new(
            rep_nr,
            density,
            outdegree_sd,
            outdegree_mean,
            indegree_sd,
            indegree_mean,
            supernode_centrality,
            clust_coeff,
            community_count,
            opinionsd,
            opchange_delta_mean
        )
    end

end

Base.:+(runeval1::RunEval, runeval2::RunEval) = RunEval(
        0,
        runeval1.density + runeval2.density,
        runeval1.outdegree_sd + runeval2.outdegree_sd,
        runeval1.outdegree_mean + runeval2.outdegree_mean,
        runeval1.indegree_sd + runeval2.indegree_sd,
        runeval1.indegree_mean + runeval2.indegree_mean,
        runeval1.supernode_centrality + runeval2.supernode_centrality,
        runeval1.clust_coeff + runeval2.clust_coeff,
        runeval1.community_count + runeval2.community_count,
        runeval1.opinionsd + runeval2.opinionsd,
        runeval1.opchange_delta_mean + runeval2.opchange_delta_mean
    )

Base.:-(runeval1::RunEval, runeval2::RunEval) = RunEval(
        runeval1.rep_nr + runeval2.rep_nr,
        abs(runeval1.density - runeval2.density),
        abs(runeval1.outdegree_sd - runeval2.outdegree_sd),
        abs(runeval1.outdegree_mean - runeval2.outdegree_mean),
        abs(runeval1.indegree_sd - runeval2.indegree_sd),
        abs(runeval1.indegree_mean - runeval2.indegree_mean),
        abs(runeval1.supernode_centrality - runeval2.supernode_centrality),
        abs(runeval1.clust_coeff - runeval2.clust_coeff),
        abs(runeval1.community_count - runeval2.community_count),
        abs(runeval1.opinionsd - runeval2.opinionsd),
        abs(runeval1.opchange_delta_mean - runeval2.opchange_delta_mean)
    )

Base.:/(runeval::RunEval, i::Int64) = RunEval(
    0,
    runeval.density / i,
    runeval.outdegree_sd / i,
    runeval.outdegree_mean / i,
    runeval.indegree_sd / i,
    runeval.indegree_mean / i,
    runeval.supernode_centrality / i,
    runeval.clust_coeff / i,
    runeval.community_count / i,
    runeval.opinionsd / i,
    runeval.opchange_delta_mean / i
)

Base.:isless(runeval1::RunEval, runeval2::RunEval) = (
    (runeval1.outdegree_sd + runeval1.outdegree_mean + runeval1.indegree_sd
    + runeval1.indegree_mean + runeval1.supernode_centrality + runeval1.clust_coeff
    + runeval1.community_count + runeval1.opinionsd + runeval1.opchange_delta_mean)
    < (runeval2.outdegree_sd + runeval2.outdegree_mean + runeval2.indegree_sd
    + runeval2.indegree_mean + runeval2.supernode_centrality + runeval2.clust_coeff
    + runeval2.community_count + runeval2.opinionsd + runeval2.opchange_delta_mean)
)

function get_prototype(current_run::Vector{Simulation})
    runevals = [RunEval(simulation) for simulation in current_run]
    rep_mean = mean(runevals)
    runeval_dist = [runeval - rep_mean for runeval in runevals]
    best_candidate = minimum(runeval_dist)
    return current_run[best_candidate.rep_nr]
end

# suppress output of include()
;


text:
import Dates

"""
    HolidayRoundingNBD()

Move to the next business day when a date falls on a holiday.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingNBD()), Date(2021,7,10), cal_mgr)
2021-07-12
julia> apply(rd"0d@WEEKEND[NBD]", Date(2021,7,10), cal_mgr)
2021-07-12
```
"""
struct HolidayRoundingNBD <: HolidayRoundingConvention end
function apply(::HolidayRoundingNBD, date::Dates.Date, calendar::Calendar)::Dates.Date
    while is_holiday(calendar, date)
        date += Dates.Day(1)
    end
    date
end
Base.show(io::IO, ::HolidayRoundingNBD) = print(io, "NBD")

"""
    HolidayRoundingPBD()

Move to the previous business day when a date falls on a holiday.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingPBD()), Date(2021,7,10), cal_mgr)
2021-07-09
julia> apply(rd"0d@WEEKEND[PBD]", Date(2021,7,10), cal_mgr)
2021-07-09
```
"""
struct HolidayRoundingPBD <: HolidayRoundingConvention end
function apply(::HolidayRoundingPBD, date::Dates.Date, calendar::Calendar)::Dates.Date
    while is_holiday(calendar, date)
        date -= Dates.Day(1)
    end
    date
end
Base.show(io::IO, ::HolidayRoundingPBD) = print(io, "PBD")

"""
    HolidayRoundingNBDSM()

Move to the next business day when a date falls on a holiday, unless the adjusted
date would be in the next month, then go the previous business date instead.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingNBDSM()), Date(2021,7,10), cal_mgr)
2021-07-12
julia> apply(rd"0d@WEEKEND[NBDSM]", Date(2021,7,10), cal_mgr)
2021-07-12
julia> apply(rd"0d@WEEKEND[NBDSM]", Date(2021,7,31), cal_mgr)
2021-07-30
```
"""
struct HolidayRoundingNBDSM <: HolidayRoundingConvention end
function apply(::HolidayRoundingNBDSM, date::Dates.Date, calendar::Calendar)::Dates.Date
    new_date = date
    while is_holiday(calendar, new_date)
        new_date += Dates.Day(1)
    end

    if Dates.month(new_date) != Dates.month(date)
        new_date = date
        while is_holiday(calendar, new_date)
            new_date -= Dates.Day(1)
        end
    end

    new_date
end
Base.show(io::IO, ::HolidayRoundingNBDSM) = print(io, "NBDSM")

"""
    HolidayRoundingPBDSM()

Move to the previous business day when a date falls on a holiday, unless the adjusted
date would be in the previous month, then go the next business date instead.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingPBDSM()), Date(2021,7,10), cal_mgr)
2021-07-09
julia> apply(rd"0d@WEEKEND[PBDSM]", Date(2021,7,10), cal_mgr)
2021-07-09
julia> apply(rd"0d@WEEKEND[NBDSM]", Date(2021,8,1), cal_mgr)
2021-08-02
```
"""
struct HolidayRoundingPBDSM <: HolidayRoundingConvention end
function apply(::HolidayRoundingPBDSM, date::Dates.Date, calendar::Calendar)::Dates.Date
    new_date = date
    while is_holiday(calendar, new_date)
        new_date -= Dates.Day(1)
    end

    if Dates.month(new_date) != Dates.month(date)
        new_date = date
        while is_holiday(calendar, new_date)
            new_date += Dates.Day(1)
        end
    end

    new_date
end
Base.show(io::IO, ::HolidayRoundingPBDSM) = print(io, "PBDSM")

"""
    HolidayRoundingNR()

No rounding, so just give back the same date even though it's a holiday. Uses
"NR" for short hand.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingNR()), Date(2021,7,10), cal_mgr)
2021-07-10
julia> apply(rd"0d@WEEKEND[NR]", Date(2021,7,10), cal_mgr)
2021-07-10
```
"""
struct HolidayRoundingNR <: HolidayRoundingConvention end
apply(::HolidayRoundingNR, date::Dates.Date, ::Calendar) = date
Base.show(io::IO, ::HolidayRoundingNR) = print(io, "NR")

"""
    HolidayRoundingNearest()

Move to the nearest business day, with the next one taking precedence in a tie.
This is commonly used for U.S. calendars which will adjust Sat to Fri and Sun to Mon
for their fixed date holidays.

#### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(RDates.CalendarAdj("WEEKEND", rd"0d", RDates.HolidayRoundingNearest()), Date(2021,7,10), cal_mgr)
2021-07-09
julia> apply(rd"0d@WEEKEND[NEAR]", Date(2021,7,10), cal_mgr)
2021-07-09
julia> apply(rd"0d@WEEKEND[NEAR]", Date(2021,7,11), cal_mgr)
2021-07-12
```
"""
struct HolidayRoundingNearest <: HolidayRoundingConvention end
function apply(::HolidayRoundingNearest, date::Dates.Date, cal::Calendar)
    is_holiday(cal, date) || return date
    count = 0
    result = nothing
    while result === nothing
        up_date = date + Dates.Day(count)
        if !is_holiday(cal, up_date)
            result = up_date
        else
            down_date = date - Dates.Day(count)
            if !is_holiday(cal, down_date)
                result = down_date
            end
        end
        count += 1
    end
    result
end
Base.show(io::IO, ::HolidayRoundingNearest) = print(io, "NEAR")

const HOLIDAY_ROUNDING_MAPPINGS = Dict(
    "NR" => HolidayRoundingNR(),
    "NBD" => HolidayRoundingNBD(),
    "PBD" => HolidayRoundingPBD(),
    "NBDSM" => HolidayRoundingNBDSM(),
    "PBDSM" => HolidayRoundingPBDSM(),
    "NEAR" => HolidayRoundingNearest()
)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    ConfidenceInterval

An estimated interval around an unknown quantity at a given confidence level
"""
struct ConfidenceInterval{T <: AbstractFloat}
    lower::T
    upper::T
    confidence::T

    function ConfidenceInterval(lower::AbstractFloat, upper::AbstractFloat, confidence::AbstractFloat)
        T = promote_type(typeof(lower), typeof(upper), typeof(confidence), Float32)
        lower <= upper || error("Inconsistent confidence interval: $lower ≰ $upper")
        return new{T}(lower, upper, confidence)
    end
end
lower(interval::ConfidenceInterval) = interval.lower
upper(interval::ConfidenceInterval) = interval.upper
confidence(interval::ConfidenceInterval) = interval.confidence
length(interval::ConfidenceInterval) = upper(interval) - lower(interval)
function in(x::T, interval::ConfidenceInterval) where T <: Real
    isinf(x) && return false
    return interval.lower <= x <= interval.upper
end
function issubset(a::ConfidenceInterval, b::ConfidenceInterval)
    return a.lower >= b.lower && a.upper <= b.upper
end
function show(io::IO, interval::ConfidenceInterval)
    intervalstr = @sprintf "Confidence interval (p = %.0f%%): [%.2f − %.2f]" 100*interval.confidence interval.lower interval.upper
    print(io, intervalstr)
    return io
end


text:
using Knet, CUDArt, ProgressMeter
using Knet: Net, params, axpy!
if !isdefined(:MNIST)
    include(Pkg.dir("Knet/examples/mnist.jl"))
    mnist = minibatch(MNIST.xtrn, MNIST.ytrn, 100)
end

@knet function model02(x; winit=Gaussian(0,.1), hidden=64, f1=:relu, o...)
    h    = wf(x; out=hidden, f=f1, winit=winit)
    return wf(h; out=10, f=:soft, winit=winit)
end

function initmodel{T}(w::CudaArray{T}, g=similar(w); X=784, Y=10, H=div(length(w),X+Y), lr=0.5)
    f = compile(:model02; hidden=H)
    setp(f, lr=lr)
    forw(f, zeros(T, X, 1))
    back(f, zeros(T, Y, 1), softloss)
    f.reg[2].out0 = CudaArray(w.ptr, (H,X), w.dev)
    f.reg[2].dif0 = CudaArray(g.ptr, (H,X), g.dev)
    f.reg[5].out0 = CudaArray(w.ptr+H*X*sizeof(T), (Y,H), w.dev)
    f.reg[5].dif0 = CudaArray(g.ptr+H*X*sizeof(T), (Y,H), g.dev)
    return f
end

function train(f, data, loss; epochs=100)
    for epoch=1:epochs
        for (x,ygold) in mnist
            ypred = forw(f, x)
            back(f, ygold, loss)
            Knet.update!(f)
        end
        println((epoch,
                 test(f, data, softloss),
                 test(f, data, zeroone)))
    end
end

function test(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in mnist
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    return sumloss/numloss
end

interpolate(alpha, w0, w1)=axpy!(1-alpha, w0, scale!(alpha, copy(w1)))

wtest(w)=test(initmodel(w), MNIST.xtrn, softloss)

function wgrad(w)
    g = similar(w)
    f = initmodel(w, g)
    forw(f, MNIST.xtrn)
    back(f, MNIST.ytrn, softloss)
    return g
end

function plotdata(w0, w1; a=-1, b=2, s=0.1, file="plot.dat")
    xvec = Float32[]
    yvec = Float32[]
    @showprogress 1 "Plotting..." for x=a:s:b
        y = wtest(interpolate(x,w0,w1))
        push!(xvec,x)
        push!(yvec,y)
    end
    writedlm(file, [xvec yvec])
end

function doit()
    global w0=scale!(0.1, randn!(CudaArray(Float32, 794*64)))
    global w1=copy(w0)
    train(initmodel(w1), mnist, softloss)
    global s = axpy!(-1, w0, copy(w1))
    global g = wgrad(w0); scale!(-vecnorm(s)/vecnorm(g), g)
    global z = copy(w0); scale!(-vecnorm(s)/vecnorm(z), z)
    global r1 = randn!(copy(w0)); scale!(vecnorm(s)/vecnorm(r1), r1)
    global r2 = randn!(copy(w0)); scale!(vecnorm(s)/vecnorm(r2), r2)
    global r3 = randn!(copy(w0)); scale!(vecnorm(s)/vecnorm(r3), r3)
end

function doplot()
    plotdata(w0, axpy!(1,s,copy(w0)); a=0, b=1, s=0.02, file="w0w1.dat")
    plotdata(w0, axpy!(1,g,copy(w0)); a=0, b=1, s=0.02, file="w0g0.dat")
    plotdata(w0, axpy!(1,z,copy(w0)); a=0, b=1, s=0.02, file="w0z0.dat")
    plotdata(w0, axpy!(1,r1,copy(w0)); a=0, b=1, s=0.02, file="w0r1.dat")
    plotdata(w0, axpy!(1,r2,copy(w0)); a=0, b=1, s=0.02, file="w0r2.dat")
    plotdata(w0, axpy!(1,r3,copy(w0)); a=0, b=1, s=0.02, file="w0r3.dat")
    plotdata(w1, axpy!(1,r1,copy(w1)); a=0, b=1, s=0.02, file="w1r1.dat")
    plotdata(w1, axpy!(1,r2,copy(w1)); a=0, b=1, s=0.02, file="w1r2.dat")
    plotdata(w1, axpy!(1,r3,copy(w1)); a=0, b=1, s=0.02, file="w1r3.dat")
    plotdata(w1, axpy!(-1,s,copy(w1)); a=0, b=1, s=0.02, file="w1w0.dat")
    plotdata(w1, axpy!(-vecnorm(s)/vecnorm(w1),w1,copy(w1)); a=0, b=1, s=0.02, file="w1z0.dat")
end

function plot2(w0, u, d=vecnorm(u); s=d/30, noise=0, file="plot2.dat")
    u = copy(u)
    if noise > 0
        r = randn!(similar(u))
        axpy!(noise*vecnorm(u)/vecnorm(r), r, u)
    end
    u = scale!(1/vecnorm(u), u)
    xvec = Float32[]
    yvec = Float32[]
    w1 = similar(w0)
    @showprogress 1 "Plotting..." for x=0:s:d
        y = wtest(axpy!(x, u, copy!(w1,w0)))
        push!(xvec,x)
        push!(yvec,y)
    end
    writedlm(file, [xvec yvec])
    return u
end

cosine(x,y)=dot(to_host(x),to_host(y))/(vecnorm(x)*vecnorm(y))


text:
/* parity of positive integers by loop */

int main () {
  int y = 17;
  while (y > 0)
    y = y - 2;
  if (y < 0) {
    printInt(0);
    return 0 ;
    }
  else {
    printInt(1);
    return 0 ;
    }
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc