text:
import LinearAlgebra.BLAS: @blasfunc, libblas
import LinearAlgebra.LAPACK: liblapack
# Level 2
## mv
### gemv

gemv!(x...) = BLAS.gemv!(x...)

for (fname, elty) in ((:dgemv_,:Float64),
                      (:sgemv_,:Float32),
                      (:zgemv_,:ComplexF64),
                      (:cgemv_,:ComplexF32))
    @eval begin
             #SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             #*     .. Scalar Arguments ..
             #      DOUBLE PRECISION ALPHA,BETA
             #      INTEGER INCX,INCY,LDA,M,N
             #      CHARACTER TRANS
             #*     .. Array Arguments ..
             #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function gemv!(trans::Char, m::Int, n::Int, alpha::($elty), A::Ptr{$elty}, stA::Int,
                        X::Ptr{$elty}, incX::Int, beta::($elty), Y::Ptr{$elty}, incY::Int)
            ccall((@blasfunc($fname), libblas), Nothing,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 trans, m, n, alpha,
                 A, stA, X, incX,
                 beta, Y, incY)
            Y
        end

        function gemv!(trans::Char, alpha::($elty), A, X, beta::($elty), Y)
            m,n = size(A,1),size(A,2)
            if trans == 'N' && (length(X) != n || length(Y) != m)
                throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'C' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("A' has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'T' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("transpose(A) has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            end
            gemv!(trans, m, n, alpha, pointer(A), max(1,stride(A,2)),
                    pointer(X), stride(X,1), beta, pointer(Y), stride(Y,1))
            Y
        end
    end
end

## Level 3

# We add some missing Ptr variants

gemm!(x...) = BLAS.gemm!(x...)

for (gemm, elty) in
        ((:dgemm_,:Float64),
         (:sgemm_,:Float32),
         (:zgemm_,:ComplexF64),
         (:cgemm_,:ComplexF32))
    @eval begin
             # SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER K,LDA,LDB,LDC,M,N
             #       CHARACTER TRANSA,TRANSB
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
        function gemm!(transA::Char, transB::Char, m::Int, n::Int, k::Int, alpha::($elty),
                        A::Ptr{$elty}, stA::Int, B::Ptr{$elty}, stB::Int,
                        beta::($elty), C::Ptr{$elty},stC::Int)
            ccall((@blasfunc($gemm), libblas), Nothing,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}),
                 transA, transB, m, n,
                 k, alpha, A, stA,
                 B, stB, beta, C,
                 stC)
            C
        end

        function gemm!(transA::Char, transB::Char, alpha::($elty), A, B, beta::($elty), C)
#           if any([stride(A,1), stride(B,1), stride(C,1)] .!= 1)
#               error("gemm!: BLAS module requires contiguous matrix columns")
#           end  # should this be checked on every call?
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            gemm!(transA, transB, m, n, ka, alpha,
                    pointer(A), max(1,stride(A,2)),
                    pointer(B), max(1,stride(B,2)), beta,
                    pointer(C), max(1,stride(C,2)))
            C
        end
    end
end


text:
## Loading data into various data structures ##

struct Require
    versions::VersionSet
    systems::Vector{Symbol}
end

struct Version
    sha1::Union{Nothing, String}
    requires::Dict{String,Require}
end

struct Package
    uuid::UUID
    url::String
    versions::Dict{VersionNumber,Version}
end

Require(versions::VersionSet) = Require(versions, Symbol[])
Require(version::VersionInterval) = Require(VersionSet([version]), Symbol[])
Version(sha1::String) = Version(sha1, Dict{String,Require}())

function load_requires(path::String)
    requires = Dict{String,Require}()
    requires["julia"] = Require(VersionInterval())
    isfile(path) || return requires
    for r in filter!(r->r isa Requirement, Reqs.read(path))
        new = haskey(requires, r.package)
        versions, systems = VersionSet(r.versions.intervals), r.system
        if haskey(requires, r.package)
            versions = versions ∩ requires[r.package].versions
            systems  = systems  ∪ requires[r.package].systems
        end
        requires[r.package] = Require(versions, Symbol.(systems))
    end
    return requires
end

function load_versions(dir::String)
    versions = Dict{VersionNumber,Version}()
    isdir(dir) || return versions
    for ver in readdir(dir)
        path = joinpath(dir, ver)
        sha1 = joinpath(path, "sha1")
        isfile(sha1) || continue
        requires = load_requires(joinpath(path, "requires"))
        versions[VersionNumber(ver)] = Version(readchomp(sha1), requires)
    end
    return versions
end

function load_packages_METADATA(
    dir::AbstractString,
    selected_packages::AbstractVector{<:AbstractString}
)
    pkgs = Dict{String,Package}()
    for pkg in selected_packages
        path = joinpath(dir, pkg)
        url = joinpath(path, "url")
        versions = joinpath(path, "versions")
        isfile(url) || continue
        pkgs[pkg] = Package(
            uuid5(uuid_package, pkg),
            readchomp(url),
            load_versions(versions),
        )
    end
    pkgs
end

@eval julia_versions() = $([VersionNumber(0,m) for m=1:7])
julia_versions(f::Function) = filter(f, julia_versions())
julia_versions(vi::VersionInterval) = julia_versions(v->v in vi)

macro clean(ex) :(x = $(esc(ex)); $(esc(:clean)) &= x; x) end

function check_packages(
    pkgs::AbstractDict{<:AbstractString,Package},
    other_registries::AbstractDict{<:AbstractString, Package}
)
    for (pkg, p) in pkgs, (ver, v) in p.versions
        if ver != thispatch(ver)
            error("$pkg version $ver is illegeal: no prereleases")
        end
        if ver == v"0.0.0"
            error("$pkg version $ver: Version must be greater than 0.0.0")
        end
        for (req, r) in v.requires
            # req == "julia" && continue
            if haskey(other_registries, req)
                req_registry = other_registries
            elseif haskey(pkgs, req)
                req_registry = pkgs
            else
                error("$pkg ver $ver: requirement $req not found")
            end
            available_req_versions = keys(req_registry[req].versions)
            if ! any(w->w in r.versions, available_req_versions)
                error(
"""
$pkg ver $ver: compat range not satisfied
$pkg at $ver requires $req version in $(r.versions)
But known versions for $req are $available_req_versions
"""
                )
            end
        end
    end
    nothing
end

ext: jl


text:
struct AsymptoticLinearTest{K<:MatrixKernel,V,S} <: CalibrationTest
    """Matrix-valued kernel."""
    kernel::K
    """Calibration error estimate."""
    estimate::V
    """Standard deviation of asymptotic normal distribution."""
    std::S
end

AsymptoticLinearTest(skce::LinearUnbiasedSKCE, data) =
    AsymptoticLinearTest(skce.kernel, data)

function AsymptoticLinearTest(kernel::MatrixKernel,
                              data::Tuple{<:AbstractMatrix{<:Real},<:AbstractVector{<:Integer}})
    # check if the number of predictions and labels is equal
    predictions, labels = CalibrationErrors.get_predictions_labels(data)

    estimate, std = _calibration_std(kernel, predictions, labels)

    AsymptoticLinearTest(kernel, estimate, std)
end

HypothesisTests.testname(::AsymptoticLinearTest) =
    "Asymptotic calibration test based on the linear unbiased SKCE estimator"

HypothesisTests.pvalue(test::AsymptoticLinearTest) = normccdf(0, test.std, test.estimate)

"""
    _calibration_std(kernel, predictions, labels)

Evaluate the linear unbiased SKCE estimator with kernel `kernel` and estimate the standard
devation of its asymptotic normal distribution around the true value for the given
`predictions` and `labels`.
"""
function _calibration_std(kernel::MatrixKernel, predictions::AbstractMatrix{<:Real},
                          labels::AbstractVector{<:Integer})
    # precalculations
    nsamples = size(predictions, 2)

    # initialize estimate and sum of squares
    m = zero(skce_result_type(LinearUnbiasedSKCE(kernel), predictions))
    S = m^2

    n, i = 0, 1
    @inbounds while i < nsamples
        # update number of summands
        n += 1

        # evaluate kernel of next two samples
        j = i + 1
        s = skce_kernel(kernel, view(predictions, :, j), labels[j], view(predictions, :, i),
                        labels[i])

        # update mean
        Δm = s - m
        m += inv(n) * Δm

        # update sum of squares
        S += Δm * (s - m)

        # update index of next sample
        i += 2
    end

    m, inv(n) * sqrt(S)
end


text:
# https://arxiv.org/abs/1511.05286
struct SegmentedLSE{T, U} <: AggregationFunction
    ρ::T
    C::U
end

Flux.@functor SegmentedLSE

SegmentedLSE(d::Int) = SegmentedLSE(randn(Float32, d), zeros(Float32, d))

r_map(ρ) = softplus.(ρ)
inv_r_map = (r) -> max.(r, 0) .+ log1p.(-exp.(-abs.(r)))

