text:
# # Arbitrary Precision
#
# COSMO allows you to solve problems with arbitrary floating-point precision, e.g. by using `BigFloat` problem data. To do this, the desired floating point type has to be consistent across the model `COSMO.Model{<: AbstractFloat}`, the input data and the (optional) settings object `COSMO.Settings{<: AbstractFloat}`.
# As an example, assume we want to solve the following quadratic program:
#
# $$
# \begin{array}{ll} \text{minimize} &  1/2 x^\top P x + q^\top x \\
# \text{subject to} &  l \leq A x \leq u
# \end{array}
# $$
# where $P = \begin{bmatrix} 4 & 1 \\ 1 & 2\end{bmatrix}$, $q = [1, 1]^\top$, $A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \\ 0 & 1\end{bmatrix}$ and $l= [1,0, 0]^\top$, $u=[1, 0.7, 0.7]^\top$. We start by creating the model with the desired precision:
using COSMO, LinearAlgebra, SparseArrays
model = COSMO.Model{BigFloat}()

#-
# Next, we define the problem data as `BigFloat` arrays and create the constraint:
q = BigFloat[1; 1.]
P = sparse(BigFloat[4. 1; 1 2])
A = BigFloat[1. 1; 1 0; 0 1]
l = BigFloat[1.; 0; 0]
u = BigFloat[1; 0.7; 0.7]

constraint = COSMO.Constraint(A, zeros(BigFloat, 3), COSMO.Box(l, u))

# Notice that the constraint type parameter is dependent on the input data. The same is true for the constraint set `Box`. Next, we define the settings
settings = COSMO.Settings{BigFloat}(verbose = true, kkt_solver = QdldlKKTSolver)

# and assemble and solve the problem:
assemble!(model, P, q, constraint, settings = settings)
result = COSMO.optimize!(model);

# Moreover, notice that when no type parameter is specified, all objects default to `Float64`:
model = COSMO.Model()

# Two limitations to arbitrary precision:
# - Since we call the LAPACK function `syevr` for eigenvalue decompositions, we currently only support solving problems with PSD constraints in `Float32` and `Float64`.
# - We suggest to use the pure Julia QDLDL linear system solver (`kkt_solver = QdldlKKTSolver`) when working with arbitrary precision types as some of the other available solvers don't support all available precisions.

#md # !!! note
#md #     `JuMP` does not currently support arbitrary precision. However, if you want to use `COSMO` directly with `MathOptInterface`, you can use: `COSMO.Optimizer{<: AbstractFloat}` as your optimizer. Again, the problem data precision of your MathOptInterface-model has to agree with the optimizer's precision.


text:
using Plots, JLD, Utilities
plotly()

N_est = 5
NGWpoints = 3
Tvals =  [  600 400 200 ]#[50 250 1250]#
N_Tind = length(Tvals)[1]
scoreRescType = ""#"FISHER-EWMA"#
useStartVal = false# true# false
save_fold = "./data/estimatesTest/compComplexity/"
##

Nmin = 10; Nspacing =  10; Nmax = 70
file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"; save_path = save_fold*file_name#
@load(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)

NvalsAll = Nvals
est_timesAll = est_times
##

Nmin = 80; Nspacing =  10; Nmax = 100
file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"; save_path = save_fold*file_name#
@load(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)

NindSelect = 1
NvalsAll = [Nvals[1];NvalsAll ]
est_timesAll = cat(2, est_times[1:end,NindSelect:NindSelect,1:end,1:end],est_timesAll)

##

GWvalsLeg = ["1 par" ; "N/2 par"; " N par"]

est_times[1,1,1:end,1:end]
#est_times[indT,indN,indGW,n]
##
h = plot(); indGW = 2
 for t = 1:3
    est_timesAll[t,1:end,indGW,1:end]
    plot!(h,NvalsAll,meanSq(est_timesAll[t,1:end,indGW,1:end],2),label = "$(Tvals[t])")
 end
 plot(h, title = GWvalsLeg[indGW])

## add more single par estimates
Nmin = 500; Nspacing =  500; Nmax = 10000; NGWpoints = 1;N_est = 1
file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"; save_path = save_fold*file_name#
save_path
@load(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)
tmp = est_times[:,:,1,1]

##
NvalsAll_GW1 = [NvalsAll;Nvals]
est_timesAll_1GW = [meanSq(est_timesAll[:,:,1,:],3) tmp]

h = plot(); indGW = 1
 nnz = est_timesAll_1GW[1,:].!=0
 for t = 1:3
    est_timesAll_1GW[t,nnz]
    plot!(h,NvalsAll_GW1[nnz],est_timesAll_1GW[t,nnz],label = "$(Tvals[t])")
 end
 plot(h, title = GWvalsLeg[indGW])

#

##

ext: jl


text:
import ApproxFun: weight, dimension

export JacobiQWeight

"""
`JacobiQWeight`
weights a basis on `ℂ\𝕀` weighted by `(z+1)^α*(z-1)^β`.
Note the inconsistency of the parameters with `JacobiQ`.
When the domain is `[a,b]` the weight is inferred by mapping to `[-1,1]`
"""
immutable JacobiQWeight{S,DD} <: WeightSpace{S,RealBasis,DD,1}
    α::Float64
    β::Float64
    space::S
    function JacobiQWeight(α::Float64,β::Float64,space::S)
        if isa(space,JacobiQWeight)
            JacobiQWeight(α+space.α,β+space.β,space.space)
        else
            new(α,β,space)
        end
    end
end

JacobiQWeight(a::Number,b::Number,d::RealUnivariateSpace)=JacobiQWeight{typeof(d),typeof(domain(d))}(Float64(a),Float64(b),d)
JacobiQWeight(a::Number,b::Number,d::IntervalDomain)=JacobiQWeight(Float64(a),Float64(b),Space(d))
JacobiQWeight(a::Number,b::Number,d::Vector)=JacobiQWeight(Float64(a),Float64(b),Space(d))
JacobiQWeight(a::Number,b::Number)=JacobiQWeight(a,b,Chebyshev())

JacobiQWeight(a::Number,b::Number,s::PiecewiseSpace) = PiecewiseSpace(JacobiQWeight(a,b,vec(s)))

spacescompatible(A::JacobiQWeight,B::JacobiQWeight)=A.α==B.α && A.β == B.β && spacescompatible(A.space,B.space)
#spacescompatible{D<:IntervalDomain}(A::JacobiQWeight,B::RealUnivariateSpace{D})=spacescompatible(A,JacobiQWeight(0,0,B))
#spacescompatible{D<:IntervalDomain}(B::RealUnivariateSpace{D},A::JacobiQWeight)=spacescompatible(A,JacobiQWeight(0,0,B))

transformtimes{JW1<:JacobiQWeight,JW2<:JacobiQWeight}(f::Fun{JW1},g::Fun{JW2})=
            Fun(coefficients(transformtimes(Fun(f.coefficients,f.space.space),
                                            Fun(g.coefficients,g.space.space))),
                             JacobiQWeight(f.space.α+g.space.α,f.space.β+g.space.β,f.space.space))
transformtimes{JW<:JacobiQWeight}(f::Fun{JW},g::Fun) = Fun(coefficients(transformtimes(Fun(f.coefficients,f.space.space),g)),f.space)
transformtimes{JW<:JacobiQWeight}(f::Fun,g::Fun{JW}) = Fun(coefficients(transformtimes(Fun(g.coefficients,g.space.space),f)),g.space)

##  α and β are opposite the convention for JacobiQ polynomials
# Here, α is the left algebraic singularity and β is the right algebraic singularity.

jacobiQweight(α,β,x)=(x+1).^α.*(x-1).^β
jacobiQweight(α,β,d::Domain)=Fun([1.],JacobiQWeight(α,β,ConstantSpace(d)))
jacobiQweight(α,β)=jacobiQweight(α,β,Interval())

weight(sp::JacobiQWeight,x)=jacobiQweight(sp.α,sp.β,tocanonical(sp,x))
dimension(sp::JacobiQWeight)=dimension(sp.space)

Base.first{JW<:JacobiQWeight}(f::Fun{JW})=space(f).α>0?zero(eltype(f)):f(first(domain(f)))
Base.last{JW<:JacobiQWeight}(f::Fun{JW})=space(f).β>0?zero(eltype(f)):f(last(domain(f)))

setdomain(sp::JacobiQWeight,d::Domain)=JacobiQWeight(sp.α,sp.β,setdomain(sp.space,d))


text:
const DEFAULT_NAMESPACE = "default"
const DEFAULT_URI = "http://localhost:8001"

struct KApi
    api::DataType
    types::Module
end

mutable struct KuberContext
    client::Swagger.Client
    apis::Dict{Symbol,Vector{KApi}}
    modelapi::Dict{Symbol,KApi}
    namespace::String

    function KuberContext()
        ctx = Swagger.Client(DEFAULT_URI)
        ctx.headers["Connection"] = "close"
        new(ctx, Dict{Symbol,Vector}(), Dict{Symbol,KApi}(), DEFAULT_NAMESPACE)
    end
end

convert(::Type{Vector{UInt8}}, s::T) where {T<:AbstractString} = collect(codeunits(s))
convert(::Type{T}, json::String) where {T<:SwaggerModel} = convert(T, JSON.parse(json))
convert(::Type{Dict{String,Any}}, model::T) where {T<:SwaggerModel} = JSON.parse(JSON.json(model))

is_json_mime(mime::T) where {T <: AbstractString} = ("*/*" == mime) || occursin(r"(?i)application/json(;.*)?", mime) || occursin(r"(?i)application/(.*)-patch\+json(;.*)?", mime)

kind_to_type(ctx::KuberContext, kind::String, version::Union{String,Nothing}=nothing) = kind_to_type(ctx, Symbol(kind), version)
function kind_to_type(ctx::KuberContext, kind::Symbol, version::Union{String,Nothing}=nothing)
    types = (version === nothing) ? (ctx.modelapi[kind]).types : api_typedefs(version)
    getfield(types, kind)
end

kuber_type(ctx::KuberContext, d) = kuber_type(ctx, Any, d)
kuber_type(ctx::KuberContext, T, data::String) = kuber_type(ctx, T, JSON.parse(data))

function kuber_type(ctx::KuberContext, T, resp::HTTP.Response)
    ctype = HTTP.header(resp, "Content-Type", "application/json")
    !is_json_mime(ctype) && return T
    kuber_type(ctx, T, String(copy(resp.body)))
end

function kuber_type(ctx::KuberContext, T, j::Dict{String,Any})
    if ("kind" in keys(j)) && !isempty(ctx.apis)
        kind = j["kind"]
        version = ("apiVersion" in keys(j)) ? j["apiVersion"] : nothing
        try
            T = kind_to_type(ctx, kind, version)
        catch ex
            @warn("Type not found.", kind, version)
        end
    end
    T
end

kuber_obj(ctx::KuberContext, data::String) = kuber_obj(ctx, JSON.parse(data))
kuber_obj(ctx::KuberContext, j::Dict{String,Any}) = convert(kind_to_type(ctx, j["kind"], get(j, "apiVersion", nothing)), j)

show(io::IO, ctx::KuberContext) = print("Kubernetes namespace ", ctx.namespace, " at ", ctx.client.root)

get_server(ctx::KuberContext) = ctx.client.root
get_ns(ctx::KuberContext) = ctx.namespace

function set_server(ctx::KuberContext, uri::String=DEFAULT_URI, reset_api_versions::Bool=false; max_tries=1, kwargs...)
    rtfn = (default,data)->kuber_type(ctx, default, data)
    ctx.client = Swagger.Client(uri; get_return_type=rtfn, kwargs...)
    ctx.client.headers["Connection"] = "close"
    reset_api_versions && set_api_versions!(ctx; max_tries=max_tries)
    ctx.client
end

set_ns(ctx::KuberContext, namespace::String) = (ctx.namespace = namespace)

camel(a) = string(uppercase(a[1])) * (a[2:end])

api_group(group) = api_group(String(group))
function api_group(group::String)
    endswith(group, ".k8s.io") && (group = group[1:end-7])
    group_parts = split(group, ".")
    join([camel(String(x)) for x in group_parts])
end

api_group_type(group_ver) = api_group_type(String(group_ver))
function api_group_type(group_ver::String)
    group, ver = occursin('/', group_ver) ? split(group_ver, "/") : ("Core", group_ver)
    group = api_group(group)
    ver = camel(ver)
    getfield(@__MODULE__, Symbol(group * ver * "Api"))
end

api_typedefs(group_ver) = api_typedefs(String(group_ver))
function api_typedefs(group_ver::String)
    group, ver = occursin('/', group_ver) ? split(group_ver, "/") : ("Core", group_ver)
    group = api_group(group)
    ver = camel(ver)
    getfield(getfield(@__MODULE__, :Typedefs), Symbol(group * ver))
end

api_method(group::Symbol, verb::String, object::String) = getfield(@__MODULE__, Symbol(verb * String(group)[1:end-3] * object))
function api_method(group_ver::String, verb::String, object::String)
    group, ver = occursin('/', group_ver) ? split(group_ver, "/") : ("Core", group_ver)
    group = api_group(group)
    ver = camel(ver)
    getfield(@__MODULE__, Symbol(verb * group * ver * object))
end

function override_pref(name, server_pref, override)
    if override !== nothing
        for (n,v) in override
            (n == name) && (return v)
        end
    end
    server_pref
end

