text:
using LinearAlgebra, ForwardDiff, Random, NLPModels

export Zoom, StrongWolfe_LineSearch

function Zoom(nlp::AbstractNLPModel, x, r, alo, ahi, η1::Float64, η2::Float64; max_iter = 20)

    f(x) = obj(nlp, x)
    ∇f(x) = grad(nlp, x)

    iter = 0
    α = 0
    while true
        a = (alo + ahi) / 2

        if f(x + a * r) > f(x) + η1 * a * dot(∇f(x), r) || f(x + a * r) ≥ f(x + alo * r)
            ahi = a
        else
            if abs(dot(∇f(x + a * r), r)) ≤ - η2 * dot(∇f(x), r)
                α = a
                break
            end
            if dot(∇f(x + a * r), r) * (ahi - alo) ≥ 0
                ahi = alo
            end
            alo = a
        end

        if iter == max_iter
            α = a
            break
        end

        iter +=1
    end

    return α, iter
end

function StrongWolfe_LineSearch(nlp::AbstractNLPModel, x, r, a1::Float64, η1::Float64, η2::Float64, p::Float64; max_iter = 100)

    amax = 10 * a1

    a0 = 0
    iter = 0

    f(x) = obj(nlp, x)
    ∇f(x) = grad(nlp, x)

    f_old = 0
    a = 0

    while true
        if f(x + a1 * r) > f(x) + η1 * a1 * dot(∇f(x),r) || (iter > 0 && f(x + a1 * r) > f_old)
            a, iter = Zoom(nlp, x, r, a0, a1, η1, η2);
            break
        end

        if abs(dot(∇f(x + a1 * r), r)) ≤ - η2 * dot(∇f(x), r)
            a = a1
            break
        end

        if dot(∇f(x + a1 * r), r) ≥ 0
        a, iter = Zoom(nlp, x, r, a1, a0, η1, η2);
            break
        end

        iter += 1
        a0 = a1
        a1 = p * a0 + (1 - p) * amax
        # a1 += p * (amax - a1)
        f_old = f(x + a0 * r)

        if iter == max_iter
            # println("LineaSearch Max iter")
            a = a1
            break
        end

    end
    return a, iter
end


text:
include("../../fancy_toys.jl")

version = v"7.6.5"

name = "CUDNN_CUDA$(cuda_version.major).$(cuda_version.minor)"

script = raw"""
cd ${WORKSPACE}/srcdir
if [[ ${target} == x86_64-linux-gnu ]]; then
    cd cuda
    find .

    install_license NVIDIA_SLA_cuDNN_Support.txt

    mv lib64/libcudnn.so* ${libdir}
    mv include/* ${prefix}/include
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    cd cuda
    find .

    install_license NVIDIA_SLA_cuDNN_Support.txt

    mv bin/cudnn64_*.dll ${libdir}
    mv include/* ${prefix}/include
elif [[ ${target} == x86_64-apple-darwin* ]]; then
    cd cuda
    find .

    install_license NVIDIA_SLA_cuDNN_Support.txt

    mv lib/libcudnn.*dylib ${libdir}
    mv include/* ${prefix}/include
fi
"""

products = [
    LibraryProduct(["libcudnn", "cudnn64_$(version.major)"], :libcudnn),
]

dependencies = [Dependency(PackageSpec(name="CUDA_jll", version=cuda_version))]

non_reg_ARGS = filter(arg -> arg != "--register", ARGS)

if @isdefined(sources_linux) && should_build_platform("x86_64-linux-gnu")
    build_tarballs(non_reg_ARGS, name, version, sources_linux, script,
                   [Linux(:x86_64)], products, dependencies)
end

if @isdefined(sources_windows) && should_build_platform("x86_64-w64-mingw32")
    build_tarballs(ARGS, name, version, sources_windows, script,
                   [Windows(:x86_64)], products, dependencies)
end

if @isdefined(sources_macos) && should_build_platform("x86_64-apple-darwin14")
    build_tarballs(ARGS, name, version, sources_macos, script,
                   [MacOS(:x86_64)], products, dependencies)
end

ext: jl

max_stars_repo_name: jbrea/Yggdrasil

max_issues_repo_name: jbrea/Yggdrasil

max_forks_repo_name: jbrea/Yggdrasil


text:
Base.@kwdef struct FEMesh{T}
    point_map::Vector{Vector{Int}}
    points::Vector{Matrix{T}}
    facets::Matrix{Matrix{Int}}
    elements::Vector{Matrix{Int}}
end


text:
module Utility

using LinearAlgebra
using Statistics, StatsBase

export sinkhorn

"""
    sinkhorn(A; r=[], c=[], maxit=1000, δ=1e-6, verbose=false)

Compute the row and column multiplicative factors that constrain marginals of `A` to unity.
Returns a boolean flag if algorithm converges within tolerance `δ` within `maxit` iterations.
If `r` is given and non-empty, row marginals will be constrained to user-supplied inputs.
If `c` is given and non-empty, column marginals will be constrained to user-supplied inputs.
"""
function sinkhorn(A; r=[], c=[], maxit=1000, δ=1e-6, verbose=false)
    if length(r) == 0
        r = size(A,2)
    end

    if length(c) == 0
        c = size(A,1)
    end

    x = ones(size(A,1))
    y = ones(size(A,2))
    for i ∈ 1:maxit
        δr = maximum(abs.(x.*(A*y)  .- r))
        δc = maximum(abs.(y.*(A'*x) .- c))

        if verbose
            @show minimum(A'*x), minimum(A*y)
            @show r, c
            @show i, δr, δc
        end

        (isnan(δr) || isnan(δc)) && return x, y, false
        (δr < δ && δc < δ) 		 && return x, y, true

        y = c ./ (A'*x)
        x = r ./ (A*y)
        if verbose
            @show mean(x), mean(y)
        end
    end

    return x, y, false
end

function sinkhorn(N, ϕ; maxit=1000, δ=1e-6, verbose=false)
    nrow = size(N,1)
    ncol = size(N,2)

    N² = N.^2
    u, v = ones(size(N,1)), ones(size(N,2))
    for i ∈ 1:maxit
        σ² = ((Diagonal(u.^2 ./ (1 .+ ϕ))*N²*Diagonal(v.^2)) .+ (Diagonal(u.*ϕ ./ (1 .+ ϕ))*N*Diagonal(v)))
        δr = maximum(abs.(sum(σ², dims=1) .- nrow))
        δc = maximum(abs.(sum(σ², dims=2) .- ncol))

        if verbose
            @show i, δr, δc
        end

        (isnan(δr) || isnan(δc)) && return u, v, false
        (δr < δ && δc < δ) 		 && return u, v, true

        A = vec(sum(Diagonal(u.^2 ./ (1 .+ ϕ)) * N², dims=1))
        B = vec(sum(Diagonal(u.*ϕ ./ (1 .+ ϕ)) * N , dims=1))

        v = @. (sqrt(B^2 + 4*A*nrow) - B) / (2*A)

        A = vec(sum(N² * Diagonal(v.^2),           dims=2))
        B = vec(sum(Diagonal(ϕ) * N * Diagonal(v), dims=2))

        u = @. (sqrt(B^2 + 4*A*ncol*(1+ϕ)) - B) / (2*A)
    end

    return u, v, false
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
const CHR_REG = r"chr[1-9]+"

check_chr_format(x) = @assert all(occursin.(CHR_REG, x))

function allele_specific_binding_snps(asb_prefix)
    bQTLs = DataFrame()
    asbdir_, prefix_ = splitdir(asb_prefix)
    asbdir = asbdir_ == "" ? "." : asbdir_
    for file in readdir(asbdir)
        if occursin(prefix_, file)
            new_bQTLs = CSV.read(joinpath(asbdir_, file), DataFrame, select=[:ID, :CHROM, :isASB])
            bQTLs = vcat(
                bQTLs,
                DataFrames.select(filter(x -> x.isASB === true && x.CHROM ∉ ("chrX", "chrY"), new_bQTLs), [:ID, :CHROM])
            )
        end
    end
    check_chr_format(bQTLs.CHROM)
    bQTLs.QTL_TYPE .= :bQTL
    return unique(bQTLs, :ID)
end

function trans_actors(path)
    eQTLs = CSV.read(path, DataFrame, select=[:ID, :CHROM])
    check_chr_format(eQTLs.CHROM)
    eQTLs.QTL_TYPE .= :eQTL
    return eQTLs
end

function add_chrfiles(snps, chr_prefix)
    chr_dir_, prefix_ = splitdir(chr_prefix)
    chr_dir = chr_dir_ == "" ? "." : chr_dir_
    required_chrs = unique(snps.CHROM)
    chr_files = DataFrame(CHROM=[], BGEN_FILE=[])
    for filename in readdir(chr_dir)
        if occursin(prefix_, filename) && endswith(filename, "bgen")
            regexp_match = match(CHR_REG, filename)
            if regexp_match !== nothing && regexp_match.match ∈ required_chrs
                filepath = joinpath(chr_dir_, filename)
                sample_filepath = string(filepath[1:end-4], "sample")
                idx_filepath = string(filepath, ".bgi")
                bgenfile = Bgen(filepath, sample_path=sample_filepath, idx_path=idx_filepath)
                push!(chr_files, [regexp_match.match, bgenfile])
            end
        end
    end
    nb_required_chrs = size(required_chrs, 1)
    nb_actual_chrs = size(chr_files, 1)
    @assert nb_required_chrs == nb_actual_chrs "Required "*
        "$nb_required_chrs chromosome files, "*
        "got: $nb_actual_chrs. Chromosome files: \n $(chr_files.CHROM)"
    return innerjoin(snps, chr_files, on=:CHROM)
end

exclude!(snps, exclusion_file::Nothing) = nothing
exclude!(snps, exclusion_file) = filter!(:ID => x -> x ∉ readlines(exclusion_file), snps)

function impute_genotypes(probabilities, variant_genotypes, threshold)
    n = size(probabilities)[2]
    t = Vector{Union{String, Missing}}(missing, n)
    for i in 1:n
        # If no allele has been annotated with sufficient confidence
        # the sample is declared as missing for this variant
        sample_gen_index = findfirst(x -> x >= threshold, probabilities[:, i])
        sample_gen_index isa Nothing || (t[i] = variant_genotypes[sample_gen_index])
    end
    return t
end

"""
Only bi-allelic SNPs are considered. We enforce that the allele₁, allele₍
genotype is encoded with the major variant coming first.
"""
function genotypes(variant::Variant, minor, major)
    alleles_ = alleles(variant)
    @assert size(alleles_, 1) == 2 "$(variant.rsid) is not bi-allelic"
    all₁, all₂ = alleles_
    @assert (all₁, all₂) == (minor, major) || (all₁, all₂) == (major, minor)
    return [all₁*all₁, major*minor, all₂*all₂]
end

function impute_genotypes(snps, threshold)
    snps_genotypes = nothing
    rsid_to_minor_major = Dict()
    for group in groupby(snps, :CHROM)
        bgenfile = first(group.BGEN_FILE)
        chr_genotypes= DataFrame(SAMPLE_ID=bgenfile.samples)
        # Iterate over variants in this chromosome
        for rsid in group.ID
            v = variant_by_rsid(bgenfile, rsid)
            minor_allele_dosage!(bgenfile, v)
            major = major_allele(v)
            minor = minor_allele(v)
            rsid_to_minor_major[rsid] = (minor=minor, major=major)
            variant_genotypes = genotypes(v, minor, major)
            probabilities = probabilities!(bgenfile, v)
            chr_genotypes[!, rsid] = impute_genotypes(probabilities, variant_genotypes, threshold)
        end
        # I think concatenating should suffice but I still join as a safety
        snps_genotypes = snps_genotypes isa Nothing ? chr_genotypes :
            innerjoin(snps_genotypes, chr_genotypes, on=:SAMPLE_ID)
    end
    return snps_genotypes, rsid_to_minor_major
end

function build_queries(eQTLs, bQTLs, genotypes, rsid_to_major_minor, threshold)
    queries = []
    n = size(genotypes, 1)
    for eqtl in eachrow(eQTLs)
        # Control is MAJOR*MAJOR
        # Treatment is MAJOR*MINOR
        eqtl_minor, eqtl_major = rsid_to_major_minor[eqtl.ID]
        eqtl_control = eqtl_major*eqtl_major
        eqtl_case = eqtl_major*eqtl_minor
        for bqtl in eachrow(bQTLs)
            bqtl_minor, bqtl_major = rsid_to_major_minor[bqtl.ID]
            bqtl_control = bqtl_major*bqtl_major
            bqtl_case = bqtl_major*bqtl_minor

            counts = DataFrames.combine(
                groupby(genotypes, [eqtl.ID, bqtl.ID]; skipmissing=true),
                x -> size(x, 1)
            )
            of_interest = filter(x-> x[1] ∈ (eqtl_control, eqtl_case) && x[2] ∈ (bqtl_control, bqtl_case), counts)
            if size(of_interest, 1) == 4 && minimum(of_interest[!, 3])/n >= threshold
                push!(queries, Dict(
                        "SNPS" => Dict(eqtl.ID => eqtl.CHROM, bqtl.ID => bqtl.CHROM),
                        "HOMOZYGOUS_MAJOR_TO_MAJOR_MINOR" => Dict(
                            eqtl.ID => eqtl_control*" -> "*eqtl_case,
                            bqtl.ID => bqtl_control*" -> "*bqtl_case
                        )
                    )
                )
            end
        end
    end
    return queries
end

function write_outputs(genotypes, queries, outdir)
    CSV.write(joinpath(outdir, "genotypes.csv"), genotypes)
    for query in queries
        query_filename = "query_"
        for (snp, _) in query["SNPS"]
            query_filename = string(query_filename, snp, "__")
        end
        query_filename = string(query_filename[1:end-2], ".toml")
        open(joinpath(outdir, query_filename), "w") do io
            TOML.print(io, query)
        end
    end
end

function read_queries(query_prefix)
    queries = Dict[]
    querydir_, prefix = splitdir(query_prefix)
    querydir = querydir_ == "" ? "." : querydir_
    for filename in readdir(querydir)
        if occursin(prefix, filename)
            push!(queries, TOML.parse(open(joinpath(querydir_, filename))))
        end
    end
    return queries
end

function snps_from_queries(queries)
    snps = DataFrame(ID=String[], CHROM=String[])
    for query in queries
        for (snp, chr) in query["SNPS"]
            push!(snps, (snp, chr))
        end
    end
    return unique(snps, :ID)
end

function maybe_filter()

end

function asb_generated_mode(parsed_args)
    # Load SNPS
    eQTLs = trans_actors(parsed_args["trans-actors"])
    bQTLs = allele_specific_binding_snps(parsed_args["asb-prefix"])
    # Build genotypes
    snps = vcat(eQTLs, bQTLs)
    exclude!(snps, parsed_args["exclude"])
    snps = add_chrfiles(snps, parsed_args["chr-prefix"])
    genotypes, rsid_to_minor_major = impute_genotypes(snps, parsed_args["call-threshold"])
    # Build queries
    eQTLs = filter(:QTL_TYPE => ==(:eQTL), snps)
    bQTLs = filter(:QTL_TYPE => ==(:bQTL), snps)
    queries = build_queries(eQTLs, bQTLs, genotypes, rsid_to_minor_major, parsed_args["minor-genotype-freq"])
    # write genotypes and queries
    write_outputs(genotypes, queries, parsed_args["outdir"])
end

function given_queries_mode(parsed_args)
    # Read provided queries
    queries = read_queries(parsed_args["query-prefix"])
    # Get snps from queries and add chromosome file information
    snps = snps_from_queries(queries)
    exclude!(snps, parsed_args["exclude"])
    snps = add_chrfiles(snps, parsed_args["chr-prefix"])
    # Get required genotypes
    genotypes, _ = impute_genotypes(snps, parsed_args["call-threshold"])
    # write genotypes and rewrite queries to outdir
    # maybe in the future check that genotypes in queries
    #are not too rare as for the asb mode
    write_outputs(genotypes, queries, parsed_args["outdir"])
end

function build_genotypes_and_queries(parsed_args)
    if parsed_args["mode"] == "asb"
        asb_generated_mode(parsed_args)
    elseif parsed_args["mode"] == "given"
        given_queries_mode(parsed_args)
    else
        throw(ArgumentError("Not implemented yet."))
    end
end


text:
n = 30
m = 20
p = 5

@testset "Testing $(MLK.unsafe_base_evaluate)" begin
    for f in base_functions
        F = (f)()
        f_agg_tmp = get(base_functions_aggregate, f, (s,x,y) -> NaN)
        f_ret_tmp = get(base_functions_return, f, s -> s)

        for T in FloatingPointTypes
            x = rand(T,p)
            y = rand(T,p)

            s = zero(T)
            for i = 1:p
                s = f_agg_tmp(s, x[i], y[i])
            end
            s = f_ret_tmp(s)

            @test isapprox(s, MLK.unsafe_base_evaluate(F, x, y))
        end
    end
end

@testset "Testing $(MLK.base_evaluate)" begin
    for f in base_functions
        F = (f)()

        for T in FloatingPointTypes
            v = rand(T,p+1)
            x = rand(T,p)
            y = rand(T,p)

            @test isapprox(MLK.base_evaluate(F, x, y),       MLK.unsafe_base_evaluate(F, x, y))
            @test isapprox(MLK.base_evaluate(F, x[1], y[1]), MLK.unsafe_base_evaluate(F, x[1:1], y[1:1]))

            @test_throws DimensionMismatch MLK.base_evaluate(F, x, v)
            @test_throws DimensionMismatch MLK.base_evaluate(F, v, x)
            @test_throws DimensionMismatch MLK.base_evaluate(F, T[], T[])
        end
    end
end

@testset "Testing $(MLK.allocate_basematrix)" begin
    for T in FloatingPointTypes
        X = rand(T,n,m)

        for layout in (Val(:row), Val(:col))
            n_tmp = layout == Val(:row) ? n : m

            K = MLK.allocate_basematrix(layout, X)
            @test size(K) == (n_tmp,n_tmp)
            @test eltype(K) == T
        end
    end
end

@testset "Testing $(MLK.checkdimensions)" begin
    for layout in (Val(:row), Val(:col))
        isrowmajor = layout == Val(:row)
        dim = isrowmajor ? 1 : 2

        X =     isrowmajor ? Array{Float64}(undef, n,p)   : Array{Float64}(undef, p,n)
        Y =     isrowmajor ? Array{Float64}(undef, m,p)   : Array{Float64}(undef, p,m)
        Y_bad = isrowmajor ? Array{Float64}(undef, m,p+1) : Array{Float64}(undef, p+1,m)

        @test MLK.checkdimensions(layout, Array{Float64}(undef, n,n), X) == n
        @test_throws DimensionMismatch MLK.checkdimensions(layout, Array{Float64}(undef, n,n+1), X)
        @test_throws DimensionMismatch MLK.checkdimensions(layout, Array{Float64}(undef, n+1,n+1), X)

        @test MLK.checkdimensions(layout, Array{Float64}(undef, n,m), X, Y) == (n,m)
        @test_throws DimensionMismatch MLK.checkdimensions(layout, Array{Float64}(undef, n,m+1), X, Y)
        @test_throws DimensionMismatch MLK.checkdimensions(layout, Array{Float64}(undef, n+1,m), X, Y)
        @test_throws DimensionMismatch MLK.checkdimensions(layout, Array{Float64}(undef, n,m), X, Y_bad)
    end
end

@testset "Testing $(MLK.squared_distance!)" begin
    for T in FloatingPointTypes
        Set_X = [rand(T,p) for i = 1:n]
        Set_Y = [rand(T,p) for i = 1:m]

        X = permutedims(hcat(Set_X...))
        Y = permutedims(hcat(Set_Y...))

        P = [LinearAlgebra.dot(x-y,x-y) for x in Set_X, y in Set_X]
        G = MLK.gramian!(Val(:row), Array{T}(undef, n,n), X, true)
        xtx = MLK.dotvectors(Val(:row), X)

        @test isapprox(MLK.squared_distance!(G, xtx, true), P)

        P = [LinearAlgebra.dot(x-y,x-y) for x in Set_X, y in Set_Y]
        G = MLK.gramian!(Val(:row), Array{T}(undef, n,m), X, Y)
        xtx = MLK.dotvectors(Val(:row), X)
        yty = MLK.dotvectors(Val(:row), Y)
        MLK.squared_distance!(G, xtx, yty)

        @test isapprox(G, P)
        @test all(G .>= 0)

        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 3,4), Array{T}(undef, 3), true)
        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 4,3), Array{T}(undef, 3), true)

        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 3,4), Array{T}(undef, 2), Array{T}(undef, 4))
        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 3,4), Array{T}(undef, 4), Array{T}(undef, 4))
        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 3,4), Array{T}(undef, 3), Array{T}(undef, 3))
        @test_throws DimensionMismatch MLK.squared_distance!(Array{T}(undef, 3,4), Array{T}(undef, 3), Array{T}(undef, 5))
    end
