text:
# generic classical discrete orthogonal polynomial, CDOP

abstract type AbstractCDOP{T,X,N} <: AbstractCOP{T,X,N} end
"""
     AbstractCDOP{T,X,N} <: AbstractCOP{T,X,N}

Following [Koepf  and Schmersau](https://arxiv.org/pdf/math/9703217.pdf), a family `y(x)=p_n(x)=k_x⋅x^n +  ...`
for  `n  ∈  {0, 1,…}, k_n ≠ 0` of polynomials is a family of classic *discrete* orthogonal polynomials if it  is  a
solution of a differential equation

`(a⋅x²+b⋅x+c) ⋅ Δ∇y + (d⋅x + e) ⋅ ∇' + λᵢ⋅ y = 0`,

where  `Δy(x) = y(x+1) - y(x)` and `∇y(x) = y(x) - y(x-1)`.

A family is characterized by the 5 coefficients: `a,b,c,d,e`.
Let `σ = (a⋅x²+b⋅x+c)`, `τ = (d⋅x + e).`

As in the classical-continuous-orthogonal-polynomial case
[`AbstractCCOP`](@ref), from these 5 values the cofficients in the
there-point recursion, and other structural equations can be
represented. These allow polynomial multiplication, integration,
differentiation, conversion, etc. to be defined generically.

[Koekoek and Swarttouw](https://arxiv.org/pdf/math/9602214.pdf)
present an encyclopedia of formula characterizing families of
orthogonal polynomials.

For example, on p29 they give  formulas for Hahn polynomials through:

`n(n+α+β+1)y(x) = B(x)y(x+1) -[B(x)+D(x)]y(x) + D(x)y(x-1)`,  with  explicit values  for  `B` and `D`. Reexpressing gives:
`BΔy(x) - D∇y(x) -λ y(x)  = 0`. From the rexpressed Eqn (4) for Koepf & Schemersau we have the identification:
`σ+τ =  B; σ=D`,  so  `τ=B-D`. From this `a,b,c,d,e` can be  gleaned.

The above, is termed the eigevalue equation (e.g. [Goertz and Offner](https://arxiv.org/pdf/1609.07291.pdf)), as it can be reexpressed as

`Δ(D(x)⋅ω(x)⋅∇yᵢ(x) = λᵢ⋅ω(x)⋅yᵢ(x)`

"""
AbstractCDOP

## Implemented  families
##                    a     b     c      d      e
##                    -------------------------------
## Charlier{μ}        0     1     0     -1      μ
## Meixner{γ,μ}       0     1     0    μ-1     γ⋅μ
## Krawchouk{p,𝐍}     0    p-1    0     1      p/𝐍
## Hahn{α,β,𝐍}        1 -(β+𝐍+1)  0   α+β+2  -𝐍⋅(α+1)

# subtypes  to keep track of number of parameters
# passed to  @registerN macros
abstract type AbstractCDOP0{T,X,N} <: AbstractCDOP{T,X,N} end
abstract type AbstractCDOP1{α,T,X,N} <: AbstractCDOP{T,X,N} end
abstract type AbstractCDOP2{α, β,T,X,N} <: AbstractCDOP{T,X,N}  end
abstract type AbstractCDOP3{α, β,γ,T,X,N} <: AbstractCDOP{T,X,N}  end

⟒(P::Type{<:AbstractCDOP1{α}}) where {α} = constructorof(P){α}
⟒(P::Type{<:AbstractCDOP2{α,β}}) where {α, β} = constructorof(P){α, β}
⟒(P::Type{<:AbstractCDOP3{α,β,γ}}) where {α, β, γ} = constructorof(P){α, β, γ}

# for conversion to base case
# \upstigma[tab]
ϛ(P::Type{<:AbstractCDOP}) = FallingFactorial

# compose with FallingFactorial
function Base.convert(::Type{Q}, p::P)  where  {Q <: AbstractCDOP,  P <: AbstractCDOP}
    T = eltype(P)
    _convert_cop(Q, _convert_cop(FallingFactorial{T},  p))
end

Δₓ(p::AbstractCDOP) = p(variable(p)+1)-p(variable(p))
∇ₓ(p::AbstractCDOP) = p(variable(p))-p(variable(p)-1)

function innerproduct(P::Type{<:AbstractCDOP}, f, g)
    a, b = extrema(P)
    (isinf(a) || isinf(b)) &&  throw(ArgumentError("Only defined for discrete polynomials  with a finite domain"))
    ω  = weight_function(P)
    sum(f(i)*g(i)*ω(i) for i in a:b)
end

## Structural Equations
## https://dlmf.nist.gov/18.22
## we have σ⋅Δ∇p + τ⋅∇p + λp =  0
## Using parameterization   of  above, τ = C, σ = A - C
function Ãn(P::Type{<:AbstractCDOP}, a,b,c,d,e, n::Int)
    one(eltype(P))
end

function B̃n(P::Type{<:AbstractCDOP}, a,b,c,d,e, n::Int)

    S = eltype(P)

    num = n*(d + 2b) * (d + a*n - a) + e * (d - 2a)
    den = (2a*n - 2a   + d )  *  (d + 2a*n)

    iszero(den) && return B̃n(P, Val(n))

    val = one(S) * num / den

    val
end

function C̃n(P::Type{<:AbstractCDOP}, a,b,c,d,e, n::Int)

    S = eltype(P)

    num = one(S)
    num *= (n - 1) * (d + a*n - a)
    num *= a*n*d - d*b - a*d + a^2*n^2 - 2a^2*n + 4c*a + a^2 + 2e*a - b^2
    num += -d*b*e + d^2*c + a*e^2

    num *= (a*n +  d - 2a) * n
    den = (d - a + 2a*n) * (d + 2a*n - 3a) * (2a*n - 2a + d)^2

    iszero(den) && return C̃n(P, Val(n))

    val = -one(S) * num  / den

    val
end