(m::SegmentedLSE)(x::MaybeMatrix, bags::AbstractBags, w=nothing) = segmented_lse_forw(x, m.C, r_map(m.ρ), bags)
function (m::SegmentedLSE)(x::AbstractMatrix, bags::AbstractBags, w::AggregationWeights, mask::AbstractVector)
    segmented_lse_forw(x .+ typemin(T) * mask', m.C, r_map(m.ρ), bags)
end

function _lse_precomp(x::AbstractMatrix, r, bags)
    M = zeros(eltype(x), length(r), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if !isempty(b)
            for i in eachindex(r)
                M[i, bi] = r[i] * x[i, first(b)]
            end
            for j in b[2:end]
                for i in eachindex(r)
                    M[i, bi] = max(M[i, bi], r[i] * x[i, j])
                end
            end
        end
    end
    M
end

function _segmented_lse_norm(x::AbstractMatrix, C, r, bags::AbstractBags, M)
    y = zeros(eltype(x), length(r), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                y[i, bi] = C[i]
            end
        else
            for j in b
                for i in eachindex(r)
                    y[i, bi] += exp.(r[i] * x[i, j] - M[i, bi])
                end
            end
            for i in eachindex(r)
                y[i, bi] = (log(y[i, bi]) - log(length(b)) + M[i, bi]) / r[i]
            end
        end
    end
    y
end

segmented_lse_forw(::Missing, C::AbstractVector, r, bags::AbstractBags) = repeat(C, 1, length(bags))
function segmented_lse_forw(x::AbstractMatrix, C, r, bags::AbstractBags)
    M = _lse_precomp(x, r, bags)
    _segmented_lse_norm(x, C, r, bags, M)
end

function segmented_lse_back(Δ, y, x, C, r, bags, M)
    dx = zero(x)
    dC = zero(C)
    dr = zero(r)
    s1 = zeros(eltype(x), length(r))
    s2 = zeros(eltype(x), length(r))

    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                dC[i] += Δ[i, bi]
            end
        else
            for i in eachindex(r)
                s1[i] = s2[i] = 0
            end
            for j in b
                for i in eachindex(r)
                    e = exp(r[i] * x[i, j] - M[i, bi])
                    s1[i] += e
                    s2[i] += x[i, j] * e
                end
            end
            for j in b
                for i in eachindex(r)
                    dx[i, j] = Δ[i, bi] * exp(r[i] * x[i, j] - M[i, bi]) / s1[i]
                end
            end
            for i in eachindex(r)
                dr[i] += Δ[i, bi] * (s2[i]/s1[i] - y[i, bi]) / r[i]
            end
        end
    end
    dx, dC, dr, nothing, nothing
end

function segmented_lse_back(Δ, ::Missing, C, bags)
    dC = zero(C)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(C)
            dC[i] += Δ[i, bi]
        end
    end
    nothing, dC, nothing, nothing, nothing
end

@adjoint function segmented_lse_forw(x::AbstractMatrix, C::AbstractVector, r::AbstractVector, bags::AbstractBags)
    M = _lse_precomp(x, r, bags)
    y = _segmented_lse_norm(x, C, r, bags, M)
    grad = Δ -> segmented_lse_back(Δ, y, x, C, r, bags, M)
    y, grad
end

@adjoint function segmented_lse_forw(x::Missing, C::AbstractVector, r::AbstractVector, bags::AbstractBags)
    y = segmented_lse_forw(x, C, r, bags)
    grad = Δ -> segmented_lse_back(Δ, x, C, bags)
    y, grad
end


text:
#

using Base.Test
using Compat
using BinDeps

Pkg.build("Cairo")  # Tests apt-get code paths
using Cairo
Pkg.build("HttpParser")  # Tests build-from-source code paths
using HttpParser
Pkg.build("GSL")
using GSL

# PR 171
@test BinDeps.lower(nothing, nothing) === nothing


text:
"""
Reverse plus
"""
function plus_rev(a::Interval, b::Interval, c::Interval)  # a = b + c
    # a = a ∩ (b + c)  # add this line for plus contractor (as opposed to reverse function)
    b_new = b ∩ (a - c)
    c_new = c ∩ (a - b)

    return a, b_new, c_new
end

plus_rev(a,b,c) = plus_rev(promote(a,b,c)...)

"""
Reverse minus
"""
function minus_rev(a::Interval, b::Interval, c::Interval)  # a = b - c

    b_new = b ∩ (a + c)
    c_new = c ∩ (b - a)

    return a, b_new, c_new
end

minus_rev(a,b,c) = minus_rev(promote(a,b,c)...)

function minus_rev(a::Interval, b::Interval)  # a = -b
    b_new = b ∩ (-a)
    return (a, b_new)
end

"""
Reverse multiplication
"""
function mul_rev(a::Interval, b::Interval, c::Interval)  # a = b * c

    # ((0.0 ∉ a) || (0.0 ∉ b)) && (c = c ∩ (a / b))
    # ((0.0 ∉ a) || (0.0 ∉ c)) && (b = b ∩ (a / c))

    # a = a ∩ (b * c)  # ?

    if 0 ∈ b
        temp = c .∩ extended_div(a, b)
        c′ = union(temp[1], temp[2])

    else
        c′ = c ∩ (a / b)
    end

    if 0 ∈ c
        temp = b .∩ extended_div(a, c)
        b′ = union(temp[1], temp[2])

    else
        b′ = b ∩ (a / c)
    end

    return a, b′, c′
end

mul_rev(a,b,c) = mul_rev(promote(a,b,c)...)

"""
Reverse division
"""
function div_rev(a::Interval, b::Interval, c::Interval)  # a = b / c

    b = b ∩ (a * c)
    c = c ∩ (b / a)

    return a, b, c
end

div_rev(a,b,c) = div_rev(promote(a,b,c)...)

"""
Reverse inverse
"""
function inv_rev(a::Interval, b::Interval)  # a = inv(b)

    b_new = b ∩ inv(a)

    return a, b_new
end

inv_rev(a,b) = inv_rev(promote(a,b)...)

"""
Reverse power
"""
function power_rev(a::Interval, b::Interval, n::Integer)  # a = b^n,  log(a) = n.log(b),  b = a^(1/n)

    if n == 2  # a = b^2
        root = √a
        b1 = b ∩ root
        b2 = b ∩ (-root)

    elseif iseven(n)
        root = a^(1//n)

        b1 = b ∩ root
        b2 = b ∩ (-root)

    elseif isodd(n)
        pos_root = (a ∩ (0..∞)) ^ (1//n)
        neg_root = -( ( (-a) ∩ (0..∞) ) ^ (1//n) )

        b1 = b ∩ pos_root
        b2 = b ∩ neg_root

    end

    b = hull(b1, b2)

    return (a, b, n)
end

function power_rev(a::Interval, b::Interval, c::Interval)  # a = b^c

    if isinteger(c)
        temp = power_rev(a, b, Int(inf(c)))  # use version with integer
        return (temp[1], temp[2], interval(temp[3]))
    end

    b_new = b ∩ ( a^(inv(c) ))
    c_new = c ∩ (log(a) / log(b))

    return a, b_new, c_new
end

power_rev(a, b, c) = power_rev(promote(a, b, c)...)

"""
Reverse square root
"""
function sqrt_rev(a::Interval, b::Interval)  # a = sqrt(b)

    b_new = b ∩ (a^2)

    return a, b_new
end

sqrt_rev(a,b) = sqrt_rev(promote(a,b)...)

# IEEE-1788 style

"""
Reverse sqr
"""
function sqr_rev(c, x)   # c = x^2;  refine x

    root = sqrt(c)

    x1 = x ∩ root
    x2 = x ∩ (-root)

    return (c, hull(x1, x2))
end

sqr_rev(c) = sqr_rev(c, -∞..∞)

"""
Reverse abs
"""
function abs_rev(y, x)   # y = abs(x); refine x

    y_new = y ∩ (0..∞)

    x1 = y_new ∩ x
    x2 = -(y_new ∩ (-x))

    return (y, hull(x1, x2))
end
#=
"""
Reverse sign
"""
function sign_rev(a::Interval, b::Interval)  # a = sqrt(b)

    (a == 1.0) && b = b ∩ (0..∞)
    (a == 0.0) && b = b ∩ (0.0..0.0)
    (a == -1.0) && b = b ∩ (-∞..0.0)

    return a, b
end
sign_rev(a,b) = sign_rev(promote(a,b)...)
=#
## IEEE-1788 versions:

"""
According to the IEEE-1788 standard:

- `∘_rev1(b, c, x)` is the subset of `x` such that `x ∘ b` is defined and in `c`;
- `∘_rev2(a, c, x)` is the subset of `x` such that `a ∘ x` is defined and in `c`

When `∘` is commutative, these agree and we write `∘_rev(b, c, x)`.
"""

function mul_rev_IEEE1788(b, c, x)   # c = b*x
    return x ∩ (c / b)
end

function pow_rev1(b, c, x)   # c = x^b
    return x ∩ c^(1/b)  # replace by 1//b
end

function pow_rev2(a, c, x)   # c = a^x
    return x ∩ (log(c) / lob(a))
end


text:
module FactorioJSON

export recipes

using JSON

const recipeLUA = """

  {
    typ : "recipe",
    name : "player-port",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 10},
      {"iron-gear-wheel", 5},
      {"iron-plate", 1 }
    },
    result : "player-port"
  },
  {
    typ : "recipe",
    name : "fast-transport-belt",
    enabled : false,
    ingredients :
    {
      {"iron-gear-wheel", 5},
      {"transport-belt", 1}
    },
    result : "fast-transport-belt",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "express-transport-belt",
    category : "crafting-with-fluid",
    enabled : false,
    ingredients :
    {
      {"iron-gear-wheel", 5},
      {"fast-transport-belt", 1},
      {type:"fluid", name:"lubricant", amount:2},
    },
    result : "express-transport-belt",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "solar-panel",
    energy_required : 10,
    enabled : false,
    ingredients :
    {
      {"steel-plate", 5},
      {"electronic-circuit", 15},
      {"copper-plate", 5}
    },
    result : "solar-panel"
  },
  {
    typ : "recipe",
    name : "assembling-machine-2",
    enabled : false,
    ingredients :
    {
      {"iron-plate", 9},
      {"electronic-circuit", 3},
      {"iron-gear-wheel", 5},
      {"assembling-machine-1", 1}
    },
    result : "assembling-machine-2",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "assembling-machine-3",
    enabled : false,
    ingredients :
    {
      {"speed-module", 4},
      {"assembling-machine-2", 2}
    },
    result : "assembling-machine-3"
  },
  {
    typ : "recipe",
    name : "car",
    enabled : false,
    ingredients :
    {
      {"engine-unit", 8},
      {"iron-plate", 20},
      {"steel-plate", 5}
    },
    result : "car"
  },
  {
    typ : "recipe",
    name : "tank",
    enabled : false,
    ingredients :
    {
      {"engine-unit", 16},
      {"steel-plate", 50},
      {"iron-gear-wheel", 15},
      {"advanced-circuit", 5}
    },
    result : "tank"
  },
  {
    typ : "recipe",
    name : "rail",
    enabled : false,
    ingredients :
    {
      {"stone", 1},
      {"iron-stick", 1},
      {"steel-plate", 1}
    },
    result : "rail",
    result_count : 2,
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "diesel-locomotive",
    enabled : false,
    ingredients :
    {
      {"engine-unit", 20},
      {"electronic-circuit", 10},
      {"steel-plate", 30}
    },
    result : "diesel-locomotive"
  },
  {
    typ : "recipe",
    name : "cargo-wagon",
    enabled : false,
    ingredients :
    {
      {"iron-gear-wheel", 10},
      {"iron-plate", 20},
      {"steel-plate", 20}
    },
    result : "cargo-wagon"
  },
  {
    typ : "recipe",
    name : "train-stop",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 5},
      {"iron-plate", 10},
      {"steel-plate", 3}
    },
    result : "train-stop"
  },
  {
    typ : "recipe",
    name : "rail-signal",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 1},
      {"iron-plate", 5}
    },
    result : "rail-signal",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "rail-chain-signal",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 1},
      {"iron-plate", 5}
    },
    result : "rail-chain-signal",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "heavy-armor",
    enabled : false,
    energy_required : 8,
    ingredients : {{ "copper-plate", 100}, {"steel-plate", 50}},
    result : "heavy-armor"
  },
  {
    typ : "recipe",
    name : "modular-armor",
    enabled : false,
    energy_required : 15,
    ingredients : {{ "advanced-circuit", 30}, {"processing-unit", 5}, {"steel-plate", 50}},
    result : "modular-armor"
  },
  {
    typ : "recipe",
    name : "power-armor",
    enabled : false,
    energy_required : 20,
    ingredients : {{ "processing-unit", 40}, {"electric-engine-unit", 20}, {"steel-plate", 40}, {"alien-artifact", 10}},
    result : "power-armor",
    requester_paste_multiplier : 1
  },
  {
    typ : "recipe",
    name : "power-armor-mk2",
    enabled : false,
    energy_required : 25,
    ingredients : {{ "effectivity-module-3", 5}, {"speed-module-3", 5}, {"processing-unit", 40}, {"steel-plate", 40}, {"alien-artifact", 50}},
    result : "power-armor-mk2",
    requester_paste_multiplier : 1
  },
  {
    typ : "recipe",
    name : "iron-chest",
    enabled : true,
    ingredients : {{"iron-plate", 8}},
    result : "iron-chest",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "steel-chest",
    enabled : false,
    ingredients : {{"steel-plate", 8}},
    result : "steel-chest",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "stone-wall",
    enabled : false,
    ingredients : {{"stone-brick", 5}},
    result : "stone-wall",
    requester_paste_multiplier : 10
  },
  {
    typ : "recipe",
    name : "gate",
    enabled : false,
    ingredients : {{"stone-wall", 1}, {"steel-plate", 2}, {"electronic-circuit", 2}},
    result : "gate"
  },
  {
    typ : "recipe",
    name : "flame-thrower",
    enabled : false,
    energy_required : 10,
    ingredients :
    {
      {"steel-plate", 5},
      {"iron-gear-wheel", 10}
    },
    result : "flame-thrower"
  },
  {
    typ : "recipe",
    name : "land-mine",
    enabled : false,
    energy_required : 5,
    ingredients :
    {
      {"steel-plate", 1},
      {"explosives", 2}
    },
    result : "land-mine",
    result_count : 4
  },
  {
    typ : "recipe",
    name : "rocket-launcher",
    enabled : false,
    energy_required : 5,
    ingredients :
    {
      {"iron-plate", 5},
      {"iron-gear-wheel", 5},
      {"electronic-circuit", 5}
    },
    result : "rocket-launcher"
  },
  {
    typ : "recipe",
    name : "shotgun",
    enabled : false,
    energy_required : 4,
    ingredients :
    {
      {"iron-plate", 15},
      {"iron-gear-wheel", 5},
      {"copper-plate", 10},
      {"wood", 5}
    },
    result : "shotgun"
  },
  {
    typ : "recipe",
    name : "combat-shotgun",
    enabled : false,
    energy_required : 8,
    ingredients :
    {
      {"steel-plate", 15},
      {"iron-gear-wheel", 5},
      {"copper-plate", 10},
      {"wood", 10}
    },
    result : "combat-shotgun"
  },
  {
    typ : "recipe",
    name : "railgun",
    enabled : false,
    energy_required : 8,
    ingredients :
    {
      {"steel-plate", 15},
      {"copper-plate", 15},
      {"electronic-circuit", 10},
      {"advanced-circuit", 5},
    },
    result : "railgun"
  },
  {
    typ : "recipe",
    name : "science-pack-1",
    energy_required : 5,
    ingredients :
    {
      {"copper-plate", 1},
      {"iron-gear-wheel", 1}
    },
    result : "science-pack-1"
  },
  {
    typ : "recipe",
    name : "science-pack-2",
    energy_required : 6,
    ingredients :
    {
      {"inserter", 1},
      {"transport-belt", 1}
    },
    result : "science-pack-2"
  },
  {
    typ : "recipe",
    name : "science-pack-3",
    enabled : false,
    energy_required : 12,
    ingredients :
    {
      {"battery", 1},
      {"advanced-circuit", 1},
      {"filter-inserter", 1},
      {"steel-plate", 1},
    },
    result : "science-pack-3"
  },
  {
    typ : "recipe",
    name : "alien-science-pack",
    enabled : false,
    energy_required : 12,
    ingredients : {{"alien-artifact", 1}
    },
    result : "alien-science-pack",
    result_count : 10
  },
  {
    typ : "recipe",
    name : "lab",
    energy_required : 5,
    ingredients :
    {
      {"electronic-circuit", 10},
      {"iron-gear-wheel", 10},
      {"transport-belt", 4}
    },
    result : "lab"
  },
  {
    typ : "recipe",
    name : "red-wire",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 1},
      {"copper-cable", 1}
    },
    result : "red-wire",
    requester_paste_multiplier : 10
  },
  {
    typ : "recipe",
    name : "green-wire",
    enabled : false,
    ingredients :
    {
      {"electronic-circuit", 1},
      {"copper-cable", 1}
    },
    result : "green-wire",
    requester_paste_multiplier : 10
  },
  {
    typ : "recipe",
    name : "underground-belt",
    enabled : false,
    energy_required : 1,
    ingredients :
    {
      {"iron-plate", 10},
      {"transport-belt", 5}
    },
    result_count : 2,
    result : "underground-belt",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "fast-underground-belt",
    enabled : false,
    ingredients :
    {
      {"iron-gear-wheel", 20},
      {"underground-belt", 2}
    },
    result_count : 2,
    result : "fast-underground-belt",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "express-underground-belt",
    category : "crafting-with-fluid",
    enabled : false,
    ingredients :
    {
      {"iron-gear-wheel", 40},
      {"fast-underground-belt", 2},
      {type:"fluid", name:"lubricant", amount:4},
    },
    result_count : 2,
    result : "express-underground-belt"
  },
  {
    typ : "recipe",
    name : "loader",
    enabled : false,
    energy_required : 1,
    ingredients :
    {
      {"inserter", 5},
      {"electronic-circuit", 5},
      {"iron-gear-wheel", 5},
      {"iron-plate", 5},
      {"transport-belt", 5}
    },
    result : "loader"
  },
  {
    typ : "recipe",
    name : "fast-loader",
    enabled : false,
    energy_required : 3,
    ingredients :
    {
      {"fast-transport-belt", 5},
      {"loader", 1}
    },
    result : "fast-loader"
  },
  {
    typ : "recipe",
    name : "express-loader",
    enabled : false,
    energy_required : 10,
    ingredients :
    {
      {"express-transport-belt", 5},
      {"fast-loader", 1}
    },
    result : "express-loader"
  },
  {
    typ : "recipe",
    name : "splitter",
    enabled : false,
    energy_required : 1,
    ingredients :
    {
      {"electronic-circuit", 5},
      {"iron-plate", 5},
      {"transport-belt", 4}
    },
    result : "splitter",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "fast-splitter",
    enabled : false,
    energy_required : 2,
    ingredients :
    {
      {"splitter", 1},
      {"iron-gear-wheel", 10},
      {"electronic-circuit", 10}
    },
    result : "fast-splitter",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "express-splitter",
    category : "crafting-with-fluid",
    enabled : false,
    energy_required : 2,
    ingredients :
    {
      {"fast-splitter", 1},
      {"iron-gear-wheel", 10},
      {"advanced-circuit", 10},
      {type:"fluid", name:"lubricant", amount:8}
    },
    result : "express-splitter"
  },
  {
    typ : "recipe",
    name : "advanced-circuit",
    enabled : false,
    energy_required : 8,
    ingredients :
    {
      {"electronic-circuit", 2},
      {"plastic-bar", 2},
      {"copper-cable", 4}
    },
    result : "advanced-circuit",
    requester_paste_multiplier : 5
  },
  {
    typ : "recipe",
    name : "processing-unit",
    category : "crafting-with-fluid",
    enabled : false,
    energy_required : 15,
    ingredients :
    {
      {"electronic-circuit", 20},
      {"advanced-circuit", 2},
      {type:"fluid", name : "sulfuric-acid", amount : 0.5}
    },
    result : "processing-unit"
  },
  {
    typ : "recipe",
    name : "logistic-robot",
    enabled : false,
    ingredients :
    {
      {"flying-robot-frame", 1},
      {"advanced-circuit", 2}
    },
    result : "logistic-robot"
  },
  {
    typ : "recipe",
    name : "construction-robot",
    enabled : false,
    ingredients :
    {
      {"flying-robot-frame", 1},
      {"electronic-circuit", 2}
    },
    result : "construction-robot"
  },
  {
    typ : "recipe",
    name : "logistic-chest-passive-provider",
    enabled : false,
    ingredients :
    {
      {"steel-chest", 1},
      {"electronic-circuit", 3},
      {"advanced-circuit", 1}
    },
    result : "logistic-chest-passive-provider",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "logistic-chest-active-provider",
    enabled : false,
    ingredients :
    {
      {"steel-chest", 1},
      {"electronic-circuit", 3},
      {"advanced-circuit", 1}
    },
    result : "logistic-chest-active-provider",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "logistic-chest-storage",
    enabled : false,
    ingredients :
    {
      {"steel-chest", 1},
      {"electronic-circuit", 3},
      {"advanced-circuit", 1}
    },
    result : "logistic-chest-storage",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "logistic-chest-requester",
    enabled : false,
    ingredients :
    {
      {"steel-chest", 1},
      {"electronic-circuit", 3},
      {"advanced-circuit", 1}
    },
    result : "logistic-chest-requester",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "rocket-silo",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 1000},
      {"concrete", 1000},
      {"pipe", 100},
      {"processing-unit", 200},
      {"electric-engine-unit", 200}
    },
    energy_required : 30,
    result : "rocket-silo",
    requester_paste_multiplier : 1
  },
  {
    typ : "recipe",
    name : "roboport",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 45},
      {"iron-gear-wheel", 45},
      {"advanced-circuit", 45}
    },
    result : "roboport",
    energy_required : 15
  },
  {
    typ : "recipe",
    name : "steel-axe",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 5},
      {"iron-stick", 2}
    },
    result : "steel-axe",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "big-electric-pole",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 5},
      {"copper-plate", 5}
    },
    result : "big-electric-pole",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "substation",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 10},
      {"advanced-circuit", 5},
      {"copper-plate", 5}
    },
    result : "substation"
  },
  {
    typ : "recipe",
    name : "medium-electric-pole",
    enabled : false,
    ingredients :
    {
      {"steel-plate", 2},
      {"copper-plate", 2}
    },
    result : "medium-electric-pole",
    requester_paste_multiplier : 4
  },
  {
    typ : "recipe",
    name : "accumulator",
    energy_required : 10,
    enabled : false,
    ingredients :
    {
      {"iron-plate", 2},
      {"battery", 5}
    },
    result : "accumulator"
  },
  {
    typ : "recipe",
    name : "steel-furnace",
    ingredients : {{"steel-plate", 8}, {"stone-brick", 10}},
    result : "steel-furnace",
    energy_required : 3,
    enabled : false
  },
  {
    typ : "recipe",
    name : "electric-furnace",
    ingredients : {{"steel-plate", 15}, {"advanced-circuit", 5}, {"stone-brick", 10}},
    result : "electric-furnace",
    energy_required : 5,
    enabled : false
  },
  {
    typ : "recipe",
    name : "beacon",
    enabled : false,
    energy_required : 15,
    ingredients :
    {
      {"electronic-circuit", 20},
      {"advanced-circuit", 20},
      {"steel-plate", 10},
      {"copper-cable", 10}
    },
    result : "beacon"
  },
  {
    typ : "recipe",
    name : "blueprint",
    energy_required : 1,
    ingredients :
    {
      {"advanced-circuit", 1}
    },
    result : "blueprint",
    enabled : false
  },
  {
    typ : "recipe",
    name : "blueprint-book",
    energy_required : 5,
    ingredients :
    {
      {"advanced-circuit", 15}
    },
    result : "blueprint-book",
    enabled : false
  },
  {
    typ : "recipe",
    name : "deconstruction-planner",
    energy_required : 1,
    ingredients :
    {
      {"advanced-circuit", 1}
    },
    result : "deconstruction-planner",
    enabled : false
  },
  {
    typ : "recipe",
    name : "pumpjack",
    energy_required : 20,
    ingredients :
    {
      {"steel-plate", 15},
      {"iron-gear-wheel", 10},
      {"electronic-circuit", 10},
      {"pipe", 10},
    },
    result : "pumpjack",
    enabled : false
  },
  {
    typ : "recipe",
    name : "oil-refinery",
    energy_required : 20,
    ingredients :
    {
      {"steel-plate", 15},
      {"iron-gear-wheel", 10},
      {"stone-brick", 10},
      {"electronic-circuit", 10},
      {"pipe", 10}
    },
    result : "oil-refinery",
    enabled : false
  },
  {
    typ : "recipe",
    name : "engine-unit",
    energy_required : 20,
    category : "advanced-crafting",
    ingredients :
    {
      {"steel-plate", 1},
      {"iron-gear-wheel", 1},
      {"pipe", 2}
    },
    result : "engine-unit",
    enabled : false
  },
  {
    typ : "recipe",
    name : "electric-engine-unit",
    category : "crafting-with-fluid",
    energy_required : 20,
    ingredients :
    {
      {"engine-unit", 1},
      {type:"fluid", name: "lubricant", amount : 2},
      {"electronic-circuit", 2}
    },
    result : "electric-engine-unit",
    enabled : false
  },
  {
    typ : "recipe",
    name : "flying-robot-frame",
    energy_required : 20,
    ingredients :
    {
      {"electric-engine-unit", 1},
      {"battery", 2},
      {"steel-plate", 1},
      {"electronic-circuit", 3}
    },
    result : "flying-robot-frame",
    enabled : false
  },
  {
    typ : "recipe",
    name : "explosives",
    energy_required : 5,
    enabled : false,
    category : "chemistry",
    ingredients :
    {
      {type:"item", name:"sulfur", amount:1},
      {type:"item", name:"coal", amount:1},
      {type:"fluid", name:"water", amount:1},
    },
    result: "explosives"
  },
  {
    typ : "recipe",
    name : "battery",
    category : "chemistry",
    energy_required : 5,
    enabled : false,
    ingredients :
    {
      {type:"fluid", name:"sulfuric-acid", amount:2},
      {"iron-plate", 1},
      {"copper-plate", 1}
    },
    result: "battery"
  },
  {
    typ : "recipe",
    name : "storage-tank",
    energy_required : 3,
    enabled : false,
    ingredients :
    {
      {"iron-plate", 20},
      {"steel-plate", 5}
    },
    result: "storage-tank"
  },
  {
    typ : "recipe",
    name : "small-pump",
    energy_required : 2,
    enabled : false,
    ingredients :
    {
      {"electric-engine-unit", 1},
      {"steel-plate", 1},
      {"pipe", 1}
    },
    result: "small-pump"
  },
  {
    typ : "recipe",
    name : "chemical-plant",
    energy_required : 10,
    enabled : false,
    ingredients :
    {
      {"steel-plate", 5},
      {"iron-gear-wheel", 5},
      {"electronic-circuit", 5},
      {"pipe", 5}
    },
    result: "chemical-plant"
  },
  {
    typ : "recipe",
    name : "small-plane",
    energy_required : 30,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"plastic-bar", 100},
      {"advanced-circuit", 200},
      {"electric-engine-unit", 20},
      {"battery", 100}
    },
    result: "small-plane"
  },
  {
    typ : "recipe",
    name : "arithmetic-combinator",
    enabled : "false",
    ingredients :
    {
      {"copper-cable", 5},
      {"electronic-circuit", 5},
    },
    result : "arithmetic-combinator"
  },
  {
    typ : "recipe",
    name : "decider-combinator",
    enabled : "false",
    ingredients :
    {
      {"copper-cable", 5},
      {"electronic-circuit", 5},
    },
    result : "decider-combinator"
  },
  {
    typ : "recipe",
    name : "constant-combinator",
    enabled : "false",
    ingredients :
    {
      {"copper-cable", 5},
      {"electronic-circuit", 2},
    },
    result : "constant-combinator"
  },
  {
    typ : "recipe",
    name : "power-switch",
    enabled : "false",
    energy_required : 2,
    ingredients :
    {
      {"iron-plate", 5},
      {"copper-cable", 5},
      {"electronic-circuit", 2},
    },
    result : "power-switch"
  },
  {
    typ : "recipe",
    name : "low-density-structure",
    energy_required : 30,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"steel-plate", 10},
      {"copper-plate", 5},
      {"plastic-bar", 5}
    },
    result: "low-density-structure"
  },
  {
    typ : "recipe",
    name : "rocket-fuel",
    energy_required : 30,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"solid-fuel", 10}
    },
    result: "rocket-fuel"
  },
  {
    typ : "recipe",
    name : "rocket-control-unit",
    energy_required : 30,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"processing-unit", 1},
      {"speed-module", 1}
    },
    result: "rocket-control-unit"
  },
  {
    typ : "recipe",
    name : "rocket-part",
    energy_required : 3,
    enabled : false,
    hidden : true,
    category : "rocket-building",
    ingredients :
    {
      {"low-density-structure", 10},
      {"rocket-fuel", 10},
      {"rocket-control-unit", 10}
    },
    result: "rocket-part"
  },
  {
    typ : "recipe",
    name : "satellite",
    energy_required : 3,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"low-density-structure", 100},
      {"solar-panel", 100},
      {"accumulator", 100},
      {"radar", 5},
      {"processing-unit", 100},
      {"rocket-fuel", 50}
    },
    result: "satellite"
  },
  {
    typ : "recipe",
    name : "concrete",
    energy_required : 10,
    enabled : false,
    category : "crafting-with-fluid",
    ingredients :
    {
      {"stone-brick", 5},
      {"iron-ore", 1},
      {type:"fluid", name:"water", amount:10}
    },
    result: "concrete",
    result_count : 10
  },
  {
    typ : "recipe",
    name : "hazard-concrete",
    energy_required : 0.25,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"concrete", 10}
    },
    result: "hazard-concrete",
    result_count : 10
  },
  {
    typ : "recipe",
    name : "landfill",
    energy_required : 0.5,
    enabled : false,
    category : "crafting",
    ingredients :
    {
      {"stone", 20}
    },
    result: "landfill",
    result_count : 1
  },
  {
    typ : "recipe",
    name : "electric-energy-interface",
    energy_required : 0.5,
    enabled : false,
    ingredients :
    {
      {"iron-plate", 2},
      {"electronic-circuit", 5}
    },
    result : "electric-energy-interface"
  }

"""

function readLuaf(fn)
	#JSON.parsefile(filename::AbstractString; dicttype=Dict, use_mmap=true)

	print(JSON.parse(recipeLUA))
	close(fid)

end

function recipes()
	print(recipeLUA)
	readLuaf("recipe.txt")
end

# staaahhppp
end

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
using Alea, FactCheck

facts("GPC") do
  germ = GPCGerm(Normal(0,1), 4)
  #=
  @show rand(germ)
  @show pdf(germ, [0.3, 0.4, 0.5, 0.6])
  @show cdf(germ, [0.3, 0.4, 0.5, 0.6])
  @show cdf(germ, [0.0, 0.0, 0.0, 0.0])
  @show cdf(germ, [Inf, 0.0, 0.0, Inf])
  =#
  context("cdf") do
    @fact cdf(germ, zeros(4)) --> 0.5^4
    @fact cdf(germ, [0.0, 0.0, 0.0, Inf]) --> 0.5 ^ 3
    @fact cdf(germ, [0.0, -Inf, 0.0, Inf]) --> 0
  end
end


text:
## Code for PCG ALGORITHM #######################################################

mutable struct PCG
    state::UInt64
    inc::UInt64

    function PCG(init_state::UInt64 = UInt64(42), init_seq::UInt64 = UInt64(54))
        pcg = new(UInt64(0), UInt64((init_seq << 1 | 1)))
        pcg_rand(pcg)
        pcg.state += init_state
        pcg_rand(pcg)
        pcg
    end
end

function pcg_rand(pcg::PCG)
    oldstate = pcg.state
    pcg.state = UInt64(oldstate * 6364136223846793005 + pcg.inc)
    xorshifted = UInt32(((oldstate >> UInt64(18)) ⊻ oldstate) >> UInt64(27) & typemax(UInt32))
    rot = oldstate >> UInt64(59)
    return UInt32((xorshifted >> rot) | (xorshifted << ((-rot) & UInt32(31))))
end

pcg_randf(pcg) = Float32(pcg_rand(pcg)/typemax(UInt32))


text:
include("SymbolicCone.jl")

mutable struct Cone
    S::SymbolicCone
    sign::Int
    Cone(S::SymbolicCone, sign::Int) =
        new(S, sign)
end

function PrintCone(C::Cone)
    PrintSymbolicCone(C.S)
    println("sign: ", C.sign)
end

max_stars_repo_name: kevseryolcu/PolyhedralOmega

max_issues_repo_name: kevseryolcu/PolyhedralOmega

max_forks_repo_name: kevseryolcu/PolyhedralOmega


text:
part1(input) = repeated_look_and_say(input, 40)
part2(input) = repeated_look_and_say(input, 50)

function repeated_look_and_say(input, n)
    x = parse.(Int, collect(readchomp(input)))
    for i = 1:n
        x = look_and_say(x)
    end
    return length(x)
end

function look_and_say(x)
    y = Int[]
    last_d = first(x)
    run_length = 0
    for d in x
        if d == last_d
            run_length += 1
        else
            if last_d > 0
                push!(y, run_length, last_d)
            end
            run_length = 1
            last_d = d
        end
    end
    push!(y, run_length, last_d)
end


text:
@testset "PLL discriminator" begin

    correlator_minus60off = EarlyPromptLateCorrelator(
        -0.5 + sqrt(3) / 2im,
        -1 + sqrt(3) * 1im,
        -0.5 + sqrt(3) / 2im
    )
    correlator_0off = EarlyPromptLateCorrelator(
        0.5 + 0.0im,
        1 + 0.0im,
        0.5 + 0.0im
    )
    correlator_plus60off = EarlyPromptLateCorrelator(
        0.5 + sqrt(3) / 2im,
        1 + sqrt(3) * 1im,
        0.5 + sqrt(3) / 2im
    )

    @test @inferred(Tracking.pll_disc(GPSL1, correlator_minus60off)) == -π / 3  #-60°
    @test @inferred(Tracking.pll_disc(GPSL1, correlator_0off)) == 0
    @test @inferred(Tracking.pll_disc(GPSL1, correlator_plus60off)) == π / 3  #+60°
end

@testset "DLL discriminator" begin

    very_early_correlator = EarlyPromptLateCorrelator(0.0 + 0.0im, 0.5 + 0.0im, 1.0 + 0.0im)
    early_correlator = EarlyPromptLateCorrelator(0.25 + 0.0im, 0.75 + 0.0im, 0.75 + 0.0im)
    prompt_correlator = EarlyPromptLateCorrelator(0.5 + 0.0im, 1 + 0.0im, 0.5 + 0.0im)
    late_correlator = EarlyPromptLateCorrelator(0.75 + 0.0im, 0.75 + 0.0im, 0.25 + 0.0im)
    very_late_correlator = EarlyPromptLateCorrelator(1.0 + 0.0im, 0.5 + 0.0im, 0.0 + 0.0im)

    sample_shift = 2
    delta = 0.25

    @test @inferred(
        Tracking.dll_disc(GPSL1, very_early_correlator, sample_shift, delta)
    ) == -0.5
    @test @inferred(
        Tracking.dll_disc(GPSL1, early_correlator, sample_shift, delta)
    ) == -0.25
    @test @inferred(
        Tracking.dll_disc(GPSL1, prompt_correlator, sample_shift, delta)
    ) == 0
    @test @inferred(
        Tracking.dll_disc(GPSL1, late_correlator, sample_shift, delta)
    ) == 0.25
    @test @inferred(
        Tracking.dll_disc(GPSL1, very_late_correlator, sample_shift, delta)
    ) == 0.5
end


text:
using CommonSubexpressions

ex1a = :( x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[2]*x[5] )
cse(ex1a)

ex1b = :( x1*x2*x3 + x1*x2*x4 + x1*x2*x5 )
cse(ex1b)

display(@macroexpand(@cse x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[2]*x[5]))
display(@macroexpand(@cse x1*x2*x3 + x1*x2*x4 + x1*x2*x5))

ex = :(
   x1*x2*x3*x4*x5 +
   x1*x2*x3*x4*x6 +
   x1*x2*x3*x4*x8 +
   x1*x2*x5*x6*x7 +
   x1*x3*x5*x6*x7 +
   x1*x5*x6*x7*x8 +
   x1*x2*x3*x4*x7 +
   x1*x4*x5*x6*x7 +
   x1*x2*x5*x8*x9 +
   x2*x3*x5*x8*x9 +
   x2*x5*x6*x8*x9 +
   x1*x2*x3*x4*x9 +
   x2*x4*x5*x8*x9 +
   x1*x5*x6*x7*x9 +
   x2*x5*x7*x8*x9 +
   x1*x3*x6*x8*x10 +
   x2*x3*x6*x8*x10 +
   x3*x5*x6*x8*x10 +
   x1*x2*x3*x4*x10 +
   x3*x4*x6*x8*x10 +
   x1*x5*x6*x7*x10 +
   x3*x6*x7*x8*x10 +
   x2*x5*x8*x9*x10 +
   x3*x6*x8*x9*x10 +
   x1*x4*x7*x9*x10 +
   x2*x4*x7*x9*x10 +
   x4*x5*x7*x9*x10 +
   x3*x4*x7*x9*x10 +
   x4*x6*x7*x9*x10 +
   x4*x7*x8*x9*x10
)

const ex_cse = cse(ex)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using DtCalc
using Test

@testset "DtCalc.jl" begin
    # Write your tests here.
end


text:
using StructJuMP

numScen = 2
m = StructuredModel(numScen)

@defVar(m, 0 <= x <= 1)
@defVar(m, 0 <= y <= 1)

@addConstraint(m, x + y == 1)
setObjective(m, :Min, x*x + y)

for i in 1:numScen
    bl = StructuredModel(parent=m)
    @defVar(bl, w >= 0)
    @addConstraint(bl, w - x - y <= 1)
    setObjective(bl, :Min, w*w + w)
end


text: include("karlgaard.jl")

max_stars_repo_name: glemieux/RelativeEstimationProject

max_issues_repo_name: glemieux/RelativeEstimationProject

max_forks_repo_name: glemieux/RelativeEstimationProject


text:
"""
    convert_JSON(file; <keyword arguments>)
    convert_JSON(dict; <keyword arguments>)

Convert a JSON `file` or a `dict` conforming to the FlexPlan WP3 API into a FlexPlan.jl dict.

Costs are scaled with the assumption that every representative year represents 10 years.

# Arguments
- `number_of_hours::Union{Int,Nothing}=nothing`: parse only the first hours of the
  file/dict.
- `number_of_scenarios::Union{Int,Nothing}=nothing`: parse only the first scenarios of the
  file/dict.
- `number_of_years::Union{Int,Nothing}=nothing`: parse only the first years of the
  file/dict.
- `cost_scale_factor::Float64=1.0`: scale factor for all costs.
- `init_data_extensions::Vector{<:Function}=Function[]`: functions to be applied to the
  target dict after its initialization. They must have exactly one argument (the target
  dict) and can modify it; the return value is unused.
- `sn_data_extensions::Vector{<:Function}=Function[]`: functions to be applied to the
  single-network dictionaries containing data for each single year, just before
  `_FP.make_multinetwork` is called. They must have exactly one argument (the single-network
  dict) and can modify it; the return value is unused.

# Extended help
Features of FlexPlan WP3 API not supported in FlexPlan.jl:
- scenario probabilities depending on year (only constant probabilities are allowed);
- number of hours depending on scenario (all scenarios must have the same number of hours);
- combined transmission and distribution networks;
- PSTs;
- curtailment of renewable generators;
- `gridModelInputFile.converters.ratedActivePowerDC`;
- `gridModelInputFile.storage.minEnergy`;
- `gridModelInputFile.storage.maxAbsRamp`;
- `gridModelInputFile.storage.maxInjRamp`;
- uniqueness of candidate components (each candidate can be reinvested at the end of its
  lifetime).
"""
function convert_JSON end

function convert_JSON(file::String; kwargs...)
    source_dict = JSON.parsefile(file)
    convert_JSON(source_dict; kwargs...)
end

function convert_JSON(source::AbstractDict;
        number_of_hours::Union{Int,Nothing} = nothing,
        number_of_scenarios::Union{Int,Nothing} = nothing,
        number_of_years::Union{Int,Nothing} = nothing,
        cost_scale_factor::Float64 = 1.0,
        init_data_extensions::Vector{<:Function} = Function[],
        sn_data_extensions::Vector{<:Function} = Function[],
    )

    # Define target dict

    target = Dict{String, Any}(
        "nw"           => Dict{String,Any}(),
        "multinetwork" => true,
        "per_unit"     => true,
    )

    # Add dimensions

    if length(unique(vcat(source["genericParameters"]["nbHours"]...))) > 1
        Memento.error(_LOGGER, "All scenarios must have the same number of hours.") # Dimensions are implemented as a multidimensional array, so the length of one dimension cannot depend on the id along another dimension.
    end
    if isnothing(number_of_hours)
        number_of_hours = first(first(source["genericParameters"]["nbHours"]))
    elseif number_of_hours > first(first(source["genericParameters"]["nbHours"]))
        Memento.error(_LOGGER, "$number_of_hours hours requested, but only " * string(first(first(source["genericParameters"]["nbHours"]))) * " found in input dict.")
    end
    _FP.add_dimension!(target, :hour, number_of_hours)

    if isnothing(number_of_scenarios)
        number_of_scenarios = source["genericParameters"]["nbScenarios"]
    elseif number_of_scenarios > source["genericParameters"]["nbScenarios"]
        Memento.error(_LOGGER, "$number_of_scenarios scenarios requested, but only " * string(source["genericParameters"]["nbScenarios"]) * " found in input dict.")
    end
    if haskey(source["genericParameters"], "scenarioProbabilities")
        if maximum(length.(unique.(source["genericParameters"]["scenarioProbabilities"]))) > 1
            Memento.warn(_LOGGER, "Only constant probabilities are supported for scenearios. Using first year probabilities for every year.")
        end
        scenario_probabilities = first(first.(source["genericParameters"]["scenarioProbabilities"]), number_of_scenarios) # The outermost `first` is needed if the user has specified a number of scenarios lower than that available.
    else
        scenario_probabilities = fill(1/number_of_scenarios,number_of_scenarios)
    end
    scenario_properties = Dict(id => Dict{String,Any}("probability"=>prob) for (id,prob) in enumerate(scenario_probabilities))
    _FP.add_dimension!(target, :scenario, scenario_properties)

    if isnothing(number_of_years)
        number_of_years = length(source["genericParameters"]["years"])
    elseif number_of_years > length(source["genericParameters"]["years"])
        Memento.error(_LOGGER, "$number_of_years years requested, but only " * string(length(source["genericParameters"]["years"])) * " found in input dict.")
    end
    year_scale_factor = 10 # Assumption: every representative year represents 10 years
    _FP.add_dimension!(target, :year, number_of_years; metadata = Dict{String,Any}("scale_factor"=>year_scale_factor))

    # Generate ID lookup dict

    lookup_acBranches = id_lookup(source["gridModelInputFile"]["acBranches"])
    lookup = Dict(
        "acBuses"           => id_lookup(source["gridModelInputFile"]["acBuses"]),
        "acBranches"        => lookup_acBranches,
        "transformers"      => id_lookup(source["gridModelInputFile"]["transformers"]; offset=length(lookup_acBranches)), # AC branches are split between `acBranches` and `transformers` dicts in JSON files
        "generators"        => id_lookup(source["gridModelInputFile"]["generators"]),
        "loads"             => id_lookup(source["gridModelInputFile"]["loads"]),
        "storage"           => id_lookup(source["gridModelInputFile"]["storage"]),
        "dcBuses"           => id_lookup(source["gridModelInputFile"]["dcBuses"]),
        "dcBranches"        => id_lookup(source["gridModelInputFile"]["dcBranches"]),
        "converters"        => id_lookup(source["gridModelInputFile"]["converters"]),
    )
    if haskey(source, "candidatesInputFile")
        lookup_cand_acBranches = id_lookup(source["candidatesInputFile"]["acBranches"], "acBranch")
        lookup["cand_acBranches"]   = lookup_cand_acBranches
        lookup["cand_transformers"] = id_lookup(source["candidatesInputFile"]["transformers"], "acBranch"; offset=length(lookup_cand_acBranches)) # AC branches are split between `acBranches` and `transformers` dicts in JSON files
        lookup["cand_storage"]      = id_lookup(source["candidatesInputFile"]["storage"], "storage")
        lookup["cand_dcBranches"]   = id_lookup(source["candidatesInputFile"]["dcBranches"], "dcBranch")
        lookup["cand_converters"]   = id_lookup(source["candidatesInputFile"]["converters"], "converter")
    end

    # Compute availability of candidates

    if haskey(source, "candidatesInputFile")
        year_scale_factor = _FP.dim_meta(target, :year, "scale_factor")
        year_lookup = Dict{Int,Int}((year,y) for (y,year) in enumerate(source["genericParameters"]["years"]))
        cand_availability = Dict{String,Any}(
            "acBranches"   => availability(source, "acBranches",    "acBranch",  year_lookup, year_scale_factor, number_of_years),
            "transformers" => availability(source, "transformers",  "acBranch",  year_lookup, year_scale_factor, number_of_years),
            "loads"        => availability(source, "flexibleLoads", "load",      year_lookup, year_scale_factor, number_of_years),
            "storage"      => availability(source, "storage",       "storage",   year_lookup, year_scale_factor, number_of_years),
            "dcBranches"   => availability(source, "dcBranches",    "dcBranch",  year_lookup, year_scale_factor, number_of_years),
            "converters"   => availability(source, "converters",    "converter", year_lookup, year_scale_factor, number_of_years),
        )
    end

    # Apply init data extensions

    for f! in init_data_extensions
        f!(target)
    end

    # Build data year by year

    for y in 1:number_of_years
        sn_data = haskey(source, "candidatesInputFile") ? nw(source, lookup, cand_availability, y) : nw(source, lookup, y)
        _FP.scale_data!(sn_data; year_idx=y, number_of_years, year_scale_factor, cost_scale_factor)
        sn_data["dim"] = target["dim"]

        # Apply single network data extensions
        for f! in sn_data_extensions
            f!(sn_data)
        end

        time_series = make_time_series(source, lookup, y, sn_data; number_of_hours, number_of_scenarios)
        year_data = _FP.make_multinetwork(sn_data, time_series; number_of_nws=number_of_hours*number_of_scenarios, nw_id_offset=number_of_hours*number_of_scenarios*(y-1))
        add_singular_data!(year_data, source, lookup, y)
        _FP.import_nws!(target, year_data)
    end

    return target
end

# Define a bijective map from existing JSON String ids to FlexPlan Int ids (generated in _id_lookup)
function id_lookup(component_vector::Vector; offset::Int=0)
    json_ids = [d["id"] for d in component_vector]
    return _id_lookup(json_ids, offset)
end

function id_lookup(component_vector::Vector, sub_key::String; offset::Int=0)
    json_ids = [d[sub_key]["id"] for d in component_vector]
    return _id_lookup(json_ids, offset)
end

function _id_lookup(json_ids::Vector, offset::Int)
    sort!(json_ids) # Sorting prevents changes due to the order of elements in JSON files
    int_ids = range(1+offset; length=length(json_ids))
    lookup = Dict{String,Int}(zip(json_ids, int_ids))
    if length(lookup) < length(json_ids)
        Memento.error(_LOGGER, "IDs must be unique (found only $(length(lookup)) unique IDs, should be $(length(json_ids))).")
    end
    return lookup
end

function availability(source::AbstractDict, comp_name::String, sub_key::String, year_lookup::AbstractDict, year_scale_factor::Int, number_of_years::Int)
    target = Dict{String,Vector{Bool}}()
    for comp in source["candidatesInputFile"][comp_name]
        id = comp[sub_key]["id"]
        investment_horizon = [year_lookup[year] for year in comp["horizons"]]
        if last(investment_horizon) - first(investment_horizon) ≥ length(investment_horizon)
            Memento.warn(_LOGGER, "Horizon of $comp_name $id is not a contiguous set.")
        end
        raw_lifetime = comp["lifetime"] ÷ year_scale_factor
        availability_horizon_start = first(investment_horizon)
        availability_horizon_end = min(last(investment_horizon)+raw_lifetime-1, number_of_years)
        target[id] = [availability_horizon_start ≤ y ≤ availability_horizon_end for y in 1:number_of_years]
    end
    return target
end


text:
using AccuracyAtTopDual
using DatasetProvider
using EvalMetrics
using Plots

# data preparation
T = Float32;
dataset = Dataset(
    DatasetProvider.Spambase;
    asmatrix = true,
    shuffle = true,
    binarize = true,
    poslabels = 1,
    seed = 123,
);
train, test = load(TrainTest(0.8), dataset);
Xtrain, ytrain = Array{T}(train[1]), train[2];
Xtest, ytest = Array{T}(test[1]), test[2];

# training
model = TopPush(; C = 1, ϑ = 1, surrogate = Quadratic, T)
kernel = KernelType(Gaussian; γ = 1, scale = true, precomputed = false, T)

h = solve!(model, Xtrain, ytrain, kernel; maxiter = 20000, seed = 123, ε = 1e-4)

# evaluation
s_train = h.solution[:train].s;
s_test = predict(model, Xtrain, ytrain, kernel, Xtest);

plot(
    rocplot(ytrain, s_train; title = "Train"),
    rocplot(ytest, s_test; title = "Test");
    xlims = (1e-5, 1),
    xscale = :log10,
    fill = false,
    size = (700, 400),
)


text:
module PIV3D

using FFTW, FileIO, LIBTIFF, Statistics

export PIV

# Shortened type names to keep function signatures short ( when possible one-liners )
include("units.jl")

F      = Float64;
I      = Int64;
A      = Array;
S      = String;
C      = Complex;
II     = Tuple{Int64,Int64};
III    = Tuple{Int64,Int64,Int64};
UnI    = Union{I,II,III}
SySySy = Tuple{Symbol,Symbol,Symbol}
DU     = DistanceUnits
TU     = TimeUnits
UNTS   = Tuple{DU,DU,DU,TU}

struct ZNCC end
struct FFT  end

include("parameters.jl")
include("pivIO.jl")
include("utils.jl")

include("PIV/gaussianSubpixel.jl")
include("PIV/signalToNoise.jl")
include("PIV/ccrZNCC.jl")
include("PIV/ccrFFT.jl")

include("postprocessing.jl")

function PIV( data1::Array{T,N}, data2::Array{T,N}, params::PIVParameters; checkParams=true, units=nothing ) where {T<:Real,N}

    checkParams && checkPIVParameters( size(data1), N, params )

    #Continuing to PIV analysis
    _PIV( data1, data2, params )
end

""" Calling 2D PIV analysis """
function _PIV( img1::A{T,2}, img2::A{T,2}, p::PIVParameters; units=nothing ) where {T<:Real}

    u, v, sn = PIV_2D( p.corr, img1, img2, p.args... );

    if units !== nothing
        ratios = changeUnits.(  p.units[1:2], units[1:2] )./changeUnits( p.units, units[4]  )
        for e in 1:length(u)
            u[e] = u[e]*ratios[1]
            v[e] = v[e]*ratios[2]
        end
    end
    return u, v, sn;
end

""" Calling 3D PIV analysis """
function _PIV( vol1::A{T,3}, vol2::A{T,3}, p::PIVParameters; units=nothing ) where {T<:Real}

    u, v, w, sn = PIV_3D( p.corr, vol1, vol2, p.args... );

    if units !== nothing
        ratios = changeUnits.( p.units, units[1:3] )./changeUnits( p.units, units[4] )
        for e in 1:length(u)
            u[e] = u[e]*ratios[1]
            v[e] = v[e]*ratios[2]
            w[e] = w[e]*ratios[3]
        end
    end

    return u, v, w, sn;
end

"""
    Automated PIV analysis between each pair of consecutive timepoints from a time-sequence data.

    This function assumes that the files follow a consistent naming scheme, where the name of each
    file contains a varying index indicating its time index. Aside from the varying index, the rest
    of the file name should be the same between all files. For example:

        result_001_monday, result_002_monday, result_003_monday, etc..

    To run a sequential PIV analysis on the previous example the following arguments should be give:

        sequencePIV( pp, "result_", "_monday", 1, 3, digits=3 )

    where pp is a PIV parameter object, "result_" and "_monday" are the invariant parts in each file name
    before and after the varying index, 1 is the first index, 3 is the last index to analyze, and digits=3
    indicates the number of digits of the varying index. 0's are prepended to reach "digits".

"""
function sequencePIV( PIVparams::PIVParameters,
                      before_id::String, after_id::String,
                      id0::Integer, idStop::Integer;
                      digits::Integer=1,
                      path::String=pwd(),
                      subsample::Int64=1,
                      typ=Float32 )

    n = idStop - id0

    U = Array{ Array{Float32,3}, 1 }( undef, n );
    V = Array{ Array{Float32,3}, 1 }( undef, n );
    W = Array{ Array{Float32,3}, 1 }( undef, n );

    if path[end] !== '/'
        path *= '/'
    end

    for idx in id0:idStop-1

        idx1  = numberDigits(  idx , maxdigits=digits );
        idx2  = numberDigits( idx+1, maxdigits=digits );
        file1 = before_id * idx1 * after_id;
        file2 = before_id * idx2 * after_id;
        vol1  = PIVload( path*file1, sample=subsample, typ=typ );
        vol2  = PIVload( path*file2, sample=subsample, typ=typ );

        println( "Performing PIV between:" )
        println( "\t$(path*file1)" )
        println( "\t$(path*file2)" )

        u, v, w = PIV( vol1, vol2, PIVparams );

        U[1+idx-id0] = u;
        V[1+idx-id0] = v;
        W[1+idx-id0] = w;
    end

    return U, V, W
end

# projectionStackPIVOFF in smbnotes
"""
    PIV3D.projectionStackPIV( PIVParams, filename )

    Computes 2D PIV between stack[idx] and stack[idx+1], for all consecutive
    projections in the stack.
"""
function projectionStackPIV( PIVParams, fn::String; threshold=-1, lo=1, hi=-1, reps=Int64[], typ=nothing )

    stack = PIVload( fn, typ=typ );
    size( stack, 3 ) < 2 && error("Stack only contains one image, no PIV possible")

    lo = ( lo < size( stack, 3 ) ) ? lo :  1;
    hi = ( hi == -1 ) ? size(stack,3)-1 : hi;

    i1 = stack[ :, :,  lo  ];
    i2 = stack[ :, :, lo+1 ];

    u, v, _  = PIV( i1, i2, PIVParams );

    U = zeros( size(u,1), size(u,2), hi - lo + 2 );
    V = zeros( size(v,1), size(v,2), hi - lo + 2 );

    U[:,:,1] .= u;
    V[:,:,1] .= v;

    for idx in lo+1:hi
        println("Analyzing $(idx-lo)/$(hi-lo)");

        i1 = stack[ :, :,  idx  ];
        i2 = stack[ :, :, idx+1 ];

        u, v, _ = PIV( i1, i2 , PIVParams );

        U[:,:,idx-lo+1] .= u;
        V[:,:,idx-lo+1] .= v;
    end

    return U, V
end

function projectionStackPIV( PIVParams, stack::Array{T,3};
                             lo=1, hi=-1, reps=Int64[] ) where {T<:Real}

    lo = ( lo < size( stack, 3 ) ) ? lo :  1;
    hi = ( hi == -1 ) ? size(stack,3)-1 : hi;

    i1 = stack[ :, :,  lo  ];
    i2 = stack[ :, :, lo+1 ];

    u, v, _  = PIV( i1, i2, PIVParams );

    U = zeros( size(u,1), size(u,2), hi - lo + 2 );
    V = zeros( size(v,1), size(v,2), hi - lo + 2 );

    U[:,:,1] .= u;
    V[:,:,1] .= v;

    for idx in lo+1:hi
        println("Analyzing $(idx-lo)/$(hi-lo)");

        i1 = stack[ :, :,  idx  ];
        i2 = stack[ :, :, idx+1 ];

        u, v, _ = PIV( i1, i2 , PIVParams );

        U[:,:,idx-lo+1] .= u;
        V[:,:,idx-lo+1] .= v;
    end

    return U, V
end

end


text:
@testset "setdiffsingleton.jl" begin
    @testset "SetDiffSingleton()" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        b = SeaPearl.IntSetVar(2, 5, "a", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)

        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        @test constraint in a.onDomainChange
        @test constraint in b.onDomainChange
        @test constraint in x.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::SetDiffSingleton)" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(1, 3, "a", trailer)
        b = SeaPearl.IntSetVar(2, 5, "b", trailer)
        x = SeaPearl.IntVar(0, 3, "x", trailer)

        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        SeaPearl.require!(a.domain, 2)

        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([2, 3, 4, 5])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([1, 3])
        @test SeaPearl.required_values(a.domain) == Set{Int}([2])
        @test length(x.domain) == 4

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        SeaPearl.saveState!(trailer)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([3])
        @test SeaPearl.required_values(b.domain) == Set{Int}([2])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([3])
        @test 0 in x.domain
        @test 1 in x.domain
        @test 3 in x.domain
        @test !(2 in x.domain)
        @test constraint.active.value
        @test ("x" => [2]) in prunedDomains

        SeaPearl.assign!(x, 3)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.required_values(b.domain) == Set{Int}([2])
        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([3])
        @test SeaPearl.required_values(a.domain) == Set{Int}([2])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}()
        @test !constraint.active.value

        SeaPearl.restoreState!(trailer)
        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([2, 3, 4, 5])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([1, 3])
        @test SeaPearl.required_values(a.domain) == Set{Int}([2])
        @test length(x.domain) == 4

        # Bug in tsptw:
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(1, 1, "a", trailer)
        b = SeaPearl.IntSetVar(1, 2, "b", trailer)
        x = SeaPearl.IntVar(2, 2, "x", trailer)
        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([1, 2])
        @test SeaPearl.required_values(b.domain) == Set{Int}()
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([1])
        @test SeaPearl.required_values(a.domain) == Set{Int}()
        @test SeaPearl.isbound(x)
        @test SeaPearl.assignedValue(x) == 2

        # more complex
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(1, 5, "a", trailer)
        b = SeaPearl.IntSetVar(1, 8, "b", trailer)
        SeaPearl.exclude!(b.domain, 6)
        SeaPearl.exclude!(b.domain, 2)
        SeaPearl.exclude!(b.domain, 3)
        SeaPearl.require!(b.domain, 1)
        SeaPearl.require!(b.domain, 4)
        x = SeaPearl.IntVar(4, 4, "x", trailer)
        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([5])
        @test SeaPearl.required_values(b.domain) == Set{Int}([1, 4])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([5])
        @test SeaPearl.required_values(a.domain) == Set{Int}([1])

        # if x is assigned
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(1, 2, "a", trailer)
        b = SeaPearl.IntSetVar(1, 1, "b", trailer)
        x = SeaPearl.IntVar(2, 2, "x", trailer)
        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.possible_not_required_values(b.domain) == Set{Int}([1])
        @test SeaPearl.required_values(b.domain) == Set{Int}()
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}([1])
        @test SeaPearl.required_values(a.domain) == Set{Int}()
        @test SeaPearl.isbound(x)
        @test SeaPearl.assignedValue(x) == 2

        # Infeasibility
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(1, 3, "a", trailer)
        b = SeaPearl.IntSetVar(2, 5, "b", trailer)
        x = SeaPearl.IntVar(3, 3, "x", trailer)

        constraint = SeaPearl.SetDiffSingleton(a, b, x, trailer)

        SeaPearl.require!(a.domain, 3)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

    end
end


text:
using STAC
using Documenter

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

makedocs(;
    modules=[STAC],
    authors="Alexander Barth <barth.alexander@gmail.com> and contributors",
    repo="https://github.com/JuliaClimate/STAC.jl/blob/{commit}{path}#{line}",
    sitename="STAC.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaClimate.github.io/STAC.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaClimate/STAC.jl",
    devbranch="main",
)


text:
mutable struct SymmetryVars
    # type of symmetry allowed for the crystal
    #  0 : only the identity element is used
    #  1 : full symmetry group is used
    #  2 : only symmorphic symmetries are allowed
    symtype::Int64
    # number of Bravais lattice point group symmetries
    nsymlat::Int64
    # Bravais lattice point group symmetries
    symlat::Array{Int64,3} #(3,3,48)
    # determinants of lattice symmetry matrices (1 or -1)
    symlatd::Vector{Int64} #(48)
    # index to inverses of the lattice symmetries
    isymlat::Vector{Int64} #(48)
    # lattice point group symmetries in Cartesian coordinates
    symlatc::Array{Float64,3} #(3,3,48)
    # tshift is .true. if atomic basis is allowed to be shifted
    tshift::Bool
    # tsyminv is .true. if the crystal has inversion symmetry
    tsyminv::Bool
    # maximum of symmetries allowed
    maxsymcrys::Int64 # 192
    # number of crystal symmetries
    nsymcrys::Int64
    # crystal symmetry translation vector in lattice and Cartesian coordinates
    vtlsymc::Array{Float64,2} # (3,maxsymcrys)
    vtcsymc::Array{Float64,2} # (3,maxsymcrys)
    # tv0symc is .true. if the translation vector is zero
    tv0symc::Vector{Bool} #(maxsymcrys)
    # spatial rotation element in lattice point group for each crystal symmetry
    lsplsymc::Vector{Int64} #(maxsymcrys)
    # global spin rotation element in lattice point group for each crystal symmetry
    lspnsymc::Vector{Int64} #(maxsymcrys)
    # equivalent atom index for each crystal symmetry
    ieqatom::Array{Float64,3} #(:,:,:)
    # eqatoms(ia,ja,is) is .true. if atoms ia and ja are equivalent
    eqatoms::Array{Bool,3} #(:,:,:)
    # number of site symmetries
    nsymsite::Vector{Int64} #(:)
    # site symmetry spatial rotation element in lattice point group
    lsplsyms::Array{Int64,2} #(:,:)
    # site symmetry global spin rotation element in lattice point group
    lspnsyms::Array{Int64,2} #(:,:)
end

function SymmetryVars()
    symtype = 1
    nsymlat = 1
    symlat = zeros(Int64,3,3,48)
    symlatd = zeros(Int64, 48)
    isymlat = zeros(Int64, 48)
    symlatc = zeros(Float64,3,3,48)
    tshift = false
    tsyminv = false
    maxsymcrys = 192
    nsymcrys = 1
    vtlsymc = zeros(Float64,3,maxsymcrys)
    vtcsymc = zeros(Float64,3,maxsymcrys)
    tv0symc = zeros(Bool,maxsymcrys)
    lsplsymc = zeros(Bool,maxsymcrys)
    lspnsymc = zeros(Int64,maxsymcrys)
    ieqatom = zeros(Float64,1,1,1)
    eqatoms = zeros(Bool,1,1,1)
    nsymsite = zeros(Int64,1)
    lsplsyms = zeros(Int64,1,1)
    lspnsyms = zeros(Int64,1,1)

    return SymmetryVars(
        symtype, nsymlat, symlat, symlatd, isymlat, symlatc,
        tshift, tsyminv, maxsymcrys, nsymcrys, vtlsymc, vtcsymc,
        tv0symc, lsplsymc, lspnsymc, ieqatom, eqatoms,
        nsymsite, lsplsyms, lspnsyms
    )
end


text:
@doc """
    write_hdf5( hdf_out::String, S::GphysData[, KWs] )

Write data to file `hdf_out` from structure `S` in a seismic HDF5 format.

    write_hdf5( hdf_out::String, W::SeisEvent[, KWs] )

Write data to file `hdf_out` from structure `W` in a seismic HDF5 format.
If the format doesn't record event header and source info, only `W.data` is
stored.

## Keywords
#### GphysData
|KW           | Type      | Default   | Meaning                               |
|:---         |:---       |:---       |:---                                   |
| add         | Bool      | false     | Add new traces to file as needed?     |
| chans       | ChanSpec  | 1:S.n     | Channels to write to file             |
| len         | Period    | Day(1)    | Length of new traces added to file    |
| ovr         | Bool      | false     | Overwrite data in existing traces?    |
| tag         | String    | ""        | Tag for trace names in ASDF volumes   |
| v           | Integer   | 0         | verbosity                             |

#### SeisEvent
|KW           | Type      | Default   | Meaning                               |
|:---         |:---       |:---       |:---                                   |
| chans       | ChanSpec  | 1:S.data.n| Channels to write to file             |
| tag         | String    | ""        | Tag for trace names in ASDF volumes   |
| v           | Integer   | 0         | verbosity                             |

## Write Methods
### Add (add = true)
This KW determines the start and end times of all data in `chans`, and
initializes new traces (filled with NaNs) of length = `len`.

#### ASDF behavior
Mode `add=true` follows these steps in this order:
1. Determine times of all data in `S[chans]` and all traces in "Waveforms/".
1. If data lie outside existing trace bounds, new traces are initialized.
1. For each segment in `S[chans]`:
  + Merge the header data in `S[chans]` into the relevant station XML.
  + Overwrite part of the relevant trace in `Waveforms/`.

Thus, unless `len` exactly matches the time boundaries of each segment in `S`,
the traces created will be intentionally larger.

### Overwrite (ovr = true)
If `ovr=true` is specified, but `add=false`, `write_hdf5` *only* overwrites
*existing* data in `hdf_out`.
* No new trace data objects are created in `hdf_out`.
* No new file is created. If `hdf_out` doesn't exist, nothing happens.
* If no traces in `hdf_out` overlap segments in `S`, `hdf_out` isn't modified.
* In ASDF format, station XML is merged in channels that are partly overwritten.

!!! warning

    `add=true`/`ovr=true` changes `:t` on file to begin at an exact sample time.

See also: `read_hdf5`
""" write_hdf5
function write_hdf5(file::String, S::GphysData;
  chans     ::ChanSpec = Int64[], # channels
  fmt       ::String    = "asdf",           # data format
  add       ::Bool      = false,            # add traces
  ovr       ::Bool      = false,            # overwrite trace data
  len       ::Period    = Day(1),           # length of added traces
  tag       ::String    = "",               # trace tag (ASDF)
  v         ::Integer   = KW.v              # verbosity
  )

  chans = mkchans(chans, S, keepirr=false)
  if fmt == "asdf"
    # write_asdf(file, S, chans, add=add, len=len, ovr=ovr, tag=tag, v=v)
    write_asdf(file, S, chans, add, "", ovr, len, tag, v)
  else
    error("Unknown file format (possibly NYI)!")
  end

# log write operation
  opts = string(", chans=", chans,
                ", add=", add,
                ", fmt=\"", fmt,
                "\", len=", len,
                ", ovr=", ovr,
                ", tag=\"", tag,
                "\", v=", v)
  for i in chans
    fwrite_note!(S, i, "write_hdf5", file, opts)
  end
  return nothing
end

function write_hdf5(file::String, C::GphysChannel;
  add       ::Bool      = false,            # add traces
  fmt       ::String    = "asdf",           # data format
  len       ::Period    = Day(1),           # length of added traces
  ovr       ::Bool      = false,            # overwrite trace data
  tag       ::String    = "",               # trace tag (ASDF)
  v         ::Integer   = KW.v              # verbosity
  )

  S = SeisData(C)
  write_hdf5(file, S, fmt=fmt, ovr=ovr, tag=tag, v=v)
  opts = string(", add=", add,
                ", fmt=\"", fmt,
                "\", len=", len,
                ", ovr=", ovr,
                ", tag=\"", tag,
                "\", v=", v)
  fwrite_note!(C, "write_hdf5", file, opts)
  return nothing
end

function write_hdf5(file::String, W::SeisEvent;
  chans     ::ChanSpec = Int64[], # channels
  fmt       ::String    = "asdf",           # data format
  tag       ::String    = "",               # trace tag (ASDF)
  v         ::Integer   = KW.v              # verbosity
  )

  S = getfield(W, :data)
  chans = mkchans(chans, S, keepirr=false)
  if fmt == "asdf"
    H = getfield(W, :hdr)
    R = getfield(W, :source)
    # write_asdf(file, S, chans, evid=H.id, tag=tag, v=v)
    write_asdf(file, S, chans, false, H.id, false, Day(1), tag, v)
    asdf_wqml(file, [H], [R], v=v)
  else
    error("Unknown file format (possibly NYI)!")
  end

  # logging
  opts = string(".data, chans=", chans,
                ", fmt=\"", fmt,
                ", tag=\"", tag,
                "\", v=", v)
  for i in chans
    fwrite_note!(W.data, i, "write_hdf5", file, opts)
  end
  return nothing
end


text:
@testset "899.orderly-queue.jl" begin
    @test orderly_queue("cba", 1) == "acb"
    @test orderly_queue("baaca", 3) == "aaabc"
end


text:
# Logging functionality

export @debug, @trace, @logging_ccall, logging_run

# I/O without libuv, for use after STDOUT is finalized
raw_print(msg::AbstractString...) =
    ccall(:write, Cssize_t, (Cint, Cstring, Csize_t), 1, join(msg), length(join(msg)))
raw_println(msg::AbstractString...) = raw_print(msg..., "\n")

# safe version of `Base.print_with_color`, switching to raw I/O before finalizers are run
# (see `atexit` in `__init_logging__`)
const after_exit = Ref{Bool}(false)
function safe_print_with_color(color::Union{Int, Symbol}, io::IO, msg::AbstractString...)
    if after_exit[]
        raw_print(msg...)
    else
        print_with_color(color, io, msg...)
    end
end

# Display a trace message. Only results in actual printing if the TRACE environment variable
# is set.
const TRACE = haskey(ENV, "TRACE")
@inline function trace(io::IO, msg...; prefix="TRACE: ", line=true)
    @static if TRACE
        safe_print_with_color(:cyan, io, prefix, chomp(string(msg...)), line ? "\n" : "")
    end
end
@inline trace(msg...; kwargs...) = trace(STDERR, msg...; kwargs...)
macro trace(args...)
    TRACE && return Expr(:call, :trace, map(esc, args)...)
end

# Display a debug message. Only results in actual printing if the TRACE or DEBUG environment
# variable is set.
const DEBUG = TRACE || haskey(ENV, "DEBUG")
@inline function debug(io::IO, msg...; prefix="DEBUG: ", line=true)
    @static if DEBUG
        safe_print_with_color(:green, io, prefix, chomp(string(msg...)), line ? "\n" : "")
    end
end
@inline debug(msg...; kwargs...) = debug(STDERR, msg...; kwargs...)
macro debug(args...)
    DEBUG && return Expr(:call, :debug, map(esc, args)...)
end

# Create an indented string from any value (instead of escaping endlines as \n)
function repr_indented(ex; prefix=" "^7, abbrev=true)
    io = IOBuffer()
    print(io, ex)
    str = String(take!(io))

    # Limit output
    if abbrev && length(str) > 256
        if isa(ex, Array)
            T = eltype(ex)
            dims = join(size(ex), " by ")
            if method_exists(zero, (T,)) && zeros(ex) == ex
                str = "$T[$dims zeros]"
            else
                str = "$T[$dims elements]"
            end
        else
            if contains(strip(str), "\n")
                str = str[1:100] * "…\n\n[snip]\n\n…" * str[end-100:end]
            else
                str = str[1:100] * "…" * str[end-100:end]
            end
        end
    end

    lines = split(strip(str), '\n')
    if length(lines) > 1
        for i = 1:length(lines)
            lines[i] = prefix * lines[i]
        end

        lines[1] = "\"\n" * lines[1]
        lines[length(lines)] = lines[length(lines)] * "\""

        return join(lines, '\n')
    else
        return str
    end
end

# ccall wrapper logging the call, its arguments, and the returned value.
# Only logs if TRACE environment variable is set.
macro logging_ccall(fun, target, rettyp, argtypes, args...)
    blk = Expr(:block)

    # print the function name & arguments
    if TRACE
        push!(blk.args, :(trace($(sprint(Base.show_unquoted,fun.args[1])*"("); line=false)))
        i = length(args)
        for arg in args
            i -= 1
            sep = (i>0 ? ", " : "")

            # TODO: we should only do this if evaluating `arg` has no side effects
            push!(blk.args, :(trace(repr_indented($(esc(arg))), $sep;
                  prefix=$(sprint(Base.show_unquoted,arg))*"=", line=false)))
        end
        push!(blk.args, :(trace(""; prefix=") =", line=false)))
    end

    # actual ccall
    @gensym ret
    push!(blk.args, quote
        $ret = ccall($(esc(target)), $(esc(rettyp)), $(esc(argtypes)), $(map(esc, args)...))
    end)

    # print results
    if TRACE
        push!(blk.args, :(trace($ret; prefix=" ")))
    end

    push!(blk.args, :($ret))

    return blk
end

function logging_run(cmd)
    if TRACE
        println(cmd)
    end
    run(cmd)
end

function __init_logging__()
    if TRACE
        trace("CUDA packages running in trace mode, this will generate a lot of additional output")
    elseif DEBUG
        debug("CUDA packages running in debug mode, this will generate additional output")
        debug("Run with TRACE=1 to enable even more output")
    end

    atexit(()->begin
        debug("Dropping down to post-finalizer I/O")
        after_exit[]=true
    end)
end


text:
using NNlib

Jax.@jaxfunc NNlib.leakyrelu(x) = Core.jax.nn.leaky_relu(x)
Jax.@jaxfunc NNlib.celu(x) = Core.jax.nn.celu(x)
Jax.@jaxfunc NNlib.relu(x) = Core.jax.nn.relu(x)
Jax.@jaxfunc NNlib.selu(x) = Core.jax.nn.selu(x)
Jax.@jaxfunc NNlib.gelu(x) = Core.jax.nn.gelu(x)
Jax.@jaxfunc NNlib.elu(x) = Core.jax.nn.elu(x)
Jax.@jaxfunc NNlib.sigmoid(x) = Core.jax.nn.sigmoid(x)
Jax.@jaxfunc NNlib.logsigmoid(x) = Core.jax.nn.log_sigmoid(x)
Jax.@jaxfunc NNlib.swish(x) = Core.jax.nn.swish(x)
Jax.@jaxfunc NNlib.softplus(x) = Core.jax.nn.softplus(x)
Jax.@jaxfunc NNlib.softsign(x) = Core.jax.nn.soft_sign(x)
Jax.@jaxfunc NNlib.hardtanh(x) = Core.jax.nn.hard_tanh(x)

NNlib.softmax(x::AbstractJaxArray; dims=1) =
  jax.nn.softmax(x, axis=_convert_dims(x, dims))

NNlib.logsoftmax(x::AbstractJaxArray; dims=1) =
    jax.nn.log_softmax(x, axis=_convert_dims(x, dims))

NNlib.batched_mul(a::AbstractJaxArray, b::AbstractJaxArray) =
  Core.jax.lax.batch_matmul(b, a)


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 3df3da8e-2bc9-44f5-bbe0-2cabdf3a2a81
begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the global environment
    Pkg.activate("../")

    using PolaronMobility
	using QuadGK
end

# ╔═╡ 69e7d9cd-7900-474d-aeba-06b38e62d756
using Gnuplot

# ╔═╡ 4df7f468-544d-11ec-06a3-0b757c39e918
# Rubrene Ordejon
# Rubrene electron-phonon matrix elements

# Ordejón, P., Boskovic, D., Panhans, M., Ortmann, F., 2017. Ab initio study of
# electron-phonon coupling in rubrene. Phys. Rev. B 96, 035202.
# https://doi.org/10.1103/PhysRevB.96.035202

# ╔═╡ 78d3ef38-fd4a-4957-bc9b-4523cb641df4
begin
	CM1inTHz = 33.356
	hbar=ħ=1.0545718E-34
	q=1.6E-19
	const Boltzmann = const kB =  1.3806504e-23
	const me=MassElectron = 9.10938188e-31
end

# ╔═╡ c4d0bfad-d837-4814-aa60-d71bd41ce28b
# Table II
# Units: cm^-1 meV | convert to THz meV
EffectiveH=[1208.9 106.8] # Holstein

# ╔═╡ 43f24329-1a03-4012-9b31-cf72d273c3d6
cm1tomeV(cm1)=1000 * 2π* ħ*cm1/CM1inTHz*1E12 / q

# ╔═╡ 139b3aed-b4f8-4705-8eff-d1e0d6ea9d12
cm1tomeV(EffectiveH[1]) # I got 149.88 meV by hand

# ╔═╡ e1573de4-450a-4ed8-a9db-003fedcc9cc3
g0 = EffectiveH[2] / cm1tomeV(EffectiveH[1]) # I got 0.71 by hand

# ╔═╡ cc9153d2-86d6-4272-a266-053165197d99
EffectiveP=[117.9 21.9]

# ╔═╡ d303b91b-567c-4705-b51e-5f32c91b3d80
gi = EffectiveP[2] / cm1tomeV(EffectiveP[1])

# ╔═╡ f192deb9-5de9-40e1-959b-a518aee7a6a3
α=gi+g0 # I got 2.21 by hand; this is essentially the Frohlich alpha

# ╔═╡ b1c36a20-9c5b-4fc2-9907-8061ed78d88f
v,w=PolaronMobility.feynmanvw(α)

# ╔═╡ 7f6aedfa-700b-47b5-8030-14d1448c4737
# We've just solved the Feynman polaron theory!

# ╔═╡ 9829e516-9240-48b2-b4f6-9eff06bf5018
EffectiveFreq=(g0*EffectiveH[1] + gi*EffectiveP[1])/α

# ╔═╡ edacbf0c-959f-4a2e-b406-afe7fdfa98bc
ω=1E12 * EffectiveFreq / CM1inTHz

# ╔═╡ 6d6d286a-2f36-4762-803b-59299d7373a1
EffectiveMass(;a=1E-9, J=0.030q)=ħ^2/(2J*a^2)
# Tight-binding model, then cos(2θ) expansion around extremum

# ╔═╡ 88f68d4c-4782-47da-ab8f-6f7ea728f23c

# ╔═╡ 79c95049-89a3-42cc-962a-e5266d8902a1
# Rubrene crystals have an orthorhombic unit cell with lattice parameters a = 1.44 nm, b = 0.718 nm, and c = 2.69 nm.
# https://www.lehigh.edu/~inlo/rubrene.html#:~:text=Rubrene%20crystals%20have%20an%20orthorhombic,the%20a%20and%20c%20directions.
EffectiveMass(a=0.718E-9, J=0.134*q)/MassElectron

# ╔═╡ 58828fe1-12c4-499d-837d-ae45049758e7
mb = 0.88 * MassElectron

# ╔═╡ ac1f76d8-1cdd-4de3-b1e1-049ba33374d2
βred = (q/1000 * cm1tomeV(EffectiveFreq) ) / (kB*300)

# ╔═╡ 691ace17-eb19-444a-9807-887156a813aa

# ╔═╡ 5f62b499-3363-4d0f-8b70-6dabe5e579ed
cm1tomeV(EffectiveH[1])

# ╔═╡ 755da976-d29f-4c82-9645-7acef0a140c5
EffectiveH[1]

# ╔═╡ 78695321-c088-4e99-987c-fca646aa7af3
# Hellwarth1999 - directly do contour integration in Feynman1962, for
        # finite temperature DC mobility
        # Hellwarth1999 Eqn (2) and (1) - These are going back to the general
        # (pre low-T limit) formulas in Feynman1962.  to evaluate these, you
        # need to do the explicit contour integration to get the polaron
        # self-energy
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)

# ╔═╡ 76794859-df83-4c35-bb51-3f34774798e0
#b=R*βred/sinh(b*βred*v/2) # This self-references b! What on Earth?
        # OK! I now understand that there is a typo in Hellwarth1999 and
        # Biaggio1997. They've introduced a spurious b on the R.H.S. compared to
        # the original, Feynman1962:
        b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)

# ╔═╡ 10adfb5e-d9ed-4f5b-bc1e-98b8a452e886
a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))