function fetch_misc_apis_versions(ctx::KuberContext; override=nothing, verbose::Bool=false, max_tries=1)
    apis = ctx.apis
    vers = @repeat max_tries try
        getAPIVersions(ApisApi(ctx.client))
    catch e
        @retry if isa(e, IOError)
            @debug("Retrying getAPIVersions")
            sleep(2)
        end
    end
    api_groups = vers.groups
    for apigrp in api_groups
        name = apigrp.name
        pref_vers_type = apigrp.preferredVersion
        pref_vers_version = override_pref(name, pref_vers_type.version, override)
        pref_vers = name * "/" * pref_vers_version
        verbose && @info("$name ($(api_group(name))) versions", supported=join(map(x->x.version, apigrp.versions), ", "), preferred=pref_vers_version)

        try
            apis[Symbol(api_group(name))] = [KApi(api_group_type(pref_vers), api_typedefs(pref_vers))]
        catch ex
            if isa(ex, UndefVarError)
                @info("unsupported $pref_vers")
                continue
            else
                rethrow()
            end
        end

        for api_vers in apigrp.versions
            try
                gt = api_group_type(api_vers.groupVersion)
                td = api_typedefs(api_vers.groupVersion)
                ka = KApi(gt, td)
                kalist = apis[Symbol(api_group(name))]
                (ka == kalist[1]) || push!(kalist, ka)
            catch
                @info("unsupported $(api_vers.groupVersion)")
            end
        end
    end
    apis
end

function fetch_core_version(ctx::KuberContext; override=nothing, verbose::Bool=false, max_tries=1)
    apis = ctx.apis
    api_vers = @repeat max_tries try
        getCoreAPIVersions(CoreApi(ctx.client))
    catch e
        @retry if isa(e, IOError)
            @debug("Retrying getCoreAPIVersions")
            sleep(2)
        end
    end
    name = "Core"
    pref_vers = override_pref(name, api_vers.versions[1], override)
    verbose && @info("Core versions", supported=join(api_vers.versions, ", "), preferred=pref_vers)
    apis[:Core] = [KApi(getfield(@__MODULE__, Symbol("Core" * camel(pref_vers) * "Api")), getfield(getfield(@__MODULE__, :Typedefs), Symbol("Core" * camel(pref_vers))))]
    for api_vers in api_vers.versions
        try
            gt = getfield(@__MODULE__, Symbol("Core" * camel(api_vers) * "Api"))
            td = getfield(getfield(@__MODULE__, :Typedefs), Symbol("Core" * camel(api_vers)))
            ka = KApi(gt, td)
            kalist = apis[:Core]
            (ka == kalist[1]) || push!(kalist, ka)
        catch
            @info("unsupported Core $api_vers")
        end
    end
    apis
end

function fetch_api_machinery(ctx::KuberContext)
    apis = ctx.apis
    apis[:Apis] = [KApi(Kuber.Kubernetes.ApisApi, Kuber.Typedefs.Apis)]
    apis[:MetaV1] = [KApi(Kuber.Kubernetes.ApisApi, Kuber.Typedefs.MetaV1)]
    apis
end

function build_model_api_map(ctx::KuberContext)
    #@info("building model - api map...")
    modelapi = ctx.modelapi
    for (apigroup,apivers) in ctx.apis
        apiver = apivers[1]
        types = apiver.types
        #@info("building model for $apiver")

        for name in names(types; all=true)
            (name in [:eval, Symbol("#eval"), :include, Symbol("#include"), Symbol(split(string(types), '.')[end])]) && continue
            modelapi[name] = apiver
        end
    end
    modelapi
end

function set_api_versions!(ctx::KuberContext; override=nothing, verbose::Bool=false, max_tries=1)
    empty!(ctx.apis)
    empty!(ctx.modelapi)

    # bootstrap: the api machinery types
    fetch_api_machinery(ctx)
    build_model_api_map(ctx)

    # fetch apis and map the types
    fetch_core_version(ctx; override=override, verbose=verbose, max_tries=max_tries)
    fetch_misc_apis_versions(ctx; override=override, verbose=verbose, max_tries=max_tries)
    build_model_api_map(ctx)

    # add custom models
    ctx.modelapi[:PodLog] = ctx.modelapi[:Pod]
    nothing
end


text:
using HSL

T = Float64
# definitions applicable to all packages
const data_map = Dict{Type, Type}(Float32 => Cfloat,
                                  Float64 => Cdouble,
                                  Complex64 => Cfloat,
                                  Complex128 => Cdouble)

rows = Int32[1, 1, 2, 2, 3, 3, 5]
cols = Int32[1, 2, 3, 5, 3, 4, 5]
vals = T[2, 3, 4, 6, 1, 5, 1]
A = sparse(rows, cols, vals) ; A = A + triu(A, 1)'

H57 = convert(SparseMatrixCSC{data_map[T],Int32}, A)
L = tril(H57)
M = ma57_coord(L.n, findnz(L)...)

b = T[8, 45, 31, 15, 17]
xexact = T[1, 2, 3, 4, 5]

ϵ = sqrt(eps(eltype(A)))
ma57_factorize(M)
x = ma57_solve(M, b)

(L, D57, s, p) = ma57_get_factors(M)
S = spdiagm(s)
P = speye(M.n)[p, :]

r0 = P * S * A * S * P' - L * D57 * L'
pi = invperm(p)
r0b = S[p,:] * A * S[:,pi] - L * D57 * L'

Sm1 = spdiagm(1.0 ./ s)