end

@testset "Testing $(MLK.basematrix!)" begin
    for T in (Float32, Float64)
        X_set = [rand(T,p) for i = 1:n]
        Y_set = [rand(T,p) for i = 1:m]

        P_tst_nn = Array{T}(undef, n, n)
        P_tst_nm = Array{T}(undef, n, m)

        for layout in (Val(:row), Val(:col))
            X = layout == Val(:row) ? permutedims(hcat(X_set...)) : hcat(X_set...)
            Y = layout == Val(:row) ? permutedims(hcat(Y_set...)) : hcat(Y_set...)

            for f in base_functions
                F = (f)()

                P = [MLK.base_evaluate(F,x,y) for x in X_set, y in X_set]
                @test isapprox(P, MLK.basematrix!(layout, P_tst_nn, F, X, true))

                P = [MLK.base_evaluate(F,x,y) for x in X_set, y in Y_set]
                @test isapprox(P, MLK.basematrix!(layout, P_tst_nm, F, X, Y))
            end
        end
    end

    for layout in (Val(:row), Val(:col))
        # Test case where squared_distance! may return a negative value
        v1 = [0.2585096115890490597877260370296426117420196533203125
              0.9692536801431554938091039730352349579334259033203125
              0.4741774214537994858176261914195492863655090332031250]

        v2 = [0.2585096115824996876320085448242025449872016906738281
              0.9692536801897344567180425656260922551155090332031250
              0.4741774214882835125628446348855504766106605529785156]

        v3 = rand(Float64, 3)

        F = SquaredEuclidean()
        X = layout == Val(:row) ? permutedims(hcat(v1, v2)) : hcat(v1, v2)
        Y = layout == Val(:row) ? permutedims(hcat(v1, v2, v3)) : hcat(v1, v2, v3)

        @test all(MLK.basematrix!(layout, Array{Float64}(undef, 2,2), F, X, true) .>= 0.0)
        @test all(MLK.basematrix!(layout, Array{Float64}(undef, 2,3), F, X, Y) .>= 0.0)
    end
end


text:
# function query(key::AbstractString, fieldname="acronym")
#     rq = HTTP.request("GET", "http://api.brain-map.org/api/v2/data/SectionDataSet/query.json?criteria=products[id\$eq1],genes[$fieldname\$eq'$key']&include=genes,section_images")
# end

const planedict = Dict("coronal"=>1, "sagittal"=>2)

function query_insitu(key::AbstractString; plane="coronal", fieldname="acronym")
    rq = HTTP.request("GET", "http://api.brain-map.org/api/v2/data/SectionDataSet/query.json?criteria=products[id\$eq1],genes[$fieldname\$eq'$key']&include=section_images,alignment3d,section_images(alignment2d)")
    data = JSON.parse(String(rq.data))
    planeid = planedict[plane]
    for sectiondataset in data["msg"]
        if sectiondataset["plane_of_section_id"] == planeid
            return sectiondataset
        end
    end
    error("No data found for $fieldname $key, section plane $plane")
end

function query_insitu(key::AbstractString, x, y, z; plane="coronal", fieldname="acronym")
    sectiondataset = query_insitu(key; plane=plane, fieldname=fieldname)
    query_reference2image(sectiondataset, x, y, z)
end

default_insitu_imagename(id) = "id"*lpad(string(id), 4, "0")*".jpg"

function download_insitu_images(genekey::AbstractString, bb::BoundingBox, dirname=genekey; plane="coronal", fieldname="acronym", downsample=0)
    @assert plane == "coronal"  # not tested yet for sagittal
    sectiondataset = query_insitu(genekey; plane=plane, fieldname=fieldname)
    download_insitu_images(sectiondataset, bb, dirname; downsample=downsample)
end

function download_insitu_images(sectiondataset::Dict, bb::BoundingBox, dirname; downsample=0)
    if !isdir(dirname)
        mkpath(dirname)
    end
    ids = [simg["id"] for simg in sectiondataset["section_images"]]  # ones for this gene
    res = sectiondataset["section_images"][1]["resolution"]*μm
    xmin, ymin, zmin = map(minimum, bb.intervals)
    xmax, ymax, zmax = map(maximum, bb.intervals)
    q = query_reference2image(sectiondataset, xmin, ymin, zmin)
    qmax = query_reference2image(sectiondataset, xmin, ymax, zmax)
    w, h = qmax["x"] - q["x"], qmax["y"] - q["y"]
    @show w h (ymax-ymin)/res (zmax-zmin)/res
    filenames = String[]
    processed = Int[]
    slicepos = Vector{NTuple{3,Float64}}(0)
    Δx = 50μm
    for x in xmin:Δx:xmax
        q = query_reference2image(sectiondataset, x, ymin, zmin)
        id = q["section_image_id"]
        id ∈ ids || continue
        id ∈ processed && continue
        fn = default_insitu_imagename(id)
        push!(filenames, fn)
        sectionimage(q, w, h, joinpath(dirname, fn); downsample=downsample)
        push!(processed, id)
        # Record the real-world slice position. (The slice may not be where we asked.)
        push!(slicepos, query_image2reference(id, q["x"], q["y"]))
    end
    # Also save data we'll use later in reconstructing the 3d space
    FileIO.save(joinpath(dirname, "metadata.jld2"),
                Dict("filenames"=>filenames,
                     "sectionposition"=>slicepos,
                     "downsample"=>downsample,
                     "sectiondataset"=>sectiondataset,
                     "boundingbox"=>bb))