# αn, βn,γn
# σ⋅pn' = [αn, βn,γn] ⋅ [p_{n+1},p_n,p_{n-1}]
#αn(P::Type{<:AbstractCDOP}, n::Int) = α̃(P,n) / k1k0(P,n)
function α̃n(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = a * n
    val = one(S) *  num

    return val
end

function β̃n(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = -n *  (d  + a*n -  a) * (2*a*n*d - a*d - d*b + 2e*a - 2a^2*n +  2a^2*n^2)
    den = (2a*n - 2a +  d) * (d + 2a*n)

    iszero(den) &&  return βn(P, Val(n))

    val = (one(S) *  num)  / den

    return val
end

#γn(P::Type{<:AbstractCDOP}, n::Int) =  γ̃n(P,n) * k1k0(P,n-1)
function γ̃n(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = (n - 1)  * (d + a*n  - a)
    num *= (a*n*d - d*b - a*d + a^2*n^2 - 2a^2*n + 4c*a + a^2 + 2e*a - b^2)
    num += -d*b*e + d^2*c + a*e^2
    num *= (d + a*n -  a) *  (a*n + d - 2a) * n
    den = (d - a + 2a*n) * (d + 2a*n - 3a) * (2a*n - 2a + d)^2
    iszero(den) &&  return γn(P, Val(n))

    val = one(S) * num / den

    return val
end

function abcdeᴵ(P::Type{<:AbstractCDOP})
    a,b,c,d,e = abcde(P).a, abcde(P).b, abcde(P).c, abcde(P).d, abcde(P).e
    NamedTuple{(:a,:b,:c,:d,:e)}((a, b, c, d + 2a, d + e + a + b))
end

function ẫn(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    val = one(S)
    val /= n+1
    return val
end

function b̂̃n(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = -2*a*n * (d + a*n - a) - d*b + a*d - d^2 +2e*a
    den = (2a*n -  2a  +  d) * (d + 2a*n)
    iszero(den) && return b̂̃n(P, Val(n))

    val = one(S) * num/ den
    return val

end

function ĉ̃n(P::Type{<:AbstractCDOP}, n::Int)
    a,b,c,d,e = abcde(P)
    S = eltype(P)

    num = (n-1)*(d + a*n -a) * (a*n*d  -d*b  -a*d + a^2*n^2 -2a^2*n + 4c*a + a^2 +  2e*a - b^2)
    num += -d*b*e + d^2*c +  a*e^2
    num  *=  a*n
    den = (d - a + 2a*n)  * (d + 2a*n -3a) * (2a*n - 2a  + d)^2
    iszero(den)  &&  return ĉ̃n(P, Val(n))

    val = (one(S)  *  num) / den
    return val
end

function ⊗(::Type{<:AbstractCDOP}, p::P, q::Q) where {P <: AbstractCDOP, Q <: AbstractCDOP}

    #@assert  ⟒(P) == ⟒(Q)
    #@assert eltype(p) == eltype(q)
    isconstant(p)  && return q * constantterm(p)
    isconstant(q)  && return p * constantterm(q)
    assert_same_variable(p,q)

    convert(⟒(P), convert(FallingFactorial, p) * convert(FallingFactorial, q))

#    R = eltype(p)
#    convert(⟒(P){R}, convert(Polynomial, p) * convert(Polynomial, q))
end


text:
"""
    Generate proto files.

[x] Proto python
[x] Proto Julia

"""

using ProtoBuf
using PyCall

"""
    Gets proto files include folder.
"""
function get_proto_include_dir()
    py"""
    def get_python_path() -> str:
        import os
        import inspect

        return inspect.getfile(os)
    """

    python_os_path = py"get_python_path()"

    python_os_directory, _ = splitdir(python_os_path)

    protobuf_include = joinpath(python_os_directory, "site-packages", "grpc_tools", "_proto")
    return protobuf_include
end

"""
    Generates Julia gRPC files.
"""
function generate_julia_grpc(proto_dir::String, jl_out_dir::String)
    if !isdir(jl_out_dir)
        println("Generating Julia proto file ...")
        mkdir(jl_out_dir)
    end

    proto_include_dir = "$(get_proto_include_dir())"

    args = `-I=$(proto_include_dir) --proto_path=$(proto_dir) --julia_out=$(jl_out_dir) route_guide.proto helloworld.proto`
    ProtoBuf.protoc(args)

    return nothing
end

"""
    Generates Python gRPC files.

    Equivalent of running the command:
    python -m grpc_tools.protoc -I../../protos --python_out=. --grpc_python_out=. ../../protos/route_guide.proto
"""
function generate_python_grpc(proto_dir::String, py_out_dir::String)
    #if !isdir(py_out_dir)

    println("Generating Python proto file from folder: $(pwd())...")
    #mkdir(py_out_dir)

    py_command = pyimport("grpc.tools.command")

    proto_include_dir = "$(get_proto_include_dir())"

    args = `--proto_path=$(proto_include_dir) --proto_path=./test --python_out=$(py_out_dir) --grpc_python_out=$(py_out_dir) proto/route_guide.proto proto/helloworld.proto`
    py_command.protoc.main(args)

    return nothing
end

proto_dir = "$(pwd())/test/proto"
jl_out_dir = "$(pwd())/test/proto/proto_jl_out"
py_out_dir = "$(pwd())/test/python_test/mymodule"

# Julia codegen
generate_julia_grpc(proto_dir, jl_out_dir)

# Python codegen
generate_python_grpc(proto_dir, py_out_dir)

ext: jl


text:
"""
    `creat_cp_rs_struct(RS_REF)`

    Function that returns a RegStruct object RS with fields taken from RS_REF.

    # Examples
    ```julia-repl
    julia> CpelNano.creat_cp_rs_struct(rs_ref)
    ```
"""
function cp_rs_flds(rs_ref::RegStruct)::RegStruct

    # Create structure
    rs = deepcopy(rs_ref)

    # Reset fields
    rs.m = 0
    rs.nme = []
    rs.mml = []
    rs.ϕhat = []
    rs.calls = []
    rs.logZ = NaN
    rs.tm = TransferMat()
    rs.exps = Expectations()
    rs.nls_rgs = AnalysisRegions()

    # Return
    return rs

end
"""
    `pmap_two_samp_null_stats(MOD_REF,CALLS_S1,CALLS_S2,PERM,CONFIG)`

    Function that computes set of null statistics for a given data permutation.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_two_samp_null_stats(mod_ref,calls_s1,calls_s2,perm,config)
    ```
"""
function pmap_two_samp_null_stats(rs_ref::RegStruct, calls_s1::Vector{Vector{MethCallCpgGrp}}, calls_s2::Vector{Vector{MethCallCpgGrp}}, perm_ids::Vector{Int64}, config::CpelNanoConfig)::NTuple{3,Vector{Float64}}

    # Create two aux structs
    aux_rs_s1 = cp_rs_flds(rs_ref)
    aux_rs_s2 = cp_rs_flds(rs_ref)
    aux_rs_s1.m = length(calls_s1)
    aux_rs_s2.m = length(calls_s2)

    # Assign calls
    aux_calls = vcat(calls_s1, calls_s2)
    aux_rs_s1.calls = aux_calls[perm_ids]
    deleteat!(aux_calls, perm_ids)
    aux_rs_s2.calls = aux_calls

    # Init output
    tmml = fill(NaN, rs_ref.nls_rgs.num)
    tnme = fill(NaN, rs_ref.nls_rgs.num)
    tcmd = fill(NaN, rs_ref.nls_rgs.num)

    # Do estimation
    get_ϕhat!(aux_rs_s1, config)
    get_ϕhat!(aux_rs_s2, config)

    # Leave if we couldn't estimate ϕ
    (length(aux_rs_s1.ϕhat) > 0 && length(aux_rs_s2.ϕhat) > 0) || return tmml, tnme, tcmd
    config.verbose && print_log("ϕhat1=$(aux_rs_s1.ϕhat); ϕhat2=$(aux_rs_s2.ϕhat)")

    # Re-scale to per cpg site resolution
    rscle_grp_mod!(aux_rs_s1)
    rscle_grp_mod!(aux_rs_s2)

    # Statistical summaries
    get_stat_sums!(aux_rs_s1, config)
    get_stat_sums!(aux_rs_s2, config)

    # Compute stats
    tmml = aux_rs_s1.mml - aux_rs_s2.mml
    tnme = aux_rs_s1.nme - aux_rs_s2.nme
    tcmd = comp_cmd(aux_rs_s1, aux_rs_s2)
    config.verbose && print_log("tmml=$(tmml)")
    config.verbose && print_log("tnme=$(tnme)")
    config.verbose && print_log("tcmd=$(tcmd)")

    # Return tuple
    return tmml, tnme, tcmd

end
"""
    `pmap_two_samp_null_stats_test_fun(MOD_REF,CALLS_S1,CALLS_S2,PERM,CONFIG)`

    Function that computes set of null statistics for a given data permutation.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_two_samp_null_stats_test_fun(mod_ref,calls_s1,calls_s2,perm,config)
    ```
"""
function pmap_two_samp_null_stats_test_fun(rs_ref::RegStruct, calls_s1::Vector{Vector{MethCallCpgGrp}}, calls_s2::Vector{Vector{MethCallCpgGrp}}, perm_ids::Vector{Int64}, config::CpelNanoConfig)::NTuple{3,Vector{Float64}}

    # Create two aux structs
    aux_rs_s1 = cp_rs_flds(rs_ref)
    aux_rs_s2 = cp_rs_flds(rs_ref)
    aux_rs_s1.m = length(calls_s1)
    aux_rs_s2.m = length(calls_s2)

    # Assign calls
    aux_calls = vcat(calls_s1, calls_s2)
    aux_rs_s1.calls = aux_calls[perm_ids]
    deleteat!(aux_calls, perm_ids)
    aux_rs_s2.calls = aux_calls

    # Init output
    tmml = fill(NaN, rs_ref.nls_rgs.num)
    tnme = fill(NaN, rs_ref.nls_rgs.num)
    tcmd = fill(NaN, rs_ref.nls_rgs.num)

    # Do estimation
    get_ϕhat!(aux_rs_s1, config)
    get_ϕhat!(aux_rs_s2, config)

    # Leave if we couldn't estimate ϕ
    (length(aux_rs_s1.ϕhat) > 0 && length(aux_rs_s2.ϕhat) > 0) || return tmml, tnme, tcmd
    print_log("ϕhat1=$(aux_rs_s1.ϕhat); ϕhat2=$(aux_rs_s2.ϕhat)")

    # Statistical summaries
    get_stat_sums!(aux_rs_s1, config)
    get_stat_sums!(aux_rs_s2, config)

    # Compute stats
    tmml = aux_rs_s1.mml - aux_rs_s2.mml
    tnme = aux_rs_s1.nme - aux_rs_s2.nme
    tcmd = comp_cmd(aux_rs_s1, aux_rs_s2)
    print_log("μ1=$(aux_rs_s1.mml)")
    print_log("μ2=$(aux_rs_s2.mml)")
    print_log("tmml=$(tmml)")
    print_log("h1=$(aux_rs_s1.nme)")
    print_log("h2=$(aux_rs_s2.nme)")
    print_log("tnme=$(tnme)")
    print_log("tcmd=$(tcmd)")
    print_log("********************")

    # Return tuple
    return tmml, tnme, tcmd

end
"""
    `pmap_diff_two_samp_comp(MODS_CHR_S1,MODS_CHR_S2,CONFIG)`

    Function that performs differential methylation analysis between two samples given two
    CpelNano models.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_diff_two_samp_comp(mods_chr_s1,mods_chr_s2,config)
    ```
"""
function pmap_diff_two_samp_comp(mod_s1::RegStruct, mod_s2::RegStruct, nano_s1::String, nano_s2::String, fa_rec::FASTA.Record, config::CpelNanoConfig)::RegStatTestStruct

    print_log("Working on: $(mod_s1.chr)-$(mod_s1.chrst)-$(mod_s1.chrend)")

    # Init output struct
    test_struct = RegStatTestStruct(mod_s1)

    # Observed statistics
    tmml_obs = mod_s1.mml - mod_s2.mml
    tnme_obs = mod_s1.nme - mod_s2.nme
    tcmd_obs = comp_cmd(mod_s1, mod_s2)

    # Init p-values
    tmml_pvals = fill(NaN, length(tmml_obs))
    tnme_pvals = fill(NaN, length(tnme_obs))
    tcmd_pvals = fill(NaN, length(tcmd_obs))

    # Compute p-values if need be
    if config.pval_comp

        # Get group calls
        get_calls!(mod_s1, nano_s1, config)
        get_calls!(mod_s2, nano_s2, config)

        # Number of reads
        n_rds_s1 = length(mod_s1.calls)
        n_rds_s2 = length(mod_s2.calls)

        # Compute number of possible randomizations
        L = binomial(BigInt(n_rds_s1 + n_rds_s2), BigInt(n_rds_s1))

        # Create iteratable object with all combinations
        comb_iter = combinations(1:(n_rds_s1 + n_rds_s2), n_rds_s1)

        # If enough permutations
        if L > 20

            #  Check if exact test
            exact = L < config.LMAX_TWO_SAMP

            # Go back to group scale
            get_grp_info!(mod_s1, fa_rec, config.min_grp_dist)
            get_grp_info!(mod_s2, fa_rec, config.min_grp_dist)

            # Store calls in aux vars
            calls_s1 = deepcopy(mod_s1.calls)
            calls_s2 = deepcopy(mod_s2.calls)

            # Get read combinations to use
            if exact
                # Use all sample assignments
                comb_iter_used = comb_iter
            else
                # Sample assignments to group 1
                comb_iter_used = fill(Vector{Int64}(), config.LMAX_TWO_SAMP)
                @inbounds for i = 1:config.LMAX_TWO_SAMP
                    comb_iter_used[i] = sort!(StatsBase.sample(1:(n_rds_s1 + n_rds_s2), n_rds_s1;replace=false))
                end
            end

            ## Null statistics

            # Compute null statistics
            pmap_out = pmap(perm -> pmap_two_samp_null_stats(mod_s1, calls_s1, calls_s2, perm, config), comb_iter_used)

            # Distribute statistics
            tmml_perms = [x[1] for x in pmap_out]
            tnme_perms = [x[2] for x in pmap_out]
            tcmd_perms = [x[3] for x in pmap_out]

            ## P-value computation

            # Loop over analysis regions
            @inbounds for k in 1:length(tmml_obs)

                # Check if data
                isnan(tmml_obs[k]) && continue

                # Get permutation stats from k-th analysis region
                tmml_perms_k = [perm[k] for perm in tmml_perms]
                tnme_perms_k = [perm[k] for perm in tnme_perms]
                tcmd_perms_k = [perm[k] for perm in tcmd_perms]

                # Clean NaNs in null stats
                tmml_perms_k = tmml_perms_k[.!isnan.(tmml_perms_k)]
                tnme_perms_k = tnme_perms_k[.!isnan.(tnme_perms_k)]
                tcmd_perms_k = tcmd_perms_k[.!isnan.(tcmd_perms_k)]

                # Check enough null stats after filtering NaNs
                length(tmml_perms_k) > 20 || continue

                # Get number of permutation stats equal or above observed
                tmml_pval_k = sum(abs.(tmml_perms_k) .>= abs(tmml_obs[k]))
                tnme_pval_k = sum(abs.(tnme_perms_k) .>= abs(tnme_obs[k]))
                tcmd_pval_k = sum(tcmd_perms_k .>= tcmd_obs[k])

                # Compute p-values
                tmml_pvals[k] = exact ? tmml_pval_k / length(tmml_perms_k) : (1.0 + tmml_pval_k) / (1.0 + length(tmml_perms_k))
                tnme_pvals[k] = exact ? tnme_pval_k / length(tnme_perms_k) : (1.0 + tnme_pval_k) / (1.0 + length(tnme_perms_k))
                tcmd_pvals[k] = exact ? tcmd_pval_k / length(tcmd_perms_k) : (1.0 + tcmd_pval_k) / (1.0 + length(tcmd_perms_k))

            end

        end
    end

    # Fill return object
    @inbounds for k in 1:length(tmml_obs)
        test_struct.tests.tmml_test[k] = (tmml_obs[k], tmml_pvals[k])
        test_struct.tests.tnme_test[k] = (tnme_obs[k], tnme_pvals[k])
        test_struct.tests.tcmd_test[k] = (tcmd_obs[k], tcmd_pvals[k])
    end

    # Return test struct
    return test_struct

end
"""
    `pmap_diff_two_samp_comp_test(MODS_CHR_S1,MODS_CHR_S2,CONFIG)`

    Function that performs differential methylation analysis between two samples given two
    CpelNano models. This function is for testing purposes.

    # Examples
    ```julia-repl
    julia> CpelNano.pmap_diff_two_samp_comp_test(mods_chr_s1,mods_chr_s2,config)
    ```
"""
function pmap_diff_two_samp_comp_test(mod_s1::RegStruct, mod_s2::RegStruct, config::CpelNanoConfig)::RegStatTestStruct

    # Init output struct
    test_struct = RegStatTestStruct(mod_s1)

    # Observed statistics
    tmml_obs = mod_s1.mml - mod_s2.mml
    tnme_obs = mod_s1.nme - mod_s2.nme
    tcmd_obs = comp_cmd(mod_s1, mod_s2)

    # Init p-values
    tmml_pvals = fill(NaN, length(tmml_obs))
    tnme_pvals = fill(NaN, length(tnme_obs))
    tcmd_pvals = fill(NaN, length(tcmd_obs))

    # Compute p-values if need be
    if config.pval_comp

        # Number of reads
        n_rds_s1 = length(mod_s1.calls)
        n_rds_s2 = length(mod_s2.calls)

        # Compute number of possible randomizations
        L = binomial(BigInt(n_rds_s1 + n_rds_s2), BigInt(n_rds_s1))

        # Create iteratable object with all combinations
        comb_iter = combinations(1:(n_rds_s1 + n_rds_s2), n_rds_s1)

        # If enough permutations
        if L > 20

            #  Check if exact test
            exact = L < config.LMAX_TWO_SAMP

            # # Go back to group scale
            # get_grp_info!(mod_s1, fa_rec, config.min_grp_dist)
            # get_grp_info!(mod_s2, fa_rec, config.min_grp_dist)

            # Store calls in aux vars
            calls_s1 = deepcopy(mod_s1.calls)
            calls_s2 = deepcopy(mod_s2.calls)

            # Get read combinations to use
            if exact
                # Use all sample assignments
                comb_iter_used = comb_iter
            else
                # Sample assignments to group 1
                comb_iter_used = fill(Vector{Int64}(), config.LMAX_TWO_SAMP)
                @inbounds for i = 1:config.LMAX_TWO_SAMP
                    comb_iter_used[i] = sort!(StatsBase.sample(1:(n_rds_s1 + n_rds_s2), n_rds_s1;replace=false))
                    print_log("$(comb_iter_used[i])")
                end
            end

            ## Null statistics

            # Compute null statistics
            pmap_out = pmap(perm -> pmap_two_samp_null_stats_test_fun(mod_s1, calls_s1, calls_s2, perm, config), comb_iter_used)

            # Distribute statistics
            tmml_perms = [x[1] for x in pmap_out]
            tnme_perms = [x[2] for x in pmap_out]
            tcmd_perms = [x[3] for x in pmap_out]

            ## P-value computation

            # Loop over analysis regions
            @inbounds for k in 1:length(tmml_obs)

                # Check if data
                isnan(tmml_obs[k]) && continue

                # Get permutation stats from k-th analysis region
                tmml_perms_k = [perm[k] for perm in tmml_perms]
                tnme_perms_k = [perm[k] for perm in tnme_perms]
                tcmd_perms_k = [perm[k] for perm in tcmd_perms]

                # Clean NaNs in null stats
                tmml_perms_k = tmml_perms_k[.!isnan.(tmml_perms_k)]
                tnme_perms_k = tnme_perms_k[.!isnan.(tnme_perms_k)]
                tcmd_perms_k = tcmd_perms_k[.!isnan.(tcmd_perms_k)]

                # Check enough null stats after filtering NaNs
                length(tmml_perms_k) > 20 || continue

                # Get number of permutation stats equal or above observed
                tmml_pval_k = sum(abs.(tmml_perms_k) .>= abs(tmml_obs[k]))
                tnme_pval_k = sum(abs.(tnme_perms_k) .>= abs(tnme_obs[k]))
                tcmd_pval_k = sum(tcmd_perms_k .>= tcmd_obs[k])

                # Compute p-values
                tmml_pvals[k] = exact ? tmml_pval_k / length(tmml_perms_k) : (1.0 + tmml_pval_k) / (1.0 + length(tmml_perms_k))
                tnme_pvals[k] = exact ? tnme_pval_k / length(tnme_perms_k) : (1.0 + tnme_pval_k) / (1.0 + length(tnme_perms_k))
                tcmd_pvals[k] = exact ? tcmd_pval_k / length(tcmd_perms_k) : (1.0 + tcmd_pval_k) / (1.0 + length(tcmd_perms_k))

            end

        end
    end

    # Fill return object
    @inbounds for k in 1:length(tmml_obs)
        test_struct.tests.tmml_test[k] = (tmml_obs[k], tmml_pvals[k])
        test_struct.tests.tnme_test[k] = (tnme_obs[k], tnme_pvals[k])
        test_struct.tests.tcmd_test[k] = (tcmd_obs[k], tcmd_pvals[k])
    end

    # Return test struct
    return test_struct

end
"""
    `diff_two_samp_comp(MODEL_FILE_G1,MODEL_FILE_G2,FASTA,CONFIG)`

    Function that performs differential methylation analysis between two samples given two
    CpelNano model files.

    # Examples
    ```julia-repl
    julia> CpelNano.diff_two_samp_comp(mods_path_s1,mods_path_s2,fasta,config)
    ```
"""
function diff_two_samp_comp(mods_path_s1::String, mods_path_s2::String, nano_s1::String, nano_s2::String, fasta::String, config::CpelNanoConfig)::Nothing

    # Find chromosomes
    chr_names, chr_sizes = get_genome_info(fasta)

    # Create output directory if not existant
    isdir(config.out_dir) || mkdir(config.out_dir)

    # Get output file names
    config.out_diff_files = OutputDiffFiles(config.out_dir, config.out_prefix)

    # Check if output file exists
    check_diff_output_exists(config) && return nothing

    # Loop over chromosomes
    for (i, chr) in enumerate(chr_names)

        print_log("Processing chromosome $(chr)")

        ## Get fasta record

        fa_rec = get_chr_fa_rec(chr, fasta)

        ## Read in models

        mods_chr_s1 = read_mod_file_chr(mods_path_s1, chr, fa_rec, config)
        mods_chr_s2 = read_mod_file_chr(mods_path_s2, chr, fa_rec, config)

        ## Find regions with at least one sample per group

        comm_regs = find_comm_regs([mods_chr_s1], [mods_chr_s2])

        # If no data in chromosome continue
        length(comm_regs) > 0 || continue
        config.verbose && print_log("Comm regions: $(comm_regs)")

        ## Process each region
        pmap_out = map(x -> pmap_diff_two_samp_comp(mods_chr_s1[x], mods_chr_s2[x], nano_s1, nano_s2, fa_rec, config), comm_regs)

        ## Write chromosome output
        write_diff_output(pmap_out, config)

    end

    # Perform multiple hypothesis testing
    if config.pval_comp
        print_log("Multiple hypothesis testing")
        mult_hyp_corr(config)
    end

    # Return nothing
    return nothing

end


text:
struct Sprite
    x::Integer
    y::Integer

    width::Integer
    height::Integer

    offsetX::Integer
    offsetY::Integer
    realWidth::Integer
    realHeight::Integer

    surface::Cairo.CairoSurface
    filename::String
end

function loadSprites(metaFn::String, spritesFn::String)
    surface = open(Cairo.read_from_png, spritesFn)
    fh = open(metaFn)

    res = Dict{String, Sprite}()

    # Mostly useless information
    read(fh, Int32)
    Maple.readString(fh)
    read(fh, Int32)

    count = read(fh, Int16)
    for i in 1:count
        dataFile = Maple.readString(fh)
        sprites = read(fh, Int16)

        for j in 1:sprites
            path = replace(Maple.readString(fh), "\\", "/")

            res[path] = Sprite(
                read(fh, Int16),
                read(fh, Int16),

                read(fh, Int16),
                read(fh, Int16),

                read(fh, Int16),
                read(fh, Int16),
                read(fh, Int16),
                read(fh, Int16),

                surface,
                spritesFn
            )
        end
    end

    close(fh)

    return res
end

max_stars_repo_name: npdmfixup/Fuji

max_issues_repo_name: npdmfixup/Fuji

max_forks_repo_name: npdmfixup/Fuji


text:
include("Functions.jl")
include("Algorithms.jl")
module MMA
  using Algorithms
  using Functions
export k_means, gram_schmidt, nearest_neighbor, avg, std,
       rms, angle, corrcoef, dist
end

max_stars_repo_name: reesepathak/mma-functions

max_issues_repo_name: reesepathak/mma-functions

max_forks_repo_name: reesepathak/mma-functions


text:
"""
        mabac(decisionMat, weights, fns)

Apply MABAC (Multi-Attributive Border Approximation area Comparison) for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n alternatives and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
mabac() applies the MABAC method to rank n alternatives subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::MABACResult`: MABACResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [2 1 4 7 6 6 7 3000;
       4 1 5 6 7 7 6 3500;
       3 2 6 6 5 6 8 4000;
       5 1 5 7 6 7 7 3000;
       4 2 5 6 7 7 6 3000;
       3 2 6 6 6 6 6 3500]
6×8
Array{Int64,2}:
 2  1  4  7  6  6  7  3000
 4  1  5  6  7  7  6  3500
 3  2  6  6  5  6  8  4000
 5  1  5  7  6  7  7  3000
 4  2  5  6  7  7  6  3000
 3  2  6  6  6  6  6  3500

julia> df = makeDecisionMatrix(decmat)
6×8 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6     Crt7     Crt8
     │ Float64  Float64  Float64  Float64  Float64  Float64  Float64  Float64
─────┼────────────────────────────────────────────────────────────────────────
   1 │     2.0      1.0      4.0      7.0      6.0      6.0      7.0   3000.0
   2 │     4.0      1.0      5.0      6.0      7.0      7.0      6.0   3500.0
   3 │     3.0      2.0      6.0      6.0      5.0      6.0      8.0   4000.0
   4 │     5.0      1.0      5.0      7.0      6.0      7.0      7.0   3000.0
   5 │     4.0      2.0      5.0      6.0      7.0      7.0      6.0   3000.0
   6 │     3.0      2.0      6.0      6.0      6.0      6.0      6.0   3500.0

julia> weights = [0.293, 0.427, 0.067, 0.027, 0.053, 0.027, 0.053, 0.053];

julia> fns = [maximum, maximum, maximum, maximum, maximum, maximum, maximum, minimum];

julia> result = mabac(df, weights, fns);

julia> result.scores
6-element Array{Float64,1}:
 -0.3113160790692055
 -0.10898274573587217
  0.2003505875974611
  0.0421839209307945
  0.3445172542641278
  0.2003505875974611

julia> result.bestIndex
5
```

# References

Pamučar, D., & Ćirović, G. (2015). The selection of transport and handling resources in logistics centers using Multi-Attributive Border Approximation area Comparison (MABAC). Expert Systems with Applications, 42(6), 3016–3028. doi:10.1016/j.eswa.2014.11.057

Ulutaş, A. (2019). Entropi ve MABAC yöntemleri ile personel seçimi. OPUS–International Journal of Society Researches, 13(19), 1552-1573. DOI: 10.26466/opus.580456
"""
function mabac(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::MABACResult

    row, col = size(decisionMat)

    w = unitize(weights)

    colMax = colmaxs(decisionMat)
    colMin = colmins(decisionMat)

    A = similar(decisionMat)

    for i in 1:row
        for j in 1:col
            if fns[j] == maximum
                @inbounds A[i, j] = (decisionMat[i, j] - colMin[j]) / (colMax[j] - colMin[j])
            elseif fns[j] == minimum
                @inbounds A[i, j] = (decisionMat[i, j] - colMax[j]) / (colMin[j] - colMax[j])
            end
        end
    end

    wA = w * (A .+ 1)

    g = zeros(Float64, col)

    for i in 1:col
        g[i] = geomean(wA[:, i])
    end

    Q = wA .- g'

    scores = zeros(Float64, row)
    for i in 1:row
        scores[i] = sum(Q[i, :])
    end

    rankings = sortperm(scores)

    bestIndex = rankings |> last

    result = MABACResult(
        decisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
"""
    regret(P::MultivariateDistribution, MC::Vector{MultivariateDistribution}, x::AbstractVector)

The regret value of the probability distribution `P` w.r.t the model class `MC`
in the point `x`.
"""
regret(P, MC, x::AbstractVector) = -logpdf(P, x) - minimum(map(p->-logpdf(p,x), MC))
regret(P, MC, X::AbstractMatrix) = -logpdf(P, X) .- min.((-logpdf(p, X) for p in MC)...)

"""
    regretmax(P::MultivariateDistribution, MC::Vector{MultivariateDistribution}, X::AbstractMatrix)

The maximal regret of the probability distribution `P` w.r.t the model class `MC`
for all points of `X`.
"""
regretmax(P, MC, X) = maximum(regret(P,MC,X))

"""
    refinedmdl(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) -> (T, T, T)

Calculate a refined MDL value for the metaclustering `meta` w.r.t. to the model class 'mcr'
for all points of `X`. Returns tuple of MDL, likelihood, and complexity values.
"""
function refinedmdl(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) where {T <: AbstractFloat}
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    cs = counts(mcr)
    ps = map(i->cs[i]./T(sum(cs[i])), meta)
    logpms = mixlogpdf(logps, ps, meta)
    return _refinedmdl(logpms)
end

function _refinedmdl(logpds::AbstractMatrix{T}) where {T <: AbstractFloat}
    n,k = size(logpds)
    minll = minimum(-logpds, dims=2)
    NLL = minll |> sum
    Rₘₐₓ = maximum(-logpds .- minll, dims=1) |> maximum # REG_max(p, M)
    return NLL + Rₘₐₓ, NLL, Rₘₐₓ
end

"""
    nml(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) -> (T, T, T)

Calculate a normalized message length for the metaclustering `meta` w.r.t. to the model class 'mcr'
for all points of `X`. Returns tuple of NML, likelihood, and complexity values.
"""
function nml(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) where {T <: AbstractFloat}
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    cs = counts(mcr)
    ps = map(i->cs[i]./T(sum(cs[i])), meta)
    logpms = mixlogpdf(logps, ps, meta)
    return _nml(logpms)
end

function _nml(logpds::AbstractMatrix)
    n = size(logpds,1)
    minll = minimum(-logpds, dims=2)
    NLL = sum(minll)
    COMP = logsumexp(-minll)
    return NLL + n*COMP, NLL, COMP
end

"""
    nll(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) -> (T, T, T)

Calculate a negative log-likelyhood for the metaclustering `meta` w.r.t. to the model class 'mcr'
for all points of `X`.
"""
function nll(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix{T}) where {T <: AbstractFloat}
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    cs = counts(mcr)
    ps = map(i->cs[i]./T(sum(cs[i])), meta)
    logpms = mixlogpdf(logps, ps, meta)
    minll = minimum(-logpms, dims=2)
    NLL = minll |> sum
    return NLL, NLL, zero(T)
end

"""
    mdldiff(P′::AbstractMixtureModel, MC::Vector{MultivariateDistribution}, X::AbstractMatrix) -> Float64

Returns MDL difference between original and merged P′ = {Pᵢ ∪ Pⱼ} clusterings.
"""
function mdldiff(meta::Vector{Vector{Int}}, mcr::ModelClusteringResult, X::AbstractMatrix; kwargs...)
    logps = hcat((logpdf(p, X) for p in models(mcr))...)
    return mdldiff(meta, logps, assignments(mcr); kwargs...)
end

function mdldiff(meta::Vector{Vector{Int}}, logps::AbstractMatrix{T},
                 assign::Vector{Int}; kwargs...) where {T <: AbstractFloat}
    n = size(logps,2)
    cs = zeros(Int, n)
    for i in assign
        cs[i] += 1
    end
    sz = map(i->sum(cs[i]), meta)
    ps = map(i->cs[i]./sum(cs[i]), meta)
    mrgidxs = vcat((findall(i-> i ∈ m, assign) for m in meta)...)
    vlogps = view(logps, mrgidxs, :)
    logpms = mixlogpdf(vlogps, ps, meta)
    return _mdldiff(logpms, vlogps, meta, sz)
end

function _mdldiff(logP′::AbstractMatrix{T}, logPs::AbstractMatrix{T},
                  meta::Vector{Vector{Int}}, sz::Vector{Int}) where {T <: AbstractFloat}
    n = sum(sz)
    NLL = maximum(logP′, dims=2) |> sum
    pv = exp.(logP′)
    s1 = max.((minimum(view(-logPs,:,i), dims=2) for i in meta)...)
    s2 = -log.(sz[2]./view(pv,:,1) .+ sz[1]./view(pv,:,2))
    return NLL + (maximum(s1 .- s2) - log(n))*n  # normalize cost
end

max_stars_repo_name: wildart/PLMC

max_issues_repo_name: wildart/PLMC

max_forks_repo_name: wildart/PLMC


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "qlopt"
version = v"0.1.0"

# Collection of sources required to build qlopt
sources = [
    "https://github.com/jgoldfar/qlopt.git" =>
    "87c839112d76a7957dfcbdc09491acdce92cdba4",

    "https://dl.bintray.com/boostorg/release/1.69.0/source/boost_1_69_0.tar.gz" =>
    "9a2c2819310839ea373f42d69e733c339b4e9a19deab6bfec448281554aa4dbb",

    "http://bitbucket.org/eigen/eigen/get/3.3.7.tar.gz" =>
    "7e84ef87a07702b54ab3306e77cea474f56a40afa1c0ab245bb11725d006d0da",

]

# Bash recipe for building across all platforms
script = raw"""
# Boost build
cd $WORKSPACE/srcdir/boost_1_69_0/
./bootstrap.sh --with-libraries="system,math,test" --prefix=$prefix
./b2 -d0 -q --with-math --with-test --with-system install

# Eigen build
cd $WORKSPACE/srcdir/eigen-eigen-323c052e1731/
mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make install

# QLopt build
cd $WORKSPACE/srcdir/qlopt/
make all CXXFLAGS_EXTRA="-I${prefix}/include"
# Install script. TODO: Add install target
mkdir -p $prefix/bin
mv bin/* $prefix/bin
mkdir -p $prefix/lib
mv build/*.a $prefix/lib/
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:x86_64, libc=:musl)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "lotkaVolterraExample", :lotkaVolterraExample),
    LibraryProduct(prefix, "libboost_prg_exec_monitor", Symbol("libboost_prg_exec_monitor.so")),
    LibraryProduct(prefix, "libboost_math_c99l", Symbol("libboost_math_c99l.so")),
    LibraryProduct(prefix, "libboost_math_tr1f", Symbol("libboost_math_tr1f.so")),
    LibraryProduct(prefix, "libboost_math_tr1l", Symbol("libboost_math_tr1l.so")),
    LibraryProduct(prefix, "libboost_math_tr1", Symbol("libboost_math_tr1.so")),
    LibraryProduct(prefix, "libboost_math_c99f", Symbol("libboost_math_c99f.so")),
    LibraryProduct(prefix, "libboost_chrono", Symbol("libboost_chrono.so")),
    ExecutableProduct(prefix, "lotkaVolterraWithType2RegsExample", :lotkaVolterraWithType2RegsExample),
    LibraryProduct(prefix, "libboost_math_c99", Symbol("ibboost_math_c99.so")),
    LibraryProduct(prefix, "libboost_timer", Symbol("libboost_timer.so")),
    LibraryProduct(prefix, "libboost_unit_test_framework", Symbol("libboost_unit_test_framework.so")),
    LibraryProduct(prefix, "libboost_system", Symbol("libboost_system.so")),
    ExecutableProduct(prefix, "ThreeStepMetabolicNetworkExample", :ThreeStepMetabolicNetworkExample)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jgoldfar/qloptBuilder

max_issues_repo_name: jgoldfar/qloptBuilder

max_forks_repo_name: jgoldfar/qloptBuilder


text:
println("Init Python deps")
import PyCall: pyimport

const PIP_PACKAGES = ["PyCIM","cimpy"]

sys = pyimport("sys")
subprocess = pyimport("subprocess")
subprocess.check_call([sys.executable, "-m", "pip", "install", "--user", "--upgrade", "--force-reinstall", PIP_PACKAGES...])


text:
#!/usr/bin/env julia

using LinearResponseVariationalBayes

include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "test", "test_hessian_reparam.jl"))
include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "test", "test_exponential_families.jl"))
include(joinpath(Pkg.dir("LinearResponseVariationalBayes"), "test", "test_sub_models.jl"))


text:
#
# impl.jl -
#
# Low level implementation to directly call the C functions of Nüvü Camēras
# SDK.
#

module Impl
include("types.jl")
include("constants.jl")
include("ccalls.jl")
end


text:
# entry point for calculating node coordinates and metrics for either
# linear or curvilinear meshes
# also functions used by both linear and curvilienar calulcation

# files containing the code that do the calculations
include("metrics_linear.jl")
include("metrics_curvilinear.jl")

#------------------------------------------------------------------------------
# entry point for both linear and curvilinear metric calculation
"""
  This function gets all the node coordinates, face coordinates, normal
  vectors, and metrics terms for both linear and curvilinear elements

  For linear meshes, dxidx and jac are calculated at the volume nodes and
  interpolated to the face.  For curvilinear meshes the face normals are
  calculated directly, and dxidx and jac are not calculated (the arrays have
  dimensions of all zeros).  This is because it is not clear how to define
  jac and dxidx at the face such that they satisfy the metric invarients.
  All code should eventually use mesh.nrm_face rather than dxidx_face and
  jac_face (and the corresponding fields for boundaryfaces and shared faces),
  but until that is done, linear meshes will still calculate the face
  quantities but curvilinear meshes will not

  All functions called by this function use smart allocators, so it is
  efficient to call this function multiple times (for example, after
  updating the mesh coordinates)

  This function is collective mesh.comm.

  Fields populated by this function:
    Linear meshes:
      vert_coords
      coords
      coords_bndry
      coords_interface
      coords_sharedface
      dxidx
      dxidx_face
      dxidx_sharedface
      dxidx_bndry
      jac
      jac_face
      jac_sharedface
      jac_bndry
      nrm_bndry
      nrm_face
      nrm_sharedface
      remote_metrics
      min_el_size
      volume

    Curvilinear meshes:
      vert_coords
      coords
      coords_bndry
      coords_interface
      coords_sharedface
      dxidx
      jac
      nrm_bndry
      nrm_face
      nrm_sharedface
      min_el_size
      remote_metrics
      volume

"""
function getAllCoordinatesAndMetrics(mesh, sbp, opts; verify=true)

  if opts["use_linear_metrics"]
#   if mesh.coord_order == 1
    @assert mesh.coord_order == 1
    getCoordinates(mesh, sbp)  # store coordinates of all nodes into array
    getMetrics(mesh, sbp)

    if mesh.isInterpolated
      interpolateCoordinatesAndMetrics(mesh)
    end

    getFaceNormals(mesh, sbp)
  else  # curvilinear
    # 3rd order and above are not supported yet
    @assert mesh.coord_order <= 2

    getMeshCoordinates(mesh, sbp)
    getFaceCoordinatesAndNormals(mesh, sbp)
    getCurvilinearCoordinatesAndMetrics(mesh, sbp)

  end

  # make sure the mapping jacobian is > 0
  if verify
    checkMapping(mesh)
  end

  # calculate things that depend on the above
  mesh.min_el_size = getMinElementSize(mesh)
  mesh.volume = calcVolumeIntegral(mesh, sbp)

  # send metric information in parallel
  exchangeMetricInfo(mesh, sbp)

  return nothing
end

"""
  This function recalculates the volume node coordinates and metrics from
  the vertex coordinates in mesh.vert_coords.  This function populates all
  the same fields as [`getAllCoordinatesAndMetics`](@ref) except for
  `mesh.vert_coords`.

  This function does not support opts["use_linear_metrics"].

  This function is collective on mesh.comm

  **Inputs**

   * mesh: the mesh
   * sbp: the SBP operator
   * opts: options dictionary

  **Keyword Arguments**

   * verify: if true, check the mapping jacobian is valid and throw an exception
             if it is not, default true
"""
function recalcCoordinatesAndMetrics(mesh, sbp, opts; verify=true)

  @assert !opts["use_linear_metrics"]
  @assert mesh.coord_order <= 2

  getFaceCoordinatesAndNormals(mesh, sbp)
  getCurvilinearCoordinatesAndMetrics(mesh, sbp)

  # make sure the mapping jacobian is > 0
  if verify
    checkMapping(mesh)
  end

  # calculate things that depend on the above
  mesh.min_el_size = getMinElementSize(mesh)
  mesh.volume = calcVolumeIntegral(mesh, sbp)

  # send metric information in parallel
  exchangeMetricInfo(mesh, sbp)

  return nothing
end

"""
  TODO: this is the reverse mode of `recalcCoordinatesAndMetrics`

  Reverse mode of getAllCoordinatesAndMetrics.  Back propigates
  mesh.nrm_*_bar, mesh.dxidx_bar, mesh.jac_bar to mesh.vert_coords_bar.
  Also does communication to get the contribution from `mesh.remote_metrics_bar`

  Users should generally call zeroBarArrays() before calling this function
  a second time, to zero out all intermediate arrays.

  This function also recalculates the face normal vectors and coordinates,
  overwriting the relevent arrays in the mesh object.

  This function does not yet support the reverse mode of mesh.remote_metrics.
"""
function getAllCoordinatesAndMetrics_rev(mesh, sbp, opts)

  # parallel communication to get contributions from mesh.remote_metrics_bar
  exchangeMetricInfo_rev(mesh, sbp)

  if opts["use_linear_metrics"]
#   if mesh.coord_order == 1
    @assert mesh.coord_order == 1
    interpolateMapping_rev(mesh)
    getVertCoords_rev(mesh, sbp)
  else
    @assert mesh.coord_order <= 2

    getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
    getFaceCoordinatesAndNormals_rev(mesh, sbp)
  end

  return nothing
end

"""
  This method outputs the result in vector form, rather than the 3D array form
  of `mesh.vert_coords_bar`.  In parallel, for coordinates that live on shared
  entities, the entity on the owning process will get the value.
  The value will be zero on non-owning processes.

  **Inputs**

   * mesh
   * sbp
   * opts

  **Inputs/Outputs**

   * xvec_bar: vector, length `mesh.dim*mesh.coord_numNodes` to overwrite
               with the result of back propigating the metrics

  **Keyword Arguments**

   * parallel: if true (default), do parallel communication for shared entities
               (as described above).  Otherwise, do only local back-propigation.
"""
function getAllCoordinatesAndMetrics_rev(mesh, sbp, opts, xvec_bar::AbstractVector; parallel::Bool=true)
# output 1D vector

  #TODO: it would be better for parallel efficiency to do:
  #  1. reverse mode for elements on parallel boundaries
  #  2. start communication for coords3DTo1D
  #  3. reverse mode for elements on interior
  #  4. finish coords3DTo1D

  getAllCoordinatesAndMetrics_rev(mesh, sbp, opts)
  coords3DTo1D(mesh, mesh.vert_coords_bar, xvec_bar, parallel=parallel)

  return nothing
end

#------------------------------------------------------------------------------
# functions common to linear and curvilinear
"""
  This function allocates the arrays of normal vector for interfaces, boundary
  faces, and shared faces.

  For interfaces, the normal is calculated for elementL
"""
function allocateNormals(mesh::PumiMeshDG{Tmsh}, sbp) where Tmsh

  dim = mesh.dim
  numfacenodes = mesh.numNodesPerFace

  if !isFieldDefined(mesh, :nrm_bndry, :nrm_face, :nrm_sharedface)
    mesh.nrm_bndry = Array{Tmsh}(dim, numfacenodes, mesh.numBoundaryFaces)
    mesh.nrm_face = Array{Tmsh}(mesh.dim, numfacenodes, mesh.numInterfaces)
    mesh.nrm_sharedface = Array{Array{Tmsh, 3}}(mesh.npeers)

    for i=1:mesh.npeers
      mesh.nrm_sharedface[i] = Array{Tmsh}(mesh.dim, numfacenodes, mesh.peer_face_counts[i])
    end

    # adjoint parts
    mesh.nrm_bndry_bar = zeros(mesh.nrm_bndry)
    mesh.nrm_face_bar = zeros(mesh.nrm_face)
    mesh.nrm_sharedface_bar = Array{Array{Tmsh, 3}}(mesh.npeers)
    for i=1:mesh.npeers
      mesh.nrm_sharedface_bar[i] = zeros(mesh.nrm_sharedface[i])
    end

  else
    fill!(mesh.nrm_bndry, 0.0)
    fill!(mesh.nrm_face, 0.0)
    for i=1:mesh.npeers
      fill!(mesh.nrm_sharedface[i], 0.0)
    end
#=
    fill!(mesh.nrm_bndry_bar, 0.0)
    fill!(mesh.nrm_face_bar, 0.0)
    for i=1:mesh.npeers
      fill!(mesh.nrm_sharedface[i], 0.0)
    end
=#
  end

  return nothing
end

"""
  Send and recieve the metric information for elements on parallel boundaries.
  This function is collective

  **Inputs**

   * mesh: a DG mesh. mesh.remote_metrics, mesh.remote_metrics_bar are
           populated by this function
"""
function exchangeMetricInfo(mesh::PumiMeshDG{Tmsh}, sbp) where Tmsh

  remote_metrics = Array{RemoteMetrics{Tmsh}}(mesh.npeers)  # receive buffers
  remote_metrics_bar = Array{RemoteMetrics{Tmsh}}(mesh.npeers)
  local_metrics = Array{RemoteMetrics{Tmsh}}(mesh.npeers)  # send buffers
  send_reqs = Array{MPI.Request}(4, mesh.npeers)
  recv_reqs = Array{MPI.Request}(4, mesh.npeers)

  # allocate the arrays and post the MPI receives
  for i=1:mesh.npeers
    remote_metrics[i] = RemoteMetrics(mesh, i, islocal=false)
    remote_metrics_bar[i] = RemoteMetrics(mesh, i, islocal=false)
    local_metrics[i] = RemoteMetrics(mesh, i, islocal=true)

    obj = remote_metrics[i]
    peernum = mesh.peer_parts[i]
    recv_reqs[1, i] = MPI.Irecv!(obj.vert_coords, peernum, 1, mesh.comm)
    recv_reqs[2, i] = MPI.Irecv!(obj.coords, peernum, 2, mesh.comm)
    recv_reqs[3, i] = MPI.Irecv!(obj.dxidx, peernum, 3, mesh.comm)
    recv_reqs[4, i] = MPI.Irecv!(obj.jac, peernum, 4, mesh.comm)
  end

  # extract the local values and send them
  for i=1:mesh.npeers
    obj = local_metrics[i]
    peernum = mesh.peer_parts[i]
    getLocalMetrics(mesh, obj)

    send_reqs[1, i] = MPI.Isend(obj.vert_coords, peernum, 1, mesh.comm)
    send_reqs[2, i] = MPI.Isend(obj.coords, peernum, 2, mesh.comm)
    send_reqs[3, i] = MPI.Isend(obj.dxidx, peernum, 3, mesh.comm)
    send_reqs[4, i] = MPI.Isend(obj.jac, peernum, 4, mesh.comm)
  end

  # wait for all communication to finish
  # this is not strictly necessary, the real requirement is that the
  # communication finishes before the mesh constructor exits
  for i=1:length(send_reqs)
    MPI.Wait!(send_reqs[i])
    MPI.Wait!(recv_reqs[i])
  end

  mesh.remote_metrics = remote_metrics
  mesh.remote_metrics_bar = remote_metrics_bar

  return nothing
end

"""
  Reverse mode of `exchangeMetricInfo`.
"""
function exchangeMetricInfo_rev(mesh::PumiMeshDG{Tmsh}, sbp) where Tmsh

  local_metrics_bar = Array{RemoteMetrics{Tmsh}}(mesh.npeers)  # send buffers
  send_reqs = Array{MPI.Request}(4, mesh.npeers)
  recv_reqs = Array{MPI.Request}(4, mesh.npeers)

  # allocate the arrays and post the MPI receives
  for i=1:mesh.npeers
    local_metrics_bar[i] = RemoteMetrics(mesh, i, islocal=true)

    obj = local_metrics_bar[i]
    peernum = mesh.peer_parts[i]
    recv_reqs[1, i] = MPI.Irecv!(obj.vert_coords, peernum, 5, mesh.comm)
    recv_reqs[2, i] = MPI.Irecv!(obj.coords, peernum, 6, mesh.comm)
    recv_reqs[3, i] = MPI.Irecv!(obj.dxidx, peernum, 7, mesh.comm)
    recv_reqs[4, i] = MPI.Irecv!(obj.jac, peernum, 8, mesh.comm)
  end

  # extract the local values and send them
  for i=1:mesh.npeers
    obj = mesh.remote_metrics_bar[i]
    peernum = mesh.peer_parts[i]

    send_reqs[1, i] = MPI.Isend(obj.vert_coords, peernum, 5, mesh.comm)
    send_reqs[2, i] = MPI.Isend(obj.coords, peernum, 6, mesh.comm)
    send_reqs[3, i] = MPI.Isend(obj.dxidx, peernum, 7, mesh.comm)
    send_reqs[4, i] = MPI.Isend(obj.jac, peernum, 8, mesh.comm)
  end

  # wait for communications to finish and unpack the receive buffer
  wait_leaders = Array{MPI.Request}(mesh.npeers)
  for i=1:mesh.npeers
    # wait on the dxidx array because it is probably largest/slowest
    wait_leaders[i] = recv_reqs[3, i]
  end

  for i=1:mesh.npeers
    idx, state = MPI.Waitany!(wait_leaders)
    # wait for other communications on the same object
    MPI.Wait!(recv_reqs[1, idx])
    MPI.Wait!(recv_reqs[2, idx])
    MPI.Wait!(recv_reqs[4, idx])

    getLocalMetrics_rev(mesh, local_metrics_bar[idx])
  end

  # wait for sends to complete before exiting
  for i=1:mesh.npeers
    for j=1:4
      MPI.Wait!(send_reqs[j, i])
    end
  end

  return nothing
end

"""
  Copy the metric values from the main mesh arrays into the RemoteMetrics
  arrays.

  **Inputs**

   * mesh: a mesh object

  **Inputs/Outputs**

   * obj: a RemoteMetrics object
"""
function getLocalMetrics(mesh, obj::RemoteMetrics)

  @assert obj.islocal

  numEl = size(obj.vert_coords, 3)
  elnums = mesh.local_element_lists[obj.peer_idx]
  for i=1:numEl
    elnum = elnums[i]

    # vert_coords
    for j=1:size(obj.vert_coords, 2)
      for k=1:size(obj.vert_coords, 1)
        obj.vert_coords[k, j, i] = mesh.vert_coords[k, j, elnum]
      end
    end

    # coords
    for j=1:size(obj.coords, 2)
      for k=1:size(obj.coords, 1)
        obj.coords[k, j, i] = mesh.coords[k, j, elnum]
      end
    end

    # dxidx
    for j=1:mesh.numNodesPerElement
      for d2=1:mesh.dim
        for d1=1:mesh.dim
          obj.dxidx[d1, d2, j, i] = mesh.dxidx[d1, d2, j, elnum]
        end
      end
    end

    # jac
    for j=1:mesh.numNodesPerElement
      obj.jac[j, i] = mesh.jac[j, elnum]
    end

  end  # end loop i

  return nothing
end

"""
  Reverse mode of `getLocalMetrics`.

  **Inputs**

   * mesh
   * obj: a `RemoteMetrics` object.  This object should already have received
          the bar metrics from the remote process.
"""
function getLocalMetrics_rev(mesh, obj::RemoteMetrics)

  @assert obj.islocal

  numEl = size(obj.vert_coords, 3)
  elnums = mesh.local_element_lists[obj.peer_idx]
  for i=1:numEl
    elnum = elnums[i]

    # vert_coords
    for j=1:size(obj.vert_coords, 2)
      for k=1:size(obj.vert_coords, 1)
        mesh.vert_coords_bar[k, j, elnum] += obj.vert_coords[k, j, i]
      end
    end

    # coords
    # TODO: support mesh.coords_bar
    #for j=1:size(obj.coords, 2)
    #  for k=1:size(obj.coords, 1)
    #    mesh.coords_bar[k, j, elnum] += obj.coords[k, j, i]
    #  end
    #end

    # dxidx
    for j=1:mesh.numNodesPerElement
      for d2=1:mesh.dim
        for d1=1:mesh.dim
          mesh.dxidx_bar[d1, d2, j, elnum] += obj.dxidx[d1, d2, j, i]
        end
      end
    end

    # jac
    for j=1:mesh.numNodesPerElement
      mesh.jac_bar[j, elnum] += obj.jac[j, i]
    end

  end  # end loop i

  return nothing
end


text:
using AdaOPS
using Test

using POMDPs
using POMDPModels
using POMDPSimulators
using Random
using POMDPModelTools
using ParticleFilters
using BeliefUpdaters
using StaticArrays
using POMDPPolicies
using Plots
ENV["GKSwstype"] = "100"
theme(:mute)

# include("baby_sanity_check.jl")
include("independent_bounds.jl")

pomdp = BabyPOMDP()
pomdp.discount = 1.0
p = solve(AdaOPSSolver(tree_in_info=true), pomdp)
Random.seed!(p, 1)

K = 10
b0 = initialstate(pomdp)
o = false
tval = 7.0
a, info = action_info(p, b0)
tree = info[:tree]
tree.obs[2] = o
b = WPFBelief([rand(b0)], [1.0], 2, 1, tree, o)

pol = FeedWhenCrying()
rng = MersenneTwister(2)

# AbstractParticleBelief interface
@testset "WPFBelief" begin
    @test n_particles(b) == 1
    s = particle(b,1)
    @test rand(rng, b) == s
    @test pdf(b, rand(rng, b)) == 1
    sup = support(b)
    @test length(sup) == 1
    @test first(sup) == s
    @test mode(b) == s
    @test mean(b) == s
    @test first(particles(b)) == s
    @test first(weights(b)) == 1.0
    @test first(weighted_particles(b)) == (s => 1.0)
    @test weight_sum(b) == 1.0
    @test weight(b, 1) == 1.0
    @test currentobs(b) == o
    @test history(b)[end].o == o
    @test initialize_belief(KMarkovUpdater(2), b)[end] == o
end

# Light Dark Test
pomdp = LightDark1D()
random = solve(RandomSolver(), pomdp)
Base.convert(::Type{SVector{1,Float64}}, s::LightDark1DState) = SVector{1,Float64}(s.y)
grid = StateGrid(range(-10, stop=15, length=26))
bds = IndependentBounds(FORollout(random), pomdp.correct_r)
solver = AdaOPSSolver(bounds=bds,
                    grid=grid,
                    m_min=30,
                    delta=1.0,
                    tree_in_info=true
                    )
planner = solve(solver, pomdp)
action(planner, initialstate(pomdp))
hr = HistoryRecorder(max_steps=50)
@time hist = simulate(hr, pomdp, planner)
hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# BabyPOMDP Test
Base.convert(::Type{SVector{1,Float64}}, s::Bool) = SVector{1,Float64}(s)
# Type stability
pomdp = BabyPOMDP()
bds = (pomdp, b)->(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = AdaOPSSolver(bounds=bds,
                      rng=MersenneTwister(4),
                      m_min=100,
                      tree_in_info=true
                     )
p = solve(solver, pomdp)

b0 = initialstate(pomdp)
D, Depth = @inferred AdaOPS.build_tree(p, b0)
@inferred action_info(p, b0)
a, info = action_info(p, b0)
info_analysis(info)
@inferred AdaOPS.explore!(D, 1, p)
Δu, Δl = @inferred AdaOPS.expand!(D, D.b, p)
@inferred AdaOPS.backup!(D, 1, p, Δu, Δl)
@inferred AdaOPS.next_best(D, 1, p)
@inferred AdaOPS.excess_uncertainty(D, 1, AdaOPS.solver(p).xi, p)
@inferred action(p, b0)

pomdp = BabyPOMDP()

# constant bounds
bds = IndependentBounds(reward(pomdp, true, false)/(1-discount(pomdp)), 0.0)
solver = AdaOPSSolver(bounds=bds, m_min=200, tree_in_info=true, num_b=10_000)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=20)
@time hist = simulate(hr, pomdp, planner)
hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# SemiPORollout lower bound
bds = IndependentBounds(SemiPORollout(FeedWhenCrying()), 0.0)
solver = AdaOPSSolver(bounds=bds, m_min=200, tree_in_info=true, num_b=10_000)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=20)
@time hist = simulate(hr, pomdp, planner)
hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# PO policy lower bound
bds = IndependentBounds(PORollout(FeedWhenCrying(), PreviousObservationUpdater()), 0.0)
solver = AdaOPSSolver(bounds=bds, m_min=200, tree_in_info=true, num_b=10_000)
planner = solve(solver, pomdp)
hr = HistoryRecorder(max_steps=20)
@time hist = simulate(hr, pomdp, planner)
hist_analysis(hist)
println("Discounted reward is $(discounted_reward(hist))")

# from README:
using POMDPs, POMDPModels, POMDPSimulators, AdaOPS

pomdp = TigerPOMDP()

solver = AdaOPSSolver(bounds=IndependentBounds(-20.0, 0.0))
planner = solve(solver, pomdp)

for (s, a, o) in stepthrough(pomdp, planner, "s,a,o", max_steps=10)
    println("State was $s,")
    println("action $a was taken,")
    println("and observation $o was received.\n")
end


text:
# Line plots

function PlotLine(label_id, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, y::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, x, y, count, offset, stride)
end

function PlotLine(label_id, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, y::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotLine(label_id, Float64.(x), Float64.(y), count, offset, stride)
end

function PlotLine(label_id, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer, xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, x, count, xscale, x0, offset, stride)
end

function PlotLine(label_id, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer, xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotLine(label_id, Float64.(x), count, xscale, x0, offset, stride)
end

function PlotLine(x::AbstractArray{T}, y::AbstractArray{T};
                  count::Integer = min(length(x), length(y)), offset::Integer = 0,
                  stride::Integer = 1, label_id::String = "") where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotLine(x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real, T2<:Real}
    PlotLine(promote(x,y)...; kwargs...)
end

function PlotLine(y::AbstractArray{T}; label_id::String="", count::Integer=length(y),
                  xscale::Real = 1.0, x0::Real = 0.0, offset::Integer=0,
                  stride::Integer=1) where {T<:ImPlotData}
    LibCImPlot.PlotLine(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotLine(x::UnitRange{<:Integer}, y::AbstractArray{T}; xscale::Real = 1.0,
                  x0::Real = 0.0, label_id::String="") where {T<:ImPlotData}

    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride::Cint = sizeof(T)
    LibCImPlot.PlotLine(label_id, y, count,  xscale, x0, offset, stride)
end

function PlotLine(x::StepRange, y::AbstractArray{T}; xscale::Real = 1.0, x0::Real = 0.0,
                  label_id::String= "") where {T<:ImPlotData}

    x.stop < 1 && throw("Range out of bounds")
    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(x.step * sizeof(T))
    LibCImPlot.PlotLine(label_id, y, count, xscale, x0, offset, stride)
end

PlotLine(x::UnitRange{<:Integer}, y::AbstractArray{<:Real}; kwargs...) = PlotLine(x, Float64.(y); kwargs...)
PlotLine(x::StepRange, y::AbstractArray{<:Real}; kwargs...) = PlotLine(x, Float64.(y); kwargs...)

# xfield, yfield should be propertynames of eltype(structvec)
function PlotLine(
    structvec::Vector{T}, xfield::Symbol, yfield::Symbol;
    count::Integer = length(structvec), offset::Integer = 0,
    stride::Integer = 1, label_id::String = ""
) where T

    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = (pointer(structvec, 1) + x_offset) |> Ptr{Tx}
    y_ptr = (pointer(structvec, 1) + y_offset) |> Ptr{Ty}

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own = false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own = false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    LibCImPlot.PlotLine(label_id, x, y, count, offset, stride)
end

PlotLineG(label_id, getter, data, count, offset = 0) =
LibCImPlot.PlotLineG(label_id, getter, data, count, offset)


text:
using LogNumbers
using Test

@testset "LogNumber.jl" begin
    a = LogNumber(exp(1))
    @test a.sign == 1
    @test a.exponent ≈ 1
    a = LogNumber(-exp(1))
    @test a.sign == -1
    @test a.exponent ≈ 1

end

const test_pairs = (
    (0.99, 0.876),
    (10.99, 0.01),
    (-0.98, 0.87),
    (0.88,-0.99),
    (-1.88, -0.98)
)

@testset "+" begin
    for (x, y) in test_pairs
        a = LogNumber(x)
        b = LogNumber(y)
        @test a+b isa LogNumber
        @test to_number(a+b) ≈ x + y
    end
end

@testset "-" begin
    for (x, y) in test_pairs
        a = LogNumber(x)
        b = LogNumber(y)
        @test a-b isa LogNumber
        @test to_number(a-b) ≈ x - y
    end
end

@testset "*" begin
    for (x, y) in test_pairs
        a = LogNumber(x)
        b = LogNumber(y)
        @test a*b isa LogNumber
        @test to_number(a*b) ≈ x * y
    end
end

@testset "/" begin
    for (x, y) in test_pairs
        a = LogNumber(x)
        b = LogNumber(y)
        @test a/b isa LogNumber
        @test to_number(a/b) ≈ x / y
    end
end


text:
const GENS =  64 # number of conventional generators
const PEAKERS =  80 # number of conventional generators
const STORG_UNITS = 8 # number of storage units
const N_ZONES = 2
const M_ZONES = 2
const BLOCKS = 7

const INITIAL_DAY = 1
const FINAL_DAY = 60
# Running time for the 1st Week-Ahead UC run.
# INIT_HR_FUCR=0 means the FUCR's optimal outcomes are ready at 00:00
const INIT_HR_FUCR = 6
# Running time for the 2nd WA unit commitment run.
# INIT_HR_SUCR=17 means the SUCR's optimal outcomes are ready at 17:00
const INIT_HR_SUCR = 17
# HRS_FUCR = 168-INIT_HR_FUCR, and INIT_HR_FUCR runs from 0 to 23;
# INIT_HR_FUCR=6 ==> HRS_FUCR =162
const HRS_FUCR = 162
# HRS_SUCR = 168-INIT_HR_SUCR, and INIT_HR_SUCR runs from 17 to 23;
# INIT_HR_SUCR=17 ==> HRS_FUCR =168
const HRS_SUCR = 151

const LOAD_SHED_PENALTY = 3000; # Load-Shedding Penalty
const OVERGEN_PENALTY = 4000; # Over-generation Penalty
const LOAD_SHED_MAX = 100; # Load-Shedding Penalty
const OVERGEN_MAX = 100; # Over-generation Penalty
const VIOLATION_PENALTY = 500;
const VIOLATION_MAX = 10;
#Solver's optimality gap that serves as the optimization termination criteria
const SOLVER_EPGAP =0.005;

const FILE_ACCESS_OVER = "w+";
const FILE_ACCESS_APPEND = "a+";

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using Test
using LogicCircuits
import .LogicCircuits.IO:
   load_cnf, load_dnf
using Pkg.Artifacts

include("../helper/ValidateSdd.jl")

@testset "CNF file parser tests" begin

   cnfs = [
            ("easy","C17_mince",32,92,45)
            ("easy","majority_mince",32,132,61)
            ("easy","b1_mince",8,169,84)
            ("easy","cm152a_mince",2048,127,62)
            ("iscas89","s208.1.scan",262144,1942,927)
          ]

   for (suite, name, count, size, nodes) in cnfs

      cnf = zoo_cnf("$suite/$name.cnf")
      vtree = zoo_vtree("$suite/$name.min.vtree");

      mgr = SddMgr(TrimSddMgr, vtree)
      # cnfΔ = @time compile_cnf(mgr, cnf)
      cnfΔ = node2dag(compile_cnf(mgr, cnf))

      # println("node count ($suite/$name): ", num_nodes(cnfΔ))
      # println("edge count ($suite/$name): ", num_edges(cnfΔ))
      # println("Final SDD model count: ", model_count(cnfΔ))

      validate(cnfΔ)

      @test model_count(cnfΔ) == count
      @test sdd_size(cnfΔ) == size
      @test sdd_num_nodes(cnfΔ) == nodes

   end

end


text:
# Cartesian categories
quote
struct Product{N, A, T} <: Arrow{A, T}
    factors::NTuple{N, Arrow{A, C} where C}
    # TODO: checks
    Product(arrows...) = begin
        N = length(arrows)
        if N == 0
            # TODO: how to best define this?
            Identity{Nothing}()
        else
            src = Nothing
            T = []
            # This is some weird logic to let you
            # take products where some arrow have source A and
            # some have Nothing.  This is done by replacing each
            # arrow α :: Nothing --> T with α' = α o Terminal{A}
            # This is pretty important to make things work nicely with
            # constants.
            for a in arrows
                push!(T, target(a))
                s = source(a)
                if s != Nothing
                    if src == Nothing
                        src = s
                    elseif src != s
                        error("Incompatible source types $(src), $s")
                    end
                end
            end

            if src != Nothing
                arrows = tuple([source(a) == Nothing ? Cat.compose(a, Terminal{src}()) : a for a in arrows]...)
            end

            new{N, src, Tuple{T...}}(arrows)
        end
    end
end

struct Proj{A<:Tuple, T} <: Arrow{A, T}
    m::Int64
end

function Cat.compose(g::Proj{A, T}, f::Product{N, R, A}) where {A,R,N,T}
    f.factors[g.m]
end

Base.getindex(m::Product{N, A, T}, k::Int64) where {N, A, T} = begin
    if 0 < k <= N
        m.factors[k]
    else
        error("")
    end
end

Base.getindex(m::Arrow{A, B}, k::Int64) where {A, B<:Tuple} = begin
    Cat.compose(Proj{B, B.parameters[k]}(k), m)
end
end # quote


text:
module PkgMining

using TOML
using ProgressMeter
using DataFrames

export pkg_mining

function pkg_mining()
    general = joinpath(homedir(), ".julia", "registries", "General")
    registry = TOML.parsefile(joinpath(general, "Registry.toml"))

    df = DataFrame(
        name = String[],
        uuid = String[],
        owner = String[],
        firstversion = String[],
        lastversion = String[],
        repo = String[],
        sha1 = String[]
    )

    @showprogress "Mining $(length(registry["packages"])) packages: " for (uuid, info) in
                                                                        registry["packages"]
        pkg = TOML.parsefile(joinpath(general, info["path"], "Package.toml"))
        versions = TOML.parsefile(joinpath(general, info["path"], "Versions.toml"))
        firstversion, lastversion = string.(extrema(VersionNumber.(keys(versions))))
        sha1 = versions[lastversion]["git-tree-sha1"]
        owner = split(pkg["repo"], '/')[end-1]
        push!(
            df,
            (
                info["name"],
                uuid,
                owner,
                firstversion,
                lastversion,
                pkg["repo"],
                sha1
            )
        )

    end

    return df
end

end


text:
import ApproxFun: dotu

# This solves as a boundary value provblem

stieltjesbackward(S::Space,z::Number) = JacobiZ(S,z)\[stieltjesmoment(S,0,z)]

# This solves via forward substitution
function forwardsubstitution!(ret,B,n,μ1,μ2)
    if n≥1
        ret[1]=μ1
    end
    if n≥2
        ret[2]=μ2
    end
    if n≥3
        for k=2:n-1
            ret[k+1]=-(B[k,k-1]*ret[k-1]+B[k,k]*ret[k])/B[k,k+1]
        end
    end
    ret
end

forwardsubstitution(R,n,μ1,μ2)=forwardsubstitution!(Array(promote_type(eltype(R),typeof(μ1),typeof(μ2)),n),R,n,μ1,μ2)

stieltjesforward(sp::Space,n,z,s...)=forwardsubstitution(JacobiZ(sp,z),n,
                        stieltjesmoment(sp,0,z,s...),stieltjesmoment(sp,1,z,s...))

function stieltjesintervalrecurrence(S,f::AbstractVector,z)
    tol=1./floor(Int,sqrt(length(f)))
    if (abs(real(z))≤1.+tol) && (abs(imag(z))≤tol)
        cfs = stieltjesforward(S,length(f),z)
        dotu(cfs,f)
    else
        cfs = stieltjesbackward(S,z)
        dotu(cfs,f)
    end
end
stieltjesintervalrecurrence(S,f::AbstractVector,z::AbstractArray) = reshape(promote_type(eltype(f),eltype(z))[ stieltjesintervalrecurrence(S,f,z[i]) for i in eachindex(z) ], size(z))

function stieltjes{D<:Interval}(S::PolynomialSpace{D},f,z::Number)
    if domain(S)==Interval()
        #TODO: check tolerance
        stieltjesintervalrecurrence(Legendre(),coefficients(f,S,Legendre()),z)
    else
        stieltjes(setdomain(S,Interval()),f,mobius(S,z))
    end
end

function stieltjes{D<:Interval}(S::PolynomialSpace{D},f,z::Number,s::Bool)
    @assert domain(S)==Interval()

   cfs=stieltjesforward(Legendre(),length(f),z,s)
   dotu(cfs,coefficients(f,S,Legendre()))
end

# Sum over all inverses of fromcanonical, see [Olver,2014]
function stieltjes{SS,L<:Line}(S::Space{SS,L},f,z,s...)
    if domain(S)==Line()
        # TODO: rename tocanonical
        stieltjes(setcanonicaldomain(S),f,tocanonical(S,z),s...) +
            stieltjes(setcanonicaldomain(S),f,(-1-sqrt(1+4z.^2))./(2z))
    else
        stieltjes(setdomain(S,Line()),f,mappoint(domain(S),Line(),z),s...)
    end
end


text:
using IterativeSolvers, KrylovKit, LinearAlgebra

# In this file, we regroud a way to provide linear solver for the Package

abstract type AbstractLinearSolver end

# The function linsolve(J, x) must return whether the solve was successfull and how many steps were required for the solve.

# the following function can be used to cache some factorization, see DefaultLS() case for example
function (ls::AbstractLinearSolver)(J, rhs1, rhs2; kwargs...)
	sol1, flag1, it1 = ls(J, rhs1; kwargs...)
	sol2, flag2, it2 = ls(J, rhs2; kwargs...)
	return sol1, sol2, flag1 & flag2, (it1, it2)
end

####################################################################################################
# The following functions are used for the Continuation of Hopf points and the computation of Floquet multipliers

# this function returns a₀ * I + a₁ .* J and ensures that we don't do unnecessary computations like 0*I + 1*J
function _axpy(J, a₀, a₁)
	if a₀ == 0
		if a₁ == 1
			return J
		else
			return a₁ .* J
		end
	elseif a₀ == 1
		if a₁ == 1
			return I + J
		else
			return I + a₁ .* J
		end
	else
		return a₀ * I + a₁ .* J
	end
end

function _axpy_op(J, v::AbstractArray, a₀, a₁)
	if a₀ == 0
		if a₁ == 1
			return apply(J, v)
		else
			return a₁ .* apply(J, v)
		end
	elseif a₀ == 1
		if a₁ == 1
			return v .+ apply(J, v)
		else
			return v .+ a₁ .* apply(J, v)
		end
	else
		return a₀ .* v .+ a₁ .* apply(J, v)
	end
end

####################################################################################################
# Solvers for default \ operator (backslash)
####################################################################################################
"""
The struct `Default` is used to  provide the backslash operator to our Package
"""
struct DefaultLS <: AbstractLinearSolver end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the options a₀, a₁ are only used for the Hopf Newton / Continuation
function (l::DefaultLS)(J, rhs; a₀ = 0, a₁ = 1, kwargs...)
	return _axpy(J, a₀, a₁) \ rhs, true, 1
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# with multiple RHS. We can cache the factorization in this case
# the options a₀, a₁ are only used for the Hopf Newton / Continuation
function (l::DefaultLS)(J, rhs1, rhs2; a₀ = 0, a₁ = 1, kwargs...)
	Jfact = lu(_axpy(J, a₀, a₁))
	return Jfact \ rhs1, Jfact \ rhs2, true, (1, 1)
end
####################################################################################################
# Solvers for IterativeSolvers
####################################################################################################
@with_kw mutable struct GMRESIterativeSolvers{T, Tl, Tr} <: AbstractLinearSolver
	tol::T = 1e-4							# tolerance for solver
	restart::Int64 = 200					# number of restarts
	maxiter::Int64 = 100
	N::Int64 = 0							# dimension of the problem
	verbose::Bool = false					# display information during iterations
	log::Bool = true						# record information
	initially_zero::Bool = true				# start with zero guess
	Pl::Tl = IterativeSolvers.Identity()	# left preconditioner
	Pr::Tr = IterativeSolvers.Identity()	# right preconditioner
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the optional shift is only used for the Hopf Newton / Continuation
function (l::GMRESIterativeSolvers{T, Tl, Tr})(J, rhs; a₀ = 0, a₁ = 1, kwargs...) where {T, Ts, Tl, Tr}
	# no need to use fancy axpy! here because IterativeSolvers "only" handles AbstractArray
	J_map = v -> _axpy_op(J, v, a₀, a₁)
	Jmap = LinearMap{T}(J_map, l.N, l.N ; ismutating = false)
	res = IterativeSolvers.gmres(Jmap, rhs; tol = l.tol, log = l.log, verbose = l.verbose, restart = l.restart, maxiter = l.maxiter, initially_zero = l.initially_zero, Pl = l.Pl, Pr = l.Pr, kwargs...)
	(res[2].iters >= l.maxiter) && (@warn "IterativeSolvers.gmres iterated maxIter = $(res[2].iters) times without achieving the desired tolerance.\n")
	return res[1], length(res) > 1, res[2].iters
end

@with_kw mutable struct GMRESIterativeSolvers!{T, Tl, Tr} <: AbstractLinearSolver
	tol::T = 1e-4							# tolerance for solver
	restart::Int64 = 200					# number of restarts
	maxiter::Int64 = 100
	N::Int64 = 0							# dimension of the problem
	verbose::Bool = false					# display information during iterations
	log::Bool = true						# record information
	initially_zero::Bool = true				# start with zero guess
	Pl::Tl = IterativeSolvers.Identity()	# left preconditioner
	Pr::Tr = IterativeSolvers.Identity()	# right preconditioner
end

function (l::GMRESIterativeSolvers!{T, Tl, Tr})(J, rhs; kwargs...) where {T, Ts, Tl, Tr}
	# no need to use fancy axpy! here because IterativeSolvers "only" handles AbstractArray
	Jmap = LinearMap{T}((o, v) -> J(o, v), l.N, l.N ; ismutating = true)
	res = IterativeSolvers.gmres(Jmap, rhs; tol = l.tol, log = l.log, verbose = l.verbose, restart = l.restart, maxiter = l.maxiter, initially_zero = l.initially_zero, Pl = l.Pl, Pr = l.Pr, kwargs...)
	(res[2].iters >= l.maxiter) && (@warn "IterativeSolvers.gmres iterated maxIter = $(res[2].iters) times without achieving the desired tolerance.\n")
	return res[1], length(res) > 1, res[2].iters
end
####################################################################################################
# Solvers for KrylovKit
####################################################################################################
@with_kw mutable struct GMRESKrylovKit{T, Tl} <: AbstractLinearSolver
	dim::Int64 = KrylovDefaults.krylovdim # Krylov Dimension
	atol::T  = KrylovDefaults.tol		  # absolute tolerance for solver
	rtol::T  = KrylovDefaults.tol		  # relative tolerance for solver
	restart::Int64 = 200				  # number of restarts
	maxiter::Int64 = KrylovDefaults.maxiter
	verbose::Int64 = 0
	issymmetric::Bool = false				# if the linear map is symmetric, only meaningful if T<:Real
	ishermitian::Bool = false 				# if the linear map is hermitian
	isposdef::Bool    = false 				# if the linear map is positive definite
	Pl::Tl			  = nothing				# left preconditioner
end

# this function is used to solve (a₀ * I + a₁ * J) * x = rhs
# the optional shift is only used for the Hopf Newton / Continuation
function (l::GMRESKrylovKit{T, Tl})(J, rhs; a₀ = 0, a₁ = 1, kwargs...) where {T, Tl}
	if Tl == Nothing
		res, info = KrylovKit.linsolve(J, rhs, a₀, a₁; rtol = l.rtol, verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, atol = l.atol, issymmetric = l.issymmetric, ishermitian = l.ishermitian, isposdef = l.isposdef, kwargs...)
		info.converged == 0 && (@warn "KrylovKit.linsolve solver did not converge")
	else # use preconditioner
		res, info = KrylovKit.linsolve(x -> (out = apply(J, x); ldiv!(l.Pl, out)), ldiv!(l.Pl, copy(rhs)), a₀, a₁; rtol = l.rtol, verbosity = l.verbose, krylovdim = l.dim, maxiter = l.maxiter, atol = l.atol, issymmetric = l.issymmetric, ishermitian = l.ishermitian, isposdef = l.isposdef, kwargs...)
		info.converged == 0 && (@warn "KrylovKit.linsolve solver did not converge")
	end
	return res, true, info.numops
end


text:
function bisect_right(data::Array{Int64}, item::Int64)::Int64
    low = 0
    high = length(data)
    while low < high
        middle = i32::from(((low + high) / 2))
        if item < data[middle]
            high = middle
        else

            low = (middle + 1)
        end
    end
    return low
end

function bin_it(limits::Array{Int64}, data::Array{Int64})::Array{Int64}
    bins = [0]
    for x in limits
        bins.push(0)
    end
    for d in data
        bins[bisect_right(limits, d)] += 1
    end
    return bins
end

function main()
    limits = [23, 37, 43, 53, 67, 83]
    data = [
        95,
        21,
        94,
        12,
        99,
        4,
        70,
        75,
        83,
        93,
        52,
        80,
        57,
        5,
        53,
        86,
        65,
        17,
        92,
        83,
        71,
        61,
        54,
        58,
        47,
        16,
        8,
        9,
        32,
        84,
        7,
        87,
        46,
        19,
        30,
        37,
        96,
        6,
        98,
        40,
        79,
        97,
        45,
        64,
        60,
        29,
        49,
        36,
        43,
        55,
    ]
    @assert(bin_it(limits, data) == [11, 4, 2, 6, 9, 5, 13])
end


text:
# Type definition
struct Constant{T} <: FunctionSeries{T}
    coeffs::Base.RefValue{T}
    Constant(x::T) where T = new{T}(Ref(x))
end

# Basic properties
domain(f::Constant) = (-Inf,+Inf)
period(f::Constant) = 0

coefficients(f::Constant) = f.coeffs

Base.:(==)(f1::Constant, f2::Constant) = (f1[] == f2[])

# Indexing, getting and setting coefficients
Base.eachindex(f::Constant) = 0:0
Base.getindex(f::Constant) = f.coeffs[]
Base.getindex(f::Constant, i::Integer) = iszero(i) ? getindex(f) : throw(BoundsError(f, i))
Base.setindex!(f::Constant, v) = f.coeffs[] = v
Base.getindex(f::Constant, v, i::Integer) =
    iszero(i) ? setindex!(f, v) : throw(BoundsError(f, i))

# Use as function
(f::Constant)(x) = f[]

# Change number of coefficients
truncate!(f::Constant; kwargs...) = f

# special purpose constructor
Base.similar(f::Constant, ::Type{T} = scalartype(f)) where {T} =
    Constant(zero(T)*f[])

Base.zero(f::Constant) = Constant(zero(f[]))
Base.one(f::Constant) = Constant(one(f[]))

# Arithmetic (out of place)
Base.copy(f::Constant) = Constant(copy(f[]))

Base.:-(f::Constant) = Constant(-f[])
Base.:+(f::Constant) = Constant(+f[])

Base.:*(f::Constant, a::Const) = Constant(f[]*a)
Base.:*(a::Const, f::Constant) = Constant(a*f[])
Base.:/(f::Constant, a::Const) = Constant(f[]/a)
Base.:\(a::Const, f::Constant) = Constant(a\f[])

Base.:+(f1::Constant, f2::Constant) = Constant(f1[]+f2[])
Base.:-(f1::Constant, f2::Constant) = Constant(f1[]-f2[])
Base.:*(f1::Constant, f2::Constant) = Constant(f1[]*f2[])
Base.:/(f1::Constant, f2::Constant) = Constant(f1[]/f2[])
Base.:\(f1::Constant, f2::Constant) = Constant(f1[]\f2[])

truncmul(f1::Constant, f2::Constant) = Constant(f1[]*f2[])

Base.conj(f::Constant) = Constant(conj(f[]))
Base.adjoint(f::Constant) = Constant(adjoint(f[]))
Base.transpose(f::Constant) = Constant(transpose(f[]))

LinearAlgebra.tr(f::Constant) = Constant(tr(f[]))

Base.real(f::Constant) = Constant(real(f[]))
Base.imag(f::Constant) = Constant(imag(f[]))

# Arithmetic (in place / mutating methods)
function Base.copy!(fdst::Constant, fsrc::Constant)
    if eltype(fdst) <: Number
        fdst[] = fsrc[]
    else
        copy!(fdst[], fsrc[])
    end
    return fdst
end
function LinearAlgebra.rmul!(f::Constant, α::Number)
    if eltype(f) <: Number || !(α isa Number)
        f[] = f[]*α
    else
        rmul!(f[], α)
    end
    return f
end
function LinearAlgebra.lmul!(α::Number, f::Constant)
    if eltype(f) <: Number || !(α isa Number)
        f[] = α*f[]
    else
        lmul!(α, f[])
    end
    return f
end

LinearAlgebra.axpy!(α::Number, fx::Constant, fy::Constant) =
    truncadd!(fy, fx, α)
LinearAlgebra.axpby!(α::Number, fx::Constant, β::Number, fy::Constant) =
    truncadd!(fy, fx, α, β)
LinearAlgebra.mul!(fy::Constant, s::Number, fx::Constant, α = true, β = false) =
    truncmul!(fy, s, fx, α, β)
LinearAlgebra.mul!(fy::Constant, fx::Constant, s::Number, α = true, β = false) =
    truncmul!(fy, fx, s, α, β)
LinearAlgebra.mul!(f::Constant, f1::Constant, f2::Constant,
                    α = true, β = false) = truncmul!(f, f1, f2, α, β)

function truncadd!(fy::Constant, fx::Constant, α = true, β = true; kwargs...)
    if eltype(fy) <: Number
        fy[] = α*fx[] + β*fy[]
    else
        if iszero(β)
            mul!(fy[], α, fx[])
        elseif isone(β)
            axpy!(α, fx[], fy[])
        else
            axpby!(α, fx[], β, fy[])
        end
    end
    return fy
end

truncmul!(fdst::Constant, α₁::Number, fsrc::Constant, α₂ = true, β = false;
            kwargs...) = truncadd!(fdst, fsrc, α₁ * α₂, β; kwargs...)
truncmul!(fdst::Constant, fsrc::Constant, α₁::Number, α₂ = true, β = false;
            kwargs...) = truncadd!(fdst, fsrc, α₁ * α₂, β; kwargs...)

function truncmul!(f::Constant, f1::Constant, f2::Constant,
                    α = true, β = false;
                    Kmax::Integer = 0, tol::Real = 0)

    if eltype(f) <: Number
        f[] = α * f1[] + β * f2[]
    else
        if eltype(f1) <: Number
            axpby!(α*f1[], f2[], β, f[])
        elseif eltype(f2) <: Number
            axpby!(α*f2[], f1[], β, f[])
        else
            mul!(f[], f1[], f2[], α, β)
        end
    end
    return f
end

# Inner product and norm
localdot(f1::Constant, f2::Constant) = Constant(dot(f1[], f2[]))

LinearAlgebra.dot(f1::Constant, f2::Constant) = dot(f1[], f2[])
LinearAlgebra.norm(f::Constant) = norm(f[])

differentiate(f::Constant) = zero(f)
integrate(f::Constant, (a,b)::Tuple{Real,Real}) = f[]*(b-a)

# Fit constant: take average over N points
fit(f, ::Type{Constant}, (a,b)::Tuple{Real,Real}; numpoints = 5) =
    Constant(sum(f, range(a, b; length = numpoints))/numpoints)

# Inverse and square root
Base.inv(f::Constant) = Constant(inv(f[]))
Base.sqrt(f::Constant) = Constant(sqrt(f[]))


text:
# Original file is a part of CoolProp project http://www.coolprop.org/coolprop/wrappers/Julia/index.html#julia

module CoolProp

export F2K, K2F, HAPropsSI, PropsSI, PhaseSI, get_global_param_string, get_param_index, get_input_pair_index, AbstractState_factory, AbstractState_free, AbstractState_set_fractions, AbstractState_update, AbstractState_keyed_output

@windows_only const CoolPropLib = abspath(joinpath(@__FILE__,"..","..","lib","CoolProp.dll"));
@linux_only const CoolPropLib = abspath(joinpath(@__FILE__,"..","..","lib","libCoolProp.so"));

errcode = Ref{Clong}(0)

const buffer_length = 255
message_buffer = Array(UInt8, buffer_length)

function F2K(TF::Number)
  return ccall( (:F2K, CoolPropLib), Cdouble, (Cdouble,), TF)
end

function K2F(TK::Number)
  return ccall( (:K2F, CoolPropLib), Cdouble, (Cdouble,), TK)
end

function HAPropsSI(Output::AbstractString, Name1::AbstractString, Value1::Number, Name2::AbstractString, Value2::Number, Name3::AbstractString, Value3::Number)
  val = ccall( (:HAPropsSI, CoolPropLib), Cdouble, (Ptr{UInt8},Ptr{UInt8},Float64,Ptr{UInt8},Float64,Ptr{UInt8},Float64), Output,Name1,Value1,Name2,Value2,Name3,Value3)
  if val == Inf
    error("CoolProp: ", get_global_param_string("errstring"))
  end
  return val
end

function PropsSI(Output::AbstractString, Name1::AbstractString, Value1::Number, Name2::AbstractString, Value2::Number, Fluid::AbstractString)
  val = ccall( (:PropsSI, CoolPropLib), Cdouble, (Ptr{UInt8},Ptr{UInt8},Float64,Ptr{UInt8},Float64,Ptr{UInt8}), Output,Name1,Value1,Name2,Value2,Fluid)
  if val == Inf
    error("CoolProp: ", get_global_param_string("errstring"))
  end
  return val
end

function PropsSI(FluidName::AbstractString, Output::AbstractString)
  val = ccall( (:Props1SI, CoolPropLib), Cdouble, (Ptr{UInt8},Ptr{UInt8}), FluidName,Output)
  if val == Inf
    error("CoolProp: ", get_global_param_string("errstring"))
  end
  return val
end

function PhaseSI(Name1::AbstractString, Value1::Number, Name2::AbstractString, Value2::Number, Fluid::AbstractString)
  val = ccall( (:PhaseSI, CoolPropLib), Int32, (Ptr{UInt8},Float64,Ptr{UInt8},Float64,Ptr{UInt8}, Ptr{UInt8}, Int), Name1,Value1,Name2,Value2,Fluid,message_buffer::Array{UInt8,1},buffer_length)
  return bytestring(convert(Ptr{UInt8}, pointer(message_buffer::Array{UInt8,1})))
end

# This function returns the output string in pre-allocated char buffer.  If buffer is not large enough, no copy is made
function get_global_param_string(Key::AbstractString)
  val = ccall( (:get_global_param_string, CoolPropLib), Clong, (Ptr{UInt8},Ptr{UInt8},Int), Key,message_buffer::Array{UInt8,1},buffer_length)
  return bytestring(convert(Ptr{UInt8}, pointer(message_buffer::Array{UInt8,1})))
end

# Get the index for a parameter "T", "P", etc.
# returns the index as a long.
function get_param_index(Param::AbstractString)
  val = ccall( (:get_param_index, CoolPropLib), Clong, (Ptr{UInt8},), Param)
  if val == -1
    error("CoolProp: Unknown parameter: ", Param)
  end
  return val
end

# Get the index for an input pair for AbstractState.update function
# returns the index as a long.
function get_input_pair_index(Param::AbstractString)
  val = ccall( (:get_input_pair_index, CoolPropLib), Clong, (Ptr{UInt8},), Param)
  if val == -1
    error("CoolProp: Unknown input pair: ", Param)
  end
  return val
end

# ---------------------------------
#        Low-level access
# ---------------------------------

# Generate an AbstractState instance, return an integer handle to the state class generated to be used in the other low-level accessor functions
# param backend The backend you will use, "HEOS", "REFPROP", etc.
# param fluids '&' delimited list of fluids
# return A handle to the state class generated
function AbstractState_factory(backend::AbstractString, fluids::AbstractString)
  AbstractState = ccall( (:AbstractState_factory, CoolPropLib), Clong, (Ptr{UInt8},Ptr{UInt8},Ref{Clong},Ptr{UInt8},Clong), backend,fluids,errcode::RefValue{Clong,1},message_buffer::Array{UInt8,1},buffer_length)
  if errcode[] != 0
    if errcode[] == 1
      error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer))))
    elseif errcode[] == 2
      error("CoolProp: message buffer too small")
    else # == 3
      error("CoolProp: unknown error")
    end
  end
  return AbstractState
end

# Release a state class generated by the low-level interface wrapper
# param handle The integer handle for the state class stored in memory
function AbstractState_free(handle::Clong)
  ccall( (:AbstractState_free, CoolPropLib), Void, (Clong,Ref{Clong},Ptr{UInt8},Clong), handle,errcode::RefValue{Clong,1},message_buffer::Array{UInt8,1},buffer_length)
  if errcode[] != 0
    if errcode[] == 1
      error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer))))
    elseif errcode[] == 2
      error("CoolProp: message buffer too small")
    else # == 3
      error("CoolProp: unknown error")
    end
  end
  return nothing
end

# Set the fractions (mole, mass, volume) for the AbstractState
# param handle The integer handle for the state class stored in memory
# param fractions The array of fractions
function AbstractState_set_fractions(handle::Clong,fractions::Array)
  ccall( (:AbstractState_set_fractions, CoolPropLib), Void, (Clong,Ptr{Cdouble},Clong,Ref{Clong},Ptr{UInt8},Clong), handle,fractions,length(fractions),errcode::RefValue{Clong,1},message_buffer::Array{UInt8,1},buffer_length)
  if errcode[] != 0
    if errcode[] == 1
      error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer))))
    elseif errcode[] == 2
      error("CoolProp: message buffer too small")
    else # == 3
      error("CoolProp: unknown error")
    end
  end
  return nothing
end

# Update the state of the AbstractState
# param handle The integer handle for the state class stored in memory
# param input_pair The integer value for the input pair obtained from get_input_pair_index(Param::AbstractString)
# param value1 The first input value
# param value2 The second input value
function AbstractState_update(handle::Clong,input_pair::Clong,value1::Number,value2::Number)
  ccall( (:AbstractState_update, CoolPropLib), Void, (Clong,Clong,Cdouble,Cdouble,Ref{Clong},Ptr{UInt8},Clong), handle,input_pair,value1,value2,errcode::RefValue{Clong,1},message_buffer::Array{UInt8,1},buffer_length)
  if errcode[] != 0
    if errcode[] == 1
      error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer))))
    elseif errcode[] == 2
      error("CoolProp: message buffer too small")
    else # == 3
      error("CoolProp: unknown error")
    end
  end
  return nothing
end

# Get an output value from the AbstractState using an integer value for the desired output value
# param handle The integer handle for the state class stored in memory
# param param The integer value for the parameter you want
function AbstractState_keyed_output(handle::Clong, param::Clong)
  output = ccall( (:AbstractState_keyed_output, CoolPropLib), Cdouble, (Clong,Clong,Ref{Clong},Ptr{UInt8},Clong), handle,param,errcode::RefValue{Clong,1},message_buffer::Array{UInt8,1},buffer_length)
  if errcode[] != 0
    if errcode[] == 1
      error("CoolProp: ", bytestring(convert(Ptr{UInt8}, pointer(message_buffer))))
    elseif errcode[] == 2
      error("CoolProp: message buffer too small")
    else # == 3
      error("CoolProp: unknown error")
    end
  elseif output == -Inf
    error("CoolProp: no correct state has been set with AbstractState_update")
  end
  return output
end

end #module


text:
module LocallyNameless

using Reexport
@reexport using ..LambdaCalculus

import Base: show
import ..LambdaCalculus: boundvartype, freevartype, reify

export Index,
    Term,
    BVar,
    FVar,
    Abs,
    App

abstract type Term <: AbstractTerm end

struct BVar <: Term
    index::Index
    BVar(index) = index ≥ 1 ? new(index) : error("index must be at least 1")
end

struct FVar <: Term
    name::Symbol
end

struct Abs <: Term
    # boundname::Union{Symbol, Nothing}
    body::Term
end

struct App <: Term
    car::Term
    cdr::Term
end

show(io::IO, t::Abs) = print(io, "(λ.", t.body, ")")
show(io::IO, t::App) = print(io, "(", t.car, " ", t.cdr, ")")
show(io::IO, t::FVar) = print(io, "⌊", t.name, "⌋")
show(io::IO, t::BVar) = print(io, "⌈", t.index, "⌉")

reify(f::FVar) = :(FVar($(f.name)))
reify(b::BVar) = :(BVar($(b.index)))
reify(t::Abs) = :(Abs($(reify(t.body))))
reify(t::App) = :(App($(reify(t.car)), $(reify(t.cdr))))

boundvartype(::Type{<:Term}) = BVar
freevartype(::Type{<:Term}) = FVar

include("syntactic.jl")
include("meta.jl")

end # module LocallyNameless


text:
using Distributions, Random, Statistics, Plots; pyplot()
Random.seed!(2)

n, N, alpha = 10, 10^7, 0.05
mActual = 0.75
dist0, dist1 = Uniform(0,1), Uniform(0,mActual)

ts(sample) = maximum(sample) - minimum(sample)

empiricalDistUnderH0 = [ts(rand(dist0,n)) for _ in 1:N]
rejectionValue = quantile(empiricalDistUnderH0,alpha)

sample = rand(dist1,n)
testStat = ts(sample)
pValue = sum(empiricalDistUnderH0 .<= testStat)/N

if testStat > rejectionValue
    print("Didn't reject: ", round(testStat,digits=4))
    print(" > ", round(rejectionValue,digits=4))
else
    print("Reject: ", round(testStat,digits=4))
    print(" <= ", round(rejectionValue,digits=4))
end
println("\np-value = $(round(pValue,digits=4))")

stephist(empiricalDistUnderH0, bins=100, c=:blue, normed=true, label="")
plot!([testStat, testStat], [0,4], c=:red, label="Observed test statistic")
plot!([rejectionValue, rejectionValue], [0,4], c=:black, ls=:dash,
	label="Critical value boundary", legend=:topleft, ylims=(0,4),
    	xlabel = "x", ylabel = "Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
This function should works for Vector{Float}, Vector{DateTime}, etc...
"""
function resample_nearest(dtv::AbstractVector, r::StepRange)
    n = length(r)
    idx_vec = Vector{Int}(undef, n)
    i = 1

    while r[i] <= dtv[1]
        if i > n
            return idx_vec
        end
        idx_vec[i] = 1
        i += 1
    end
    for idx in 2:length(dtv)
        while true
            if i > n
                return idx_vec
            end
            if r[i] > dtv[idx]
                break
            end
            idx_vec[i] = r[i] - dtv[idx-1] <= dtv[idx] - r[i] ? idx - 1 : idx
            i += 1
        end
    end
    idx_vec[i:n] .= length(dtv)
    return idx_vec
end

abstract type ResampleMethod end

struct Nearest <: ResampleMethod end

function resample(::Nearest, ddf::DateDataFrame, r::StepRange{DateTime})
    idx_vec = resample_nearest(ddf.timestamp, r)
    return DateDataFrame(r, ddf.df[idx_vec, :])
end

_get_delta(step::Period, probe_coef) = Millisecond(step) * probe_coef
_get_delta(step::Period, probe_coef::Int) = step * probe_coef

"""
probe_coef is used to handle following format to encode [0, 1/24], [1/24, 2/24], ...:

nextfloat(0)
prevfloat(1/24)
nextfloat(1/24)
prevfloat(2/24)
...
"""
function resample(::Nearest, ddf::DateDataFrame, step::Period; probe_coef=0, drop=true)
    delta = _get_delta(step, probe_coef)

    bt = round(ddf.timestamp[1], typeof(step))
    et = round(ddf.timestamp[end], typeof(step))

    if drop
        while bt + delta < ddf.timestamp[1]
            bt += step
        end
        while et + delta > ddf.timestamp[end]
            et -= step
        end
    end

    r = bt:step:et
    test_r = (bt+delta):step:(et+delta)

    idx_vec = resample_nearest(ddf.timestamp, test_r)

    return DateDataFrame(r, ddf.df[idx_vec, :])
end


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

mutable struct PrivateEndpointConnectionProvisioningState <: SwaggerModel

    function PrivateEndpointConnectionProvisioningState(;)
        o = new()
        o
    end
end # type PrivateEndpointConnectionProvisioningState

const _property_map_PrivateEndpointConnectionProvisioningState = Dict{Symbol,Symbol}()
const _property_types_PrivateEndpointConnectionProvisioningState = Dict{Symbol,String}()
Base.propertynames(::Type{ PrivateEndpointConnectionProvisioningState }) = collect(keys(_property_map_PrivateEndpointConnectionProvisioningState))
Swagger.property_type(::Type{ PrivateEndpointConnectionProvisioningState }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PrivateEndpointConnectionProvisioningState[name]))}
Swagger.field_name(::Type{ PrivateEndpointConnectionProvisioningState }, property_name::Symbol) =  _property_map_PrivateEndpointConnectionProvisioningState[property_name]

function check_required(o::PrivateEndpointConnectionProvisioningState)
    true
end

function validate_property(::Type{ PrivateEndpointConnectionProvisioningState }, name::Symbol, val)
end


text:
using Pkg
Pkg.activate(".")
using LTMSim
using SimpleHypergraphs
using Random
using Plots
using DataFrames
using TableView

#metaV .= 1
#metaE .= 1
h, actV, actE, metaV, metaE = generateHypergraph();
bisect(h, metaV, metaE)

res = DataFrame(id=Int[],meanV=Float64[],tsssize=Int[])
for id in 1:1000
    h = randomH(8,6)
    metaV = randMetaV(h)
    metaE = proportionalMetaE(h,0.5)
    r = bisect(h,metaV,metaE)
    push!(res, Dict(
        :id => id,
        :meanV => mean(metaV),
        :tsssize => r,
    ))
end
using Plots
Plots.gr()

Plots.scatter(res.meanV .+ (rand(length(res.tsssize))*0.07) ,res.tsssize .+ (rand(length(res.tsssize)).*0.5), markersize=1)

using PyPlot
pygui(:qt)
plot(1:5,[1,4,1,5,1])
gcf()

using PyCall
using PyPlot
pygui(:qt)
hnx = pyimport("hypernetx")

function plotH(h, metaV, metaE, actV, actE)
    H = hnx.Hypergraph()

    edges = [ hnx.Entity("E$e [$(metaE[e])] $(actE[e] ? "X" : "_")",
        elements=["V$v [$(metaV[v])]$(actV[v] ? "X" : "_")"  for v in keys(h.he2v[e])])
        for e in 1:nhe(h)  ]
    H.add_edges_from(edges)
    clf()
    vs = ["V$v [$(metaV[v])]$(actV[v] ? "X" : "_")"  for v in 1:nhv(h)]
    Random.seed!(0)
    hnx.draw(H, pos=Dict(
        vs .=> [(3,0),(0,3),(4,3),(0,0),(3,1.5)] ) )
    gcf()
end

h, actV, actE, metaV, metaE = generateHypergraph();
plotH(h,metaV,metaE, actV, actE)
simulate!(h,actV, actE, metaV, metaE;max_step=1)
plotH(h,metaV,metaE, actV, actE)

bisect(h, metaV, metaE)

#H.incidence_dict

res = DataFrame(id=Int[],meanV=Float64[],propMetaV=Float64[],tsssize=Int[])
for id in 1:1000
    h = randomH(30,15)
    for propMetaV in 0.1:0.1:0.9
        metaV = proportionalMetaV(h,propMetaV)
        metaE = proportionalMetaE(h,0.5)
        r = bisect(h,metaV,metaE)
        push!(res, Dict(
            :id => id,
            #:meanV => mean(something.(h,0)),
            :meanV => std(length.(h.he2v)),
            :propMetaV => propMetaV,
            :tsssize => r,
        ))
    end
end
using Plots
Plots.gr()

Plots.scatter(res.meanV .+ (rand(length(res.tsssize))*0.07) ,res.tsssize .+ (rand(length(res.tsssize)).*0.5), markersize=1)

unique(res.tsssize)

matplotlib_cm = pyimport("matplotlib.cm")
matplotlib_colors = pyimport("matplotlib.colors")
cmap = matplotlib_cm.get_cmap("cool");

res.meanV

getColor(x) = RGB(cmap(x)[1:3]...)

using TableView
TableView.showtable(res)

maxmv = maximum(res.meanV)
minmv = minimum(res.meanV)

Plots.scatter(
    res.propMetaV.+ rand(nrow(res)).*0.06,
    res.tsssize.+ rand(nrow(res)).*0.4,
    markersize=1,
    color=getColor.((res.meanV.-minmv)./(maxmv-minmv) ))

using Statistics
cor(Matrix(res[:,[:meanV,:propMetaV,:tsssize]]))

data_to_plot = Vector{Array{Any, 1}}()

for propMetaV in unique(res, :propMetaV).propMetaV
    series = filter(x->x.propMetaV == propMetaV, res).tsssize

    push!(data_to_plot, series)
end

clf()
fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(111)
ax.boxplot(data_to_plot)

xlabels = unique(res, :propMetaV).propMetaV
xticks(1:length(xlabels), xlabels, rotation=0)

ylabel("TSS size")
xlabel("Vertices thresholds")

gcf()

plt.tight_layout(.5)