x1 = ma57_solve(M, b, job = :LS)
r1 = (L*S * x1) - (P*S*b)
r1b = (L * x1) - (P*S*b)
x2 = ma57_solve(M, b, job = :LPS)
r2 = (L'*P*Sm1 * x2) - (Sm1*b)
r2b = (L'*P*Sm1 * x2) - b

H97 = convert(SparseMatrixCSC{data_map[T],Int64}, A)
PLD = Ma97(H97, print_level=-1)
ma97_factorize(PLD)

ret = ma97_inquire(PLD)
vD = ret[2]
pivot = ret[1]
# Take care of zero valued pivots
z(a)= a==0
pivot[find(z,pivot)] = collect(1:length(pivot))[find(z,pivot)]

vD1 = copy(vD[1:1,:])'[:,1:1]
vD2 = copy(vD[2:2,1:end-1])'[:,1:1]

D97 = spdiagm((vD1,vD2,vD2),[0 -1 1]) # don't forget inquire returns D^-1


text:
using SparseArrays
using LinearAlgebra

"""Load an LDPC matrix from a text file in alist format."""
function load_alist(file_path::AbstractString; check_redundant=false,)
    if file_extension(file_path) != ".alist"
        @warn "load_alist called on file with extension '$(file_extension(file_path))', expected '.alist'"
    end

    file = open(file_path, "r")
    nVN, nCN = space_sep_ints(readline(file))
    dmax_VN, dmax_CN = space_sep_ints(readline(file))
    var_node_degs = space_sep_ints(readline(file))
    check_node_degs = space_sep_ints(readline(file))
    remaining_lines = readlines(file)
    close(file)

    if length(remaining_lines) != nVN + nCN
        error("Number of lines in $file_path is inconcistent with stated matrix size.")
    end

    if dmax_CN != maximum(check_node_degs)
        error("Alist file $file_path claims: max. CN degree=$dmax_CN but contents give $(maximum(check_node_degs)).")
    end

    if dmax_VN != maximum(var_node_degs)
        error("Alist file $file_path claims: max. VN degree=$dmax_CN but contents give $(maximum(var_node_degs)).")
    end

    # parity check matrix
    H = spzeros(Int8, nCN, nVN)

    # fill the matrix
    for col_ind in 1:nVN
        rows = space_sep_ints(remaining_lines[col_ind])

        if check_redundant && length(rows) != var_node_degs[col_ind]
            error("Variable node degree in $file_path inconcistent with below data for VN $col_ind.")
        end

        for row_ind in rows
            H[row_ind, col_ind] = 1
        end
    end

    # the second half of the alist file is redundant. Check that it is consistent.
    if check_redundant
        entry_counter = 0
        for row_ind in 1:nCN
            cols = space_sep_ints(remaining_lines[nVN + row_ind])

            check_node_degree = length(cols)
            if check_node_degree != check_node_degs[row_ind]
                error("Check node degree in $file_path inconcistent with below data for CN $row_ind.")
            end

            entry_counter += check_node_degree
            for col_ind in cols
                if H[row_ind, col_ind] != 1
                    error("VN and CN specifications in $file_path disagree on matrix entry ($row_ind, $col_ind).")
                end
            end
        end

        if entry_counter != sum(H)
            error("VN and CN specification in $file_path are inconsistent.")
        end
    end

    return H
end

"""
    function save_to_alist(matrix::AbstractArray{Int8,2}, out_file_path::String)

Save LDPC matrix to file in alist format. For details about the format, see:
https://aff3ct.readthedocs.io/en/latest/user/simulation/parameters/codec/ldpc/decoder.html#dec-h-path-image-required-argument
http://www.inference.org.uk/mackay/codes/alist.html
todo test this carefully
"""
function save_to_alist(matrix::AbstractArray{Int8,2}, out_file_path::String)

    (the_M, the_N) = size(matrix)

    variable_node_degrees = get_variable_node_degrees(matrix)
    check_node_degrees = get_check_node_degrees(matrix)

    # write data as specified by the alist format
    lines = String[]
    # -- Part 1 --
    # 'the_N' is the total number of variable nodes and 'the_M' is the total number of check nodes
    push!(lines, "$the_N $the_M")

    # 'dmax_VN' is the highest variable node degree and 'dmax_CN' is the highest check node degree
    push!(lines, "$(maximum(variable_node_degrees)) $(maximum(check_node_degrees))")

    # list of the degrees for each variable node
    push!(lines, join(["$deg" for deg in variable_node_degrees], " "))

    # list of the degrees for each check node
    push!(lines, join(["$deg" for deg in check_node_degrees], " "))

    # -- Part 2 --
    # each following line describes the check nodes connected to a variable node, the first
    # check node index is '1' (and not '0')
    # variable node '1'
    """
        Get indices of elements equal to one in a matrix.
        :param matrix: 2d numpy array
        :return: Array of strings, one string with indices for each line in the matrix.
    """
    function get_node_indices(matrix::AbstractArray{Int8,2})
        # the first check node index is '1' (and not '0')
        degrees = [findall(row .== 1) for row in eachrow(matrix)]
        return [join(string.(index_list), " ") for index_list in degrees]
    end
    append!(lines, get_node_indices(transpose(matrix)))

    # -- Part 3 --
    # each following line describes the variables nodes connected to a check node, the first
    # variable node index is '1' (and not '0')
    # check node '1'
    append!(lines, get_node_indices(matrix))

    open(out_file_path, "w+") do file
        for line in lines
            println(file, line)
        end
    end

    return nothing
end

# helper methods for testing the parity check matrix
function get_variable_node_degrees(matrix::AbstractArray{Int8,2})
    @assert(length(size(matrix)) == 2, "Matrix required. Wrong number of dimensions")
    return [sum(row) for row in eachcol(matrix)]
end

function get_check_node_degrees(matrix::AbstractArray{Int8,2})
    @assert(length(size(matrix)) == 2, "Matrix required. Wrong number of dimensions")
    return [sum(row) for row in eachrow(matrix)]
end


text:
module TestAudioIONodes

using FactCheck
using AudioIO
import AudioIO: AudioSample, AudioBuf, AudioRenderer, AudioNode
import AudioIO: DeviceInfo, render

include("testhelpers.jl")

# A TestNode just renders out 1:buf_size each frame
type TestRenderer <: AudioRenderer
    buf::AudioBuf
    TestRenderer(buf_size::Integer) = new(AudioSample[1:buf_size])
end

typealias TestNode AudioNode{TestRenderer}
TestNode(buf_size) = TestNode(TestRenderer(buf_size))

function render(node::TestRenderer,
                device_input::AudioBuf,
                info::DeviceInfo)
    return node.buf
end

test_info = DeviceInfo(44100, 512)
dev_input = zeros(AudioSample, test_info.buf_size)

facts("Validating TestNode allocation") do
    # first validate that the TestNode doesn't allocate so it doesn't mess up our
    # other tests
    test = TestNode(test_info.buf_size)
    # JIT
    render(test, dev_input, test_info)
    @fact (@allocated render(test, dev_input, test_info)) => 16
end

#### AudioMixer Tests ####

# TODO: there should be a setup/teardown mechanism and some way to isolate
# tests

facts("AudioMixer") do
    context("0 Input Mixer") do
        mix = AudioMixer()
        render_output = render(mix, dev_input, test_info)
        @fact render_output => AudioSample[]
        @fact (@allocated render(mix, dev_input, test_info)) => 48
    end

    context("1 Input Mixer") do
        testnode = TestNode(test_info.buf_size)
        mix = AudioMixer([testnode])
        render_output = render(mix, dev_input, test_info)
        @fact render_output => AudioSample[1:test_info.buf_size]
        @fact (@allocated render(mix, dev_input, test_info)) => 64
    end

    context("2 Input Mixer") do
        test1 = TestNode(test_info.buf_size)
        test2 = TestNode(test_info.buf_size)
        mix = AudioMixer([test1, test2])
        render_output = render(mix, dev_input, test_info)
        # make sure the two inputs are being added together
        @fact render_output => 2 * AudioSample[1:test_info.buf_size]
        @fact (@allocated render(mix, dev_input, test_info)) => 96
        # now we'll stop one of the inputs and make sure it gets removed
        stop(test1)
        render_output = render(mix, dev_input, test_info)
        # make sure the two inputs are being added together
        @fact render_output => AudioSample[1:test_info.buf_size]

        stop(mix)
        render_output = render(mix, dev_input, test_info)
        @fact render_output => AudioSample[]
    end
end

MSE_THRESH = 1e-7

facts("SinOSC") do
    freq = 440
    # note that this range includes the end, which is why there are
    # sample_rate+1 samples
    t = linspace(0, 1, int(test_info.sample_rate+1))
    test_vect = convert(AudioBuf, sin(2pi * t * freq))
    context("Fixed Frequency") do
        osc = SinOsc(freq)
        render_output = render(osc, dev_input, test_info)
        @fact mse(render_output, test_vect[1:test_info.buf_size]) =>
                lessthan(MSE_THRESH)
        render_output = render(osc, dev_input, test_info)
        @fact mse(render_output,
                test_vect[test_info.buf_size+1:2*test_info.buf_size]) =>
                lessthan(MSE_THRESH)
        @fact (@allocated render(osc, dev_input, test_info)) => 64
        stop(osc)
        render_output = render(osc, dev_input, test_info)
        @fact render_output => AudioSample[]
    end

    context("Testing SinOsc with signal input") do
        t = linspace(0, 1, int(test_info.sample_rate+1))
        f = 440 .- t .* (440-110)
        dt = 1 / test_info.sample_rate
        # NOTE - this treats the phase as constant at each sample, which isn't strictly
        # true. Unfortunately doing this correctly requires knowing more about the
        # modulating signal and doing the real integral
        phase = cumsum(2pi * dt .* f)
        unshift!(phase, 0)
        expected = convert(AudioBuf, sin(phase))

        freq = LinRamp(440, 110, 1)
        osc = SinOsc(freq)
        render_output = render(osc, dev_input, test_info)
        @fact mse(render_output, expected[1:test_info.buf_size]) =>
                lessthan(MSE_THRESH)
        render_output = render(osc, dev_input, test_info)
        @fact mse(render_output,
                expected[test_info.buf_size+1:2*test_info.buf_size]) =>
                lessthan(MSE_THRESH)
        # give a bigger budget here because we're rendering 2 nodes
        @fact (@allocated render(osc, dev_input, test_info)) => 160
    end
end

facts("AudioInput") do
    node = AudioInput()
    test_data = rand(AudioSample, test_info.buf_size)
    render_output = render(node, test_data, test_info)
    @fact render_output => test_data
end

facts("ArrayPlayer") do
    context("playing long sample") do
        v = rand(AudioSample, 44100)
        player = ArrayPlayer(v)
        render_output = render(player, dev_input, test_info)
        @fact render_output => v[1:test_info.buf_size]
        render_output = render(player, dev_input, test_info)
        @fact render_output => v[(test_info.buf_size + 1) : (2*test_info.buf_size)]
        @fact (@allocated render(player, dev_input, test_info)) => 192
        stop(player)
        render_output = render(player, dev_input, test_info)
        @fact render_output => AudioSample[]
    end

    context("testing end of vector") do
        # give a vector just a bit larger than 1 buffer size
        v = rand(AudioSample, test_info.buf_size + 1)
        player = ArrayPlayer(v)
        render(player, dev_input, test_info)
        render_output = render(player, dev_input, test_info)
        @fact render_output => v[test_info.buf_size+1:end]
    end
end

facts("Gain") do
    context("Constant Gain") do
        gained = TestNode(test_info.buf_size) * 0.75
        render_output = render(gained, dev_input, test_info)
        @fact render_output => 0.75 * AudioSample[1:test_info.buf_size]
        @fact (@allocated render(gained, dev_input, test_info)) => 32
    end
    context("Gain by a Signal") do
        gained = TestNode(test_info.buf_size) * TestNode(test_info.buf_size)
        render_output = render(gained, dev_input, test_info)
        @fact render_output => AudioSample[1:test_info.buf_size] .* AudioSample[1:test_info.buf_size]
        @fact (@allocated render(gained, dev_input, test_info)) => 48
    end
end

facts("LinRamp") do
    ramp = LinRamp(0.25, 0.80, 1)
    expected = convert(AudioBuf, linspace(0.25, 0.80, int(test_info.sample_rate+1)))
    render_output = render(ramp, dev_input, test_info)
    @fact mse(render_output, expected[1:test_info.buf_size]) =>
            lessthan(MSE_THRESH)
    render_output = render(ramp, dev_input, test_info)
    @fact mse(render_output,
            expected[(test_info.buf_size+1):(2*test_info.buf_size)]) =>
            lessthan(MSE_THRESH)
    @fact (@allocated render(ramp, dev_input, test_info)) => 64
end

facts("Offset") do
    offs = TestNode(test_info.buf_size) + 0.5
    render_output = render(offs, dev_input, test_info)
    @fact render_output => 0.5 + AudioSample[1:test_info.buf_size]
    @fact (@allocated render(offs, dev_input, test_info)) => 32
end

end # module TestAudioIONodes


text:
"""
    findall(in(interval), x::AbstractRange)

Return all indices `i` for which `x[i] ∈ interval`, specialized for
the case where `x` is a range, which enables constant-time complexity.

# Examples

```jldoctest
julia> x = range(0,stop=3,length=10)
0.0:0.3333333333333333:3.0

julia> collect(x)'
1×10 adjoint(::Vector{Float64}) with eltype Float64:
 0.0  0.333333  0.666667  1.0  1.33333  1.66667  2.0  2.33333  2.66667  3.0

julia> findall(in(1..6), x)
4:10
```

It also works for decreasing ranges:
```jldoctest
julia> y = 8:-0.5:0
8.0:-0.5:0.0

julia> collect(y)'
1×17 adjoint(::Vector{Float64}) with eltype Float64:
 8.0  7.5  7.0  6.5  6.0  5.5  5.0  4.5  4.0  3.5  3.0  2.5  2.0  1.5  1.0  0.5  0.0

julia> findall(in(1..6), y)
5:15

julia> findall(in(Interval{:open,:closed}(1,6)), y) # (1,6], does not include 1
5:14
```
"""
function Base.findall(interval_d::Base.Fix2{typeof(in),Interval{L,R,T}}, x::AbstractRange)  where {L,R,T}
    isempty(x) && return 1:0

    interval = interval_d.x
    il, ir = firstindex(x), lastindex(x)
    δx = step(x)
    a,b = if δx < zero(δx)
        rev = findall(in(interval), reverse(x))
        isempty(rev) && return rev

        a = (il+ir)-last(rev)
        b = (il+ir)-first(rev)

        a,b
    else
        lx, rx = first(x), last(x)
        l = max(leftendpoint(interval), lx - oneunit(δx))
        r = min(rightendpoint(interval), rx + oneunit(δx))

        (l > rx || r < lx) && return 1:0

        a = il + max(0, round(Int, cld(l-lx, δx)))
        a += (a ≤ ir && (x[a] == l && L == :open || x[a] < l))

        b = min(ir, round(Int, cld(r-lx, δx)) + il)
        b -= (b ≥ il && (x[b] == r && R == :open || x[b] > r))

        a,b
    end
    # Reversing a range could change sign of values close to zero (cf
    # sign of the smallest element in x and reverse(x), where x =
    # range(BigFloat(-0.5),stop=BigFloat(1.0),length=10)), or more
    # generally push elements in or out of the interval (as can cld),
    # so we need to check once again.
    a += +(a < ir && x[a] ∉ interval) - (il < a && x[a-1] ∈ interval)
    b += -(il < b && x[b] ∉ interval) + (b < ir && x[b+1] ∈ interval)

    a:b
end

# We overload Base._findin to avoid an ambiguity that arises with
# Base.findall(interval_d::Base.Fix2{typeof(in),Interval{L,R,T}}, x::AbstractArray)
function Base._findin(a::Union{AbstractArray, Tuple}, b::Interval)
    ind  = Vector{eltype(keys(a))}()
    @inbounds for (i,ai) in pairs(a)
        ai in b && push!(ind, i)
    end
    ind
end


text:
struct MathLinkError <: Exception
    err::Error
    msg::String
end
function MathLinkError(link::Link)
    err = ccall((:MLError, mlib), Error, (CLink,), link)
    msg = unsafe_string(ccall((:MLErrorMessage, mlib), Cstring, (CLink,), link))
    clearerror(link)
    newpacket(link)
    MathLinkError(err, msg)
end

clearerror(link::Link) =
    ccall((:MLClearError, mlib), Cint, (CLink,), link)
endpacket(link::Link) =
    ccall((:MLEndPacket, mlib), Cint, (CLink,), link)
newpacket(link::Link) =
    ccall((:MLNewPacket, mlib), Cint, (CLink,), link)

nextpacket(link::Link) =
    ccall((:MLNextPacket, mlib), Packet, (CLink,), link)

macro wschk(expr)
    link = expr.args[5] # first argument
    quote
        if $(esc(expr)) == 0
            throw(MathLinkError($(esc(link))))
        end
    end
end

function getargcount(link::Link)
    r = Ref{Cint}()
    # int WSGetArgCount(WSLINK link,int *n)
    @wschk ccall((:MLGetArgCount, mlib), Cint,
                 (CLink, Ptr{Cint}),
                 link, r)
    r[]
end

function putargcount(link::Link, n::Integer)
    # int WSPutArgCount(WSLINK link,int n)
    @wschk ccall((:MLPutArgCount, mlib), Cint,
                 (CLink, Cint),
                 link, n)
    nothing
end
function putfunction(link::Link, w::WSymbol, nargs::Integer)
    # int MLPutUTF8Function(( MLLINK l , const unsigned char * s , int v , int n )
    @wschk ccall((:MLPutUTF8Function, mlib), Cint,
                 (CLink, Ptr{UInt8}, Cint, Cint),
                 link, w.name, sizeof(w.name), nargs)
    nothing
end
function getfunction(link::Link)
    r_str = Ref{Ptr{UInt8}}()
    r_bytes = Ref{Cint}()
    r_nargs = Ref{Cint}()
    # int MLGetUTF8Function( MLLINK l , const unsigned char ** s , int * v , int * n )
    @wschk ccall((:MLGetUTF8Function, mlib), Cint,
                 (CLink, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}),
                 link, r_str, r_bytes, r_nargs)
    str = unsafe_string(r_str[], r_bytes[])
    # void MLReleaseUTF8Symbol(MLLINK link,const unsigned char *s,int len)
    ccall((:MLReleaseUTF8Symbol, mlib), Cvoid,
          (CLink, Ptr{UInt8}, Cint),
          link, r_str[], r_bytes[])
    return WSymbol(str), r_nargs[]
end

# WSymbol
function put(link::Link, w::WSymbol)
    # int MLPutUTF8Symbol(MLLINK link,const unsigned char *s,int len)
    @wschk ccall((:MLPutUTF8Symbol, mlib), Cint,
                 (CLink, Ptr{UInt8}, Cint),
                 link, w.name, sizeof(w.name))
    nothing
end
function get(link::Link, ::Type{WSymbol})
    r_str = Ref{Ptr{UInt8}}()
    r_bytes = Ref{Cint}()
    r_chars = Ref{Cint}()
    # int MLGetUTF8Symbol(MLLINK link,const unsigned char **s,int *b,int *c)
    @wschk ccall((:MLGetUTF8Symbol, mlib), Cint,
                 (CLink, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}),
                 link, r_str, r_bytes, r_chars)
    str = unsafe_string(r_str[], r_bytes[])
    # void MLReleaseUTF8Symbol(MLLINK link,const unsigned char *s,int len)
    ccall((:MLReleaseUTF8Symbol, mlib), Cvoid,
          (CLink, Ptr{UInt8}, Cint),
          link, r_str[], r_bytes[])
    return WSymbol(str)
end

# WReal/WInteger
function put(link::Link, w::WReal)
    # int MLPutRealNumberAsUTF8String(MLLINK l, const unsigned char *s, int n)
    @wschk ccall((:MLPutRealNumberAsUTF8String, mlib), Cint,
                 (CLink, Ptr{UInt8}, Cint),
                 link, w.value, sizeof(w.value))
    nothing
end
function get(link::Link, ::Type{W}) where W <: Union{WReal,WInteger}
    r_str = Ref{Ptr{UInt8}}()
    r_bytes = Ref{Cint}()
    r_chars = Ref{Cint}()
    # int MLGetNumberAsUTF8String(MLLINK l, const unsigned char **s, int *v, int *c)
    @wschk ccall((:MLGetNumberAsUTF8String, mlib), Cint,
                 (CLink, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}),
                 link, r_str, r_bytes, r_chars)
    str = unsafe_string(r_str[], r_bytes[])
    # void MLReleaseUTF8String(MLLINK link,const unsigned char *s,int len)
    ccall((:MLReleaseUTF8String, mlib), Cvoid,
          (CLink, Ptr{UInt8}, Cint),
          link, r_str[], r_bytes[])
    return W(str)
end

# String
function put(link::Link, str::Union{String,SubString})
    # int MLPutUTF8String(MLLINK link,const unsigned char *s,int len)
    @wschk ccall((:MLPutUTF8String, mlib), Cint,
                 (CLink, Ptr{UInt8}, Cint),
                 link, str, sizeof(str))
    nothing
end
function get(link::Link, ::Type{String})
    r_str = Ref{Ptr{UInt8}}()
    r_bytes = Ref{Cint}()
    r_chars = Ref{Cint}()
    # int MLGetUTF8String(MLLINK link,const unsigned char **s,int *b,int *c)
    @wschk ccall((:MLGetUTF8String, mlib), Cint,
                 (CLink, Ptr{Ptr{UInt8}}, Ptr{Cint}, Ptr{Cint}),
                 link, r_str, r_bytes, r_chars)
    str = unsafe_string(r_str[], r_bytes[])
    # void MLReleaseUTF8String(MLLINK link,const unsigned char *s,int len)
    ccall((:MLReleaseUTF8String, mlib), Cvoid,
          (CLink, Ptr{UInt8}, Cint),
          link, r_str[], r_bytes[])
    return str
end