end

    # showall(filenames)
    # showall(slicepos)
    # sz = imagesize(joinpath(dirname, filenames[1]))
    # for i = 2:length(filenames)
    #     sz = map(max, sz, imagesize(joinpath(dirname, filenames[i])))
    # end
    # img2d = load(joinpath(dirname, first(filenames)))
    # img = similar(img2d, sz[2], sz[1], length(filenames))
    # fill!(img, zero(eltype(img)))
    # for i = 1:length(filenames)
    #     fn = joinpath(dirname, filenames[i])
    #     if imagesize(fn) == sz
    #         img[:,:,i] = load(fn)
    #     end
    # end
    # resp = median(diff(slicepos))
    # FileIO.save(joinpath(dirname, "merged.nrrd"), AxisArray(img, (:P, :I, :R), (resp, res, res)))

# function download_insitu_images(genekey::AbstractString, x, y, dirname=genekey; plane="coronal", fieldname="acronym", downsample=0)
#     if !isdir(dirname)
#         mkpath(dirname)
#     end
#     sectiondataset = query_insitu(key; plane=plane, fieldname=fieldname)

#     sections = sectiondataset["section_images"]
#     cd(dirname) do
#         @showprogress 1 "Downloading images for $key: " for i = 1:length(sections)
#             id = sections[i]["id"]
#             download("http://api.brain-map.org/api/v2/image_download/$id?downsample=$downsample", default_insitu_imagename(id))
#         end
#         open("metadata.json", "w") do io
#             write(io, JSON.json(sectiondataset))
#         end
#     end
#     # TODO: how to handle downsample in the metadata? relevant for the alignment
#     nothing
# end

function assemble_insitu_images(dirname)
    cd(dirname) do
        sectiondataset = open("metadata.json") do io
            JSON.parse(io)
        end
        imgsinfo = sectiondataset["section_images"]
        aff3 = parse_3daffine(sectiondataset["alignment3d"])
        n = length(imgsinfo)
        local inds, ds
        img = Array{RGB{N0f8}}(map(length, inds)..., n)
        for (i, info) in enumerate(imgsinfo)
            id = info["id"]
            img2d = load(default_insitu_imagename(id))
            aff2 = parse_2daffine(info["alignment2d"])
            if i == 1
                ds = downsampling(img2d, info)
                inds = (1:500, 1:500)  # fixme
            end
            aff2s = scaletform(aff2, ds)
            img[:,:,i] = warp(img2d, inv(aff2s), inds)
        end
        img
    end
end

function downsampling(img, info::Dict)
    h, w = info["height"], info["width"]
    ih, iw = size(img)
    ph, pw = round(Int, log2(h/ih)), round(Int, log2(w/iw))
    @assert ph == pw
    ph
end

function scaletform(tform2, ds)
    f = 2.0^ds
    M = diagm([f,f,1])
    tfm = M\([tform2; [0 0 1]] * M)
    LinearMap(SMatrix{2, 2, Float64}(tfm[1:2,1:2])) ∘ Translation(tfm[1:2,3]...)
end

function parse_2daffine(dct::Dict)
    # section-to-volume = "tsv"
    # http://api.brain-map.org/doc/Alignment2d.html
    tsv = zeros(2,3)
    for j = 1:2, i = 1:2
        tsv[i,j] = dct["tsv_0"*string(sub2ind((2,2), j, i)-1)]  # row major
    end
    tsv[1,3] = dct["tsv_04"]
    tsv[2,3] = dct["tsv_05"]
    tsv
end

function parse_3daffine(dct::Dict)
    tvr = zeros(3,4)
    for j = 1:4, i = 1:3
        tvr[i,j] = dct["trv_"*lpad(string(sub2ind((3,4), i, j)-1), 2, "0")]
    end
    tvr
end


text:
function simple_estimator(model::XY, T::Real, Js::AbstractArray, _=nothing)
    nsites = numsites(model)
    nbonds = numbonds(model)
    D = dim(model)
    invN = 1.0/nsites
    beta = 1.0/T

    M = zeros(2)
    E = 0.0
    U1 = zeros(2)
    U2 = zeros(2)

    @inbounds for s in 1:nsites
        M[1] += cospi(2model.spins[s])
        M[2] += sinpi(2model.spins[s])
    end

    @inbounds for b in bonds(model)
        i, j = source(b), target(b)
        dir = bonddirection(b)
        dt = mod(model.spins[j] - model.spins[i] + 2.0, 1.0)
        dt = ifelse(0.5 <= dt, dt - 1.0, dt)
        E -= cospi(2dt) * Js[bondtype(b)]

        for d in 1:D
            U1[d] += cospi(2dt) * dir[d]^2
            U2[d] += sinpi(2dt) * dir[d]
        end
    end
    for d in 1:2
        M[d] *= invN
        U1[d] -= beta * U2[d]^2
        U1[d] *= invN
    end
    E /= nsites

    res = Measurement()
    for (i,c) in enumerate(["x","y"])
        res["Magnetization $c"] = M[i]
        res["|Magnetization $c|"] = abs(M[i])
        res["Magnetization $c^2"] = M[i]^2
        res["Magnetization $c^4"] = M[i]^4
        res["Helicity Modulus $c"] = U1[i]
    end
    M2 = sum(abs2, M)
    res["|Magnetization|"] = sqrt(M2)
    res["|Magnetization|^2"] = M2
    res["|Magnetization|^4"] = M2^2
    res["Energy"] = E
    res["Energy^2"] = E^2

    return res
end


text:
using Test
using SumOfSquares
using DynamicPolynomials

function quartic_ideal_test(optimizer, config::MOIT.TestConfig,
                            degree::Union{Nothing, Int}, remainder::Bool)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @polyvar x
    # Set {-1, 0, 1}
    K = @set x^3 == x
    p = (x^2 - 1)^2
    cref = @constraint(model, p in SOSCone(), domain = K, maxdegree=degree, newton_of_remainder = remainder)
    optimize!(model)

    if remainder && (degree === nothing || degree < 4)
        @test termination_status(model) in [MOI.INFEASIBLE, MOI.INFEASIBLE_OR_UNBOUNDED]
    else
        @test termination_status(model) == MOI.OPTIMAL
        @test primal_status(model) == MOI.FEASIBLE_POINT
        μ = dual(cref)
        @test monomial(moments(μ)[1]) == x^4
        @test monomial(moments(μ)[2]) == x^2
        @test monomial(moments(μ)[3]) == 1
        μ = moments(cref)
        @test monomial(moments(μ)[1]) == x^2
        @test monomial(moments(μ)[2]) == x^1
        @test monomial(moments(μ)[3]) == 1
        @test moment_matrix(cref).basis.monomials == [x^2, x, 1]
    end
end
quartic_ideal_test(optimizer, config)   = quartic_ideal_test(optimizer, config, nothing, false)
sd_tests["quartic_ideal"] = quartic_ideal_test
quartic_ideal_4_test(optimizer, config) = quartic_ideal_test(optimizer, config, 4, false)
sd_tests["quartic_ideal_4"] = quartic_ideal_4_test
quartic_ideal_rem_test(optimizer, config)   = quartic_ideal_test(optimizer, config, nothing, true)
sd_tests["quartic_ideal_rem"] = quartic_ideal_rem_test
quartic_ideal_2_rem_test(optimizer, config) = quartic_ideal_test(optimizer, config, 2, true)
sd_tests["quartic_ideal_2_rem"] = quartic_ideal_2_rem_test
quartic_ideal_4_rem_test(optimizer, config) = quartic_ideal_test(optimizer, config, 4, true)
sd_tests["quartic_ideal_4_rem"] = quartic_ideal_4_rem_test


text:
@testset "Construction" begin end

@testset "Conversion" begin end

@testset "Evaluation" begin end

@testset "Arithmetic" begin end

@testset "Elementwise operations" begin end

@testset "Divrem" begin end

@testset "Comparisons" begin end

@testset "Derivatives and ntegrals" begin end

@testset "roots" begin end

@testset "fitting" begin end


text:
# =========================================================================== #
# Sample Container
# =========================================================================== #

abstract type AbstractSampleContainer
end

struct MatrixVectorSampleContainer{T} <: AbstractSampleContainer
    # each column is a new sample
    x::Matrix{T}
    y::Vector{T}
    function MatrixVectorSampleContainer{T}(x::Matrix{T}, y::Vector{T}) where {T}
        n_col = size(x, 2)
        L     = length(y)
        if n_col != L
            error("Invalid dimensions")
        end

        return new(x, y)
    end
end

function MatrixVectorSampleContainer(x::Matrix{T}, y::Vector{T}) where {T}
    return MatrixVectorSampleContainer{T}(x, y)
end

function Base.length(sample_cont::MatrixVectorSampleContainer)
    return size(sample_cont.x, 2)
end

function n_samples(sample_cont::MatrixVectorSampleContainer)
    return size(sample_cont.x, 2)
end

function n_dim(sample_cont::MatrixVectorSampleContainer)
    return size(sample_cont.x, 1)
end

function get_sample(sample_cont::MatrixVectorSampleContainer, idx::Int)
    nd = n_dim(sample_cont)
    x  = view(sample_cont.x, 1:nd, idx)
    y  = sample_cont.y[idx]

    return (x, y)
end

function split(sample_cont::MatrixVectorSampleContainer, idx::Int)

    ns = n_samples(sample_cont)

    x1  = sample_cont.x[:, 1:idx]
    y1  = sample_cont.y[1:idx]
    sc1 = MatrixVectorSampleContainer(x1, y1)

    x2  = sample_cont.x[:, (idx+1):ns]
    y2  = sample_cont.y[(idx+1):ns]
    sc2 = MatrixVectorSampleContainer(x2, y2)

    return (sc1, sc2)
end

function Random.shuffle!(sample_cont::MatrixVectorSampleContainer)
    N = length(sample_cont)
    xx = copy(sample_cont.x)
    yy = copy(sample_cont.y)
    idx = randperm(N)
    for ii = 1:N
        sample_cont.x[:,ii] = xx[:,idx[ii]]
        sample_cont.y[ii]   = yy[idx[ii]]
    end
end

struct MatrixMatrixSampleContainer{T} <: AbstractSampleContainer
    # each column is a new sample
    x::Matrix{T}
    y::Matrix{T}
    function MatrixMatrixSampleContainer{T}(x::Matrix{T}, y::Matrix{T}) where {T}
        n_col = size(x, 2)
        L     = size(y, 2)
        if n_col != L
            error("Invalid dimensions")
        end

        return new(x, y)
    end
end

function MatrixMatrixSampleContainer(x::Matrix{T}, y::Matrix{T}) where {T}
    return MatrixMatrixSampleContainer{T}(x, y)
end

function Base.length(sample_cont::MatrixMatrixSampleContainer)
    return size(sample_cont.x, 2)
end

function n_samples(sample_cont::MatrixMatrixSampleContainer)
    return size(sample_cont.x, 2)
end

function n_dim_x(sample_cont::MatrixMatrixSampleContainer)
    return size(sample_cont.x, 1)
end

function n_dim_y(sample_cont::MatrixMatrixSampleContainer)
    return size(sample_cont.y, 1)
end

function get_sample(sample_cont::MatrixMatrixSampleContainer, idx::Int)
    ndx = n_dim_x(sample_cont)
    ndy = n_dim_y(sample_cont)

    x  = view(sample_cont.x, 1:ndx, idx)
    y  = view(sample_cont.y, 1:ndy, idx)

    return (x, y)
end

function split(sample_cont::MatrixMatrixSampleContainer, idx::Int)
    ns = n_samples(sample_cont)

    x1  = sample_cont.x[:, 1:idx]
    y1  = sample_cont.y[:, 1:idx]
    sc1 = MatrixMatrixSampleContainer(x1, y1)

    x2  = sample_cont.x[:, (idx+1):ns]
    y2  = sample_cont.y[:, (idx+1):ns]
    sc2 = MatrixMatrixSampleContainer(x2, y2)

    return (sc1, sc2)
end

function Random.shuffle!(sample_cont::MatrixMatrixSampleContainer)
    N = length(sample_cont)
    xx = copy(sample_cont.x)
    yy = copy(sample_cont.y)
    idx = randperm(N)
    for ii = 1:N
        sample_cont.x[:, ii] = xx[:, idx[ii]]
        sample_cont.y[:, ii] = yy[:, idx[ii]]
    end
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#! format: off
############################### Reserve Variables #########################################