# ╔═╡ 1e2a33b0-624c-41a5-8c94-5f1f51fbe835
k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)

# ╔═╡ 4a36c74a-d44d-4a92-9d24-35ec22d8ae57
K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

# ╔═╡ 6f49c66c-82b8-428c-8420-4646a6394a9e
#Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K

# ╔═╡ 7255cc2c-4680-4cde-93da-1ba2154e2d22
μ=RHS^-1 * (q)/(ω*mb) #, "m^2/Vs"

# ╔═╡ e6fe559e-46f1-4336-9513-7306c32f2515
μ*100^2 , "cm^2/Vs"

# ╔═╡ ed7cb9bb-4b83-4ed6-a8fd-9936eb02643b
function FHIPmob(T; v=v, w=w, ω=ω)
	#βred = (q/1000 * cm1tomeV(EffectiveP[1]) ) / (kB*T)
	βred = ħ*ω / (kB*T)
	R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
    b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)
	a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
	k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
	K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)
	#Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
	RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
	μ=RHS^-1 * (q)/(ω*mb) #, "m^2/Vs"
	μ=μ*100^2 # convert to cm^2/Vs, from SI
end

# ╔═╡ 084ac351-17e6-4916-b0b7-12aa593984b2
Trange=90:10:400

# ╔═╡ 322962f3-c1f5-47d1-882b-f4efed0c2fa0
mobs=[FHIPmob(T) for T=Trange]