for (f, T) in [
    (:Integer8, Int8)
    (:Integer16, Int16)
    (:Integer32, Int32)
    (:Integer64, Int64)
    (:Real32, Float32)
    (:Real64, Float64)
]
    @eval begin
        function put(link::Link, x::$T)
            # note slightly bizarre handling of Float32
            @wschk ccall(($(QuoteNode(Symbol(:MLPut, f))), mlib), Cint, (CLink, $(f == :Real32 ? Float64 : T)), link, x)
            nothing
        end
        function get(link::Link, ::Type{$T})
            r = Ref{$T}()
            @wschk ccall(($(QuoteNode(Symbol(:MLGet, f))), mlib), Cint, (CLink, Ptr{$T}), link, r)
            r[]
        end
    end
end

# Get fns
gettype(link::Link) =
    ccall((:MLGetType, mlib), Token, (CLink,), link)

getrawtype(link::Link) =
    ccall((:MLGetRawType, mlib), Token, (CLink,), link)

puttype(link::Link, t::Token) =
    @wschk ccall((:MLPutType, mlib), Cint, (CLink, Token), link, t)

getnext(link::Link) =
    ccall((:MLGetNext, mlib), Token, (CLink,), link)

getnextraw(link::Link) =
    ccall((:MLGetNextRaw, mlib), Token, (CLink,), link)


text:
using Compose, Colors

# const padding = 0.5mm

function sierpinski(n::Int)
    if n == 0
        compose(context(), polygon([(1,1), (0,1), (1/2, 0)]))
    else
        t = sierpinski(n - 1)
        compose(context(),
                (context(1/4,   0, 1/2, 1/2), t),
                (context(  0, 1/2, 1/2, 1/2), t),
                (context(1/2, 1/2, 1/2, 1/2), t))
        # compose(context(),
        #         (context(0.25w + padding, padding,        0.5w - 2*padding, 0.5h - 2*padding), t),
        #         (context(padding,         0.5h + padding, 0.5w - 2*padding, 0.5h - 2*padding), t),
        #         (context(0.5w + padding,  0.5h + padding, 0.5w - 2*padding, 0.5h - 2*padding), t))
    end
end

img = SVG("sierpinski.svg", 4inch, 4(√3/2)inch)
draw(img, compose(sierpinski(2), linewidth(0.1mm), fill(nothing), stroke(colorant"black")))

img = SVG("sierpinski.svg", 4inch, 4(√3/2)inch)
@time draw(img, compose(sierpinski(8), linewidth(0.1mm), fill(nothing), stroke(colorant"black")))

#img = SVG("sierpinski.svg", 4inch, 4(√3/2)inch)
#@profile draw(img, compose(sierpinski(10), linewidth(0.1mm), fill(nothing), stroke(colorant"black")))
#Profile.print()


text:
__precompile__(true)
module mdl

using JuliaParser.Diagnostics
using JuliaParser.Lexer
using JuliaParser.Lexer: @tok, TokenStream, eof, eof_token, make_token, peekchar
using JuliaParser.Tokens: ¬, √
using Render
using Render: ParseState, modifycoord!, withem!, withpm!

export mdl_execute, mdl_parser

COMMENT = "//"
peekcomment(ts::TokenStream) = Render.peekcomment(ts, COMMENT)

function next_token{T}(ts::TokenStream{T})
    ts.ateof && return eof_token(ts)
    tmp = Lexer.skipws(ts, false)
    eof(ts) && return eof_token(ts)
    ts.isspace = tmp

    while ! eof(ts)
        peekcomment(ts) && skipcomment(ts, COMMENT)
        c = peekchar(ts)
        if eof(c)
            ts.ateof = true
            return eof_token(T)
        elseif c === ' ' || c === '\t'
            Lexer.skip(ts, 1)
            continue
        elseif c == '\r'
            Lexer.skip(ts, 1)
            peekchar(ts) === '\n' || throw(Lexer.here(ts),
                                           "'\\r' not followed by '\\n' is invalid")
            continue
        elseif Lexer.isnewline(c)
            return @tok Lexer.readchar(ts)
        elseif isdigit(c)
            return Lexer.read_number(ts, false, false)
        elseif c === '-'
            Lexer.skip(ts, 1)
            d = peekchar(ts)
            return isdigit(d) ? Lexer.read_number(ts, d === '.', true) : @tok :(-)
        elseif c === '.'
            Lexer.skip(ts, 1)
            return isdigit(peekchar(ts)) ? Lexer.read_number(ts, true, false) : @tok :(.)
        elseif Lexer.is_identifier_start_char(c)
            return @tok Lexer.accum_julia_symbol(ts, c)
        else
            @assert Lexer.readchar(ts) === c
            throw(Diagnostics.diag(Lexer.here(ts), "invalid character \"$c\""))
        end
    end

    ts.ateof = true
    return eof_token(ts)
end

# simple substitution
const F64 = Float64

# MDL commands
mdlpush(ps::ParseState) = push!(ps.coords, ps.coords[end])
mdlpop(ps::ParseState)  = length(ps.coords) > 1 && pop!(ps.coords)

knobworks(ps, k, thing) = k === nothing ? thing : thing * ps.symtab[k][ps.cframe]

mdlmove(ps::ParseState, x::F64, y::F64, z::F64, knob = nothing) = begin
    modifycoord!(ps, mktranslate(knobworks(ps, knob, [x, y, z])))
end

mdlscale(ps::ParseState, x::F64, y::F64, z::F64, knob = nothing) = begin
    modifycoord!(ps, mkscale(knobworks(ps, knob, [x, y, z])))
end

mdlrotate(ps::ParseState, dir::Symbol, d::F64, knob = nothing) = begin
    modifycoord!(ps, mkrotate(knobworks(ps, knob, d), dir))
end

# optional not at end
mdlsphere(ps::ParseState, cons, x::F64, y::F64, z::F64, r::F64, coord = nothing) = begin
    withpm!(ps, addsphere!, [x, y, z], r, SPHSTEPS)
end
mdlsphere(ps::ParseState, x::F64, y::F64, z::F64, r::F64, coord = nothing) = mdlsphere(ps, nothing, x, y, z, r, coord)

mdltorus(ps::ParseState, cons, x::F64, y::F64, z::F64, r1::F64, r2::F64, coord = nothing) = begin
    withpm!(ps, addtorus!, [x, y, z], r1, r2, TORUSTEPS)
end
mdltorus(ps::ParseState, x::F64, y::F64, z::F64, r1::F64, r2::F64, coord = nothing) = mdltorus(ps, nothing, x, y, z, r1, r2, coord)

mdlbox(ps::ParseState, cons, x1::F64, y1::F64, z1::F64, x2::F64, y2::F64, z2::F64, coord = nothing) = begin
    withpm!(ps, addbox!, [x1, y1, z1], x2, y2, z2)
end
mdlbox(ps::ParseState, x1::F64, y1::F64, z1::F64, x2::F64, y2::F64, z2::F64, coord = nothing) = mdlbox(ps, nothing, x1, y1, z1, x2, y2, z2, coord)

# there should be a cleaner way to do this
mdlline(ps::ParseState, cons, x1::F64, y1::F64, z1::F64, coord1::Symbol, x2::F64, y2::F64, z2::F64, coord2 = nothing) = begin
    withem!(ps, addedge!, [x1, y1, z1], [x2, y2, z2])
end
mdlline(ps::ParseState, x1::F64, y1::F64, z1::F64, coord1::Symbol, x2::F64, y2::F64, z2::F64, coord2 = nothing) = mdlline(ps, nothing, x1, y1, z1, coord1, x2, y2, z2, coord2)
mdlline(ps::ParseState, cons::Symbol, x1::F64, y1::F64, z1::F64, x2::F64, y2::F64, z2::F64, coord2 = nothing) = mdlline(ps, cons, x1, y1, z1, nothing, x2, y2, z2, coord2)
mdlline(ps::ParseState, x1::F64, y1::F64, z1::F64, x2::F64, y2::F64, z2::F64, coord2 = nothing) = mdlline(ps, nothing, x1, y1, z1, nothing, x2, y2, z2, coord2)

mdlbasename(ps::ParseState, bn...) = ps.basename = string(bn...)
mdlframes(ps::ParseState, f::F64) = ps.nframes = Int(f)

mdlvary(ps::ParseState, knob::Symbol, sf, ef, sv::F64, ev::F64) = begin
    sf = Int(sf)
    ef = Int(ef)
    ps.symtab[knob] = get(ps.symtab, knob, Render.Knob(sf, ef, sv, ev, ps.nframes))
    ps.symtab[knob].frdata[sf:ef - 1] = [sv + (ev - sv) / (ef - sf) * i for i in 1:(ef - sf)]
end

mdlsave(ps::ParseState, fn...) = Render.display(ps, `convert - $(string(fn...))`)

mdldisplay(ps::ParseState) = begin
    if ps.nframes > 1
        run(`animate -delay $(Render.ANIM_DELAY) "$tmp/$(ps.basename)*"`)
    else
        Render.display(ps, `display`)
    end
end

function mdl_parser(ts::TokenStream, ps::ParseState)
    Render.skipws_and_comments(ts, COMMENT)
    eof(ts) && return

    args = Vector{Any}()
    t = Render.require_token(ts, next_token)
    command = nothing

    try
        command = getfield(mdl, Symbol(:mdl, ¬t))
    catch UndefVarError
        error("$(¬t) is not a valid MDL command")
    end

    while ! eof(ts)
        x = next_token(ts)
        push!(args, ¬x isa Number ? float(¬x) : ¬x)
    end

    push!(ps.commands, Render.GCommand(command, args))
end

function mdl_execute(ps::ParseState)
    # preprocess "frames", "basename", "vary"
    i = 1
    while i <= length(ps.commands)
        c = ps.commands[i]
        if c.func === mdlframes || c.func === mdlbasename
            c.func(ps, c.args...)
            deleteat!(ps.commands, i)
        else
            i += 1
        end
    end

    i = 1
    while i <= length(ps.commands)
        c = ps.commands[i]
        if c.func === mdlvary
            mdlvary(ps, c.args...)
            deleteat!(ps.commands, i)
        else
            i += 1
        end
    end

    mktempdir() do tmp
        reference = ps
        @sync @parallel for current in 1:reference.nframes
            if reference.nframes > 1
                println("Processing frame ", current)
                ps = deepcopy(reference)
            end

            ps.cframe = current
            for c in ps.commands
                c.func(ps, c.args...)
            end

            ps.nframes > 1 && mdlsave(ps, joinpath(tmp, string(ps.basename, lpad(current, 4, 0), ".png")))
        end
    end
end
end

max_stars_repo_name: AnAverageHuman/jlgraphics

max_issues_repo_name: AnAverageHuman/jlgraphics

max_forks_repo_name: AnAverageHuman/jlgraphics


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Cache

import Base.LibGit2, ..Dir, ...Pkg.PkgError
using ..Types

rewrite_url_to = "https"

const GITHUB_REGEX =
    r"^(?:git@|git://|https://(?:[\w\.\+\-]+@)?)github.com[:/](([^/].+)/(.+?))(?:\.git)?$"i

path(pkg::AbstractString) = abspath(".cache", pkg)

function mkcachedir()
    cache = joinpath(realpath("."), ".cache")
    if isdir(cache)
        return
    end

    @static if is_unix()
        if Dir.isversioned(pwd())
            rootcache = joinpath(realpath(".."), ".cache")
            if !isdir(rootcache)
                mkdir(rootcache)
            end
            symlink(rootcache, cache)
            return
        end
    end
    mkdir(cache)
end

function prefetch(pkg::AbstractString, url::AbstractString, sha1s::Vector)
    isdir(".cache") || mkcachedir()

    cache = path(pkg)
    normalized_url = normalize_url(url)

    repo = if isdir(cache)
        LibGit2.GitRepo(cache) # open repo, free it at the end
    else
        info("Cloning cache of $pkg from $normalized_url")
        try
            # clone repo, free it at the end
            LibGit2.clone(normalized_url, cache, isbare = true, remote_cb = LibGit2.mirror_cb())
        catch err
            errmsg = if isa(err, LibGit2.Error.GitError)
                "Cannot clone $pkg from $normalized_url. $(err.msg)"
            elseif isa(err, InterruptException)
                "Package `$pkg` prefetching was interrupted."
            else
                "Unknown error: $err"
            end
            isdir(cache) && rm(cache, recursive=true)
            throw(PkgError(errmsg))
        end
    end
    try
        LibGit2.set_remote_url(repo, normalized_url)
        in_cache = BitArray(map(sha1->LibGit2.iscommit(sha1, repo), sha1s))
        if !all(in_cache)
            info("Updating cache of $pkg...")
            LibGit2.fetch(repo)
            in_cache = BitArray(map(sha1->LibGit2.iscommit(sha1, repo), sha1s))
        end
        sha1s[.!in_cache]
    finally
        close(repo) # closing repo opened/created above
    end
end
prefetch(pkg::AbstractString, url::AbstractString, sha1::AbstractString...) =
    prefetch(pkg, url, AbstractString[sha1...])

function setprotocol!(proto::AbstractString)
    global rewrite_url_to

    if length(proto) == 0
        rewrite_url_to = nothing
    else
        rewrite_url_to = proto
    end
end

function normalize_url(url::AbstractString)
    global rewrite_url_to

    m = match(GITHUB_REGEX,url)
    (m === nothing || rewrite_url_to === nothing) ?
        url : "$rewrite_url_to://github.com/$(m.captures[1]).git"
end

end # module


text:
############################################################################################
### Types
############################################################################################

mutable struct EnsembleParameterTable{C} <: AbstractParameterTable
    tables::C
end

############################################################################################
### Constructors
############################################################################################

# constuct an empty table
function EnsembleParameterTable(::Nothing)
    tables = Dict{Symbol, ParameterTable}()
    return EnsembleParameterTable(tables)