get_variable_multiplier(_, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = NaN
############################### ActivePowerReserveVariable, Reserve #########################################

get_variable_binary(::ActivePowerReserveVariable, ::Type{<:PSY.Reserve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ActivePowerReserveVariable, ::PSY.Reserve, d::PSY.Component, _) = PSY.get_max_active_power(d)
get_variable_upper_bound(::ActivePowerReserveVariable, ::PSY.Reserve, d::PSY.Storage, _) =  PSY.get_output_active_power_limits(d).max
get_variable_lower_bound(::ActivePowerReserveVariable, ::PSY.Reserve, ::PSY.Component, _) = 0.0

############################### ActivePowerReserveVariable, ReserveNonSpinning #########################################

get_variable_binary(::ActivePowerReserveVariable, ::Type{<:PSY.ReserveNonSpinning}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ActivePowerReserveVariable, ::PSY.ReserveNonSpinning, d::PSY.Component, _) = PSY.get_max_active_power(d)
get_variable_upper_bound(::ActivePowerReserveVariable, ::PSY.ReserveNonSpinning, d::PSY.Storage, _) =  PSY.get_output_active_power_limits(d).max
get_variable_lower_bound(::ActivePowerReserveVariable, ::PSY.ReserveNonSpinning, ::PSY.Component, _) = 0.0

############################### ServiceRequirementVariable, ReserveDemandCurve ################################

get_variable_binary(::ServiceRequirementVariable, ::Type{<:PSY.ReserveDemandCurve}, ::AbstractReservesFormulation) = false
get_variable_upper_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, d::PSY.Component, ::AbstractReservesFormulation) = PSY.get_max_active_power(d)
get_variable_lower_bound(::ServiceRequirementVariable, ::PSY.ReserveDemandCurve, ::PSY.Component, ::AbstractReservesFormulation) = 0.0

get_multiplier_value(::RequirementTimeSeriesParameter, d::PSY.Reserve, ::AbstractReservesFormulation) = PSY.get_requirement(d)
get_multiplier_value(::RequirementTimeSeriesParameter, d::PSY.ReserveNonSpinning, ::AbstractReservesFormulation) = PSY.get_requirement(d)

get_parameter_multiplier(::VariableValueParameter, d::Type{<:PSY.AbstractReserve}, ::AbstractReservesFormulation) = 1.0
get_initial_parameter_value(::VariableValueParameter, d::Type{<:PSY.AbstractReserve}, ::AbstractReservesFormulation) = 1.0

objective_function_multiplier(::ServiceRequirementVariable, ::StepwiseCostReserve) = 1.0
sos_status(::PSY.ReserveDemandCurve, ::StepwiseCostReserve)=SOSStatusVariable.NO_VARIABLE
uses_compact_power(::PSY.ReserveDemandCurve, ::StepwiseCostReserve)=false
#! format: on

function get_default_time_series_names(
    ::Type{<:PSY.Reserve},
    ::Type{T},
) where {T <: Union{RangeReserve, RampReserve}}
    return Dict{Type{<:TimeSeriesParameter}, String}(
        RequirementTimeSeriesParameter => "requirement",
    )
end

function get_default_time_series_names(
    ::Type{<:PSY.ReserveNonSpinning},
    ::Type{NonSpinningReserve},
)
    return Dict{Type{<:TimeSeriesParameter}, String}(
        RequirementTimeSeriesParameter => "requirement",
    )
end

function get_default_time_series_names(
    ::Type{<:PSY.Service},
    ::Type{<:AbstractServiceFormulation},
)
    return Dict{Type{<:TimeSeriesParameter}, String}()
end

function get_default_attributes(::Type{<:PSY.Service}, ::Type{<:AbstractServiceFormulation})
    return Dict{String, Any}()
end

################################## Reserve Requirement Constraint ##########################
function add_constraints!(
    container::OptimizationContainer,
    T::Type{RequirementConstraint},
    service::SR,
    ::U,
    model::ServiceModel{SR, V},
) where {
    SR <: PSY.AbstractReserve,
    V <: AbstractReservesFormulation,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    parameters = built_for_recurrent_solves(container)
    initial_time = get_initial_time(container)

    time_steps = get_time_steps(container)
    service_name = PSY.get_name(service)
    # TODO: Add a method for services that handles this better
    constraint = add_constraints_container!(
        container,
        T(),
        SR,
        [service_name],
        time_steps;
        meta=service_name,
    )
    reserve_variable =
        get_variable(container, ActivePowerReserveVariable(), SR, service_name)
    use_slacks = get_use_slacks(model)

    ts_vector = get_time_series(container, service, "requirement")

    use_slacks && (slack_vars = reserve_slacks(container, service))
    requirement = PSY.get_requirement(service)
    jump_model = get_jump_model(container)
    if parameters
        container =
            get_parameter(container, RequirementTimeSeriesParameter(), SR, service_name)
        param = get_parameter_array(container)
        for t in time_steps
            if use_slacks
                resource_expression = sum(reserve_variable[:, t]) + slack_vars[t]
            else
                resource_expression = sum(reserve_variable[:, t])
            end
            constraint[service_name, t] = JuMP.@constraint(
                jump_model,
                resource_expression >= param[service_name, t] * requirement
            )
        end
    else
        for t in time_steps
            if use_slacks
                resource_expression = sum(reserve_variable[:, t]) + slack_vars[t]
            else
                resource_expression = sum(reserve_variable[:, t])
            end
            constraint[service_name, t] = JuMP.@constraint(
                jump_model,
                resource_expression >= ts_vector[t] * requirement
            )
        end
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{RequirementConstraint},
    service::SR,
    ::U,
    model::ServiceModel{SR, V},
) where {
    SR <: PSY.StaticReserve,
    V <: AbstractReservesFormulation,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    time_steps = get_time_steps(container)
    service_name = PSY.get_name(service)
    # TODO: The constraint addition is still not clean enough
    constraint = add_constraints_container!(
        container,
        T(),
        SR,
        [service_name],
        time_steps;
        meta=service_name,
    )
    reserve_variable =
        get_variable(container, ActivePowerReserveVariable(), SR, service_name)
    use_slacks = get_use_slacks(model)
    use_slacks && (slack_vars = reserve_slacks(container, service))

    requirement = PSY.get_requirement(service)
    jump_model = get_jump_model(container)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[service_name, t] =
            JuMP.@constraint(jump_model, resource_expression >= requirement)
    end

    return
end

function objective_function!(
    container::OptimizationContainer,
    service::SR,
    ::ServiceModel{SR, T},
) where {SR <: PSY.AbstractReserve, T <: AbstractReservesFormulation}
    add_proportional_cost!(container, ActivePowerReserveVariable(), service, T())
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{RequirementConstraint},
    service::SR,
    ::U,
    ::ServiceModel{SR, StepwiseCostReserve},
) where {
    SR <: PSY.ReserveDemandCurve,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    time_steps = get_time_steps(container)
    service_name = PSY.get_name(service)
    constraint = add_constraints_container!(
        container,
        T(),
        SR,
        [service_name],
        time_steps;
        meta=service_name,
    )
    reserve_variable =
        get_variable(container, ActivePowerReserveVariable(), SR, service_name)
    requirement_variable = get_variable(container, ServiceRequirementVariable(), SR)
    jump_model = get_jump_model(container)
    for t in time_steps
        constraint[service_name, t] = JuMP.@constraint(
            jump_model,
            sum(reserve_variable[:, t]) >= requirement_variable[service_name, t]
        )
    end

    return
end

_get_ramp_limits(::PSY.Component) = nothing
_get_ramp_limits(d::PSY.ThermalGen) = PSY.get_ramp_limits(d)
_get_ramp_limits(d::PSY.HydroGen) = PSY.get_ramp_limits(d)

function _get_ramp_constraint_contributing_devices(
    service::PSY.Reserve,
    contributing_devices::Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
) where {D <: PSY.Component}
    time_frame = PSY.get_time_frame(service)
    filtered_device = Vector{D}()
    for d in contributing_devices
        ramp_limits = _get_ramp_limits(d)
        if ramp_limits !== nothing
            p_lims = PSY.get_active_power_limits(d)
            max_rate = abs(p_lims.min - p_lims.max) / time_frame
            if (ramp_limits.up >= max_rate) & (ramp_limits.down >= max_rate)
                @debug "Generator $(name) has a nonbinding ramp limits. Constraints Skipped"
                continue
            else
                push!(filtered_device, d)
            end
        end
    end
    return filtered_device
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{RampConstraint},
    service::SR,
    contributing_devices::Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    ::ServiceModel{SR, V},
) where {
    SR <: PSY.Reserve{PSY.ReserveUp},
    V <: AbstractReservesFormulation,
    D <: PSY.Component,
}
    ramp_devices = _get_ramp_constraint_contributing_devices(service, contributing_devices)
    service_name = PSY.get_name(service)
    if !isempty(ramp_devices)
        jump_model = get_jump_model(container)
        time_steps = get_time_steps(container)
        time_frame = PSY.get_time_frame(service)
        variable = get_variable(container, ActivePowerReserveVariable(), SR, service_name)
        set_name = [PSY.get_name(d) for d in ramp_devices]
        con_up = add_constraints_container!(
            container,
            T(),
            SR,
            set_name,
            time_steps;
            meta=service_name,
        )
        for d in ramp_devices, t in time_steps
            name = PSY.get_name(d)
            ramp_limits = PSY.get_ramp_limits(d)
            con_up[name, t] = JuMP.@constraint(
                jump_model,
                variable[name, t] <= ramp_limits.up * time_frame
            )
        end
    else
        @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation"
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{RampConstraint},
    service::SR,
    contributing_devices::Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    ::ServiceModel{SR, V},
) where {
    SR <: PSY.Reserve{PSY.ReserveDown},
    V <: AbstractReservesFormulation,
    D <: PSY.Component,
}
    ramp_devices = _get_ramp_constraint_contributing_devices(service, contributing_devices)
    service_name = PSY.get_name(service)
    if !isempty(ramp_devices)
        jump_model = get_jump_model(container)
        time_steps = get_time_steps(container)
        time_frame = PSY.get_time_frame(service)
        variable = get_variable(container, ActivePowerReserveVariable(), SR, service_name)
        set_name = [PSY.get_name(d) for d in ramp_devices]
        con_down = add_constraints_container!(
            container,
            T(),
            SR,
            set_name,
            time_steps;
            meta=service_name,
        )
        for d in ramp_devices, t in time_steps
            name = PSY.get_name(d)
            ramp_limits = PSY.get_ramp_limits(d)
            con_down[name, t] = JuMP.@constraint(
                jump_model,
                variable[name, t] <= ramp_limits.down * time_frame
            )
        end
    else
        @warn "Data doesn't contain contributing devices with ramp limits for service $service_name, consider adjusting your formulation"
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    T::Type{ReservePowerConstraint},
    service::SR,
    contributing_devices::U,
    ::ServiceModel{SR, V},
) where {
    SR <: PSY.VariableReserveNonSpinning,
    V <: AbstractReservesFormulation,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
} where {D <: PSY.Component}
    time_steps = get_time_steps(container)
    resolution = get_resolution(container)
    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        @warn("Not all formulations support under 1-minute resolutions. Exercise caution.")
        minutes_per_period = Dates.value(Dates.Second(resolution)) / 60
    end
    service_name = PSY.get_name(service)
    cons = add_constraints_container!(
        container,
        T(),
        SR,
        [PSY.get_name(d) for d in contributing_devices],
        time_steps;
        meta=service_name,
    )
    var_r = get_variable(container, ActivePowerReserveVariable(), SR, service_name)
    reserve_response_time = PSY.get_time_frame(service)
    jump_model = get_jump_model(container)
    for d in contributing_devices
        component_type = typeof(d)
        name = PSY.get_name(d)
        varstatus = get_variable(container, OnVariable(), component_type)
        startup_time = PSY.get_time_limits(d).up
        ramp_limits = _get_ramp_limits(d)
        if reserve_response_time > startup_time
            reserve_limit =
                PSY.get_active_power_limits(d).min +
                (reserve_response_time - startup_time) * minutes_per_period * ramp_limits.up
        else
            reserve_limit = 0
        end
        for t in time_steps
            cons[name, t] = JuMP.@constraint(
                jump_model,
                var_r[name, t] <= (1 - varstatus[name, t]) * reserve_limit
            )
        end
    end
    return
end

function objective_function!(
    container::OptimizationContainer,
    service::PSY.ReserveDemandCurve{T},
    ::ServiceModel{PSY.ReserveDemandCurve{T}, SR},
) where {T <: PSY.ReserveDirection, SR <: StepwiseCostReserve}
    add_variable_cost!(container, ServiceRequirementVariable(), service, SR())
    return
end

function add_services_to_device_model!(template)
    service_models = get_service_models(template)
    devices_template = get_device_models(template)
    for (service_key, service_model) in service_models
        S = get_component_type(service_model)
        (S <: PSY.AGC || S <: PSY.StaticReserveGroup) && continue
        contributing_devices = get_contributing_devices(service_model)
        isempty(contributing_devices) && continue
        modify_device_model!(devices_template, service_model, contributing_devices)
    end
    return
end

function modify_device_model!(
    devices_template::Dict{Symbol, DeviceModel},
    service_model::ServiceModel{<:PSY.Reserve, <:AbstractReservesFormulation},
    contributing_devices::Vector{<:PSY.Component},
)
    for dt in unique(typeof.(contributing_devices))
        for device_model in values(devices_template)
            # add message here when it exists
            get_component_type(device_model) != dt && continue
            service_model in device_model.services && continue
            push!(device_model.services, service_model)
        end
    end

    return
end

function modify_device_model!(
    devices_template::Dict{Symbol, DeviceModel},
    service_model::ServiceModel{<:PSY.ReserveNonSpinning, <:AbstractReservesFormulation},
    contributing_devices::Vector{<:PSY.Component},
)
    return
end


text:
include("etcdserver.jl")

lock_prefix(etcd::EtcdServer,
            lock::String) = "http://$(etcd.ip):$(etcd.port)/mod/$(etcd.version)/lock/$(lock)"