# ╔═╡ 5e996c83-88dc-4569-a435-4d00684de2af
FHIPmob(300)

# ╔═╡ 19b69d34-367e-4794-abb1-3f76e1356561
begin
	@gp Trange mobs " w lp title 'Rubrene - b axis'"
	@gp :- "set yrange [0:]"
	@gp :- "set xlabel 'Temperature (K)'"
	@gp :- "set ylabel 'Mobility (cm^2/Vs)'"
end

# ╔═╡ bbe9a49b-6495-4c13-91d0-3a4dbd0a2800
begin
	@gp Trange mobs " w lp title 'Rubrene - b axis'"
	#@gp :- "set yrange [0:]"
	@gp :- "set xlabel 'Temperature (K)'"
	@gp :- "set ylabel 'Mobility (cm^2/Vs)'"
	@gp :- "set logscale y"
end

# ╔═╡ 444d7454-48a9-44dd-87a3-f4cf1a495072
@printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)

# ╔═╡ b0f72692-d76c-44f1-91aa-2e9b7b882d1b
# Jan's COF
Jperylenecoff=0.204

# ╔═╡ 0c043a99-18bf-4b75-95fa-ce8aea8f8bef
relevantmode=200 #cm^-1

# ╔═╡ 06ccb204-694d-4c64-9ecd-742c009fb972
perylenecofmass=EffectiveMass(a=1.0E-9, J=Jperylenecoff*q)