end

############################################################################################
### Convert to other types
############################################################################################

import Base.Dict

function Dict(partable::EnsembleParameterTable)
    return partable.tables
end

#= function DataFrame(
        partable::ParameterTable;
        columns = nothing)
    if isnothing(columns) columns = keys(partable.columns) end
    out = DataFrame([key => partable.columns[key] for key in columns])
    return DataFrame(out)
end =#

############################################################################################
### get parameter table from RAMMatrices
############################################################################################

function EnsembleParameterTable(args...; groups)
    partable = EnsembleParameterTable(nothing)

    for (group, ram_matrices) in zip(groups, args)
        push!(partable.tables, group => ParameterTable(ram_matrices))
    end

    return partable
end

############################################################################################
### Pretty Printing
############################################################################################

function Base.show(io::IO, partable::EnsembleParameterTable)
    print(io, "EnsembleParameterTable with groups: ")
    for key in keys(partable.tables) print(io, "|", key, "|") end
    print(io, "\n")
    for key in keys(partable.tables)
        print("\n")
        print(io, key, ": \n")
        print(io, partable.tables[key])
    end
end

############################################################################################
### Additional Methods
############################################################################################

# Sorting ----------------------------------------------------------------------------------

# todo

# add a row --------------------------------------------------------------------------------

# do we really need this?
import Base.push!

function push!(partable::EnsembleParameterTable, d::AbstractDict, group)
    push!(partable.tables[group], d)
end

push!(partable::EnsembleParameterTable, d::Nothing, group) = nothing

# get group --------------------------------------------------------------------------------

get_group(partable::EnsembleParameterTable, group) = get_group(partable.tables, group)

############################################################################################
### Update Partable from Fitted Model
############################################################################################

# update generic ---------------------------------------------------------------------------
function update_partable!(partable::EnsembleParameterTable, model_identifier::AbstractDict, vec, column)
    for k in keys(partable.tables)
        update_partable!(partable.tables[k], model_identifier, vec, column)
    end
    return partable
end


text:
export FeasibilityFormNLS

"""Converts a nonlinear least-squares problem with residual ``F(x)`` to a nonlinear
optimization problem with constraints ``F(x) = r`` and objective ``\\tfrac{1}{2}\\|r\\|^2``.
In other words, converts
```math
\\begin{aligned}
       \\min_x \\quad & \\tfrac{1}{2}\\|F(x)\\|^2 \\\\
\\mathrm{s.t.} \\quad & c_L ≤ c(x) ≤ c_U \\\\
                      &   ℓ ≤   x  ≤ u
\\end{aligned}
```
to
```math
\\begin{aligned}
   \\min_{x,r} \\quad & \\tfrac{1}{2}\\|r\\|^2 \\\\
\\mathrm{s.t.} \\quad & F(x) - r = 0 \\\\
                      & c_L ≤ c(x) ≤ c_U \\\\
                      &   ℓ ≤   x  ≤ u
\\end{aligned}
```
If you rather have the first problem, the `nls` model already works as an NLPModel of
that format.
"""
mutable struct FeasibilityFormNLS{M <: AbstractNLSModel} <: AbstractNLSModel
  meta :: NLPModelMeta
  nls_meta :: NLSMeta
  internal :: M
  counters :: NLSCounters
end

show_header(io :: IO, nls :: FeasibilityFormNLS) = println(io, "FeasibilityFormNLS - Nonlinear least-squares moving the residual to constraints")

"""
    FeasibilityFormNLS(nls)
Converts a nonlinear least-squares problem with residual `F(x)` to a nonlinear
optimization problem with constraints `F(x) = r` and objective `¹/₂‖r‖²`.
"""
function FeasibilityFormNLS(nls :: AbstractNLSModel; name="$(nls.meta.name)-ffnls")
  nequ = nls.nls_meta.nequ
  meta = nls.meta
  nvar = meta.nvar + nequ
  ncon = meta.ncon + nequ
  nnzh = nls.nls_meta.nnzh + nequ + (meta.ncon == 0 ? 0 : meta.nnzh) # Some indexes can be repeated
  meta = NLPModelMeta(nvar, x0=[meta.x0; zeros(nequ)],
                      lvar=[meta.lvar; fill(-Inf, nequ)],
                      uvar=[meta.uvar; fill( Inf, nequ)],
                      ncon=ncon,
                      lcon=[zeros(nequ); meta.lcon],
                      ucon=[zeros(nequ); meta.ucon],
                      y0=[zeros(nequ); meta.y0],
                      lin=[nls.nls_meta.lin; meta.lin .+ nequ],
                      nln=[nls.nls_meta.nln; meta.nln .+ nequ],
                      nnzj=meta.nnzj + nls.nls_meta.nnzj + nequ,
                      nnzh=nnzh,
                      name=name
                     )
  nls_meta = NLSMeta(nequ, nvar, x0=[meta.x0; zeros(nequ)], nnzj=nequ, nnzh=0, lin=1:nequ, nln=Int[])

  nlp = FeasibilityFormNLS{typeof(nls)}(meta, nls_meta, nls, NLSCounters())
  finalizer(nlp -> finalize(nlp.internal), nlp)

  return nlp
end

function FeasibilityFormNLS(nls :: FeasibilityResidual)
  meta = nls.nlp.meta
  nequ = meta.ncon
  nvar = meta.nvar + nequ
  ncon = meta.ncon
  nnzj = meta.nnzj + nequ
  nnzh = meta.nnzh + nequ
  meta = NLPModelMeta(nvar, x0=[meta.x0; zeros(nequ)],
                      lvar=[meta.lvar; fill(-Inf, nequ)],
                      uvar=[meta.uvar; fill( Inf, nequ)],
                      ncon=ncon,
                      lcon=meta.lcon,
                      ucon=meta.ucon,
                      y0=meta.y0,
                      lin=meta.lin,
                      nln=meta.nln,
                      nnzj=nnzj, nnzh=nnzh
                     )
  nls_meta = NLSMeta(nequ, nvar, x0=[meta.x0; zeros(nequ)], nnzj=nequ, nnzh=0)

  nlp = FeasibilityFormNLS{FeasibilityResidual}(meta, nls_meta, nls, NLSCounters())
  finalizer(nlp -> finalize(nlp.internal), nlp)

  return nlp
end

function FeasibilityFormNLS(nls :: LLSModel; name="$(nls.meta.name)-ffnls")
  nequ = nls.nls_meta.nequ
  meta = nls.meta
  nvar = meta.nvar + nequ
  ncon = meta.ncon + nequ
  nnzh = nls.nls_meta.nnzh + nequ + (meta.ncon == 0 ? 0 : meta.nnzh) # Some indexes can be repeated
  meta = NLPModelMeta(nvar, x0=[meta.x0; zeros(nequ)],
                      lvar=[meta.lvar; fill(-Inf, nequ)],
                      uvar=[meta.uvar; fill( Inf, nequ)],
                      ncon=ncon,
                      lcon=[zeros(nequ); meta.lcon],
                      ucon=[zeros(nequ); meta.ucon],
                      y0=[zeros(nequ); meta.y0],
                      #lin=[nls.nls_meta.lin; meta.lin .+ nequ],
                      nln=[nls.nls_meta.nln; meta.nln .+ nequ],
                      nnzj=meta.nnzj + nls.nls_meta.nnzj + nequ,
                      nnzh=nnzh,
                      name=name
                     )
  nls_meta = NLSMeta(nequ, nvar, x0=[meta.x0; zeros(nequ)], nnzj=nequ, nnzh=0, lin=1:nequ, nln=Int[])

  nlp = FeasibilityFormNLS{typeof(nls)}(meta, nls_meta, nls, NLSCounters())
  finalizer(nlp -> finalize(nlp.internal), nlp)

  return nlp
end

function obj(nlp :: FeasibilityFormNLS, x :: AbstractVector)
  @lencheck nlp.meta.nvar x
  increment!(nlp, :neval_obj)
  n = nlp.internal.meta.nvar
  r = @view x[n+1:end]
  return dot(r, r) / 2
end

function grad!(nlp :: FeasibilityFormNLS, x :: AbstractVector, g :: AbstractVector)
  @lencheck nlp.meta.nvar x g
  increment!(nlp, :neval_grad)
  n = nlp.internal.meta.nvar
  g[1:n] .= 0.0
  g[n+1:end] .= @view x[n+1:end]
  return g
end

function objgrad!(nlp :: FeasibilityFormNLS, x :: Array{Float64}, g :: Array{Float64})
  @lencheck nlp.meta.nvar x g
  increment!(nlp, :neval_obj)
  increment!(nlp, :neval_grad)
  n = nlp.internal.meta.nvar
  r = @view x[n+1:end]
  f = dot(r, r) / 2
  g[1:n] .= 0.0
  g[n+1:end] .= @view x[n+1:end]
  return f, g
end

function cons!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, c :: AbstractVector)
  @lencheck nlp.meta.nvar xr
  @lencheck nlp.meta.ncon c
  increment!(nlp, :neval_cons)
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  r = @view xr[n+1:end]
  residual!(nlp.internal, x, @view c[1:ne])
  c[1:ne] .-= r
  if m > 0
    cons!(nlp.internal, x, @view c[ne+1:end])
  end
  return c
end

function jac_structure!(nlp :: FeasibilityFormNLS, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  @lencheck nlp.meta.nnzj rows cols
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  nnzjF = nlp.internal.nls_meta.nnzj
  @views jac_structure_residual!(nlp.internal, rows[1:nnzjF], cols[1:nnzjF])
  if m > 0
    idx = nnzjF .+ (1:nlp.internal.meta.nnzj)
    @views jac_structure!(nlp.internal, rows[idx], cols[idx])
    rows[idx] .+= ne
  end
  rows[end-ne+1:end] .= 1:ne
  cols[end-ne+1:end] .= n .+ (1:ne)
  return rows, cols
end

function jac_coord!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, vals :: AbstractVector)
  @lencheck nlp.meta.nvar xr
  @lencheck nlp.meta.nnzj vals
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  nnzjF = nlp.internal.nls_meta.nnzj
  nnzjc = m > 0 ? nlp.internal.meta.nnzj : 0
  I = 1:nnzjF
  @views jac_coord_residual!(nlp.internal, x, vals[I])
  if m > 0
    I = nnzjF+1:nnzjF+nnzjc
    @views jac_coord!(nlp.internal, x, vals[I])
  end
  vals[nnzjF+nnzjc+1:nnzjF+nnzjc+ne] .= -1
  return vals
end

function jprod!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, v :: AbstractVector, jv :: AbstractVector)
  @lencheck nlp.meta.nvar xr v
  @lencheck nlp.meta.ncon jv
  increment!(nlp, :neval_jprod)
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  @views jprod_residual!(nlp.internal, x, v[1:n], jv[1:ne])
  @views jv[1:ne] .-= v[n+1:end]
  if m > 0
    @views jprod!(nlp.internal, x, v[1:n], jv[ne+1:end])
  end
  return jv
end

function jtprod!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, v :: AbstractVector, jtv :: AbstractVector)
  @lencheck nlp.meta.nvar xr jtv
  @lencheck nlp.meta.ncon v
  increment!(nlp, :neval_jtprod)
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  @views jtprod_residual!(nlp.internal, x, v[1:ne], jtv[1:n])
  if m > 0
    @views jtv[1:n] .+= jtprod(nlp.internal, x, v[ne+1:end])
  end
  @views jtv[n+1:end] .= -v[1:ne]
  return jtv
end