# blocking until lock is available
# TODO should add an async lock
function lock_acquire(etcd::EtcdServer,lock::String;ttl::Union(Int,Nothing)=nothing)
    if is(ttl,nothing)
        warn("Must specify an integer value for lock ttl")
    else
        etcd_request(:post,lock_prefix(etcd,lock),{"ttl"=>ttl}) |>
        check_etcd_error
    end
end

function lock_renew(etcd::EtcdServer,lock::String;
                    index::Union(Int,Nothing)=nothing,
                    ttl::Union(Int,Nothing)=nothing)

    if is(index,nothing) && is(ttl,nothing)
        warn("Must specify an integer value for lock index and ttl")
    else
        etcd_request(:put,lock_prefix(etcd,lock),{"index"=>index,
                                                          "ttl"=>ttl}) |>
        check_etcd_error
    end
end

function lock_release(etcd::EtcdServer,lock::String;
                      index::Union(Int,Nothing)=nothing)
    if is(index,nothing)
        warn("Must specify an integer value for lock index")
    else
        etcd_request(:delete,lock_prefix(etcd,lock),{"index"=>index}) |>
        check_etcd_error
    end
end

function lock_retrieve(etcd::EtcdServer,lock::String)
    etcd_request(:get,lock_prefix(etcd,lock),{"field"=>"index"}) |>
    check_etcd_error
end


text:
# This file was generated, do not modify it. # hide
DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree

dcrm = machine(DecisionTreeRegressor(), X, y)

fit!(dcrm, rows=train)
pred_dcrm = predict(dcrm, rows=test);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module mytest
    export greet
    greet() = print("Hello World!")

end # module


text:
function readinput(filename)
    map(value -> parse(Int, value), readlines(filename))
end

function sumexist(numbers, number)
    i = 1
    for x in numbers[1:end-1]
        i += 1
        for y in numbers[i:end]
            if x + y == number
                return true
            end
        end
    end
    return false
end

function findinvalid(values, preamble = 25)
    queue = []
    for value in values[1:preamble]
        push!(queue, value)
    end
    for value in values[preamble+1:end]
        if !sumexist(queue, value)
            return value
        end
        pop = popfirst!(queue)
        push!(queue, value)
    end
    return -1
end

function findsum(values, invalid)
    n = 1
    for x in values[1:end-1]
        group = []
        n += 1
        push!(group, x)
        for y in values[n:end]
            push!(group, y)
            if sum(group) == P1
                return deepcopy(group)
            end
        end
    end
    return []
end

values = readinput("d09.in")
P1 = findinvalid(values)
sumgroup = findsum(values, P1)
P2 = minimum(sumgroup) + maximum(sumgroup)

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
# Included from build.jl and ../test/test_build.jl

using VersionParsing
import Conda, Libdl

pyvar(python::AbstractString, mod::AbstractString, var::AbstractString) = chomp(read(pythonenv(`$python -c "import $mod; print($mod.$var)"`), String))

pyconfigvar(python::AbstractString, var::AbstractString) = pyvar(python, "distutils.sysconfig", "get_config_var('$var')")
pyconfigvar(python, var, default) = let v = pyconfigvar(python, var)
    v == "None" ? default : v
end

pysys(python::AbstractString, var::AbstractString) = pyvar(python, "sys", var)

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

# print out extra info to help with remote debugging
const PYCALL_DEBUG_BUILD = "yes" == get(ENV, "PYCALL_DEBUG_BUILD", "no")

function exec_find_libpython(python::AbstractString, options)
    # Do not inline `@__DIR__` into the backticks to expand correctly.
    # See: https://github.com/JuliaLang/julia/issues/26323
    script = joinpath(@__DIR__, "find_libpython.py")
    cmd = `$python $script $options`
    if PYCALL_DEBUG_BUILD
        cmd = `$cmd --verbose`
    end
    return readlines(pythonenv(cmd))
end

function show_dlopen_error(lib, e)
    if PYCALL_DEBUG_BUILD
        println(stderr, "dlopen($lib) ==> ", e)
        # Using STDERR since find_libpython.py prints debugging
        # messages to STDERR too.
    end
end

# return libpython name, libpython pointer
function find_libpython(python::AbstractString; _dlopen = Libdl.dlopen)
    dlopen_flags = Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL

    libpaths = exec_find_libpython(python, `--list-all`)
    for lib in libpaths
        try
            return (_dlopen(lib, dlopen_flags), lib)
        catch e
            show_dlopen_error(lib, e)
        end
    end

    # Try all candidate libpython names and let Libdl find the path.
    # We do this *last* because the libpython in the system
    # library path might be the wrong one if multiple python
    # versions are installed (we prefer the one in LIBDIR):
    libs = exec_find_libpython(python, `--candidate-names`)
    for lib in libs
        lib = splitext(lib)[1]
        try
            libpython = _dlopen(lib, dlopen_flags)
            # Store the fullpath to libpython in deps.jl.  This makes
            # it easier for users to investigate Python setup
            # PyCall.jl trying to use.  It also helps PyJulia to
            # compare libpython.
            return (libpython, Libdl.dlpath(libpython))
        catch e
            show_dlopen_error(lib, e)
        end
    end

    v = pyconfigvar(python, "VERSION", "unknown")
    error("""
        Couldn't find libpython; check your PYTHON environment variable.

        The python executable we tried was $python (= version $v).
        Re-building with
            ENV["PYCALL_DEBUG_BUILD"] = "yes"
        may provide extra information for why it failed.
        """)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 41e660a2-0d51-4a05-aacd-a78728d77b1c
md"""
### Variables and Some Basic Types

In Julia, you can define a variable similar to how you define it in Python, e.g we can define a `x` using `=` (assignment)
"""

# ╔═╡ c4f6b254-07e4-4ed7-8870-452818cb6e67
x = 1

# ╔═╡ 3b1e55f8-1616-462f-8996-6c49b92b1f0a
md"""every variable has a type, you can check it using `typeof`"""

# ╔═╡ 8742c25c-fcea-4799-a91c-04094f5d0e0b
typeof(x)

# ╔═╡ 662685dd-96e5-4a18-ac48-9c3056926a5c
md"By default Julia displays the output of the last operation. (You can suppress the output by adding `;` (a semicolon) at the end.)"

# ╔═╡ 905f2ff2-29fe-4c65-b955-f1a6675ff965
md"""
### Functions

In Julia, you can also define short-form, one-line functions using `=` (assignment) similar to how you write things mathematically.
"""

# ╔═╡ 8569be06-1f4a-47b2-aa64-b65027a3d2b8
f(x) = 2x

# ╔═╡ ee81eb70-168b-4b98-8a89-bcb8d2a5eb38
md"Typing the function's name gives information about the function. To call it we must use parentheses:"

# ╔═╡ b6bdc0df-54a2-4c9a-a4e2-e3aee8fdcd05
f

# ╔═╡ a72f47f4-0c75-4433-8afd-468609d0bbaf
f(2)

# ╔═╡ 34f81e8d-7714-40ea-b546-2383d7813eed
md"For longer functions we use the following syntax with the `function` keyword and `end`:"

# ╔═╡ b5b567da-ec59-4a9e-a70e-a226a828eb50
function g(x, y)
	z = x + y
	return z^2
end

# ╔═╡ 04f8a8d5-1bbb-42df-a0d0-ed80aa0ce957
md"""
### Control Flows

In Julia, there are `for`, `if` and `while`, they look like the following
"""

# ╔═╡ 10cd61db-09b9-4876-a3c9-2114578da639
begin
	s = 0
	for i in 1:10
		s += 1
	end
end

# ╔═╡ a1aa126f-15f3-44a2-8acb-ebccd8f16cfd
md"here `begin...end` is used to wrap a code block as Pluto required, you don't have to write it outside Pluto, we can now check the value of `s` by typing it again"

# ╔═╡ 9c002ca1-28ac-451c-ada2-8e55b90d3e64
s

# ╔═╡ 6154997e-a121-430e-bafe-c794412d9478
md"Here, `1:10` is a **range** representing the numbers from 1 to 10:"

# ╔═╡ 2ae4588f-6b88-47b0-928a-36e6df4bd755
typeof(1:10)

# ╔═╡ 201c7fd9-5723-46ce-84fa-7cb400a02f7b
md"""
the if else statement looks like the following
"""

# ╔═╡ 766d92dd-3df7-43f5-853c-43c4194aa842
if s < 10
	# do something
elseif 10 < s < 13
	# do something
else
	# do something
end

# ╔═╡ c7a414ca-a016-4959-8e06-080907ff244c
md"""
### Matrix and Array

Julia carries its own `Array` type, if you use Python, it is similar to `numpy.array` in Python except:

1. index starts from 1
2. the multi-dimensional index is column-wise

You can also have list comprehension:
"""

# ╔═╡ 382f7095-a3fe-4e99-9a2f-1b745f41830f
[i for i in 1:10]

# ╔═╡ e19c3b51-8405-4c47-9e7d-df809d281b9b
md"it works for multi-dimensional case too:"

# ╔═╡ d18227df-56d5-4879-ac26-9f45d35ba144
[(i, j) for i in 1:10, j in 1:5]

# ╔═╡ 4537e253-4954-4d2c-bc33-18aee3a2bdf5
md"most functions follow the same convention as numpy or MATLAB, e.g you can create a random matrix using:"

# ╔═╡ 301a8694-0b01-4515-ab40-d0f8564cfc78
rand(5, 5)

# ╔═╡ 749bf936-eefc-4dd5-b495-939d82f496a0
md"""
if you have question about using a function, you can always type question mark `?` in your REPL following the function name

```julia
julia> ?rand
```
"""

# ╔═╡ 0571129d-bdbf-4ff4-b607-dacb6d6d3911
md"""
### Package Manager & Environments

Julia carries its own package manager, you can use it as a normal package:

```julia
julia> using Pkg
```

to install a pacakge, you can use

```julia
julia> Pkg.add("Yao")
```

to remove a pacakge, you can use

```julia
julia> Pkg.rm("Yao")
```

All Julia program **runs inside an environment**, it is the global environment by default. It is usually recommended to run your notebook in a local environment, so we won't hit **any version conflicts** between different packages.
"""

# ╔═╡ Cell order:
# ╟─41e660a2-0d51-4a05-aacd-a78728d77b1c
# ╠═c4f6b254-07e4-4ed7-8870-452818cb6e67
# ╟─3b1e55f8-1616-462f-8996-6c49b92b1f0a
# ╠═8742c25c-fcea-4799-a91c-04094f5d0e0b
# ╟─662685dd-96e5-4a18-ac48-9c3056926a5c
# ╟─905f2ff2-29fe-4c65-b955-f1a6675ff965
# ╠═8569be06-1f4a-47b2-aa64-b65027a3d2b8
# ╟─ee81eb70-168b-4b98-8a89-bcb8d2a5eb38
# ╠═b6bdc0df-54a2-4c9a-a4e2-e3aee8fdcd05
# ╠═a72f47f4-0c75-4433-8afd-468609d0bbaf
# ╟─34f81e8d-7714-40ea-b546-2383d7813eed
# ╠═b5b567da-ec59-4a9e-a70e-a226a828eb50
# ╟─04f8a8d5-1bbb-42df-a0d0-ed80aa0ce957
# ╠═10cd61db-09b9-4876-a3c9-2114578da639
# ╟─a1aa126f-15f3-44a2-8acb-ebccd8f16cfd
# ╠═9c002ca1-28ac-451c-ada2-8e55b90d3e64
# ╟─6154997e-a121-430e-bafe-c794412d9478
# ╠═2ae4588f-6b88-47b0-928a-36e6df4bd755
# ╟─201c7fd9-5723-46ce-84fa-7cb400a02f7b
# ╠═766d92dd-3df7-43f5-853c-43c4194aa842
# ╟─c7a414ca-a016-4959-8e06-080907ff244c
# ╠═382f7095-a3fe-4e99-9a2f-1b745f41830f
# ╟─e19c3b51-8405-4c47-9e7d-df809d281b9b
# ╠═d18227df-56d5-4879-ac26-9f45d35ba144
# ╟─4537e253-4954-4d2c-bc33-18aee3a2bdf5
# ╠═301a8694-0b01-4515-ab40-d0f8564cfc78
# ╟─749bf936-eefc-4dd5-b495-939d82f496a0
# ╟─0571129d-bdbf-4ff4-b607-dacb6d6d3911


text:
using Documenter, ModelPredictiveControl

makedocs(
    modules=[ModelPredictiveControl],
    format = Documenter.HTML(),
    sitename="ModelPredictiveControl.jl",

    pages=[
        "Home" => "index.md",
        "Problems" => Any["linearmpc.md"],
        "Tools" => Any["tools.md"],
        "Examples" => Any["Examples/QTP.md",
                          "Examples/mass-string-damper.md"],
        "Functions" => Any[
                        "functions.md"],
        "Style guide" => Any["style_guide.md"],

    ] )


text:
using SpecialFunctions: gamma
using PyCall
np = pyimport("numpy")
plt = pyimport("matplotlib.pyplot")
tri = pyimport("matplotlib.tri")

function tri_area(xy, pair)
    """The area of the triangle formed by point xy and another pair or points."""
    return 0.5 * np.linalg.norm(np.cross(pair[1,:] .- xy, pair[2,:] .- xy))
end

function xy2bc(xy, pair_corners; tol=1.e-4, AREA=0.5*1*0.75^0.5)
    """Converts 2D Cartesian coordinates to barycentric."""
    coords = np.array([tri_area(xy, p) for p in pair_corners]) / AREA
    return np.clip(coords, tol, 1.0 - tol)