# ╔═╡ 2ba19614-f2ef-4dca-9f61-b95908b7c81c
perylenecofmass/MassElectron

# ╔═╡ 27ee4adf-b7d9-41e0-b0a0-f414422b7b73
perylenecofg0 = 128 / cm1tomeV(relevantmode)

# ╔═╡ ed6d3a1a-3fcd-4ef1-a57f-a394895d4b29
pv,pw= feynmanvw(perylenecofg0)

# ╔═╡ fbee75ec-96f0-41b6-8c14-a724872476fb
FHIPmob(300; v=pv, w=pw, ω=relevantmode/CM1inTHz*1E12*2π)

# ╔═╡ 7d7d3ff6-30d4-477a-8f03-e4c08ac8a608
md"""
# Stuff below here for later
"""

# ╔═╡ ba64c0a2-8cd4-421a-b494-a70f1f92ca0f
EffectiveHGirlando=[1277 99]

# ╔═╡ d30b10f4-2c62-4596-aa73-e4f2796b0541
EffectivePGirlando=[77   20]

# ╔═╡ 04b0a526-b808-4c61-a565-6c60e12be2b7
# Table IV
# ωp      ωpg0  ωpgi
# cm^-1    meV   meV
PerModeOrdejon=
[57.8    -1.7   0.85
 59.6     1.4  -0.83
 89.0     1.6  -4.8
 107.3   -0.14  2.8
 139.1   -2.3  -3.7
 639.1   -7.5   1.0
 1011.2  -3.6  -0.04
 1344.7  19.8   0.04
 1593.3 -42.0  -0.12
] .* [1/CM1inTHz 1 1]