function hess_structure!(nlp :: FeasibilityFormNLS, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  @lencheck nlp.meta.nnzh rows cols
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  nnzhF = nlp.internal.nls_meta.nnzh
  nnzhc = m > 0 ? nlp.internal.meta.nnzh : 0
  I = 1:nnzhF
  @views hess_structure_residual!(nlp.internal, rows[I], cols[I])
  if m > 0
    I = nnzhF+1:nnzhF+nnzhc
    @views hess_structure!(nlp.internal, rows[I], cols[I])
  end
  I = nnzhF+nnzhc+1:nnzhF+nnzhc+ne
  rows[I] .= n+1:n+ne
  cols[I] .= n+1:n+ne
  return rows, cols
end

function hess_structure!(nlp :: FeasibilityFormNLS{LLSModel}, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  @lencheck nlp.meta.nnzh rows cols
  n, ne = nlp.internal.meta.nvar, nlp.internal.nls_meta.nequ
  rows .= n+1:n+ne
  cols .= n+1:n+ne
  return rows, cols
end

function hess_coord!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, y :: AbstractVector, vals :: AbstractVector;
                     obj_weight :: Float64=1.0)
  @lencheck nlp.meta.nvar xr
  @lencheck nlp.meta.ncon y
  @lencheck nlp.meta.nnzh vals
  increment!(nlp, :neval_hess)
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  nnzhF = nlp.internal.nls_meta.nnzh
  nnzhc = m > 0 ? nlp.internal.meta.nnzh : 0
  x = @view xr[1:n]
  y1 = @view y[1:ne]
  y2 = @view y[ne+1:ne+m]
  I = 1:nnzhF
  @views hess_coord_residual!(nlp.internal, x, y1, vals[I])
  if m > 0
    I = nnzhF+1:nnzhF+nnzhc
    @views hess_coord!(nlp.internal, x, y2, vals[I], obj_weight=0.0)
  end
  vals[nnzhF+nnzhc+1:nnzhF+nnzhc+ne] .= obj_weight
  return vals
end

function hess_coord!(nlp :: FeasibilityFormNLS{LLSModel}, xr :: AbstractVector, y :: AbstractVector, vals :: AbstractVector;
                     obj_weight :: Float64=1.0)
  @lencheck nlp.meta.nvar xr
  @lencheck nlp.meta.ncon y
  @lencheck nlp.meta.nnzh vals
  increment!(nlp, :neval_hess)
  vals .= obj_weight
  return vals
end

function hess(nlp :: FeasibilityFormNLS, xr :: AbstractVector; obj_weight :: Float64=1.0)
  @lencheck nlp.meta.nvar xr
  increment!(nlp, :neval_hess)
  n, ne = nlp.internal.meta.nvar, nlp.internal.nls_meta.nequ
  return [spzeros(n, n + ne); spzeros(ne, n) obj_weight * I]
end

function hess(nlp :: FeasibilityFormNLS, xr :: AbstractVector, y :: AbstractVector;
              obj_weight :: Float64=1.0)
  @lencheck nlp.meta.nvar xr
  @lencheck nlp.meta.ncon y
  increment!(nlp, :neval_hess)
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  @views Hx = m > 0 ? hess(nlp.internal, x, y[ne+1:end], obj_weight=0.0) : spzeros(n, n)
  Hx += hess_residual(nlp.internal, x, @view y[1:ne])
  return [Hx spzeros(n, ne); spzeros(ne, n) obj_weight * I]
end

function hprod!(nlp :: FeasibilityFormNLS, xr :: AbstractVector, y :: AbstractVector, v :: AbstractVector, hv :: AbstractVector;
    obj_weight :: Float64=1.0)
  @lencheck nlp.meta.nvar xr v hv
  @lencheck nlp.meta.ncon y
  n, m, ne = nlp.internal.meta.nvar, nlp.internal.meta.ncon, nlp.internal.nls_meta.nequ
  x = @view xr[1:n]
  if m > 0
    @views hprod!(nlp.internal, x, y[ne+1:end], v[1:n], hv[1:n], obj_weight=0.0)
  else
    fill!(hv, 0.0)
  end
  for i = 1:ne
    @views hv[1:n] .+= hprod_residual(nlp.internal, x, i, v[1:n]) * y[i]
  end
  @views hv[n+1:end] .= obj_weight * v[n+1:end]
  return hv
end

function residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, Fx :: AbstractVector)
  @lencheck nlp.meta.nvar x
  @lencheck nlp.nls_meta.nequ Fx
  increment!(nlp, :neval_residual)
  n = nlp.internal.meta.nvar
  Fx .= @view x[n+1:end]
  return Fx
end

function jac_structure_residual!(nlp :: FeasibilityFormNLS, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  @lencheck nlp.nls_meta.nnzj rows cols
  n, ne = nlp.internal.meta.nvar, nlp.internal.nls_meta.nequ
  rows .= 1:ne
  cols .= n .+ (1:ne)
  return rows, cols
end

function jac_coord_residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, vals :: AbstractVector)
  @lencheck nlp.meta.nvar x
  @lencheck nlp.nls_meta.nnzj vals
  increment!(nlp, :neval_jac_residual)
  vals[1:nlp.nls_meta.nnzj] .= 1
  return vals
end

function jprod_residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  @lencheck nlp.meta.nvar x v
  @lencheck nlp.nls_meta.nequ Jv
  increment!(nlp, :neval_jprod_residual)
  n = nlp.internal.meta.nvar
  Jv .= @view v[n+1:end]
  return Jv
end

function jtprod_residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  @lencheck nlp.meta.nvar x Jtv
  @lencheck nlp.nls_meta.nequ v
  increment!(nlp, :neval_jtprod_residual)
  n, ne = nlp.internal.meta.nvar, nlp.internal.nls_meta.nequ
  Jtv[1:n] .= 0.0
  Jtv[n+1:end] .= v
  return Jtv
end

function hess_structure_residual!(nlp :: FeasibilityFormNLS, rows :: AbstractVector{<: Integer}, cols :: AbstractVector{<: Integer})
  @lencheck nlp.nls_meta.nnzh rows cols
  return rows, cols  # Hessian of residual is zero; do not change rows and cols
end

function hess_coord_residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, v :: AbstractVector, vals :: AbstractVector)
  @lencheck nlp.meta.nvar x
  @lencheck nlp.nls_meta.nequ v
  @lencheck nlp.nls_meta.nnzh vals
  increment!(nlp, :neval_hess_residual)
  return vals
end

function jth_hess_residual(nlp :: FeasibilityFormNLS, x :: AbstractVector, i :: Int)
  @lencheck nlp.meta.nvar x
  increment!(nlp, :neval_jhess_residual)
  n = nlp.meta.nvar
  return spzeros(n, n)
end

function hprod_residual!(nlp :: FeasibilityFormNLS, x :: AbstractVector, i :: Int, v :: AbstractVector, Hiv :: AbstractVector)
  @lencheck nlp.meta.nvar x v Hiv
  increment!(nlp, :neval_hprod_residual)
  fill!(Hiv, 0.0)
  return Hiv
end


text:
module TestCopulas

using Copulas
using Base.Test

epsTol = 1e-12

###############################
## CDVine consistency checks ##
###############################

## load CDVine test results
##-------------------------

dataPath = joinpath(Pkg.dir("Copulas"),
                    "test/data/r_cop_funcs_results.csv")

testResults = readcsv(dataPath, header = true)[1]
resDict = [testResults[ii, 1] => testResults[ii, 2]
           for ii=1:size(testResults, 1)]

function tName(func, cop, copInd, ptInd)
    return string(func, "_", cop, "_", copInd, "_", ptInd)
end

## tests using metaprogramming
##----------------------------

testPoints = [(0.3, 0.4),
              (0.2, 0.8)]

funcs = [:pdf,
         :cdf,
         :hfun]

cops = []

for ptInd=1:length(testPoints)
    for copInd=1:size(cops, 1)
        for f in funcs
            copName = cops[copInd, 1]
            cop = cops[copInd, 2]
            u1, u2 = testPoints[ptInd]
            mod(copInd, 2) == 0 ? ind = 2 : ind = 1
            testName = tName(string(f), copName, ind, ptInd)
#            display(testName)
            eval(:(val = (Copulas.$f)(cop, u1, u2)))
            expOut = resDict[testName]
            @test_approx_eq_eps val expOut epsTol
        end
    end
end

## v functions
for ptInd=1:length(testPoints)
    for copInd=1:size(cops, 1)
        copName = cops[copInd, 1]
        cop = cops[copInd, 2]
        u1, u2 = testPoints[ptInd]
        mod(copInd, 2) == 0 ? ind = 2 : ind = 1
        testName = tName("vfun", copName, ind, ptInd)
        #            display(testName)
        val = Copulas.vfun(cop, u2, u1)
        expOut = resDict[testName]
        @test_approx_eq_eps val expOut epsTol
    end
end

########################
## consistency checks ##
########################

epsTol = 1e-7

uLow = 0.001
uHigh = 0.999
testGrid = [uLow, 0.1:0.1:0.9, uHigh]

cops = [
        Copulas.AsymFGMPC(0.2),
        Copulas.AsymFGMPC(0.9)
        ]

## pdf, cdf, h- and v-functions
##-----------------------------

for cop in cops
    show(cop)
    println("")
    println(" *  cdf tests:")
    Copulas.testCdf(cop, epsTol, testGrid)
    println(" *  h/v tests:")
    Copulas.testHandV(cop, epsTol, testGrid)
    println(" *  pdf tests:")
    Copulas.testPdf(cop, testGrid)
end

## inverse h- and v-functions
##---------------------------

cops = []

testPoints = [(0.3, 0.7),
              (0.00002, 0.00098),
              (0.1, 0.998),
              ## (0.8, 0.2132323543),
              (0.8, 0.2132),
              ## (0.2349248, 0.983475),
              (0.999, 0.999),
              (0.0001, 0.0001)]

println("-------------------------------")
println("inverse h- and v-function tests")
println("-------------------------------")

for pts in testPoints
    u1, u2 = pts
    for cop in cops
        copShown = false

        q = Copulas.hfun(cop, u1, u2)
        u1_hat = Copulas.hinv(cop, q, [u2])
        if abs(u1 - u1_hat)[1] > epsTol
            # show copula
            if !copShown
                show(cop)
                println("")
                copShown = true
            end

            # at least both values should be mapped to same hfun value
            println("u2=$u2: $u1 -> $u1_hat")
            if !(q == Copulas.hfun(cop, u1_hat, u2))
                println("different h-function values")
                h1 = Copulas.hfun(cop, u1, u2)
                h2 = Copulas.hfun(cop, u1_hat, u2)
                println("  h($u1) = $h1")
                println("  h($u1_hat) = $h2")
            end
        end

        q = Copulas.vfun(cop, u1, u2)
        u1_hat = Copulas.vinv(cop, q, [u2])
        if abs(u1 - u1_hat)[1] > epsTol
            # show copula
            if !copShown
                show(cop)
                println("")
                copShown = true
            end

            # at least both values should be mapped to same hfun value
            println("u2=$u2: $u1 -> $u1_hat")
            if !(q == Copulas.vfun(cop, u1_hat, u2))
                println("different v-function values")
                h1 = Copulas.vfun(cop, u1, u2)
                h2 = Copulas.vfun(cop, u1_hat, u2)
                println("  v($u1) = $h1")
                println("  v($u1_hat) = $h2")
            end
        end

        ## q = Copulas.hinv(cop, u1, u2)
        ## u1_hat = Copulas.hfun(cop, q, [u2])
        ## hdiff2 = u1 - u1_hat
        ## ## @test_approx_eq_eps u1 u1_hat epsTol
        ## q = Copulas.vfun(cop, u2, u1)
        ## u2_hat = Copulas.vinv(cop, q, [u1])
        ## vdiff1 = u2 - u2_hat
        ## q = Copulas.vfun(cop, u2, u1)
        ## u2_hat = Copulas.vinv(cop, q, [u1])
        ## vdiff2 = u2 - u2_hat
        ## ## @test_approx_eq_eps u2 u2_hat epsTol
        ## maxDiff = max(abs(hdiff1[1]),
        ##               abs(hdiff2[1]),
        ##               abs(vdiff1[1]),
        ##               abs(vdiff2[1]))
        ## if abs(maxDiff) > 1e-10
        ##     show(cop)
        ##     println("")
        ##     println("diffh1: $hdiff1")
        ##     println("diffh2: $hdiff2")
        ##     println("diffv1: $vdiff1")
        ##     println("diffv2: $vdiff2")
        ## end
    end
end

######################
## parameter bounds ##
######################

@test_throws Exception Copulas.AsymFGMPC(4.3)
@test_throws Exception Copulas.ClaytonPC(-4.)

## params function
##----------------

@test (0.8,) == Copulas.params(Copulas.AsymFGMPC(0.8))
@test (4.3,) == Copulas.params(Copulas.ClaytonPC(4.3))

end


text:
# # SCOPE Model

## load packages
using Land.CanopyLayers
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

# ## Initialization
# Besides the individual functions to initialize parameters for `CanopyLayers`,
#     a general function is provided to initialize all the parameters directly.
angles, can, can_opt, can_rad, in_rad, leaves, rt_con, rt_dim, soil, wls =
    initialize_rt_module(FT; nLayer=20, LAI=3);
#------------------------------------------------------------------------------

# ## Steps
## 1. Update canopy optical properties (required)
canopy_geometry!(can, angles, can_opt, rt_con);
## 2. Update scattering coefficients (required)
canopy_matrices!(leaves, can_opt);
## 3. Simulate short wave simulation (required)
short_wave!(can, can_opt, can_rad, in_rad, soil, rt_con);
## 4. Update integrated radiation fluxes (required for photosynthesis)
canopy_fluxes!(can, can_opt, can_rad, in_rad, soil, leaves, wls, rt_con);
## 5. Update SIF related spectrum (required for SIF)
SIF_fluxes!(leaves, can_opt, can_rad, can, soil, wls, rt_con, rt_dim);
## 6. Update thermo fluxes (required for leaf energy budget)
thermal_fluxes!(leaves, can_opt, can_rad, can, soil, [FT(400.0)], wls);

_fig,_axes = create_canvas("SIF example"; ncol=2);
_ax1,_ax2 = _axes;
_ax1.plot(wls.WL , can_rad.alb_obs, "k-");
_ax2.plot(wls.WLF, can_rad.SIF_obs, "k-");
set_xlabels!(_axes, ["Wave length (nm)" for i in 1:2], fontsize=12);
set_ylabels!(_axes, ["Albedo", "obs SIF (mW m⁻² nm⁻¹ sr⁻¹)"], fontsize=12);
_fig.set_tight_layout(true);
_fig
#------------------------------------------------------------------------------

# ## Change FQE
for leaf in leaves
    leaf.Cx = 0.5;
    leaf.fqe = 0.004;
    fluspect!(leaf, wls);
end
can.Ω = 0.48
canopy_geometry!(can, angles, can_opt, rt_con);
canopy_matrices!(leaves, can_opt);
short_wave!(can, can_opt, can_rad, in_rad, soil, rt_con);
canopy_fluxes!(can, can_opt, can_rad, in_rad, soil, leaves, wls, rt_con);
SIF_fluxes!(leaves, can_opt, can_rad, can, soil, wls, rt_con, rt_dim);
thermal_fluxes!(leaves, can_opt, can_rad, can, soil, [FT(400.0)], wls);

_fig,_axes = create_canvas("SIF example FQE"; ncol=2);
_ax1,_ax2 = _axes;
_ax1.plot(wls.WL , can_rad.alb_obs, "k-");
_ax2.plot(wls.WLF, can_rad.SIF_obs, "k-");
set_xlabels!(_axes, ["Wave length (nm)" for i in 1:2], fontsize=12);
set_ylabels!(_axes, ["Albedo", "obs SIF (mW m⁻² nm⁻¹ sr⁻¹)"], fontsize=12);
_fig.set_tight_layout(true);
_fig
#------------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#=
lbfgs.jl

    Quasi-Newton type update step routines using the limited-memory
    Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm, that uses approximations
    to the Hessian using past approximations as well as the gradient.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/01/31