end

function pdf_Dir(θ, α)
    """Plots pdf of Dirichlet distribution"""
    return prod(gamma.(α))/gamma(prod(α)) * prod(θ.^α)
end

function pdf_contours_simplex(concentration_params; nlevels=200, subdiv=8)
    """Generate pdf contours."""

    corners = np.array([[0, 0], [1, 0], [0.5, 0.75^0.5]])
    AREA = 0.5 * 1 * 0.75^0.5
    triangle = tri.Triangulation(corners[:, 1], corners[:, 2])

    pair_corners = [corners[np.roll(range(1, length=3), -i)[2:end],:] for i in range(0, length=3)]

    refiner = tri.UniformTriRefiner(triangle)
    trimesh = refiner.refine_triangulation(subdiv=subdiv)
    pvals = [pdf_Dir(xy2bc(xy, pair_corners, AREA=AREA), concentration_params) for xy in zip(trimesh.x, trimesh.y)]

    return trimesh, pvals
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
module ShaleDrillingLikelihood_check_psi

using ShaleDrillingLikelihood
using Test
using SparseArrays
using Random
using StatsFuns

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: DCDPEmax,
    ObservationDrill,
    DCDPTmpVars,
    flow!,
    end_ex0,
    reward,
    revenue,
    _dψ1dρ,
    _dψ1dθρ,
    _ρ,
    dp1space,
    sprime,
    SimulationDraws,
    SimulationDraw,
    _psi1, _psi2, _ψ1, _ψ2, _ψ,
    end_ex0,
    end_lrn,
    _D,
    ssprime,
    _dmax,
    end_inf,
    sprime,
    end_ex1,
    UnconstrainedProblem,
    vw_cost,
    vw_extend,
    vw_revenue,
    _sgnext

println("print to keep from blowing up")

@testset "Selection of psi in DDM" begin

    # set up coefs
    θρ = 0.0
    αψ = 1.0 # 0.33
    αg = 0.56

    θ_drill_u = [-6.5, -0.85, -2.8, αg, αψ, θρ]
    θ_drill_c = vcat(θ_drill_u[1:3], θρ)

    Random.seed!(1234)
    f = DrillReward(DrillingRevenue(Constrained(),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant())
    fu = UnconstrainedProblem(f)
    theta = randn(_nparm(f))
    # @show theta

    @test vw_cost(fu, θ_drill_u) == θ_drill_u[1:1]
    @test vw_extend(fu, θ_drill_u) == θ_drill_u[2:2]
    @test vw_revenue(fu, θ_drill_u) == θ_drill_u[3:end]

    nψ, dmx, nz =  13, 3, 11

    # psi
    ψs = range(-4.5; stop=4.5, length=nψ)

    # z transition
    zs = ( range(-1.5; stop=1.5, length=nz), )
    nz = prod(length.(zs))
    ztrans = spdiagm(-1 => fill(1.0/3, nz-1), 0 => fill(1.0/3, nz), 1 => fill(1.0/3, nz-1) )
    ztrans[1,1] = ztrans[1,2] = ztrans[end,end-1] = ztrans[end,end] = 0.5

    # wp = LeasedProblemContsDrill(dmx,4,5,3,2)
    # wp = LeasedProblem(dmx,4,5,3,2)
    wp = LeasedProblem(8, 8, 30, 20, 8)

    # ddm object
    ddm = DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, false)
    u,v = (0.3, -1.5)
    θρ = 0.5
    rho = _ρ(θρ)

    _ichar = (4.7, 0.2)
    _z = (0.5,)

    @test rho == logistic(θρ)
    sim = SimulationDraw(u,v, θρ)
    @test _psi1(sim) == _ψ1(sim)
    @test _psi2(sim) == _ψ2(sim)
    @test u == _psi2(sim) == _ψ2(u,v)
    @test (rho*u + sqrt(1-rho^2)*v) == _psi1(sim) == _ψ1(sim)
    @test _psi1(SimulationDraw(u,v,Inf)) == u
    @test _psi1(SimulationDraw(u,v,-Inf)) == v

    for sidx in 1:length(wp)
        D = _D(wp,sidx)
        obs = ObservationDrill(ddm, _ichar, _z, 0, sidx)

        if sidx <= end_ex0(wp)+1
            @test D == 0
            @test _ψ(ddm, sidx, sim) == _ψ1(sim)
            @test _ψ(obs,sim) == _ψ1(sim)
            @test _ψ(obs,sim) != u
            @test _ψ(obs,sim) != v

            if sidx == end_ex1(wp)
                @test _sgnext(wp,sidx) == true
                @test _sgnext(wp,sidx,0) == true
                @test _sgnext(wp,sidx,1) == false
            else
                @test _sgnext(wp,sidx) == false
                @test _sgnext(wp,sidx,0) == false
                @test _sgnext(wp,sidx,1) == false
            end

            if sidx <= end_ex0(wp)
                @test _dmax(wp,sidx) > 0
                @test ssprime(wp, sidx, 1) > end_lrn(wp)
                @test _D(wp,ssprime(wp, sidx, 1)) == 1
            else
                @test _dmax(wp,sidx) == 0
            end
        else
            @test D > 0
            if sidx > end_lrn(wp)
                @test _ψ(obs,sim) == _ψ2(sim)
                @test _ψ(obs,sim) == u
                @test _ψ(ddm, sidx, sim) == _ψ2(sim)
                @test _sgnext(wp,sidx) == false
                if sidx < end_inf(wp)
                    sp = ssprime(wp,sidx,1)
                    @test _dmax(wp,sidx) > 0
                    @test end_inf(wp) >= sp > sidx > end_lrn(wp)
                    @test sp == sprime(wp,sidx,1)
                    @test 1 <= _D(wp,sidx)
                    @test _D(wp,sidx)+1 == _D(wp,sp)
                else
                    @test _dmax(wp,sidx) == 0
                end
            end
        end

    end

end # testset

end # module


text:
# Particle: variant S

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test

include("particle_test_utils.jl")

@bench @testset "NLEIGS: Particle variant S" begin
    verbose = displaylevel

    nep, Σ, Ξ, v, nodes, xmin, xmax = particle_init(2)

    # solve nlep
    lambda, X, res, solution_info = nleigs(nep, Σ, Ξ=Ξ, logger=verbose > 0 ? 1 : 0, maxdgr=50, minit=120, maxit=200, v=v, nodes=nodes, static=true, return_details=verbose > 1)

    verify_lambdas(2, nep, lambda, X)

    if verbose > 1
        include("nleigs_residual_plot.jl")
        approx_Σ = [xmin-im*1e-10, xmin+im*1e-10, xmax+im*1e-10, xmax-im*1e-10]
        nleigs_residual_plot("Particle: variant S", solution_info, approx_Σ)
    end
end


text:
# Returns the gradient of x, computed by finited differences
#
# Except for the first and last element, all gradients are computed by central
# finite differences
function fdgrad(x::AbstractVector{T}, dt::Real) where T
    n = length(x)
    GType = typeof((zero(T) - zero(T)) / dt)
    g = Array{GType}(undef, n)
    if n == 1
        g[1] = zero(GType)
    else
        g[1] = (x[2] - x[1]) / dt
        for i = 2:(n-1)
            g[i] = (x[i+1] - x[i-1]) / 2dt
        end
        g[n] = (x[n] - x[n-1]) / dt
    end
    return g
end


text:
"""
The EUROCONTROL Demand Data Repository 'DDR2' provides air traffic management
(ATM) actors with the most accurate picture of past and future pan-European air
traffic demand (from several years ahead until the day before operation), as
well as environment data, analysis reports and tools that can read and process
the data.

All this information is management by DDR service, a EUROCONTROL cross-unit
activity, and can be accessed from the Demand Data Repository 2 'DDR2' web
portal. Access to the DDR2 web portal is restricted. Access conditions apply.

DDR future traffic can be forecast thanks to the knowledge of past traffic and
several thousands of flight intentions provided by airlines and airports that
are collected, stored, analysed and treated on a daily basis.

DDR traffic forecast supports strategic, seasonal and pre-tactical planning,
and also special events or major ATM evolution projects.

Finally, DDR provides a refined analysis of past demand to facilitate
post-operations analysis and to identify best practices for future operations.

Functionality
DDR2 gives access to:

Past traffic data - from August 2012 till now, traffic demand, last filed
flight plan traffic trajectories as well as actual trajectories are provided
for past analysis;
Past and Pre-OPS (one AIRAC in advance) environment data - they can be
downloaded and are used internally for processing future traffic trajectories.
    They contain all information necessary to analyse and process sector loads,
    capacity bottlenecks, re-routing options, etc;
Strategic traffic forecast - this covers the planning phase, from 18 months to
7 days before operations. It is used for medium- to-short-term capacity
planning and seasonal planning. Users can themselves generate, with several 4D
trajectory processing options, and download this type of forecast directly via
the DDR2 web portal;
Pre-tactical traffic forecast - it focusses on the planning phase, from 6 days
to 1 day before operations. Network pre-tactical planning is supported by the
NM PREDICT system and can be accessed via the DDR2 portal;
NEST and SAAM tools - they can be downloaded from DDR2 Web portal and are
compatible with DDR data. These tools analyse and process a great deal of information for the purpose of facilitating airspace design and capacity planning in Europe.

Users
The DDR addresses the needs of a wide range of users such as:

air navigation service providers (ANSPs), who use it to prepare and optimise
their capacity plans;
airlines, who rely on it to detect flight efficiency improvement opportunities,
by visualising and comparing flight plan trajectories for any period of time in
    the past;
airspace management actors, for airspace management and coordination of the
available airspace;
airports, with the aim of integrating their local plans with the Network
Operations Plan;
the NM at central/FAB/local level.

See EUROCONTROL NEST Manual Section 9.7 for fileformat descriptions

"""
module DDR2import
    include("util.jl")
    include("T5.jl")
    include("SO6.jl")
    include("Exp2.jl")
    include("Allftplus.jl")
    include("Ase.jl")
    include("Are.jl")
    include("Gar.jl")
    include("Gsl.jl")
    include("Frp.jl")
    include("Sls.jl")
    include("Sid.jl")
    include("Star.jl")
    include("Spc.jl")
    include("Mot.jl")
    include("Narp.jl")
    include("Nnpt.jl")
    include("Ntfv.jl")
    include("Cost.jl")
    include("Crco.jl")
    include("Ur.jl")
    include("Arp.jl")
    include("Country.jl")
    include("Operator.jl")
    include("Routes.jl")
    include("Runway.jl")
    include("For.jl")
end # module

#TODO Add tests for last element in test files.


text:
module NaturalStrings

export NaturalChar, NaturalString

import Base: codepoint, read, write, ncodeunits, codeunnit, isvalid

struct NaturalChar{T} <: AbstractChar
    parent::T
end

codepoint(c::NaturalChar{<:AbstractChar}) = codepoint(c.parent)
read(io::IO, c::NaturalChar{<:AbstractChar}) = read(io, c.parent)
write(io::IO, c::NaturalChar{<:AbstractChar}) = write(io, c.parent)
Char(c::NaturalChar{<:AbstractChar}) = Char(c.parent)

==(a::NaturalChar, b::NaturalChar) = (a.parent == b.parent)

# Even though "01" == "1", `isless` should provide a total ordering.

isless(a::NaturalChar{T}, b::NaturalChar{T}) where {T<:AbstractChar}
    a = a.parent
    b = b.parent
    isless(a,b) # TODO
end

isless(a::NaturalChar{<:AbstractChar}, b::NaturalChar{<:AbstractString}) = false
isless(a::NaturalChar{<:AbstractString}, b::NaturalChar{<:AbstractChar}) = true

struct NaturalString{T} <: AbstractString
    parent::T
end

ncodeunits(s::NaturalString) = ncodeunits(s.parent)
codeunit(s::NaturalString) = codeunit(s.parent)

# It is assumed that digits require exactly one code-unit. (This seems to be the case in all common encodings.)
function isvalid(s::NaturalString, i::Integer)
    @boundscheck 1 <= i <= ncodeunits(s) || return false
    i == 1 && return true
    isvalid(s.parent, i) || return false
    isdigit(s.parent[i]) || return true
    @inbounds isvalid(s.parent, i-1) || return true
    return !isdigit(@inbounds s.parent[i-1])
end

end # module


text:
#=
function x0_volume(model::EoS,z; phase = :unknown)
    if phase == :unknown || is_liquid(phase)
        if typeof(model)<:SAFTgammaMie
            x0 = [log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)/0.8)]
        elseif typeof(model)<:SAFT
            x0 = [log10(π/6*N_A*sum(z[i]*model.params.segment[i]*model.params.sigma[i]^3 for i in model.components)/0.8)]
        elseif typeof(model)<:Cubic
            x0 = [log10(sum(z[i]*z[j]*model.params.b[union(i,j)] for i in model.components for j in model.components)/0.8)]
        end
    elseif is_vapour(phase)
        if typeof(model)<:SAFTgammaMie
            x0 = [log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)/1e-2)]
        elseif typeof(model)<:SAFT
            x0 = [log10(π/6*N_A*sum(z[i]*model.params.segment[i]*model.params.sigma[i]^3 for i in model.components)/1e-2)]
        elseif typeof(model)<:Cubic
            x0 = [log10(sum(z[i]*z[j]*model.params.b[union(i,j)] for i in model.components for j in model.components)/1e-2)]
        end
    elseif is_supercritical(phase)
        if typeof(model)<:SAFTgammaMie
            x0 = [log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)/0.5)]
        elseif typeof(model)<:SAFT
            x0 = [log10(π/6*N_A*sum(z[i]*model.params.segment[i]*model.params.sigma[i]^3 for i in model.components)/0.5)]
        elseif typeof(model)<:Cubic
            x0 = [log10(sum(z[i]*z[j]*model.params.b[union(i,j)] for i in model.components for j in model.components)/0.5)]
        end
    end

    return x0