# ╔═╡ b38b707f-6b97-424a-a22a-cf2fdebcfb31
PerModeGirlando =
[37.4    -0.9   3.4
 66.6     1.6  -6.6
 86.7    -0.6  -9.3
 106.3    0    -4.4
 125.1    1.4  -4.7
 631.2  -10.8   1.3
 1002.3  24.6   0
 1348.6  49.9   0
 1593.8 -45.6   1.6
] .* [1/CM1inTHz 1 1]

# ╔═╡ Cell order:
# ╠═4df7f468-544d-11ec-06a3-0b757c39e918
# ╠═3df3da8e-2bc9-44f5-bbe0-2cabdf3a2a81
# ╠═78d3ef38-fd4a-4957-bc9b-4523cb641df4
# ╠═c4d0bfad-d837-4814-aa60-d71bd41ce28b
# ╠═43f24329-1a03-4012-9b31-cf72d273c3d6
# ╠═139b3aed-b4f8-4705-8eff-d1e0d6ea9d12
# ╠═e1573de4-450a-4ed8-a9db-003fedcc9cc3
# ╠═cc9153d2-86d6-4272-a266-053165197d99
# ╠═d303b91b-567c-4705-b51e-5f32c91b3d80
# ╠═f192deb9-5de9-40e1-959b-a518aee7a6a3
# ╠═b1c36a20-9c5b-4fc2-9907-8061ed78d88f
# ╠═7f6aedfa-700b-47b5-8030-14d1448c4737
# ╠═9829e516-9240-48b2-b4f6-9eff06bf5018
# ╠═edacbf0c-959f-4a2e-b406-afe7fdfa98bc
# ╠═6d6d286a-2f36-4762-803b-59299d7373a1
# ╠═88f68d4c-4782-47da-ab8f-6f7ea728f23c
# ╠═79c95049-89a3-42cc-962a-e5266d8902a1
# ╠═58828fe1-12c4-499d-837d-ae45049758e7
# ╠═ac1f76d8-1cdd-4de3-b1e1-049ba33374d2
# ╠═691ace17-eb19-444a-9807-887156a813aa
# ╠═5f62b499-3363-4d0f-8b70-6dabe5e579ed
# ╠═755da976-d29f-4c82-9645-7acef0a140c5
# ╠═78695321-c088-4e99-987c-fca646aa7af3
# ╠═76794859-df83-4c35-bb51-3f34774798e0
# ╠═10adfb5e-d9ed-4f5b-bc1e-98b8a452e886
# ╠═1e2a33b0-624c-41a5-8c94-5f1f51fbe835
# ╠═4a36c74a-d44d-4a92-9d24-35ec22d8ae57
# ╠═6f49c66c-82b8-428c-8420-4646a6394a9e
# ╠═7255cc2c-4680-4cde-93da-1ba2154e2d22
# ╠═e6fe559e-46f1-4336-9513-7306c32f2515
# ╠═ed7cb9bb-4b83-4ed6-a8fd-9936eb02643b
# ╠═084ac351-17e6-4916-b0b7-12aa593984b2
# ╠═322962f3-c1f5-47d1-882b-f4efed0c2fa0
# ╠═5e996c83-88dc-4569-a435-4d00684de2af
# ╠═69e7d9cd-7900-474d-aeba-06b38e62d756
# ╠═19b69d34-367e-4794-abb1-3f76e1356561
# ╠═bbe9a49b-6495-4c13-91d0-3a4dbd0a2800
# ╠═444d7454-48a9-44dd-87a3-f4cf1a495072
# ╠═b0f72692-d76c-44f1-91aa-2e9b7b882d1b
# ╠═0c043a99-18bf-4b75-95fa-ce8aea8f8bef
# ╠═06ccb204-694d-4c64-9ecd-742c009fb972
# ╠═2ba19614-f2ef-4dca-9f61-b95908b7c81c
# ╠═27ee4adf-b7d9-41e0-b0a0-f414422b7b73
# ╠═ed6d3a1a-3fcd-4ef1-a57f-a394895d4b29
# ╠═fbee75ec-96f0-41b6-8c14-a724872476fb
# ╠═7d7d3ff6-30d4-477a-8f03-e4c08ac8a608
# ╠═ba64c0a2-8cd4-421a-b494-a70f1f92ca0f
# ╠═d30b10f4-2c62-4596-aa73-e4f2796b0541
# ╠═04b0a526-b808-4c61-a565-6c60e12be2b7
# ╠═b38b707f-6b97-424a-a22a-cf2fdebcfb31