=#

# Method
Base.@kwdef mutable struct LBFGS{T}
    m::T = 10           # memory depth
    pseudo_iter::T = 1  # pseudo iteration counter
end
# limited-memory BFGS state
mutable struct LBFGSState{Tv, Tf, Tm}
    x::Tv       # current state
    x_prev::Tv  # previous state
    ∇f_prev::Tv # previous gradient of f(x)
    f_prev::Tf  # previous f(x)
    s::Tv       # change in state
    y::Tv       # change in gradient
    s_mem::Tm   # memory change in states
    y_mem::Tm   # memory change in gradients
    r::Tv       # search direction
    ρ::Tv       # two-loop ρ
    α::Tv       # buffer two-loop α
    q::Tv       # buffer two-loop q
end

"""
    twoloop!(state, m, pseudo_iter)

Compute descent direction using the two-loop recursion.

#### Arguments
  - `state::LBFGSState`     : state variables
  - `m::Integer`            : memory depth
  - `psuedo_iter::Integer`  : pseudo iteration counter
"""
function twoloop!(state::LBFGSState, m::Integer, pseudo_iter::Integer)
    # Upper and lower bounds
    lower= max(pseudo_iter - m, 1) # project on non-negative orthant
    upper= pseudo_iter - 1

    # Initalize q at previous gradient
    copyto!(state.q, state.∇f_prev)

    # Backward pass to update q
    @inbounds for i in upper:-1:lower
        idx= mod1(i,m)  # convert index to [1,m] range

        # Store views
        s_i= view(state.s_mem, :, idx)
        y_i= view(state.y_mem, :, idx)

        state.α[idx]= state.ρ[idx] * dot(s_i, state.q)
        state.q.-= state.α[idx] .* y_i
    end

    # Initialize search direction r
    if pseudo_iter > 1
        idx= mod1(upper,m)  # convert index to [1,m] range

        # Store views
        s_i= view(state.s_mem, :, idx)
        y_i= view(state.y_mem, :, idx)

        γ= dot(s_i, y_i) * inv(sum(abs2, y_i))
        state.r.= γ .* state.q
    else
        state.r.= state.q
    end

    # Forward pass to update search direction r
    @inbounds for i in lower:upper
        idx= mod1(i,m)  # convert index to [1,m] range

        # Store views
        s_i= view(state.s_mem, :, idx)
        y_i= view(state.y_mem, :, idx)

        β= state.ρ[idx] * dot(y_i, state.r)
        state.r.+= s_i .* (state.α[idx] - β)
    end

    # Negate search direction
    rmul!(state.r, eltype(state.r)(-1))

    return nothing
end

"""
    update_state!(state, method, ls, f)

Update state using the limited-memory BFGS inverse Hessian approximation for the
search direction with an inexact backtracking linesearch, storing the result in
`x`.

#### Arguments
  - `state::LBFGSState` : state variables
  - `ls::BackTrack`     : line search parameters
  - `ls::BackTrack`     : line search parameters
  - `f::Function`       : ``f(x)``
"""
function update_state!(state::LBFGSState, method::LBFGS, ls::BackTrack,
                        f::Function)
    # Store current state
    state.x_prev.= state.x

    # Search direction
    twoloop!(state, method.m, method.pseudo_iter)

    # Backtracking line Search
    state.f_prev= backtrack!(ls, state.x, state.x_prev, state.r, state.f_prev,
                            state.∇f_prev, f)

    return nothing
end

"""
    update_memory!(state, method)

Update memory of the limited-memory BFGS algorithm to construct the inverse
Hessian approximation for the search direction.

#### Arguments
  - `state::LBFGSState` : state variables
  - `method::LBFGS`     : method info
"""
function update_memory!(state::LBFGSState, method::LBFGS)
    # Inverse curvature condition
    ρ= inv(dot(state.y, state.s))

    if isinf(ρ)
        # Erase memory, restart
        method.pseudo_iter= 1
        # Exit
        return nothing
    end

    idx= mod1(method.pseudo_iter,method.m)  # convert index to [1,m] range
    # Update memory
    state.s_mem[:,idx].= state.s
    state.y_mem[:,idx].= state.y
    state.ρ[idx]= ρ

    # Update pseudo iteration
    method.pseudo_iter+= 1

    return nothing
end

"""
    lbfgs!(method, state, ls, f, ∇f!)

Update state and inverse Hessian approximation using the limited-memory BFGS
algorithm.

#### Arguments
  - `method::LBFGS`     : method info
  - `state::LBFGSState` : state variables
  - `ls::BackTrack`     : line search parameters
  - `f::Function`       : ``f(x)``
  - `∇f!::Function`     : gradient of `f`
"""
function lbfgs!(method::LBFGS, state::LBFGSState, ls::BackTrack, f::Function,
                ∇f!::Function)
    # Update state
    update_state!(state, method, ls, f)

    # Change in state
    state.s.= ls.α .* state.r

    # Store previous gradient
    state.y.= state.∇f_prev
    # Update gradient
    ∇f!(state.∇f_prev, state.x)
    # Change in gradient
    state.y.= state.∇f_prev .- state.y

    # Update memory
    update_memory!(state, method)

    return nothing
end


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

function random(d::UniformDistribution, n::Int)
    return rand(n) .* (d.b - d.a) .+ d.a
end

function random(d::NormalDistribution, n::Int)
    z = randn(n)
    return d.sigma .* z .+ d.mu
end

function random(d::CauchyDistribution, n::Int)
    return d.c .* tan.( pi.*( rand(n) .- 0.5) ) .+ d.mu
end

function random(d::MultiNormalDistribution, n::Int)
    n_dim = length(d.mu)
    z = randn(n_dim, n)
    return d.L' * z .+ d.mu
end

function rejection_sampling_1d(pdf_fcn, x_min, x_max, n_rand)

    # select a majorizing constant
    n_for_y_max = max( 100, convert(Int, ceil(sqrt(n_rand)) ) )
    xs          = collect( range(x_min, x_max, n_for_y_max) )
    pdf_s       = pdf_fcn.(xs)
    max_pdf     = maximum(pdf_s)

    # coordinate distribution
    x_distr     = UniformDistribution(x_min, x_max)

    n_batch     = div(n_rand, 10)
    x_rnb_total = zeros(n_rand)
    left_idx    = 1
    right_idx   = 0

    while right_idx < n_rand

        # check if the sampled point is below the curve
        x_rnb       = random( x_distr, n_batch )
        pdf_rnb     = max_pdf * rand(n_batch)
        pdf_values  = pdf_fcn.(x_rnb)
        bool        = pdf_rnb .<= pdf_values

        # select the coordinates
        x_rnb = x_rnb[bool]

        # push back
        L = length(x_rnb)
        right_idx += L
        if right_idx > n_rand
            valid_idx = L - (right_idx - n_rand)
            right_idx = n_rand
            x_rnb_total[ left_idx:right_idx ] = x_rnb[1:valid_idx]
        else
            x_rnb_total[ left_idx:right_idx ] = x_rnb
        end
        left_idx = right_idx + 1
    end

    return x_rnb_total

end

function rejection_sampling_2d(pdf_fcn, x_min, x_max, y_min, y_max, n_rand)

    # select a majorizing constant
    n_for_y_max = max( 50, convert(Int, ceil(sqrt(n_rand)) ) )
    xs          = collect( range(x_min, x_max, n_for_y_max) )
    ys          = collect( range(y_min, y_max, n_for_y_max) )

    pdf_s = zeros(n_for_y_max, n_for_y_max)
    for ii = 1:n_for_y_max
        x_ = xs[ii]
        for jj = 1:n_for_y_max
            y_ = ys[jj]
            pdf_s[ii, jj] = pdf_fcn( (x_, y_) )
        end
    end
    max_pdf     = maximum(pdf_s)

    # coordinate distribution
    x_distr     = UniformDistribution(x_min, x_max)
    y_distr     = UniformDistribution(y_min, y_max)

    n_batch     = div(n_rand, 10)

    x_rnb_total = zeros(n_rand)
    y_rnb_total = zeros(n_rand)

    left_idx    = 1
    right_idx   = 0

    while right_idx < n_rand

        # check if the sampled point is below the curve
        x_rnb       = random( x_distr, n_batch )
        y_rnb       = random( y_distr, n_batch )

        pdf_rnb     = max_pdf * rand(n_batch)

        pdf_values = zeros(n_batch)
        for ii = 1:n_batch
            pdf_values[ii] = pdf_fcn( (x_rnb[ii], y_rnb[ii]) )
        end

        bool = pdf_rnb .<= pdf_values

        # select the coordinates
        x_rnb = x_rnb[bool]
        y_rnb = y_rnb[bool]

        # push back
        L = length(x_rnb)
        right_idx += L
        if right_idx > n_rand
            valid_idx = L - (right_idx - n_rand)
            right_idx = n_rand
            x_rnb_total[ left_idx:right_idx ] = x_rnb[1:valid_idx]
            y_rnb_total[ left_idx:right_idx ] = y_rnb[1:valid_idx]
        else
            x_rnb_total[ left_idx:right_idx ] = x_rnb
            y_rnb_total[ left_idx:right_idx ] = y_rnb
        end
        left_idx = right_idx + 1
    end

    return (x_rnb_total, y_rnb_total)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
export markDiff, getQdiff!, getNdiff!
import StatsBase: mean

#### Structures #####
# Define supertypes of differentiable blocks.
const CphaseGate{N, T} = ControlBlock{N, <:ShiftGate{T}, <:Any}
const DiffBlock{N, T} = Union{RotationGate{N, T, <:Any}, CphaseGate{N, T}}

# Define the struct of a differentiable gate `QDiff{GT, N}`.
"""
    QDiff{GT, N} <: TagBlock{GT, N}
Differentiable gate. If you want to mark or return a differentiable block, use `markDiff` or `markDiff!` instead.
\n
\nFields:
\n`block::GT`: Sub-block(s) of the differentiable block.
\n`grad::Float64`: Gradient(s) of the differentiable block.
"""
mutable struct QDiff{GT, N} <: TagBlock{GT, N}
    block::GT
    grad::Float64
    mat
    QDiff(blk::DiffBlock{N}) where {N} = new{typeof(blk), N}(blk, 0, mat(blk))
end

#### Functions #####
# Make Marks of differentiable gates when printing circuit layout.
function YaoBlocks.print_annotation(io::IO, df::QDiff)
    printstyled(io, "[̂∂(MPSC)] "; bold=true, color=:yellow)
end