end
=#

#=x0_volume=#
# function x0_volume(model::SAFTgammaMie,z; phase = :unknown)
#     if phase == :unknown || is_liquid(phase)
#             x0 = [log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)/0.8)]
#     elseif is_vapour(phase)
#             x0 = [log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)/1e-2)]
#     elseif is_supercritical(phase)
#             x0 = [/0.5)]
#     end
#     return x0
# end

# function x0_volume(model::LJSAFT,z; phase = :unknown)
#     if phase == :unknown || is_liquid(phase)
#         x0 = [log10(π/6*sum(z[i]*model.params.segment[i]*model.params.b[i] for i in model.components)/0.8)]
#     elseif is_vapour(phase)
#         x0 = [log10(π/6*sum(z[i]*model.params.segment[i]*model.params.b[i] for i in model.components)/1e-2)]
#     elseif is_supercritical(phase)
#         x0 = [/0.5)]
#     end
#     return x0
# end

#=
lb_volume:

SAFTgammaMie:  log10(π/6*N_A*sum(z[i]*sum(model.group_multiplicities[i][k]*model.params.segment[k]*model.params.shapefactor[k]*model.params.sigma[k]^3 for k in @groups(i)) for i in @comps)
LJSAFT: log10(π/6*sum(z[i]*model.params.segment[i]*model.params.b[i] for i in @comps)
=#

"""
    x0_volume_liquid(model,T,z)

Returns an initial guess to the liquid volume, dependent on temperature and composition. by default is 1.25 times the [lower bound volume](@ref lb_volume).
"""
function x0_volume_liquid(model,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*1.25
end

function x0_volume_liquid(model::SAFTVRMieModel,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*1.5
end

function x0_volume_liquid(model::SAFTgammaMieModel,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*2.0
end
"""
    x0_volume_gas(model,p,T,z)

Returns an initial guess to the gas volume, depending of pressure, temperature and composition. by default uses a [virial aproximation](@ref volume_virial)
"""
function x0_volume_gas(model,p,T,z)
    return volume_virial(model,p,T,z)
end

function x0_volume_sc(model,p,T,z)
    v_sc = lb_volume(model,z)
    return v_sc*2
end

"""
    x0_volume(model::EoSModel,p,T,z; phase = :unknown)

Returns an initial guess of the volume at a pressure, temperature, composition and suggested phase.

If the suggested phase is `:unkwown` or `:liquid`, calls [`x0_volume_liquid`](@ref).

If the suggested phase is `:gas`, calls [`x0_volume_gas`](@ref).

"""
function x0_volume(model::EoSModel,p,T,z; phase = :unknown)
    phase = Symbol(phase)
    if phase === :unknown || is_liquid(phase)
        return x0_volume_liquid(model,T,z)
    elseif is_vapour(phase)
        return x0_volume_gas(model,p,T,z)
    elseif is_supercritical(phase)
     else
        error("unreachable state on x0_volume")
    end
end

#=x0_sat_pure=#

# function x0_sat_pure(model::SAFTVRQMie)
#     x0    = [log10(π/6*N_A*model.params.segment[model.components[1]]*model.params.sigma[model.components[1]]^3/0.2),
#     log10(π/6*N_A*model.params.segment[model.components[1]]*model.params.sigma[model.components[1]]^3/1e-3)]
# end

# function x0_sat_pure(model::LJSAFT)
#     x0    = [log10(π/6*model.params.segment[model.components[1]]*model.params.b[model.components[1]]/0.5),
#     log10(π/6*model.params.segment[model.components[1]]*model.params.b[model.components[1]]/1e-3)]
# end

##=lb_volume=#
#
#lb_volume(model::LJSAFT,z; phase = :unknown) = [log10(π/6*sum(z[i]*model.params.segment[i]*model.params.b[i] for i in model.components)/1)]

"""
    lb_volume(model::EoSModel,z=SA[1.0])

Returns the lower bound volume.

It has different meanings depending on the Equation of State, but symbolizes the minimum allowable volume at a certain composition:

- SAFT EoS: the packing volume
- Cubic EoS, covolume (b) parameter

On empiric equations of state, the value is chosen to match the volume of the conditions at maximum pressure and minimum temperature
, but the equation itself normally can be evaluated at lower volumes.

On SAFT and Cubic EoS, volumes lower than `lb_volume` will likely error.

The lower bound volume is used for guesses of liquid volumes at a certain pressure, saturated liquid volumes and critical volumes.
"""
function lb_volume end

function lb_volume(model::SAFTModel, z = SA[1.0])
    seg = model.params.segment.values
    σᵢᵢ = model.params.sigma.diagvalues
    val = π/6*N_A*sum(z[i]*seg[i]*σᵢᵢ[i]^3 for i in 1:length(z))
    return val
end

function lb_volume(model::CubicModel,z = SA[1.0])
    n = sum(z)
    V = 1e-5
    T = 0.
    invn = one(n)/n
    b = model.params.b.values
    c = @f(translation,model.translation)
    b̄ = dot(z,Symmetric(b),z)*invn*invn
    c̄ = dot(z,c)*invn
    return b̄-c̄
end

function lb_volume(model::ActivityModel,z = SA[1.0])
    b = maximum([model.puremodel[i].params.b.values[1,1] for i ∈ @comps])
    return b
end

function lb_volume(model::CPAModel,z = SA[1.0])
    n = sum(z)
    invn = one(n)/n
    b = model.params.b.values
    b̄ = dot(z,Symmetric(b),z)*invn*invn
    return b̄

end

function lb_volume(model::SAFTgammaMieModel, z = SA[1.0])
    vk  = model.groups.n_flattenedgroups
    seg = model.params.segment.values
    S   = model.params.shapefactor.values
    σᵢᵢ = model.params.sigma.diagvalues
    val = π/6*N_A*sum(z[i]*sum(vk[i][k]*seg[k]*S[k]*σᵢᵢ[k]^3 for k in @groups(i)) for i in @comps)
    return val
end

function lb_volume(model::BACKSAFTModel, z = SA[1.0])
    seg = model.params.segment.values
    σᵢᵢ = model.params.sigma.diagvalues
    α   = model.params.alpha.values
    val = π/6*N_A*sum(z[i]*α[i]*seg[i]*σᵢᵢ[i]^3 for i in 1:length(z))
    return val
end

function lb_volume(model::LJSAFTModel, z = SA[1.0])
    seg = model.params.segment.values
    b = model.params.b.diagvalues
    val = π/6*sum(z[i]*seg[i]*b[i] for i in 1:length(z))
    return val
end

#=scale_sat_pure=#

# function scale_sat_pure(model::SAFTgammaMie)
#     m̄  = sum(model.group_multiplicities[model.components[1]][k]*model.params.segment[k]*model.params.shapefactor[k] for k in @groups(model.components[1]))
#         σ̄3 = sum(model.group_multiplicities[model.components[1]][k]*model.group_multiplicities[model.components[1]][l]*
#                  model.params.segment[k]*model.params.segment[l]*
#                  model.params.shapefactor[k]*model.params.shapefactor[l]*
#                  model.params.sigma[union(k,l)]^3 for k in @groups(model.components[1]) for l in @groups(model.components[1]))/m̄^2
#         ϵ̄  = T_crit_pure(model)
#         p_scale    = σ̄3*N_A/R̄/ϵ̄
#         μ_scale    = 1/R̄/ϵ̄
#     return p_scale,μ_scale
# end

# function scale_sat_pure(model::LJSAFT)
#     p_scale    = model.params.b[model.components[1]]/N_A/R̄/model.params.T[model.components[1]]
#     μ_scale    = 1/R̄/model.params.T[model.components[1]]
#     return p_scale,μ_scale
# end

"""
    x0_sat_pure(model::EoSModel,T,z=SA[1.0])

Returns a 2-element vector corresponding to `[log10(Vₗ),log10(Vᵥ)]`, where Vₗ and Vᵥ are the liquid and vapor initial guesses.
Used in [`sat_pure`](@ref).
"""
function x0_sat_pure(model,T,z=SA[1.0])
    #=theory as follows
    #given T = Teos:
    #calculate B(Teos)
    PB = Peos = -2B

    veos = volume_compress(model,PB,T)
    with a (P,V,T,B) pair, change to
    (P,V,a,b)
=#
    B = second_virial_coefficient(model,T,SA[1.0])
    lb_v = lb_volume(model,SA[1.0])*one(T)
    _0 = zero(B)
    #virial volume below lower bound volume.
    #that means that we are way over the critical point
    if -2B < lb_v
        _nan = _0/_0
        return [_nan,_nan]
    end
    p = -0.25*R̄*T/B
    vl = x0_volume(model,p,T,z,phase=:l)
    vl = volume_compress(model,p,T,SA[1.0],V0=vl)
    #=the basis is that p = RT/v-b - a/v2
    we have a (p,v,T) pair
    and B = 2nd virial coefficient = b-a/RT
    with that, we solve for a and b
    as a and b are vdW, Pc and Tc can be calculated
    with Tc and Pc, we use [1] to calculate vl0 and vv0
    with Tc, we can also know in what regime we are.
    in near critical pressures, we use directly vv0 = -2B
    and vl0 = 4*lb_v

    [1]
    DOI: 10.1007/s10910-007-9272-4
    Journal of Mathematical Chemistry, Vol. 43, No. 4, May 2008 (© 2007)
    The van der Waals equation: analytical and approximate solutions
    =#
    γ = p*vl*vl/(R̄*T)
    _c = vl*vl + B*vl - γ*vl
    _b = γ - B - vl
    Δ = _b*_b - 4*_c
    if isnan(vl) | (Δ < 0)
        #fails on two ocassions:
        #near critical point, or too low.
        #old strategy
        x0l = 4*lb_v
        x0v = -2*B + 2*lb_v
        return [log10(x0l),log10(x0v)]
    end
    Δsqrt = sqrt(Δ)
    b1 = 0.5*(-_b + Δsqrt)
    b2 = 0.5*(-_b - Δsqrt)
    if b1 < 0
        b = b2
    elseif b1 > vl
        b = b2
    else
        b = b1
    end
    a = -R̄*T*(B-b)
    Vc = 3*b
    Ωa =  27/64
    Ωb =  1/8
    ar = a/Ωa
    br = b/Ωb
    Tc = ar/br/R̄
    Pc = ar/(br*br)
    Tr = T/Tc
    #Tr(vdW approx) > Tr(model)
    if Tr >= 1
        _nan = _0/_0
        return [_nan,_nan]
    end
    # if b1/b2 < -0.95, then T is near Tc.
    #if b<lb_v then we are in trouble
    #critical regime or something horribly wrong happened
    if (b1/b2 < -0.95) | (b<lb_v) | (Tr>0.99)
        x0l = 4*lb_v
        x0v = -2*B #gas volume as high as possible
        return [log10(x0l),log10(x0v)]
    end
    Vl0,Vv0 = vdw_x0_xat_pure(T,Tc,Pc,Vc)
    x0l = min(Vl0,vl)
    x0v = min(1e4*one(Vv0),Vv0) #cutoff volume
    return [log10(x0l),log10(x0v)]
end

function vdw_x0_xat_pure(T,T_c,P_c,V_c)
    Tr = T/T_c
    Trm1 = 1.0-Tr
    Trmid = sqrt(Trm1)
    if Tr >= 0.7
        c_l = 1.0+2.0*Trmid + 0.4*Trm1 - 0.52*Trmid*Trm1 +0.115*Trm1*Trm1 #Eq. 29
    else
        c_l = 1.5*(1+sqrt(1-(32/27)*Tr)) #Eq. 32
    end

    if Tr >= 0.46
        #Eq. 30, valid in 0.46 < Tr < 1
        c_v = 1.0-2.0*Trmid + 0.4*Trm1 + 0.52*Trmid*Trm1 +0.207*Trm1*Trm1
    elseif Tr <= 0.33
        #Eq. 33, valid in 0 < Tr < 0.33
        c_v = (3*c_l/(ℯ*(3-c_l)))*exp(-(1/(1-c_l/3)))

    else
        #Eq. 31 valid in 0.25 < Tr < 1
        mean_c = 1.0 + 0.4*Trm1 + 0.161*Trm1*Trm1
        c_v = 2*mean_c - c_l
    end
    #volumes predicted by vdW
    Vl0 = (1/c_l)*V_c
    Vv0 = (1/c_v)*V_c
    return (Vl0,Vv0)
end
function scale_sat_pure(model::EoSModel,z=SA[1.0])
    p    = 1/p_scale(model,z)
    μ    = 1/R̄/T_scale(model,z)
    return p,μ
end

#=x0_crit_pure=#
# x0_crit_pure(model::LJSAFT) = [1.5, log10(π/6*model.params.segment[model.components[1]]*model.params.b[model.components[1]]/0.3)]

"""
    x0_crit_pure(model::SAFTModel,z=SA[1.0])

Returns a 2-element vector corresponding to
    `[k,log10(Vc0)]`, where `k` is `Tc0/T_scale(model,z)`
"""
function x0_crit_pure end

function x0_crit_pure(model::SAFTModel,z=SA[1.0])
    lb_v = lb_volume(model,z)
    [2.0, log10(lb_v/0.3)]
end

function x0_crit_pure(model::sCKSAFTModel,z=SA[1.0])
    lb_v = lb_volume(model,z)
    res = [5.0, log10(lb_v/0.3)]
    #@show res
    return res
end

function x0_crit_pure(model::CubicModel,z=SA[1.0])
    lb_v = lb_volume(model,z)
    [1.0, log10(lb_v/0.3)]
end

function x0_crit_pure(model::CPAModel,z=SA[1.0])
    lb_v = lb_volume(model,z)
    [2.0, log10(lb_v/0.3)]
end

function x0_crit_pure(model::EoSModel,z=SA[1.0])
    lb_v = lb_volume(model,z)
    [1.5, log10(lb_v/0.3)]
end

# x0_crit_pure(model::Cubic) = [1.0, log10(model.params.b[model.components[1]]/0.3)]
# x0_crit_pure(model::CPA) = [2, log10(model.params.b[model.components[1]]/0.3)]

#=T_crit_pure=#
# function T_crit_pure(model::SAFTgammaMie)
#     m̄ = sum(model.group_multiplicities[model.components[1]][k]*model.params.segment[k]*model.params.shapefactor[k] for k in @groups(model.components[1]))
#     return sum(model.group_multiplicities[model.components[1]][k]*model.group_multiplicities[model.components[1]][l]*
#                model.params.segment[k]*model.params.segment[l]*
#                model.params.shapefactor[k]*model.params.shapefactor[l]*
#                model.params.epsilon[union(k,l)] for k in @groups(model.components[1]) for l in @groups(model.components[1]))/m̄^2
# end
# T_crit_pure(model::LJSAFT) = model.params.T[model.components[1]]

# T_crit_pure(model::Cubic) = model.params.a[model.components[1]]/model.params.b[model.components[1]]/8.314*8/27

#=
 temperature scaling factor,
on critical based EoS, is a function of critical temperature
on SAFT EoS, is a function of ϵ
=#

"""
    T_scale(model::EoS,z=SA[1.0])

Represents a temperature scaling factor.

On any EoS based on Critical parameters (Cubic or Empiric EoS), the temperature scaling factor is chosen to be the critical temperature.

On SAFT or other molecular EoS, the temperature scaling factor is chosen to be a function of the potential depth ϵ.

Used as scaling factors in [`saturation_pressure`](@ref) and as input for solving [`crit_pure`](@ref)
"""
function T_scale end

function T_scale(model::SAFTModel,z=SA[1.0])
    ϵ = model.params.epsilon.diagvalues
    return prod(ϵ[i]^z[i] for i in 1:length(z))^(1/sum(z))
end

function T_scale(model::SAFTgammaMieModel,z=SA[1.0])
    return T_scale(model.vrmodel,z)
end

function T_scale(model::LJSAFTModel,z=SA[1.0])
    T̃ = model.params.T_tilde.diagvalues
    return prod(T̃[i]^z[i] for i in 1:length(z))^(1/sum(z))
end

#dont use αa, just a, to avoid temperature dependence
function T_scale(model::CubicModel,z=SA[1.0])
    n = sum(z)
    invn2 = one(n)/(n*n)
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = dot(z, Symmetric(_a), z)*invn2/Ωa
    b = dot(z, Symmetric(_b), z)*invn2/Ωb
    return a/b/R̄
end

function T_scale(model::ActivityModel,z=SA[1.0])
    n = sum(z)
    invn2 = one(n)/(n*n)
    Ωa,Ωb = ab_consts(model.puremodel[1])
    _a = [model.puremodel[i].params.a.values[1,1] for i in @comps]
    _b = [model.puremodel[i].params.b.values[1,1] for i in @comps]
    a = dot(z, _a)*invn2/Ωa
    b = dot(z, _b)*invn2/Ωb
    return a/b/R̄
end

function T_scale(model::CPAModel,z=SA[1.0])
    n = sum(z)
    invn2 = one(n)/(n*n)
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = dot(z, Symmetric(_a), z)*invn2/Ωa
    b = dot(z, Symmetric(_b), z)*invn2/Ωb
    return a/b/R̄

end

"""
    p_scale(model::SAFTModel,z=SA[1.0])

Represents a pressure scaling factor

On any EoS based on Critical parameters (Cubic or
Empiric EoS), the pressure scaling factor is
chosen to be a function of the critical pressure.

On SAFT or other molecular EoS, the temperature
scaling factor is chosen to a function of ∑(zᵢ*ϵᵢ*(σᵢᵢ)³)

Used as scaling factors in [`saturation_pressure`](@ref) and as input for solving [`crit_pure`](@ref)

"""
function p_scale(model::SAFTModel,z=SA[1.0])
    ϵ = model.params.epsilon.diagvalues
    σᵢᵢ = model.params.sigma.diagvalues
    val =  sum(z[i]*σᵢᵢ[i]^3/ϵ[i] for i in 1:length(z))*N_A/R̄
    return 1/val
end

function p_scale(model::LJSAFTModel,z=SA[1.0])
    T̃ = model.params.T_tilde.diagvalues
    b = model.params.b.diagvalues
    val =  sum(z[i]*b[i]/T̃[i] for i in 1:length(z))/R̄
    return 1/val
end

function p_scale(model::SAFTgammaMieModel,z=SA[1.0])
    V = zero(first(z))
    T = zero(first(z))
    σ̄3 = @f(σ3x)
    ϵ̄ = T_scale(model,z)
    val    = σ̄3*N_A/R̄/ϵ̄
    return 1/val
end

function p_scale(model::CubicModel,z=SA[1.0])
    n = sum(z)
    invn2 = (1/n)^2
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = invn2*dot(z, Symmetric(_a), z)/Ωa
    b = invn2*dot(z, Symmetric(_b), z)/Ωb
    return a/ (b^2) # Pc mean
end

function p_scale(model::ActivityModel,z=SA[1.0])
    n = sum(z)
    invn2 = (1/n)^2
    Ωa,Ωb = ab_consts(model.puremodel[1])
    _a = [model.puremodel[i].params.a.values[1,1] for i in @comps]
    _b = [model.puremodel[i].params.b.values[1,1] for i in @comps]
    a = dot(z, _a)*invn2/Ωa
    b = dot(z, _b)*invn2/Ωb
    return a/ (b^2) # Pc mean
end

function p_scale(model::CPAModel,z=SA[1.0])
    n = sum(z)
    invn2 = (1/n)^2
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = invn2*dot(z, Symmetric(_a), z)/Ωa
    b = invn2*dot(z, Symmetric(_b), z)/Ωb
    return a/ (b^2) # Pc mean
end

#=
the following methods are fallbacks,
that require just the definition of T_scale,p_scale and lb_volume
respectively. if possible, each eos should define those
=#
function T_scales(model,z)
    n = length(z)
    x = zeros(n)
    res = zeros(n)
    for i = 1:n
        x[i] = 1.0
        res[i] = T_scale(model,x)
        x[i] = 0.0
    end
    return res
end

# function lb_volumes(model,z)
#     n = length(z)
#     x = zeros(n)
#     res = zeros(n)
#     for i = 1:n
#         x[i] = 1.0
#         res[i] = lb_volume(model,x)
#         x[i] = 0.0
#     end
#     return res
# end

# function p_scales(model,z)
#     n = length(z)
#     x = zeros(n)
#     res = zeros(n)
#     for i = 1:n
#         x[i] = 1.0
#         res[i] = p_scale(model,x)
#         x[i] = 0.0
#     end
# end


text:
using LightGraphs

"""
Read a graph in DIMACS format. Assumes file is correctly formatted.
"""
function read_graph(filename::AbstractString)::Graph
	f = open(filename)
	g = nothing
	for line in eachline(f)
		pieces = split(line)			#split strings into pieces separated by space (default delimiter)
		if pieces[1] == "p"
			nvertices = parse(Int, pieces[3])		#converts string to int values
			g = Graph(nvertices)
		elseif pieces[1] == "e"
			if g == nothing
				error("clq file does not define number of vertices in header")
			end
			add_edge!(g, parse(Int, pieces[2]), parse(Int, pieces[3]))
		end
	end
	return g
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
abstract type AbstractMemoryBuffer{T} <: AbstractVector{T} end

abstract type AbstractMutableBuffer{T} <: AbstractMemoryBuffer{T} end

abstract type AbstractImmutableBuffer{T} <: AbstractMemoryBuffer{T} end

mutable struct StaticBuffer{T,N} <: AbstractMutableBuffer{T}
    data::NTuple{N,T}
end

struct MutableBuffer{T} <: AbstractMutableBuffer{T}
    data::StaticBuffer{T}
    length::Int

    MutableBuffer{T}(n::Int) where {T} = new{T}(StaticBuffer{T,n}(), n)
end

struct StaticImmutableBuffer{T,N} <: AbstractImmutableBuffer{T}
    data::NTuple{N,T}
end

struct ImmutableBuffer{T} <: AbstractImmutableBuffer{T}
    data::StaticImmutableBuffer{T}
    length::Int

    ImmutableBuffer(data::StaticImmutableBuffer{T}, length::Int) where {T} = new{T}(data, length)
end

@inline function Base.unsafe_convert(::Type{Ptr{T}}, d::StaticBuffer) where {T}
    return Base.unsafe_convert(Ptr{T}, Base.pointer_from_objref(d.data.data))
end
@inline function Base.unsafe_convert(::Type{Ptr{T}}, d::MutableBuffer) where {T}
    return Base.unsafe_convert(Ptr{T}, Base.pointer_from_objref(d.data.data))
end

@propagate_inbounds function Base.getindex(m::AbstractMutableBuffer{T}, i::Int) where {T}
    @boundscheck checkbounds(m, i)
    GC.@preserve m x = unsafe_load(pointer(m), i)
    return x
end
@propagate_inbounds function Base.getindex(m::StaticImmutableBuffer{T}, i::Int) where {T}
    @boundscheck checkbounds(m, i)
    return @inbounds(getfield(m.data, i))
end
@propagate_inbounds function Base.getindex(m::ImmutableBuffer{T}, i::Int) where {T}
    @boundscheck checkbounds(m, i)
    return @inbounds(getfield(m.data.data, i))
end
@propagate_inbounds function Base.setindex!(m::AbstractMutableBuffer{T}, x, i::Int) where {T}
    @boundscheck checkbounds(m, i)
    GC.@preserve m unsafe_store!(pointer(m), convert(T, x), i)
end

ArrayInterface.known_length(::Type{StaticImmutableBuffer{T,N}}) where {T,N} = N
ArrayInterface.known_length(::Type{StaticBuffer{T,N}}) where {T,N} = N
ArrayInterface.size(m::MutableBuffer) = (static_length(m),)

""" grow_beg! """
function grow_beg!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    return unsafe_grow_beg!(x, n)
end
unsafe_grow_beg!(x::Vector, n) = Base._growbeg!(x, n)
unsafe_grow_beg!(x::AbstractVector, n) = unsafe_grow_beg!(parent(x), n)

###
### grow_end!
###
""" grow_end! """
function grow_end!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    unsafe_grow_end!(x, n)
    return x
end
unsafe_grow_end!(x::Vector, n) = Base._growend!(x, n)
unsafe_grow_end!(x::AbstractVector, n) = unsafe_grow_end!(parent(x), n)

""" grow_to! """
function grow_to!(x, n::Integer)
    len = length(x)
    if len < n
        unsafe_grow_to!(x, n)
        return x
    elseif len == n
        return x
    else
        throw(ArgumentError("new length must be ≥ than length of collection, got length $(length(x))."))
    end
end
unsafe_grow_to!(x, n::Integer) = unsafe_grow_end!(x, n - length(x))

""" shrink_beg! """
function shrink_beg!(x, n::Integer)
    n > length(x) && throw(ArgumentError("new length cannot be < 0"))
    return unsafe_shrink_beg!(x, n)
end
unsafe_shrink_beg!(x::Vector, n) = Base._deletebeg!(x, n)
unsafe_shrink_beg!(x::AbstractVector, n) = unsafe_shrink_beg!(parent(x), n)

""" shrink_end! """
function shrink_end!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    return unsafe_shrink_end!(x, n)
end
unsafe_shrink_end!(x::Vector, n) = Base._deleteend!(x, n)
unsafe_shrink_end!(x::AbstractVector, n) = unsafe_shrink_end!(parent(x), n)

""" shrink_to! """
function shrink_to!(x, n::Integer)
    len = length(x)
    if len > n
        return unsafe_shrink_to!(x, n)
    elseif len == n
        return x
    else
        throw(ArgumentError("new length must be ≤ than length of collection, got length $(length(x))."))
    end
end

unsafe_shrink_to!(x, n) = unsafe_shrink_end!(x, length(x) - n)