text:
# Cases must include `case_name` and `scm_suffix`. If `case_name`=="LES_driven_SCM", must also include `y_dir`.
struct ScmTest end

function get_reference_config(::ScmTest)
    config = Dict()

    # amip4K data: April
    cfsite_numbers = (13, 21, 23)
    les_kwargs = (forcing_model = "HadGEM2-A", month = 4, experiment = "amip4K")

    ref_dirs = [get_cfsite_les_dir(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]
    suffixes = [get_gcm_les_uuid(cfsite_number; les_kwargs...) for cfsite_number in cfsite_numbers]
    n_repeat = length(ref_dirs)
    config["case_name"] = repeat(["LES_driven_SCM"], n_repeat)
    config["y_dir"] = ref_dirs
    config["scm_suffix"] = suffixes
    return config
end


text:
img = zeros(Gray, 50, 50);

# create a cell with a background half a bright around it
centers = [Point(20, 20), Point(20, 20)]
radii = [15, 5]
# the local background will be twice as dim as the "cell"
colors = [Gray(0.5), Gray(1.0)]

ImageDraw.draw!(img, [ImageDraw.CirclePointRadius(c, r) for (c,r) in zip(centers, radii)], colors)
labels = Images.label_components(img .== 1.0)
localities = SegmentationTools.get_localities(labels, [1], dist=(0, 10))

@test median(Float64.(img[localities[1]])) == 0.5

img = zeros(Gray, 50, 50);
img[10, 10] = Gray(1.0);
img[20, 10] = Gray(1.0);
img[30, 10] = Gray(1.0);

labels = Images.label_components(img .== 1.0)
labels[labels .== 2] .= 0
ids = [1, 3]

msg = "Do not remove values from `labels`, only from `ids`. Calc might be wrong!"
@test_logs (:warn, msg) get_localities(labels, ids; dist=(5, 15))

labels = Images.label_components(img .== 1.0)

ids = [1]
localities = get_localities(labels, ids, dist=(5,50))

@test all(keys(localities) .== ids)

locality = Set(localities[1])

# all of the points should be missing from the locality
@test length(intersect(findall(img .> 0), locality)) == 0


text:
using ProjAssistant
@quickactivate

@time begin
    import SparseArrays
    import Base.Threads: @threads, threadid, SpinLock

    import SimTools
    const SimT = SimTools

    using Serialization

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006

    const iJR = ChN.iJR904
    const Nd = ChN.NanchenData # experimental data
    const Bd = ChN.BegData    # cost data

    import Chemostat
    import Chemostat.LP: MathProgBase
    const Ch = ChN.Chemostat
    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP
    const ChEP = Ch.MaxEntEP
    const ChSU = Ch.SimulationUtils

    using ArgParse
    using Statistics
    using Serialization
    using Base.Threads
    using ExtractMacro
end

## ----------------------------------------------------------------------------
# arg settings
ARGSET = ArgParse.ArgParseSettings()
@ArgParse.add_arg_table! ARGSET begin
    "--ignore-cached"
        help = "Ingnore on disk version of data"
        action = :store_true
end
ARGS_DICT = ArgParse.parse_args(ARGSET)

if isinteractive()
    # dev values
    ignore_cached = true
else
    ignore_cached = ARGS_DICT["ignore-cached"]
end
@info("ARGS", ignore_cached)

## ----------------------------------------------------------------------------
# FILE GLOBALS
const WLOCK = ReentrantLock()
const IGNORE_EXPS = []

## -------------------------------------------------------------------
# UTILS
dat_file(;kwargs...) = procdir(iJR, "maxent_ep_dat", kwargs..., ".jls")

function is_cached(;kwargs...)
    ignore_cached && return false
    thid = threadid()
    datfile = dat_file(;kwargs...)
    isfile(datfile) && lock(WLOCK) do
        @info("Cache Found",
            datfile, thid
        ); println()
    end; isfile(datfile)
end

## -------------------------------------------------------------------
# MAXENT FUNS
include("2.0.1_ME_MAX_POL.jl")
include("2.0.2_ME_MAX_POL_B0.jl")
include("2.0.3_ME_Z_EXPECTED_G_BOUNDED.jl")

## -------------------------------------------------------------------
# ME GLOBALS
sglob(iJR, :maxent, :params) do
    (;
        alpha = Inf,
        epsconv = 9e-4,
        maxiter = 3000,
        damp = 0.9,
        maxvar = 1e50,
        minvar = 1e-50,
    )
end

## -------------------------------------------------------------------
# ME_MAX_POL
let
    method = iJR.ME_MAX_POL
    model_key = "max_model"

    maxent_max_pol(method, model_key)
end

## ----------------------------------------------------------------------------
# ME_MAX_POL_B0
let
    method = iJR.ME_MAX_POL_B0
    model_key = "max_model"

    do_max_pol_b0(method, model_key)
end

## ----------------------------------------------------------------------------
# ME_Z_EXPECTED_G_BOUNDED
let
    method = iJR.ME_Z_EXPECTED_G_BOUNDED
    model_key = "fva_models"

    do_z_expected_ug_bounded(method, model_key)
end