# Define basic Yao-compatible functions for QDiff{GT, N}.
Yao.content(cb::QDiff) = cb.block
Yao.apply!(reg::AbstractRegister, db::QDiff) = apply!(reg, Yao.content(db))
Yao.mat(::Type{T}, df::QDiff) where T = Yao.mat(T, df.block)
Base.adjoint(df::QDiff) = QDiff(Yao.content(df)')

# Convert gates that are differentiable into type `QDiff{GT, N}`.
"""
    markDiff(block::AbstractBlock) -> block::AbstractBlock
Return the differentiable gate(s) `QMPS.QDiff{GT, N}` from a block or a block tree such as `ChainBlock`.
"""
## For Block trees.
function markDiff(blk::AbstractBlock)
    blks = subblocks(blk)
    isempty(blks) ? blk : chsubblocks(blk, markDiff.(blks))
end
## For differentiable blocks.
markDiff(block::DiffBlock) = QDiff(block)
## Exclude control blocks.
markDiff(block::ControlBlock) = block

"""
    getQdiff!(psifunc::Function, diffGate::Union{QMPS.QDiff, Array{QMPS.QDiff,1}}, op::AbstractBlock) -> grad::Union{Float64, Array{Float64, 1}}
Quantum Operator differentiation. When only apply `getQdiff!` to one differentiable gate, the output satisfies `grad == diffGate.grad`.
\n `psifunc = ()-> reg::ArrayReg |> c::ChainBlock`
\n `diffGate = collect_blocks(QMPS.QDiff, c)`
\n `op`: Witness Operator to measure reg.
\n
\nNOTE:
\n1) `getQdiff!` modifies only the element `grad` inside the field of `QDiff` and all other elements remain the same.
\n2) `getQdiff!` does indirectly modify the input qubit `reg` to `reg |> c` when inputting the `psifunc` as an arguement.
"""
@inline function getQdiff!(psifunc::Function, diffGate::QDiff, op::AbstractBlock)
    r1, r2 = _perturb( ()->mean( Yao.expect(op, psifunc()) ) |> real, diffGate, π/2 )
    diffGate.grad = (r2 - r1)/2
end
@inline function getQdiff!(psifunc::Function, diffGates::Array{QDiff,1}, op::AbstractBlock)
    grads = getQdiff!.(psifunc, diffGates, Ref(op))
end

"""
    getNdiff!(psifunc::Function, parGate::Union{AbstractBlock, Array{AbstractBlock,1}}, op::AbstractBlock; δ::Real=0.01) -> grad::Union{Float64, Array{Float64, 1}}
Numerical Operator differentiation. When only apply `getNdiff!!` to one differentiable gate, the output satisfies `grad == parGate.grad`.
\n `psifunc = ()-> reg::ArrayReg |> c::ChainBlock`
\n `parGate`: Paramterized block(gate) in c.
\n `op`: Witness Operator to measure reg.
\n
\nNOTE:
\n1) When applying `getNdiff!` to QMPS-DC or similar differentiable MPS circuits, the accuracy of the function can get much worse if you set `δ` to be too small (normally when < 1e-1).
\n2) `getNdiff!` indirectly modifies the input qubit `reg` to `reg |> c` when inputting the `psifunc` as an arguement.
"""
@inline function getNdiff!(psifunc::Function, parGate::AbstractBlock, op::AbstractBlock; δ::Real=0.01)
    r1, r2 = _perturb( ()->mean( Yao.expect(op, psifunc()) ) |> real, parGate, δ )
    grad = (r2 - r1) / (2δ)
end
@inline function getNdiff!(psifunc::Function, parGates::Array{<:AbstractBlock,1}, op::AbstractBlock; δ::Real=0.01)
    grads = getNdiff!.(psifunc, parGates, Ref(op), δ=δ)
end

@inline function _perturb(func, gate::AbstractBlock, δ::Real)
    dispatch!(-, gate, (δ,))
    r1 = func()
    dispatch!(+, gate, (2δ,))
    r2 = func()
    dispatch!(-, gate, (δ,))
    r1, r2
end


text:
module ColormapParams

using Colors
using Main.PNG16x16

struct ColormapParamSVG <: Main.SVG
    buf::IOBuffer
end

function ColormapParamSVG(target::Symbol)
    io = IOBuffer()
    write(io,
        """
        <svg xmlns="http://www.w3.org/2000/svg"
             xmlns:xlink="http://www.w3.org/1999/xlink"
             version="1.1"
             width="45mm" height="45mm"
             viewBox="0 0 270 270" preserveAspectRatio="none"
             shape-rendering="crispEdges" stroke="none" stroke-width="1">
        """)
    write(io,
        """
        <rect width="270" height="270" fill="#aaa" opacity="0.2" />
        """)
    if target in (:c, :s, :b)
        write_chart_csb(io, target)
    elseif target in (:w, :d)
        write_chart_wd(io, target)
    end
    write(io, "</svg>")
    ColormapParamSVG(io)
end

function write_chart_csb(io, target::Symbol)
    pp = (0.1, 0.5, 0.9)
    pc = map(_ -> 0.88, pp)
    ps = map(_ -> 0.60, pp)
    pb = map(_ -> 0.75, pp)
    markers = ("m -5,-5 v 10 h 10 v -10 z",
               "m 0,-5 l -5,9 h 10 z",
               "m 0,-4 l -4,4 4,4 4,-4 z")
    if target === :c
        col = Colors.JULIA_LOGO_COLORS.purple
        label = "contrast"
        pc = pp
    elseif target === :s
        col = Colors.JULIA_LOGO_COLORS.red
        label = "saturation"
        ps = pp
    elseif target === :b
        col = Colors.JULIA_LOGO_COLORS.green
        label = "brightness"
        pb = pp
    end
    h = round(Int, hue(convert(Luv, col)))
    n = 16
    ls = (100/2/n):(100/n):100
    cs = (175/2/n):(175/n):175
    mc(l) = Colors.find_maximum_chroma(LCHuv(l, 0.0, h))
    plane = [LCHuv(100.0 - l, min(c, mc(l)), h) for l in ls, c in cs]
    write(io,
        """
        <g transform="translate(80, 20)" fill="none" stroke-width="1.5"
           style="font-size: 18px; font-style: italic;">
            <image width="100" height="100" transform="scale(1.75, 2)" xlink:href=\"""")
    write_png_as_data(io, plane)
    write(io, "\" />\n")
    mc_coords = join(string(round(Int, mc(100 - l)), ",", 2l, " ") for l = 0:2:100)
    write(io,
        """
            <path d="m 175,0 L $mc_coords 175,200 z" fill="#aaa" />
            <path d="m 0,200 v -215 l 5,10 m -10,0 l 5,-10" stroke="currentColor" />
            <path d="m 0,200 h  185 l -10,5 m 0,-10 l 10,5" stroke="currentColor" />
            <text fill="currentColor" x="150" y="220">C*</text>
            <text fill="currentColor" x="-25" y="10" >L*</text>
        """)
    write(io, "</g>")
    for i = 1:3
        colors= sequential_palette(h, 7, c=pc[i], s=ps[i], b=pb[i], w=0, d=0)
        write_colormap(io, 25 * i - 20, pp[i], colors, markers[i])
    end
    tcolors = (sequential_palette(h, 25, c=pc[i], s=ps[i], b=pb[i], w=0, d=0) for i = 1:3)
    write_stroke(io, tcolors, markers)
    write(io,
        """
        <text fill="currentColor" y="260" style="font-size: 18px;">
            <tspan x="70" style="font-style: italic;">$target</tspan>
            <tspan x="90">- $label [0,1]</tspan>
        </text>
        """)
end

function write_chart_wd(io, target::Symbol)
    pp = 0.0:0.2:1.0
    pw = map(_ -> 0.15, pp)
    pd = map(_ -> 0.0, pp)

    if target === :w
        col = Colors.JULIA_LOGO_COLORS.blue
        cy = 40
        pw = pp
        kcol = "#ff0"
    elseif target === :d
        col = Colors.JULIA_LOGO_COLORS.green
        cy = 200
        pd = pp
        kcol = "#00f"
    end
    h = round(Int, hue(convert(Luv, col)))
    for i in eachindex(pp)
        colors = sequential_palette(h, 9, w=pw[i], d=pd[i])
        write_colormap(io, 30 * i - 10, pp[i], colors)
    end
    write(io,
        """
        <rect x="205" y="$cy" width="20" height="20" fill="$kcol" />
        """)
    write(io,
        """
        <text fill="currentColor" y="260" style="font-size: 18px;">
            <tspan x="20" style="font-style: italic;">$target</tspan>
            <tspan x="40">- strength of $(target)color [0,1]</tspan>
            <tspan x="195" y="$(cy-5)" style="font-style: italic;">$(target)color</tspan>
        </text>
        """)
end

function write_colormap(io, x, p, colors, marker="")
    y = 0
    write(io,
        """
        <g transform="translate($x, 40)" stroke-opacity="0.5">
        """)
    for col in colors
        write(io,
            """
                <rect x="0" y="$y" width="20" height="20" fill="#$(hex(col))" />
            """)
        y += 20
    end
    if !isempty(marker)
        y += 10
        write(io,
            """
                <path d="m 10,$y $marker" fill="#$(hex(last(colors)))" stroke="white" />
            """)
    end
    y += 20
    write(io,
        """
            <text x="0" y="$y" fill="currentColor" style="font-size: 15px;">$p</text>
        """)
    write(io, "</g>")
end

function write_stroke(io, tcolors, markers)
    write(io,
        """
        <g transform="translate(80, 20)" stroke="white" stroke-opacity="0.5">
        """)
    coord(c) = string(round(Int, c.c), ",", round(Int, 2 * (100 - c.l)), " ")
    for (i, colors) in enumerate(tcolors)
        coords = join(map(coord, LCHuv.(colors)))
        op = (0.4, 0.6, 1.0)[i]
        sw = (10, 4, 1)[i]
        write(io,
            """
                <path d="M $coords" fill="none" stroke-opacity="$op" stroke-width="$sw" />
            """)
    end
    for (i, colors) in enumerate(tcolors)
        mk = markers[i]
        hx = hex(last(colors))
        write(io,
            """
            <path d="M $(coord(LCHuv(first(colors))))$mk" fill="#$hx" />
            <path d="M $(coord(LCHuv(last( colors))))$mk" fill="#$hx" />
            """)
    end
    write(io, "</g>")
end

end # module


text:
# N = 10 # 2520
N = 20 # 232792560
function isPrime(n::Int)
	for i in 2:ceil(sqrt(n))
		if n%i==0 && i!=n return false end
	end
	return true
end
primes = [x for x in 1:N if isPrime(x)]
print(primes)

function factorize(n::Int)
	f = Int[]
	while n > 1
		for i in 2:n
			if n%i==0
				n /= i
				push!(f, i)
				break
			end
		end
	end
	f
end
#println(factorize(2520))

p = Dict((i, 0) for i in primes)
for i in 1:N
	s = factorize(i)
	c = Dict((i,0) for i in primes)
	for j in s
		c[j] += 1
	end
	for (k,v) in c
		p[k] = max(p[k], v)
	end
end
println(p)
prod_p = 1
for (k, v) in p
	prod_p *= k^v
end
println(N, '\t', prod_p)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
# Autogenerated wrapper script for Qt5Svg_jll for x86_64-w64-mingw32-cxx11
export libqt5svg

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Svg")
JLLWrappers.@declare_library_product(libqt5svg, "Qt5Svg.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5svg,
        "bin\\Qt5Svg.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export glscreen

function _is_alive(x::WeakRef)
    if isa(x.value, Screen)
        isopen(x.value) && return true
    end
    false
end

let screen_list = WeakRef[]
    global current_screen, add_screen, get_screens, empty_screens!
    clean!() = filter!(_is_alive, screen_list)
    function current_screen()
        clean!()
        if isempty(screen_list)
            error("No screen available. Consider creating one with the function glscreen()")
        end
        last(screen_list).value
    end
    add_screen(screen) = push!(screen_list, WeakRef(screen))
    get_screens() = (clean!(); map(x->x.value, screen_list))
    empty_screens!() = empty!(screen_list)
end

function cleanup()
    for screen in get_screens()
        destroy!(screen)
    end
    empty_screens!()
end

struct Millimeter
end
global const mm = Millimeter()

import Base: *
const pixel_per_mm = Ref(3.7)

function (*)(x::Millimeter, y::Number)
    round(Int, y * pixel_per_mm[])
end
function (*)(x::Number, y::Millimeter)
    round(Int, x * pixel_per_mm[])
end

function get_scaled_dpi(window)
    monitor = GLFW.GetPrimaryMonitor()
    props = GLWindow.MonitorProperties(monitor)

    # it seems like small displays with high dpi make mm look quite big.
    # so lets scale it a bit. 518 is a bit arbitrary, but the scale of my
    # screen on which I test everything, hence it will make you see things as I do.
    scaling = props.physicalsize[1] / 518
    min(props.dpi...) * scaling # we do not start fiddling with differently scaled xy dpi's
end

function glscreen(name = "GLVisualize";
        resolution = GLWindow.standard_screen_resolution(),
        debugging = false,
        clear = true,
        color = RGBA(1,1,1,1),
        stroke = (0f0, color),
        hidden = false,
        visible = true,
        focus = true,
        fullscreen = false
    )
    cleanup()
    screen = Screen(
        name,
        resolution = resolution, debugging = debugging,
        clear = clear, color = color, stroke = stroke,
        hidden = hidden, visible = visible, focus = focus,
        fullscreen = fullscreen
    )
    add_screen(screen)
    GLWindow.add_complex_signals!(screen) #add the drag events and such
    GLFW.MakeContextCurrent(GLWindow.nativewindow(screen))
    pixel_per_mm[] = get_scaled_dpi(screen) / 25.4
    screen
end

const timer_signal_dict = Dict{WeakRef, Dict{Int, Signal{Float64}}}()
"""
Creates a timer signal with `updates_per_second` while `window` is open.
It's reusing timer signals with the same update rate and registering the updates
with GLFW.
"""
function get_timer_signal(updates_per_second, window=current_screen())
    dict = get!(timer_signal_dict, WeakRef(window), Dict{Int, Signal{Float64}}())
    get!(dict, updates_per_second) do
        fpswhen(window.inputs[:window_open], updates_per_second)
    end
end

function fold_loop(v0, _)
    val, range, index = v0
    val = range[index]
    index += 1
    index > length(range) && (index = 1)
    (val, range, index)
end

function loop(range::Range, rate=60)
    t = get_timer_signal(rate)
    map(first, foldp(fold_loop, (first(range), range, 1), t))
end

function fold_bounce(v0, _)
    val, range, index, direction = v0
    val = range[index]
    index += direction
    if index in (length(range)+1, 0)
        direction = -direction
        index += 2direction
    end
    (val, range, index, direction)
end

function bounce(range::Range{T}, rate=60) where T
    t = get_timer_signal(rate)
    map(first, foldp(fold_bounce, (first(range), range, 1, 1), t))
end

export bounce, loop


text:
export LevenbergMarquardt_BC_Result

@doc raw"""

```julia
struct LevenbergMarquardt_BC_Result{T<:Real} <:  Abstract_BC_Solver_Result
    ...
end
```

This structure subtypes [`Abstract_BC_Solver_Result`](@ref)

"""
const LevenbergMarquardt_BC_Result = LevenbergMarquardt_Result{T} where {T<:Real}


text:
include("solvers/explore.jl")

function main_explore(args)
    println("="^70)
    println("Début de l'action explore")
    inst = Instance(args[:infile])

    sv = ExploreSolver(inst)
    itermax_default = 50*inst.nb_planes
    itermax = args[:itermax] == 0 ? itermax_default : args[:itermax]

    println("Starting solving")
    ms_start = ms() # seconde depuis le démarrage avec précision à la ms
    costs = solve(sv, itermax)
    ms_stop = ms()

    bestsol = sv.bestsol
    print_sol(bestsol)

    nb_calls = bestsol.solver.nb_calls
    nb_infeasable = bestsol.solver.nb_infeasable
    nb_sec = round(ms_stop - ms_start, digits=3)
    nb_call_per_sec = round(nb_calls/nb_sec, digits=3)
    println("Performance: ")
    println("  nb_calls=$nb_calls")
    println("  nb_infeasable=$nb_infeasable")
    println("  nb_sec=$nb_sec")
    println("  => nb_call_per_sec = $nb_call_per_sec call/sec")

    println("Fin de l'action explore")
    return costs
end

function loop_explore(args, instances::Vector{Instance})
    costs = []
    for instance in instances
        args[:infile] = "data/$instance"
        push!(costs, main_explore(args))
    end
    return costs
end

costs = main_explore(Args.args)