text:
using Transpiler
using Base.Test
import Transpiler: cli, CLMethod
using Sugar: getsource!, dependencies!

function mapkernel(f, a, b, c)
    gid = cli.get_global_id(0) + Cuint(1)
    c[gid] = f(a[gid], b[gid])
    return
end

function test_source(target, result)
    source_equal = target == result
    if source_equal
        @test true
    else
        @test false
        println("source unequal:\ntarget:\n", target)
        println("result:\n", result)
    end
end
# empty caches
Transpiler.empty_caches!()

args = (typeof(+), cli.GlobalPointer{Float32}, cli.GlobalPointer{Float32}, cli.GlobalPointer{Float32})
cl_mapkernel = CLMethod((mapkernel, args))
source = Sugar.getsource!(cl_mapkernel)
mapsource = """void mapkernel_1(Base123 f, __global float *  a, __global float *  b, __global float *  c);
void mapkernel_1(Base123 f, __global float *  a, __global float *  b, __global float *  c)
{
    uint gid;
    gid = get_global_id(0) + (uint)(1);
    float _ssavalue_0;
    _ssavalue_0 = (a)[gid - 0x00000001] + (b)[gid - 0x00000001];
    (c)[gid - 0x00000001] = _ssavalue_0;
    return;
}"""

@testset "map kernel" begin
    test_source(mapsource, source)
    deps = dependencies!(cl_mapkernel, true)
    deps_test = [
        Int64,
        UInt32,
        (+, Tuple{UInt32, UInt32}),
        (cli.get_global_id, Tuple{Int64}),
        (+, Tuple{Float32, Float32}),
        Float32,
        (-, Tuple{UInt32, UInt32}),
        typeof(+),
        cli.GlobalPointer{Float32},
        (UInt32, Tuple{Int64}),
        typeof(mapkernel),
    ]
    @test length(deps) == length(deps_test)
    for elem in deps
        @test elem.signature in deps_test
    end
end

#Broadcast
Base.@propagate_inbounds broadcast_index(arg, shape, i) = arg
Base.@propagate_inbounds function broadcast_index{T, N}(
        arg::cli.GlobalPointer{T}, shape::NTuple{N, Integer}, i
    )
    @inbounds return arg[i]
end

# The implementation of prod in base doesn't play very well with current
# transpiler. TODO figure out what Core._apply maps to!
_prod{T}(x::NTuple{1, T}) = x[1]
_prod{T}(x::NTuple{2, T}) = x[1] * x[2]

function broadcast_kernel(A, f, sz, arg1, arg2)
    i = cli.get_global_id(0) + Cuint(1)
    @inbounds if i <= _prod(sz)
        A[i] = f(
            broadcast_index(arg1, sz, i),
            broadcast_index(arg2, sz, i),
        )
    end
    return
end

args = (cli.GlobalPointer{Float32}, typeof(+), Tuple{UInt32}, cli.GlobalPointer{Float32}, Float32)
cl_mapkernel = CLMethod((broadcast_kernel, args))
source = getsource!(cl_mapkernel)
broadcastsource = """void broadcast_kernel_5(__global float *  A, Base123 f, uint sz, __global float *  arg1, float arg2);
void broadcast_kernel_5(__global float *  A, Base123 f, uint sz, __global float *  arg1, float arg2)
{
    uint i;
    i = get_global_id(0) + (uint)(1);
    if(i <= _prod_2(sz)){
        float _ssavalue_0;
        _ssavalue_0 = broadcast_index_3(arg1, sz, i) + broadcast_index_4(arg2, sz, i);
        (A)[i - 0x00000001] = _ssavalue_0;
    };
    return;
}"""

@testset "broadcast kernel" begin
    test_source(broadcastsource, source)

    deps = dependencies!(cl_mapkernel, true)
    deps_test = [
        typeof(broadcast_kernel),
        UInt32,
        Int,
        (+,Tuple{UInt32,UInt32}),
        (cli.get_global_id, Tuple{Int64}),
        (<=,Tuple{UInt32,UInt32}),
        (_prod,Tuple{Tuple{UInt32}}),
        (+,Tuple{Float32,Float32}),
        (broadcast_index,Tuple{cli.GlobalPointer{Float32},Tuple{UInt32},UInt32}),
        (broadcast_index,Tuple{Float32,Tuple{UInt32},UInt32}),
        (-,Tuple{UInt32,UInt32}),
        cli.GlobalPointer{Float32},
        typeof(+),
        Tuple{UInt32},
        Float32,
        (UInt32, Tuple{Int64}),
        Bool,
        typeof(_prod),
        typeof(broadcast_index)
    ]
    @test length(deps) == length(deps_test)
    for elem in deps
        @test elem.signature in deps_test
    end
end

@testset "for" begin
    decl = CLMethod((fortest, (Float32,)))
    source = Sugar.getsource!(decl)
    #TODO remove xxtempx4, which is unused now...
    target_source = """float fortest_6(float x);
    float fortest_6(float x)
    {
        float acc;
        long x2temp2;
        long i;
        acc = x;
        for(i = 1; i <= 5; i++){
            if(i == 1){
                acc = acc + x;
                continue;
            };
            if(i == 2){
                acc = acc - x;
                continue;
            };
            acc = acc + x * x;
        };
        return acc;
    }"""
    test_source(target_source, source)
end

function custom_index_test(x)
    x[1, 1]
end

@testset "custom getindex" begin
    source, method, name = Transpiler.kernel_source(custom_index_test, (typeof(1f0*I),))
    source_compare = """// Inbuilds
    typedef char JLBool;
    // dependencies
    // #custom_index_test
    __constant int FUNC_INST_x2custom_index_test = 0;
    typedef int x2custom_index_test; // empty type emitted as an int
    // UniformScaling{Float32}
    struct __attribute__((packed)) TYPUniformScaling_float{
        float __attribute__((aligned (4))) x4;
    };
    typedef struct TYPUniformScaling_float UniformScaling_float;

    // Base.#getindex
    __constant int FUNC_INST_Base12getindex = 0;
    typedef int Base12getindex; // empty type emitted as an int
    // Transpiler.CLIntrinsics.#cl_select
    __constant int FUNC_INST_Transpiler1CLIntrinsics12cl_select = 0;
    typedef int Transpiler1CLIntrinsics12cl_select; // empty type emitted as an int
    // (Transpiler.CLIntrinsics.cl_select, Tuple{Float32,Float32,Bool})
    float cl_select_9(float a, float b, JLBool c);
    float cl_select_9(float a, float b, JLBool c)
    {
        return select(b, a, (uint)(c));
    }
    // Symbol

    // Base.#zero
    __constant int FUNC_INST_Base12zero = 0;
    typedef int Base12zero; // empty type emitted as an int
    // Base.#oftype
    __constant int FUNC_INST_Base12oftype = 0;
    typedef int Base12oftype; // empty type emitted as an int
    // (oftype, Tuple{Float32,Int64})
    float oftype_10(float x, long c);
    float oftype_10(float x, long c)
    {
        return (float){c};
    }
    // (zero, Tuple{Float32})
    float zero_6(float x);
    float zero_6(float x)
    {
        return oftype_10(x, 0);
    }
    // (getindex, Tuple{UniformScaling{Float32},Int64,Int64})
    float getindex_7(UniformScaling_float J, long i, long j);
    float getindex_7(UniformScaling_float J, long i, long j)
    {
        return cl_select_9(J.λ, zero_6(J.λ), i == j);
    }
    // ########################
    // Main inner function
    // (custom_index_test, (UniformScaling{Float32},))
    __kernel float custom_index_test_8(UniformScaling_float x);
    float custom_index_test_8(UniformScaling_float x)
    {
        return getindex_7(x, 1, 1);
    }
    """
    test_source(source_compare, source)
end

inner(i) = Float32(i) * 77f0
function ntuple_test(::Val{N}) where N
    ntuple(inner, Val{N})
end
@testset "ntuple" begin
    source, method, name = Transpiler.kernel_source(ntuple_test, (Val{4},))
    compare_source = """// Inbuilds
    typedef char JLBool;
    // dependencies
    // #ntuple_test
    __constant int FUNC_INST_x2ntuple_test = 0;
    typedef int x2ntuple_test; // empty type emitted as an int
    // Val{4}
    typedef int Val_4; // empty type emitted as an int
    // #inner
    __constant int FUNC_INST_x2inner = 0;
    typedef int x2inner; // empty type emitted as an int
    // Type{Val{4}}
    typedef int Type5Val5466; // placeholder type instance
    __constant Type5Val5466 TYP_INST_Type5Val5466 = 0;

    // Base.#ntuple
    __constant int FUNC_INST_Base12ntuple = 0;
    typedef int Base12ntuple; // empty type emitted as an int
    // Any
    typedef int Any; // placeholder type instance
    __constant Any TYP_INST_Any = 0;

    // (inner, Tuple{Int64})
    float inner_13(long i);
    float inner_13(long i)
    {
        return (float)(i) * 77.0f;
    }
    // (ntuple, Tuple{#inner,Type{Val{4}}})
    float4 ntuple_11(x2inner f, Type5Val5466 x2unused2);
    float4 ntuple_11(x2inner f, Type5Val5466 x2unused2)
    {
        return (float4){inner_13(1), inner_13(2), inner_13(3), inner_13(4)};
    }
    // ########################
    // Main inner function
    // (ntuple_test, (Val{4},))
    __kernel float4 ntuple_test_12(Val_4 x2unused2);
    float4 ntuple_test_12(Val_4 x2unused2)
    {
        return ntuple_11(FUNC_INST_x2inner, TYP_INST_Type5Val5466);
    }
    """
    test_source(compare_source, source)
end

function testifelse(a, b)
    ifelse(a == b, a, b)
end
@testset "ifelse" begin
    source, method, name = Transpiler.kernel_source(testifelse, (Int, Int))
    testsource = """// Inbuilds
    typedef char JLBool;
    // dependencies
    // #testifelse
    __constant int FUNC_INST_x2testifelse = 0;
    typedef int x2testifelse; // empty type emitted as an int
    // Transpiler.CLIntrinsics.#cl_select
    __constant int FUNC_INST_Transpiler1CLIntrinsics12cl_select = 0;
    typedef int Transpiler1CLIntrinsics12cl_select; // empty type emitted as an int
    // (Transpiler.CLIntrinsics.cl_select, Tuple{Int64,Int64,Bool})
    long cl_select_14(long a, long b, JLBool c);
    long cl_select_14(long a, long b, JLBool c)
    {
        return select(b, a, (ulong)(c));
    }
    // ########################
    // Main inner function
    // (testifelse, (Int64, Int64))
    __kernel long testifelse_15(long a, long b);
    long testifelse_15(long a, long b)
    {
        return cl_select_14(a, b, a == b);
    }
    """
    test_source(testsource, source)
end

function testfastmath(a::Complex64)
    @fastmath exp(a)
end
@testset "fastmath" begin
    source, method, name = Transpiler.kernel_source(testfastmath, (Complex64,))
    testsource = """// Inbuilds
    typedef char JLBool;
    // dependencies
    // #testfastmath
    __constant int FUNC_INST_x2testfastmath = 0;
    typedef int x2testfastmath; // empty type emitted as an int
    // Complex{Float32}
    struct __attribute__((packed)) TYPComplex_float{
        float __attribute__((aligned (4))) re;
        float __attribute__((aligned (4))) im;
    };
    typedef struct TYPComplex_float Complex_float;

    // Base.FastMath.#exp_fast
    __constant int FUNC_INST_Base1FastMath12exp_fast = 0;
    typedef int Base1FastMath12exp_fast; // empty type emitted as an int
    // Base.FastMath.#mul_fast
    __constant int FUNC_INST_Base1FastMath12mul_fast = 0;
    typedef int Base1FastMath12mul_fast; // empty type emitted as an int
    // Any
    typedef int Any; // placeholder type instance
    __constant Any TYP_INST_Any = 0;

    // Type{Float32}
    typedef int Type5Float326; // placeholder type instance
    __constant Type5Float326 TYP_INST_Type5Float326 = 0;

    // (Complex{Float32}, Tuple{Float32,Float32})
    Complex_float x7Complex_float8_17(float re, float im);
    Complex_float x7Complex_float8_17(float re, float im)
    {
        return (Complex_float){re, im};
    }
    // Base.#real
    __constant int FUNC_INST_Base12real = 0;
    typedef int Base12real; // empty type emitted as an int
    // Symbol

    // (real, Tuple{Complex{Float32}})
    float real_16(Complex_float z);
    float real_16(Complex_float z)
    {
        return z.re;
    }
    // Base.#imag
    __constant int FUNC_INST_Base12imag = 0;
    typedef int Base12imag; // empty type emitted as an int
    // (imag, Tuple{Complex{Float32}})
    float imag_16(Complex_float z);
    float imag_16(Complex_float z)
    {
        return z.im;
    }
    // (Base.FastMath.mul_fast, Tuple{Float32,Complex{Float32}})
    Complex_float mul_fast_18(float a, Complex_float y);
    Complex_float mul_fast_18(float a, Complex_float y)
    {
        return x7Complex_float8_17(a * real_16(y), a * imag_16(y));
    }
    // Base.FastMath.#cis_fast
    __constant int FUNC_INST_Base1FastMath12cis_fast = 0;
    typedef int Base1FastMath12cis_fast; // empty type emitted as an int
    // (Base.FastMath.cis_fast, Tuple{Float32})
    Complex_float cis_fast_6(float x);
    Complex_float cis_fast_6(float x)
    {
        return x7Complex_float8_17(cos(x), sin(x));
    }
    // (Base.FastMath.exp_fast, Tuple{Complex{Float32}})
    Complex_float exp_fast_16(Complex_float x);
    Complex_float exp_fast_16(Complex_float x)
    {
        return mul_fast_18(exp(real_16(x)), cis_fast_6(imag_16(x)));
    }
    // ########################
    // Main inner function
    // (testfastmath, (Complex{Float32},))
    __kernel Complex_float testfastmath_16(Complex_float a);
    Complex_float testfastmath_16(Complex_float a)
    {
        return exp_fast_16(a);
    }
    """
    test_source(testsource, source)
end


text:
using COSMO
using Logging

DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 900000
)

"""
simulate returns normal
"""
function returns_montecarlo(Σ, r̄, numS)
    d = MvNormal(r̄, Σ)
    r = rand(d, numS)
    P = pdf(d, r)
    return r', P / sum(P)
end

## Prep data
function compute_solution(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    r = sum(r̄'w_values)
    return w_values, objective_value(model), r
end

function compute_solution_dual(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    return w_values, objective_value(model), value(model[:E])
end

function compute_solution_stoc(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    Cvar = -objective_value(model)
    r = sum(r̄'w_values)
    q1_α = value.(model[:z])
    return w_values, r, Cvar, q1_α
end

function compute_solution_stoc_2(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    r = objective_value(model)
    q1_α = value.(model[:z])
    return w_values, r, q1_α
end

ext: jl


text:
using Gadfly

n = 10
plot(x=rand(n), y=rand(n), color=rand(n),
     Scale.color_continuous(minvalue=-10, maxvalue=10))


text:
using Literate
Literate.markdown("foo.jl", "."; documenter=true)
Literate.notebook("foo.jl", "."; documenter=true)

max_stars_repo_name: christophernhill/clima-impldiff-testing

max_issues_repo_name: christophernhill/clima-impldiff-testing

max_forks_repo_name: christophernhill/clima-impldiff-testing


text:
#=
linear system solver subroutines

6x6 nonsymmetric system in (x, y, z, tau, s, kap):
A'*y + G'*z + c*tau = xrhs
-A*x + b*tau = yrhs
-G*x + h*tau - s = zrhs
-c'*x - b'*y - h'*z - kap = taurhs
(pr bar) z_k + mu*H_k*s_k = srhs_k
(du bar) mu*H_k*z_k + s_k = srhs_k
mu/(taubar^2)*tau + kap = kaprhs
=#

# calculate direction given rhs, and apply iterative refinement
function get_directions(
    stepper::Stepper{T},
    solver::Solver{T},
    min_impr_tol::T = T(0.5) # improvement tolerance for iterative refinement
    ) where {T <: Real}
    rhs = stepper.rhs
    dir = stepper.dir
    dir_temp = stepper.dir_temp
    res = stepper.temp
    syssolver = solver.syssolver
    res_norm_cutoff = solver.res_norm_cutoff
    max_ref_steps = solver.max_ref_steps

    solve_system(syssolver, solver, dir, rhs)

    iszero(max_ref_steps) && return dir

    # compute residual norm
    copyto!(dir_temp, dir.vec)
    apply_lhs(stepper, solver) # modifies res
    res.vec .-= rhs.vec
    res_norm = norm(res.vec, Inf)

    # use iterative refinement if residual norm exceeds cutoff
    if res_norm > res_norm_cutoff
        is_prev_slow = false
        prev_res_norm = res_norm

        for i in 1:max_ref_steps
            # compute refined direction
            solve_system(syssolver, solver, dir, res)
            axpby!(true, dir_temp, -1, dir.vec)

            # compute residual
            apply_lhs(stepper, solver) # modifies res
            res.vec .-= rhs.vec

            res_norm_new = norm(res.vec, Inf)
            if res_norm_new >= res_norm
                # residual has not improved
                copyto!(dir.vec, dir_temp)
                break
            end

            # residual has improved
            copyto!(dir_temp, dir.vec)
            res_norm = res_norm_new

            (res_norm < res_norm_cutoff) && break
            is_curr_slow = (res_norm > min_impr_tol * prev_res_norm)
            is_prev_slow && is_curr_slow && break # slow progress

            prev_res_norm = res_norm
            is_prev_slow = is_curr_slow
        end
    end

    @assert !isnan(res_norm) # TODO error instead
    solver.worst_dir_res = max(solver.worst_dir_res, res_norm)

    return dir
end

# calculate residual on 6x6 linear system
function apply_lhs(
    stepper::Stepper{T},
    solver::Solver{T},
    ) where {T <: Real}
    model = solver.model
    dir = stepper.dir
    res = stepper.temp
    tau_dir = dir.tau[]
    kap_dir = dir.kap[]

    # A'*y + G'*z + c*tau
    copyto!(res.x, model.c)
    mul!(res.x, model.G', dir.z, true, tau_dir)
    # -G*x + h*tau - s
    @. res.z = model.h * tau_dir - dir.s
    mul!(res.z, model.G, dir.x, -1, true)
    # -c'*x - b'*y - h'*z - kap
    res.tau[] = -dot(model.c, dir.x) - dot(model.h, dir.z) - kap_dir
    # if p = 0, ignore A, b, y
    if !iszero(model.p)
        # A'*y + G'*z + c*tau
        mul!(res.x, model.A', dir.y, true, true)
        # -A*x + b*tau
        copyto!(res.y, model.b)
        mul!(res.y, model.A, dir.x, -1, tau_dir)
        # -c'*x - b'*y - h'*z - kap
        res.tau[] -= dot(model.b, dir.y)
    end

    # s
    for (k, cone_k) in enumerate(model.cones)
        # (du bar) mu*H_k*z_k + s_k
        # (pr bar) z_k + mu*H_k*s_k
        s_res_k = res.s_views[k]
        Cones.hess_prod_slow!(s_res_k, dir.primal_views[k], cone_k)
        @. s_res_k += dir.dual_views[k]
    end

    tau = solver.point.tau[]
    res.kap[] = solver.mu / tau * tau_dir / tau + kap_dir

    return res
end

include("naive.jl")
include("naiveelim.jl")
include("symindef.jl")
include("qrchol.jl")

# reduce to 4x4 subsystem
function solve_system(
    syssolver::Union{NaiveElimSystemSolver{T}, SymIndefSystemSolver{T},
        QRCholSystemSolver{T}},
    solver::Solver{T},
    sol::Point{T},
    rhs::Point{T},
    ) where {T <: Real}
    model = solver.model

    solve_subsystem4(syssolver, solver, sol, rhs)
    tau = sol.tau[]

    # lift to get s and kap
    # s = -G*x + h*tau - zrhs
    @. sol.s = model.h * tau - rhs.z
    mul!(sol.s, model.G, sol.x, -1, true)

    # kap = -mu/taubar^2*tau + kaprhs
    taubar = solver.point.tau[]
    sol.kap[] = -solver.mu / taubar / taubar * tau + rhs.kap[]

    return sol
end

# reduce to 3x3 subsystem
function solve_subsystem4(
    syssolver::Union{SymIndefSystemSolver{T}, QRCholSystemSolver{T}},
    solver::Solver{T},
    sol::Point{T},
    rhs::Point{T},
    ) where {T <: Real}
    model = solver.model
    rhs_sub = syssolver.rhs_sub
    sol_sub = syssolver.sol_sub

    @. rhs_sub.x = rhs.x
    @. rhs_sub.y = -rhs.y
    setup_rhs3(syssolver, model, rhs, sol, rhs_sub)

    solve_subsystem3(syssolver, solver, sol_sub, rhs_sub)

    # lift to get tau
    sol_const = syssolver.sol_const
    tau_num = rhs.tau[] + rhs.kap[] + dot_obj(model, sol_sub)
    taubar = solver.point.tau[]
    tau_denom = solver.mu / taubar / taubar - dot_obj(model, sol_const)
    sol_tau = tau_num / tau_denom

    dim3 = length(sol_sub.vec)
    @. sol.vec[1:dim3] = sol_sub.vec + sol_tau * sol_const.vec
    sol.tau[] = sol_tau

    return sol
end

function setup_point_sub(
    syssolver::Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}},
    model::Models.Model{T},
    ) where {T <: Real}
    (n, p, q) = (model.n, model.p, model.q)
    dim_sub = n + p + q
    z_start = model.n + model.p

    sol_sub = syssolver.sol_sub = Point{T}()
    rhs_sub = syssolver.rhs_sub = Point{T}()
    rhs_const = syssolver.rhs_const = Point{T}()
    sol_const = syssolver.sol_const = Point{T}()
    for point_sub in (sol_sub, rhs_sub, rhs_const, sol_const)
        point_sub.vec = zeros(T, dim_sub)
        @views point_sub.x = point_sub.vec[1:n]
        @views point_sub.y = point_sub.vec[n .+ (1:p)]
        @views point_sub.z = point_sub.vec[n + p .+ (1:q)]
        point_sub.z_views = [view(point_sub.z, idxs) for idxs in model.cone_idxs]
    end
    @. rhs_const.x = -model.c
    @. rhs_const.y = model.b
    @. rhs_const.z = model.h

    return
end

dot_obj(model::Models.Model, point::Point) = dot(model.c, point.x) +
    dot(model.b, point.y) + dot(model.h, point.z)


text:
using PackageCompiler

# for PILS example
create_sysimage(
		[
		 "PyCall",
		 "FlightSims",
		 "FSimBase",
		 "FSimZoo",
		 "FSimROS",
		 "FSimPlots",
		 "Plots",
		 "UnPack",
		 "OnlineStats",
		 "DataFrames",
		 "ComponentArrays",
		 ],
		sysimage_path="sys_pils.so",
		precompile_execution_file="test/PILS/precompile_execution.jl",
	       )


text:
"""
In a graph with numbered edges, provide access to neighbors *and edge numbers*
for a specified node.

!!! note

    [`IncidenceIndex`](@ref) always includes both directions, i.e.
    for an edge `(u, v)` with edge index `ei`, it is
    `(u, ei) ∈ neighbors[v]` and `(v, ei) ∈ neighbors[u]`.
"""
struct IncidenceIndex
    idx::Vector{Int}
    pi::Vector{Tuple{Int, Int}}
end

@inline Base.getindex(c::IncidenceIndex, j::Int) =
    view(c.pi, c.idx[j]:c.idx[j+1]-1)

"""
    num_edges(::IncidenceIndex)

Actual number of (undirected) edges (**not counting** both directions).
"""
num_edges(c::IncidenceIndex) = length(c.pi) ÷ 2

num_nodes(c::IncidenceIndex) = length(c.idx) - 1

Base.length(c::IncidenceIndex) = num_nodes(c)

"""
    IncidenceIndex(n, head, tail)

Construct an index for `n` nodes of the (undirected) neighbors given
by the edges `zip(head, tail)`.
If `head[i] <= 0` then edge `i` is excluded.
"""
function IncidenceIndex(n::Int, head::Vector{Int}, tail::Vector{Int})
    @assert length(head) == length(tail)
    IncidenceIndex(n, () -> zip(head, tail))
end

"""
    IncidenceIndex(n, edges)

Same but for edges like `[(1, 2), (2, 3)]`
"""
IncidenceIndex(n::Int, edges::Vector{Tuple{Int, Int}}) =
    IncidenceIndex(n, () -> edges)

IncidenceIndex(n::Int, edges::Vector{Pair{Int, Int}}) =
    IncidenceIndex(n, () -> edges)

"""
    IncidenceIndex(n, m, iter::Function)

Provide an iterator over the edges.
If the first node of an edge is zero, skip that edge.
"""
function IncidenceIndex(n::Int, iter::Function)
    local m::Int = 0
    idx = zeros(Int, n+1)
    for (h::Int, t::Int) in iter()
        h <= 0 && continue
        m += 1
        idx[h] += 1
        idx[t] += 1
    end
    acc = 1                        # accumulate degrees ==> positions
    deg_i = 0
    deg_ii = idx[1]
    for i = 1:n
        idx[i] = acc
        acc += deg_i
        deg_i, deg_ii = deg_ii, idx[i+1]
    end
    idx[n+1] = acc
    pi = Vector{Tuple{Int,Int}}(undef, 2m)
    @assert(idx[end] + deg_i == 2m + 1,
            "idx[$(length(idx))]: $(idx[end] + deg_i) != $(2m + 1)")
    for (i, (u, v)) in enumerate(iter())
        u <= 0 && continue
        pi[idx[u+1]] = (v, i)
        idx[u+1] += 1
        pi[idx[v+1]] = (u, i)
        idx[v+1] += 1
    end
    @assert(idx[end] == 2m + 1, "$(idx[end]) vs $(2m + 1)")
    return IncidenceIndex(idx, pi)
end

import Base.==

==(a::IncidenceIndex, b::IncidenceIndex) =
    a.pi == b.pi && a.idx == b.idx

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
#===================================================================================================
  Composition Classes: Valid kernel transformations
===================================================================================================#

abstract CompositionClass{T<:AbstractFloat}

@inline eltype{T}(::CompositionClass{T}) = T

@inline iscomposable(::CompositionClass, ::Kernel) = false

@inline ismercer(::CompositionClass) = false
@inline isnegdef(::CompositionClass) = false

@inline attainszero(::CompositionClass)     = true
@inline attainspositive(::CompositionClass) = true
@inline attainsnegative(::CompositionClass) = true

function description_string(ϕ::CompositionClass, showtype::Bool = true)
    class = typeof(ϕ)
    fields = fieldnames(class)
    class_str = string(class.name.name) * (showtype ? string("{", eltype(ϕ), "}") : "")
    *(class_str, "(", join(["$field=$(getfield(ϕ,field).value)" for field in fields], ","), ")")
end

function show(io::IO, ϕ::CompositionClass)
    print(io, description_string(ϕ))
end

function convert{T<:AbstractFloat,K<:CompositionClass}(::Type{CompositionClass{T}}, ϕ::K)
    convert(K.name.primary{T}, ϕ)
end

function =={T<:CompositionClass}(ϕ1::T, ϕ2::T)
    all([getfield(ϕ1,field) == getfield(ϕ2,field) for field in fieldnames(T)])
end

#== Positive Mercer Classes ==#

abstract PositiveMercerClass{T<:AbstractFloat} <: CompositionClass{T}
@inline ismercer(::PositiveMercerClass) = true
@inline attainsnegative(::PositiveMercerClass) = false
@inline attainszero(::PositiveMercerClass) = false

doc"GammaExponentialClass(κ;α,γ) = exp(-α⋅κᵞ)"
immutable GammaExponentialClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    alpha::HyperParameter{T}
    gamma::HyperParameter{T}
    GammaExponentialClass(α::Variable{T}, γ::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open)),
        HyperParameter(γ, Interval(Bound(zero(T), :open), Bound(one(T), :closed)))
    )
end
@outer_constructor(GammaExponentialClass, (1,0.5))
@inline iscomposable(::GammaExponentialClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::GammaExponentialClass{T}, z::T) = exp(-ϕ.alpha * z^ϕ.gamma)

doc"ExponentialClass(κ;α) = exp(-α⋅κ²)"
immutable ExponentialClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    alpha::HyperParameter{T}
    ExponentialClass(α::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open))
    )
end
@outer_constructor(ExponentialClass, (1,))
@inline iscomposable(::ExponentialClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::ExponentialClass{T}, z::T) = exp(-ϕ.alpha * z)

doc"GammaRationalClass(κ;α,β,γ) = (1 + α⋅κᵞ)⁻ᵝ"
immutable GammaRationalClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    alpha::HyperParameter{T}
    beta::HyperParameter{T}
    gamma::HyperParameter{T}
    GammaRationalClass(α::Variable{T}, β::Variable{T}, γ::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open)),
        HyperParameter(β, leftbounded(zero(T), :open)),
        HyperParameter(γ, Interval(Bound(zero(T), :open), Bound(one(T), :closed)))
    )
end
@outer_constructor(GammaRationalClass, (1,1,0.5))
@inline iscomposable(::GammaRationalClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::GammaRationalClass{T}, z::T) = (1 + ϕ.alpha*z^ϕ.gamma)^(-ϕ.beta)

doc"RationalClass(κ;α,β,γ) = (1 + α⋅κ)⁻ᵝ"
immutable RationalClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    alpha::HyperParameter{T}
    beta::HyperParameter{T}
    RationalClass(α::Variable{T}, β::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open)),
        HyperParameter(β, leftbounded(zero(T), :open))
    )
end
@outer_constructor(RationalClass, (1,1))
@inline iscomposable(::RationalClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::RationalClass{T}, z::T) = (1 + ϕ.alpha*z)^(-ϕ.beta)

doc"MatérnClass(κ;ν,ρ) = 2ᵛ⁻¹(√(2ν)κ/ρ)ᵛKᵥ(√(2ν)κ/ρ)/Γ(ν)"
immutable MaternClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    nu::HyperParameter{T}
    rho::HyperParameter{T}
    MaternClass(ν::Variable{T}, ρ::Variable{T}) = new(
        HyperParameter(ν, leftbounded(zero(T), :open)),
        HyperParameter(ρ, leftbounded(zero(T), :open))
    )
end
@outer_constructor(MaternClass, (1,1))
@inline iscomposable(::MaternClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline function phi{T<:AbstractFloat}(ϕ::MaternClass{T}, z::T)
    v1 = sqrt(2ϕ.nu) * z / ϕ.rho
    v1 = v1 < eps(T) ? eps(T) : v1  # Overflow risk, z -> Inf
    2 * (v1/2)^(ϕ.nu) * besselk(ϕ.nu, v1) / gamma(ϕ.nu)
end

doc"ExponentiatedClass(κ;α) = exp(a⋅κ + c)"
immutable ExponentiatedClass{T<:AbstractFloat} <: PositiveMercerClass{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    ExponentiatedClass(a::Variable{T}, c::Variable{T}) = new(
        HyperParameter(a, leftbounded(zero(T), :open)),
        HyperParameter(c, leftbounded(zero(T), :closed))
    )
end
@outer_constructor(ExponentiatedClass, (1,0))
@inline iscomposable(::ExponentiatedClass, κ::Kernel) = ismercer(κ)
@inline phi{T<:AbstractFloat}(ϕ::ExponentiatedClass{T}, z::T) = exp(ϕ.a*z + ϕ.c)

#== Other Mercer Classes ==#

doc"PolynomialClass(κ;a,c,d) = (a⋅κ + c)ᵈ"
immutable PolynomialClass{T<:AbstractFloat,U<:Integer} <: CompositionClass{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    d::HyperParameter{U}
    PolynomialClass(a::Variable{T}, c::Variable{T}, d::Variable{U}) = new(
        HyperParameter(a, leftbounded(zero(T), :open)),
        HyperParameter(c, leftbounded(zero(T), :closed)),
        HyperParameter(d, leftbounded(one(U),  :closed))
    )
end
@outer_constructor(PolynomialClass, (1,0,3))
@inline iscomposable(::PolynomialClass, κ::Kernel) = ismercer(κ)
@inline phi{T<:AbstractFloat}(ϕ::PolynomialClass{T}, z::T) = (ϕ.a*z + ϕ.c)^ϕ.d
@inline ismercer(::PolynomialClass) = true

#== Non-Negative Negative Definite Kernel Classes ==#

abstract NonNegNegDefClass{T<:AbstractFloat} <: CompositionClass{T}
@inline isnegdef(::NonNegNegDefClass) = true
@inline attainsnegative(::NonNegNegDefClass) = false

doc"PowerClass(z;a,c,γ) = (az + c)ᵞ"
immutable PowerClass{T<:AbstractFloat} <: NonNegNegDefClass{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    gamma::HyperParameter{T}
    PowerClass(a::Variable{T}, c::Variable{T}, γ::Variable{T}) = new(
        HyperParameter(a, leftbounded(zero(T), :open)),
        HyperParameter(c, leftbounded(zero(T), :closed)),
        HyperParameter(γ, Interval(Bound(zero(T), :open), Bound(one(T), :closed)))
    )
end
@outer_constructor(PowerClass, (1,0,0.5))
@inline iscomposable(::PowerClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::PowerClass{T}, z::T) = (ϕ.a*z + ϕ.c)^(ϕ.gamma)

doc"GammmaLogClass(z;α,γ) = log(1 + α⋅zᵞ)"
immutable GammaLogClass{T<:AbstractFloat} <: NonNegNegDefClass{T}
    alpha::HyperParameter{T}
    gamma::HyperParameter{T}
    GammaLogClass(α::Variable{T}, γ::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open)),
        HyperParameter(γ, Interval(Bound(zero(T), :open), Bound(one(T), :closed)))
    )
end
@outer_constructor(GammaLogClass, (1,0.5))
@inline iscomposable(::GammaLogClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::GammaLogClass{T}, z::T) = log(ϕ.alpha*z^(ϕ.gamma) + 1)

doc"LogClass(z;α) = log(1 + α⋅z)"
immutable LogClass{T<:AbstractFloat} <: NonNegNegDefClass{T}
    alpha::HyperParameter{T}
    LogClass(α::Variable{T}) = new(
        HyperParameter(α, leftbounded(zero(T), :open))
    )
end
@outer_constructor(LogClass, (1,))
@inline iscomposable(::LogClass, κ::Kernel) = isnegdef(κ) && isnonnegative(κ)
@inline phi{T<:AbstractFloat}(ϕ::LogClass{T}, z::T) = log(ϕ.alpha*z + 1)

#== Non-Mercer, Non-Negative Definite Classes ==#

doc"SigmoidClass(κ;α,c) = tanh(a⋅κ + c)"
immutable SigmoidClass{T<:AbstractFloat} <: CompositionClass{T}
    a::HyperParameter{T}
    c::HyperParameter{T}
    SigmoidClass(a::Variable{T}, c::Variable{T}) = new(
        HyperParameter(a, leftbounded(zero(T), :open)),
        HyperParameter(c, leftbounded(zero(T), :closed))
    )
end
@outer_constructor(SigmoidClass, (1,0))
@inline iscomposable(::SigmoidClass, κ::Kernel) = ismercer(κ)
@inline phi{T<:AbstractFloat}(ϕ::SigmoidClass{T}, z::T) = tanh(ϕ.a*z + ϕ.c)

#===================================================================================================
  Kernel Composition ψ = ϕ(κ(x,y))
===================================================================================================#

doc"KernelComposition(ϕ,κ) = ϕ∘κ"
immutable KernelComposition{T<:AbstractFloat} <: StandardKernel{T}
    phi::CompositionClass{T}
    kappa::PairwiseKernel{T}
    function KernelComposition(ϕ::CompositionClass{T}, κ::PairwiseKernel{T})
        iscomposable(ϕ, κ) || error("Kernel is not composable.")
        new(ϕ, κ)
    end
end
function KernelComposition{T<:AbstractFloat}(ϕ::CompositionClass{T}, κ::PairwiseKernel{T})
    KernelComposition{T}(ϕ, κ)
end

function convert{T<:AbstractFloat}(::Type{KernelComposition{T}}, κ::KernelComposition)
    KernelComposition(convert(CompositionClass{T}, κ.phi), convert(Kernel{T}, κ.kappa))
end

function description_string(κ::KernelComposition, showtype::Bool = true)
    obj_str = string("KernelComposition", showtype ? string("{", eltype(κ), "}") : "")
    class_str = description_string(κ.phi, false)
    kernel_str = description_string(κ.kappa, false)
    string(obj_str, "(phi=", class_str, ",kappa=", kernel_str, ")")
end

function ==(ψ1::KernelComposition, ψ2::KernelComposition)
    (ψ1.phi == ψ2.phi) && (ψ1.kappa == ψ2.kappa)
end

ismercer(κ::KernelComposition) = ismercer(κ.phi)
isnegdef(κ::KernelComposition) = isnegdef(κ.phi)

attainszero(κ::KernelComposition)     = attainszero(κ.phi)
attainspositive(κ::KernelComposition) = attainspositive(κ.phi)
attainsnegative(κ::KernelComposition) = attainsnegative(κ.phi)

#== Composition Kernels ==#

doc"GaussianKernel(α) = exp(-α⋅‖x-y‖²)"
function GaussianKernel{T<:AbstractFloat}(α::Argument{T} = 1.0)
    KernelComposition(ExponentialClass(α), SquaredDistanceKernel{T}())
end
SquaredExponentialKernel = GaussianKernel
RadialBasisKernel = GaussianKernel

doc"LaplacianKernel(α) = exp(α⋅‖x-y‖)"
function LaplacianKernel{T<:AbstractFloat}(α::Argument{T} = 1.0)
    KernelComposition(GammaExponentialClass(α, convert(T, 0.5)), SquaredDistanceKernel{T}())
end

doc"PeriodicKernel(α,p) = exp(-α⋅Σⱼsin²(p(xⱼ-yⱼ)))"
function PeriodicKernel{T<:AbstractFloat}(α::Argument{T} = 1.0, p::Argument{T} = convert(T, π))
    KernelComposition(ExponentialClass(α), SineSquaredKernel(p))
end

doc"RationalQuadraticKernel(α,β) = (1 + α⋅‖x-y‖²)⁻ᵝ"
function RationalQuadraticKernel{T<:AbstractFloat}(α::Argument{T} = 1.0, β::Argument{T} = one(T))
    KernelComposition(RationalClass(α, β), SquaredDistanceKernel{T}())
end

doc"MatérnKernel(ν,θ) = 2ᵛ⁻¹(√(2ν)‖x-y‖²/θ)ᵛKᵥ(√(2ν)‖x-y‖²/θ)/Γ(ν)"
function MaternKernel{T<:AbstractFloat}(ν::Argument{T} = 1.0, θ::Argument{T} = one(T))
    KernelComposition(MaternClass(ν, θ), SquaredDistanceKernel{T}())
end
MatérnKernel = MaternKernel

doc"PolynomialKernel(a,c,d) = (a⋅xᵀy + c)ᵈ"
function PolynomialKernel{T<:AbstractFloat,U<:Integer}(
        a::Argument{T} = 1.0,
        c::Argument{T} = one(T),
        d::Argument{U} = 3
    )
    KernelComposition(PolynomialClass(a, c, d), ScalarProductKernel{T}())
end

doc"LinearKernel(α,c,d) = a⋅xᵀy + c"
function LinearKernel{T<:AbstractFloat}(a::Argument{T} = 1.0, c::Argument{T} = one(T))
    KernelComposition(PolynomialClass(a, c, 1), ScalarProductKernel{T}())
end

doc"SigmoidKernel(α,c) = tanh(a⋅xᵀy + c)"
function SigmoidKernel{T<:Real}(a::Argument{T} = 1.0, c::Argument{T} = one(T))
    KernelComposition(SigmoidClass(a, c), ScalarProductKernel{T}())
end

#== Special Compositions ==#

∘(ϕ::CompositionClass, κ::Kernel) = KernelComposition(ϕ, κ)

function ^{T<:AbstractFloat}(κ::PairwiseKernel{T}, d::Integer)
    KernelComposition(PolynomialClass(one(T), zero(T), d), κ)
end

function ^{T<:AbstractFloat}(κ::PairwiseKernel{T}, γ::T)
    KernelComposition(PowerClass(one(T), zero(T), γ), κ)
end

function exp{T<:AbstractFloat}(κ::PairwiseKernel{T})
    KernelComposition(ExponentiatedClass(one(T), zero(T)), κ)
end

function tanh{T<:AbstractFloat}(κ::PairwiseKernel{T})
    KernelComposition(SigmoidClass(one(T), zero(T)), κ)
end


text:
using Weber
using DataFrames: DataFrame, writetable

# TODO: try with mobile pre

# HOW TO USE: Run this script, while recording its sound (e.g. using your phone,
# or Audacity). Then run the script analyze_audiotiming.jl to extract the peaks
# from the recorded audio, and determine the timing accuracy of audio playback.
setup_sound(buffer_size=16)
timing = 0.2abs.(randn(Float64,100)) + 0.11
beep = sound(attenuate(ramp(tone(1000,0.05)),25))

exp = Experiment(name = "audio_timing",latency = current_sound_latency())

setup(exp) do
  addbreak(instruct("Test is ready."))

  moments = map(timing) do delta
    moment(delta) do
      play(beep)
      record("sound")
    end
  end
  addtrial(show_cross(),moments)
end

run(exp)


text:
# diagm test
x = create_var(4, 1)
d = diagm(x)
test_expr(d, (x,))

# main diagonal
x = create_var(3, 4)
d = diag(x, 0)
test_expr(d, (x,))

# off diagonals
x = create_var(3, 4)
d = diag(x, 2)
test_expr(d, (x,))

x = create_var(3, 4)
d = diag(x, -2)
test_expr(d, (x,))

info("All diag tests passed")


text:
# Autogenerated wrapper script for cif_api_jll for x86_64-apple-darwin
export libcif

using ICU_jll
using SQLite_jll
JLLWrappers.@generate_wrapper_header("cif_api")
JLLWrappers.@declare_library_product(libcif, "@rpath/libcif.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(ICU_jll, SQLite_jll)
    JLLWrappers.@init_library_product(
        libcif,
        "lib/libcif.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BinaryBuilder

# Collection of sources required to build Glib
sources = [
    "https://ftp.gnome.org/pub/gnome/sources/glib/2.54/glib-2.54.2.tar.xz" =>
    "bb89e5c5aad33169a8c7f28b45671c7899c12f74caf707737f784d7102758e6c",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd glib-2.54.2/

# Provide answers to a few configure questions automatically
cat > glib.cache <<END
glib_cv_stack_grows=no
glib_cv_uscore=no
ac_cv_path_MSGFMT=$prefix/bin/msgfmt
END

./configure LDFLAGS=-L$prefix/lib CPPFLAGS=-I$prefix/include --enable-libmount=no --cache-file=glib.cache --prefix=$prefix --host=$target

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
    BinaryProvider.Linux(:x86_64, :glibc),
    BinaryProvider.Linux(:aarch64, :glibc),
    BinaryProvider.Linux(:armv7l, :glibc),
    BinaryProvider.Linux(:powerpc64le, :glibc),
    BinaryProvider.MacOS(),
    # gettext doesn't build on windows yet :(
    #BinaryProvider.Windows(:i686),
    #BinaryProvider.Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libglib", :libglib)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # We need zlib
    "https://github.com/staticfloat/ZlibBuilder/releases/download/v1.2.11-3/build.jl",
    # We need libffi
    "https://github.com/staticfloat/libffiBuilder/releases/download/v3.2.1-0/build.jl",
    # We need gettext
    "https://github.com/staticfloat/GettextBuilder/releases/download/v0.19.8-0/build.jl",
    # We need pcre
    "https://github.com/staticfloat/PcreBuilder/releases/download/v8.41-0/build.jl",
]

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


text:
function squeezesum(A, dim)
    dropdims(sum(A, dims=dim), dims=dim)
end

function boundbelow!(x, tol=1e-10)
    for i in eachindex(x)
        if x[i] < tol
            x[i] = tol
        end
    end
end

function proportionalize(Y::Matrix, tol=1e-10)
    ϕ = Y ./ sum(Y, dims=2)
    for i in eachindex(ϕ)
        if ϕ[i] < tol
            ϕ[i] = tol
        end
    end
    ϕ ./ sum(ϕ, dims=2)
end

function clr(θ)
    x = exp.(θ)
    S = sum(x, dims=2)
    x ./ S
end

function iclr(ϕ)
    x = log.(ϕ)
    x .- mean(x, dims=2)
end

function logfgrad(θ, θ_mean, θ_var, y, m)
    e = exp.(θ)
    S = sum(e)
    resid = θ - θ_mean
    logf = dot(y, θ) - m * log(S) - 0.5sum((resid.^2) ./ θ_var)
    grad = y - m .* e ./ S - resid ./ θ_var
    logf, grad
end

function get_post(sim, data, param)
    N = data[:N]
    K = data[:K]
    L = data[:L]
    p = data[:p]
    num_blocking_factors = data[:num_blocking_factors]
    param_names = Dict{Symbol, Any}(
        :Y => ["Y[$i, $k]" for i in 1:N, k in 1:K],
        :ϕ => ["ϕ[$i, $k]" for i in 1:N, k in 1:K],
        :θ => ["θ[$i, $k]" for i in 1:N, k in 1:K],
        :θ_mean => ["θ_mean[$i, $k]" for i in 1:N, k in 1:K],
        :θ_var => ["θ_var[$k]" for k in 1:K],
        :μ => ["μ[$k]" for k in 1:K],
        :μ_var => ["μ_var"],
        :β => ["β[$k, $j]" for k in 1:K, j in 1:p],
        :β_full => ["β[$k, $j]" for k in 1:K, j in 1:p],
        :β_var => ["β_var[$j]" for j in 1:p],
        :γ => ["γ[$k, $r]" for k in 1:K, r in 1:num_blocking_factors],
        :γ_var => ["γ_var[$r]" for r in 1:num_blocking_factors],
        :ω => ["ω[$l, $j]" for j in 1:p, l in 1:L],
        :π => ["π[$j]" for j in 1:p],
        :Λ => ["Λ[$l, $k]" for l in 1:L, k in 1:K],
        :λ_var => ["λ_var[$l, $k]" for l in 1:L, k in 1:K],
        :a => ["a[$l]" for l in 1:L],
        :ψ => ["ψ[$l, $k]" for l in 1:L, k in 1:K],
        :ξ => ["ξ[$l, $k]" for l in 1:L, k in 1:K],
        :τ => ["τ[$l]" for l in 1:L],
        :ν => ["ν"],
        :g => ["g[$l, $r]" for l in 1:L, r in 1:num_blocking_factors],
        :g_var => ["g_var[$r]" for r in 1:num_blocking_factors],
        :b => ["b[$l, $j]" for l in 1:L, j in 1:p],
        :b_full => ["b[$l, $j]" for l in 1:L, j in 1:p],
        :b_var => ["b_var[$j]" for j in 1:p],
    )
    _, _, chains = size(sim)
    names = param_names[param]
    post = vcat([sim[:, vec(names), c].value for c in 1:chains]...)
    post = DataFrame(dropdims(post, dims=3))
    names!(post, [Symbol(name) for name in vec(names)])
    return post
end

function parse_config(conf)
    parsed_conf = Dict{Symbol, Any}()
    for key in keys(conf)
        param = key
        for (name, letter) in greek
            if startswith(param, name)
                param = replace(param, name => letter)
            end
        end
        parsed_conf[Symbol(param)] = conf[key]
    end
    return parsed_conf
end

greek = Dict(
    "beta" => "β",
    "epsilon" => "ϵ",
    "gamma" => "γ",
    "lambda" => "λ",
    "Lambda" => "Λ",
    "mu" => "μ",
    "nu" => "ν",
    "omega" => "ω",
    "phi" => "ϕ",
    "pi" => "π",
    "psi" => "ψ",  # must come after epsilon, or will change epsilon -> eψlon :(
    "tau" => "τ",
    "theta" => "θ",
    "xi" => "ξ"
)

latin = Dict([(v, k) for (k, v) in greek])

function load_config(filename::AbstractString)
    @assert isfile(filename)
    conf = YAML.load(open(filename, "r"))
    return parse_config(conf)
end

load_config(filenames::Vector{T}) where T <: AbstractString = merge!([load_config(f) for f in filenames]...)

max_stars_repo_name: reich-group/MIMIX

max_issues_repo_name: reich-group/MIMIX

max_forks_repo_name: reich-group/MIMIX


text:
using BandedMatrices, LazyArrays, Test

@testset "gbmm!" begin
    # test gbmm! subpieces step by step and column by column
    for n in (1,5,50), ν in (1,5,50), m in (1,5,50),
                    Al in (0,1,2,30), Au in (0,1,2,30),
                    Bl in (0,1,2,30), Bu in (0,1,2,30)
        A=brand(n,ν,Al,Au)
        B=brand(ν,m,Bl,Bu)
        α,β,T=0.123,0.456,Float64
        C=brand(Float64,n,m,A.l+B.l,A.u+B.u)
        a=pointer(A.data)
        b=pointer(B.data)
        c=pointer(C.data)
        sta=max(1,stride(A.data,2))
        stb=max(1,stride(B.data,2))
        stc=max(1,stride(C.data,2))

        sz=sizeof(T)

        mr=1:min(m,1+B.u)
        exC=(β*Matrix(C)+α*Matrix(A)*Matrix(B))
        for j=mr
            BandedMatrices.A11_Btop_Ctop_gbmv!(α,β,
                                           n,ν,m,j,
                                           sz,
                                           a,A.l,A.u,sta,
                                           b,B.l,B.u,stb,
                                           c,C.l,C.u,stc)
       end
        @test C[:,mr] ≈ exC[:,mr]

        mr=1+B.u:min(1+C.u,ν+B.u,m)
        exC=(β*Matrix(C)+α*Matrix(A)*Matrix(B))
        for j=mr
            BandedMatrices.Atop_Bmid_Ctop_gbmv!(α,β,
                                           n,ν,m,j,
                                           sz,
                                           a,A.l,A.u,sta,
                                           b,B.l,B.u,stb,
                                           c,C.l,C.u,stc)
       end
       if !isempty(mr)
           @test C[:,mr] ≈ exC[:,mr]
       end

       mr=1+C.u:min(m,ν+B.u,n+C.u)
       exC=(β*Matrix(C)+α*Matrix(A)*Matrix(B))
       for j=mr
           BandedMatrices.Amid_Bmid_Cmid_gbmv!(α,β,
                                          n,ν,m,j,
                                          sz,
                                          a,A.l,A.u,sta,
                                          b,B.l,B.u,stb,
                                          c,C.l,C.u,stc)
      end
      if !isempty(mr)
          @test C[:,mr] ≈ exC[:,mr]
      end

      mr=ν+B.u+1:min(m,n+C.u)
      exC=(β*Matrix(C)+α*Matrix(A)*Matrix(B))
      for j=mr
          BandedMatrices.Anon_Bnon_C_gbmv!(α,β,
                                         n,ν,m,j,
                                         sz,
                                         a,A.l,A.u,sta,
                                         b,B.l,B.u,stb,
                                         c,C.l,C.u,stc)
     end
     if !isempty(mr)
         @test C[:,mr] ≈ exC[:,mr]
     end
    end

    # test gbmm!

    for n in (1,5,50), ν in (1,5,50), m in (1,5,50), Al in (0,1,2,30), Au in (0,1,2,30), Bl in (0,1,2,30), Bu in (0,1,2,30)
        A=brand(n,ν,Al,Au)
        B=brand(ν,m,Bl,Bu)
        α,β,T=0.123,0.456,Float64
        C=brand(Float64,n,m,A.l+B.l,A.u+B.u)

        exC=α*Matrix(A)*Matrix(B)+β*Matrix(C)
        BandedMatrices.gbmm!('N','N', α,A,B,β,C)

        @test Matrix(exC) ≈ Matrix(C)
    end
end

@testset "Negative bands fills with zero" begin
    A = brand(10,10,2,2)
    B = brand(10,10,-2,2)
    C = BandedMatrix(Fill(NaN,10,10),(0,4))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(10,10,-2,2)
    B = brand(10,10,-2,2)
    C = BandedMatrix(Fill(NaN,10,10),(-4,4))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(10,10,-2,2)
    B = brand(10,10,2,2)
    C = BandedMatrix(Fill(NaN,10,10),(0,4))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(10,10,2,2)
    B = brand(10,10,2,-2)
    C = BandedMatrix(Fill(NaN,10,10),(4,0))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(10,10,2,-2)
    B = brand(10,10,2,-2)
    C = BandedMatrix(Fill(NaN,10,10),(4,-4))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(10,10,2,-2)
    B = brand(10,10,2,2)
    C = BandedMatrix(Fill(NaN,10,10),(4,0))
    C .= Mul(A,B)
    @test C == Matrix(A)*Matrix(B)

    A = brand(30,1,0,0)
    B = brand(1,30,17,17)
    C = BandedMatrix(Fill(NaN, 30,30), (17,17))
    C .= Mul(A,B)
    @test C == A*B
end

@testset "Not enough bands" begin
    A = BandedMatrix(Zeros(10,10), (1,1))
    A[band(0)] .= randn(10)
    B = BandedMatrix(randn(10,10), (1,1))
    C = BandedMatrix(Zeros(10,10), (1,1))

    C .= Mul(A,B)

    @test all(C .=== A*B)

    A[band(1)] .= randn(9)
    @test_throws BandError C .= Mul(A,B)

    A = BandedMatrix(randn(2,1), (2,0))
    A[1,1] = 0
    B = BandedMatrix(randn(1,1), (1,1))
    C = BandedMatrix(randn(2,1), (2,0))
    D = A * B + C
    @test BandedMatrices.gbmm!('N','N', 1.0 , A, B, 1.0, C) ≈ D
end

@testset "BandedMatrix{Int} * Vector{Vector{Int}}" begin
    A, x =  [1 2; 3 4] , [[1,2],[3,4]]
    @test BandedMatrix(A)*x == A*x
end


text:
function schooling_test(s)
    @testset "Schooling" begin
        ns = n_school(s);
        @test ns == 3

        @test ed_idx(s, :HSG) == 1
        @test ed_idx(s, :CG) == ns
        @test ed_idx(s, 2) == 2

        @test ed_label(s, ns) == "CG"
        @test endswith(ed_label(s, :CG), "CG")
        edLabelV = ed_labels(s)
        @test isequal(edLabelV[1], ed_label(s, 1))
        @test length(edLabelV) == ns

        edSuffixV = ed_suffixes(s);
        @test endswith(edSuffixV[1], "HSG")
	end
end

@testset "Schooling" begin
    for s ∈ (Schooling3{UInt8}(), )
        schooling_test(s);
    end
end

# --------------


text:
#=
countDfs:
- Julia version:
- Author: Dan
- Date: 2021-11-28
=#

push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools

pdb_filepath = "1AKE.pdb"

comp = System(pdb_filepath)
dfs = DataFrameSystem(comp)

function countDfs()
    length(filter(x -> x == "ALA", dfs.residues.name_))
end

countDfs()

elapsed = @elapsed countDfs()
@show @benchmark countDfs()
println(elapsed)


text:
struct ExtendedVector{T,A,B} <: AbstractVector{T}
  void_to_val::A
  cell_to_val::B
  oldcell_to_cell_or_void::Vector{Int}
  void_to_oldcell::Vector{Int}
  cell_to_oldcell::Vector{Int}
  function ExtendedVector(
    void_to_val::AbstractArray,
    cell_to_val::AbstractArray,
    oldcell_to_cell_or_void::Vector{Int},
    void_to_oldcell::Vector{Int},
    cell_to_oldcell::Vector{Int})

    A = typeof(void_to_val)
    B = typeof(cell_to_val)
    ai = testitem(void_to_val)
    bi = testitem(cell_to_val)
    T = eltype([ai,bi])

    new{T,A,B}(
      void_to_val,
      cell_to_val,
      oldcell_to_cell_or_void,
      void_to_oldcell,
      cell_to_oldcell)

  end
end

Base.size(a::ExtendedVector) = (length(a.oldcell_to_cell_or_void),)

Base.IndexStyle(::Type{<:ExtendedVector}) = IndexLinear()

function Base.getindex(a::ExtendedVector,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end

function array_cache(a::ExtendedVector)
  cv = array_cache(a.void_to_val)
  cc = array_cache(a.cell_to_val)
  (cv,cc)
end

@inline function getindex!(cache,a::ExtendedVector,oldcell)
  cv, cc = cache
  i = a.oldcell_to_cell_or_void[oldcell]
  if i > 0
    cell = i
    return getindex!(cc,a.cell_to_val,cell)
  else
    void = -i
    return getindex!(cv,a.void_to_val,void)
  end
end

function evaluate_field_array(a::ExtendedVector,x::AbstractArray)
  x_void = reindex(x,a.void_to_oldcell)
  x_cell = reindex(x,a.cell_to_oldcell)
  r_void = evaluate_field_array(a.void_to_val,x_void)
  r_cell = evaluate_field_array(a.cell_to_val,x_cell)
  ExtendedVector(
      r_void,
      r_cell,
      a.oldcell_to_cell_or_void,
      a.void_to_oldcell,
      a.cell_to_oldcell)
end

function evaluate_dof_array(dofs::AbstractArray,fields::ExtendedVector)
  _evaluate_dof_array_extended_fields(dofs,fields)
end

function evaluate_dof_array(dofs::AbstractArray{<:Dof},fields::ExtendedVector)
  _evaluate_dof_array_extended_fields(dofs,fields)
end

function _evaluate_dof_array_extended_fields(dofs,fields)
  dofs_void = reindex(dofs,fields.void_to_oldcell)
  dofs_cell = reindex(dofs,fields.cell_to_oldcell)
  r_void = evaluate_dof_array(dofs_void,fields.void_to_val)
  r_cell = evaluate_dof_array(dofs_cell,fields.cell_to_val)
  ExtendedVector(
      r_void,
      r_cell,
      fields.oldcell_to_cell_or_void,
      fields.void_to_oldcell,
      fields.cell_to_oldcell)
end

function reindex(a::ExtendedVector,ptrs::AbstractArray)
  _extended_reindex(a,ptrs)
end

function reindex(a::ExtendedVector,ptrs::IdentityVector)
  a
end

function reindex(a::ExtendedVector,ptrs::SkeletonPair)
  _extended_reindex(a,ptrs::SkeletonPair)
end

function apply(f::Fill,a::ExtendedVector...)
  a_void, a_cell = _split_void_non_void(a...)
  r_void = apply(f.value,a_void...)
  r_cell = apply(f.value,a_cell...)
  ExtendedVector(
      r_void,
      r_cell,
      a[1].oldcell_to_cell_or_void,
      a[1].void_to_oldcell,
      a[1].cell_to_oldcell)
end

function apply(::Type{T},f::Fill,a::ExtendedVector...) where T
  a_void, a_cell = _split_void_non_void(a...)
  r_void = apply(T,f.value,a_void...)
  r_cell = apply(T,f.value,a_cell...)
  ExtendedVector(
      r_void,
      r_cell,
      a[1].oldcell_to_cell_or_void,
      a[1].void_to_oldcell,
      a[1].cell_to_oldcell)
end

function _split_void_non_void(a::ExtendedVector...)
  @notimplementedif ! all(map(i->i.void_to_oldcell===a[1].void_to_oldcell,a))
  @notimplementedif ! all(map(i->i.cell_to_oldcell===a[1].cell_to_oldcell,a))
  a_void = map(i->i.void_to_val,a)
  a_cell = map(i->i.cell_to_val,a)
  a_void, a_cell
end

#function reindex(a::ExtendedVector,trian::Triangulation)
#  ptrs = get_cell_id(trian)
#  _extended_reindex(a,ptrs)
#end

function _extended_reindex(a,ptrs::SkeletonPair)
  left= _extended_reindex(a,ptrs.left)
  right = _extended_reindex(a,ptrs.right)
  SkeletonPair(left,right)
end

function _extended_reindex(a,ptrs::AbstractArray)
  if a.cell_to_oldcell === ptrs || a.cell_to_oldcell == ptrs
    return a.cell_to_val
  elseif a.void_to_oldcell === ptrs || a.void_to_oldcell == ptrs
    return a.void_to_val
  else
    j_to_oldcell = ptrs
    j_to_cell_or_void = a.oldcell_to_cell_or_void[j_to_oldcell]
    all_cell, all_void = _find_all_cell_all_void(j_to_cell_or_void)
    if all_cell
      return reindex(a.cell_to_val,j_to_cell_or_void)
    elseif all_void
      j_to_cell_or_void .*= -1
      return reindex(a.void_to_val,j_to_cell_or_void)
    else
      return Reindexed(a,ptrs)
    end
  end
end

function _find_all_cell_all_void(j_to_cell_or_void)
  all_cell = true
  all_void = true
  for k in j_to_cell_or_void
    all_cell = all_cell && (k>0)
    all_void = all_void && (k<0)
  end
  all_cell, all_void
end

struct VoidBasis{T,D} <: Field end

function field_cache(f::VoidBasis{T},x) where T
  Q = length(x)
  v = zeros(T,(Q,0))
  CachedArray(v)
end

@inline function evaluate_field!(cache,f::VoidBasis{T},x) where T
  Q = length(x)
  setsize!(cache,(Q,0))
  cache.array
end

function field_gradient(f::VoidBasis{T,D}) where {T,D}
  x = zero(Point{D,eltype(T)})
  G = gradient_type(T,x)
  VoidBasis{G,D}()
end

"""
"""
struct ExtendedFESpace <: SingleFieldFESpace
  space::SingleFieldFESpace
  trian::RestrictedTriangulation
end

constraint_style(::Type{ExtendedFESpace}) = Val{false}()

function get_cell_dofs(f::ExtendedFESpace)

  cell_to_val = get_cell_dofs(f.space)
  vi = testitem(cell_to_val)
  T = eltype(vi)
  void_to_val = Fill(T[],length(f.trian.void_to_oldcell))

  ExtendedVector(
    void_to_val,
    cell_to_val,
    f.trian.oldcell_to_cell,
    f.trian.void_to_oldcell,
    f.trian.cell_to_oldcell)

end

function get_cell_axes(f::ExtendedFESpace)
  _extend_cell_axes(f,get_cell_axes(f.space))
end

function get_cell_axes_with_constraints(f::ExtendedFESpace)
  _extend_cell_axes(f,get_cell_axes_with_constraints(f.space))
end

function _extend_cell_axes(f,cell_axes)
  cell_axes_void = Fill(_empty_axes(testitem(cell_axes)...),length(f.trian.void_to_oldcell))
  array = ExtendedVector(
    cell_axes_void,
    cell_axes,
    f.trian.oldcell_to_cell,
    f.trian.void_to_oldcell,
    f.trian.cell_to_oldcell)
  # TODO this is a hack
  # to circumvent a julia bug (?) that is triggered
  # by commit f8e8942c30ffa8d6756e9b68386c19d16914a3dc
  collect(array)
end

#TODO this is a little hacky
function _empty_axes(::Base.OneTo)
  (Base.OneTo(0),)
end

function _empty_axes(::Base.OneTo,::Base.OneTo)
  (Base.OneTo(1),Base.OneTo(0))
end

function _empty_axes(::BlockedUnitRange)
  (blockedrange([0]),)
end

function _empty_axes(::BlockedUnitRange,::BlockedUnitRange)
  (blockedrange([1]),blockedrange([0]))
end

function get_cell_basis(f::ExtendedFESpace)
  cell_basis = get_cell_basis(f.space)
  _extend_cell_basis(f,cell_basis)
end

function _extend_cell_basis(f,cell_basis)
  cell_to_val = get_array(cell_basis)

  xi = testitem(get_cell_coordinates(f.trian))
  vi = testitem(cell_to_val)
  Tv = field_return_type(vi,xi)
  T = eltype(Tv)
  D = num_components(eltype(xi))

  _void_to_val = Fill(VoidBasis{T,D}(),length(f.trian.void_to_oldcell))
  void_to_val = _trialize_if_needed(_void_to_val,MetaSizeStyle(cell_basis))

  array = ExtendedVector(
    void_to_val,
    cell_to_val,
    f.trian.oldcell_to_cell,
    f.trian.void_to_oldcell,
    f.trian.cell_to_oldcell)

  cm = get_cell_map(f.trian.oldtrian)
  ca = _extend_cell_axes(f,get_cell_axes(cell_basis))
  GenericCellField(array,cm,RefStyle(cell_basis),ca,MetaSizeStyle(cell_basis))
end

function _trialize_if_needed(_void_to_val,::Val{(:,)})
  _void_to_val
end

function _trialize_if_needed(_void_to_val,::Val{(1,:)})
  trialize_array_of_bases(_void_to_val)
end

function CellData.CellField(f::ExtendedFESpace,cell_vals)
  _cell_field_for_ext_fe_space(f,cell_vals,f.space)
end

function  _cell_field_for_ext_fe_space(f,cell_vals,ft::Any)
  _default_cell_field(f,cell_vals)
end

# TODO, this is a bit hacky, lincomb should also work for trial bases
function  _cell_field_for_ext_fe_space(f,cell_vals,ft::TrialFESpace)
  cell_basis = _extend_cell_basis(f,get_cell_basis(ft.space))
  cell_field = lincomb(cell_basis,cell_vals)
  cell_field
end

function get_cell_dof_basis(f::ExtendedFESpace)

  cell_to_val = get_cell_dof_basis(f.space)
  ref_trait = RefStyle(cell_to_val)
  cell_to_val = get_array(cell_to_val)

  D = num_dims(f.trian)
  T = Float64 # TODO
  void_to_val = Fill(LagrangianDofBasis(T,Point{D,T}[]),length(f.trian.void_to_oldcell))

  eb = ExtendedVector(
         void_to_val,
         cell_to_val,
         f.trian.oldcell_to_cell,
         f.trian.void_to_oldcell,
         f.trian.cell_to_oldcell)

  cell_dof_basis = GenericCellDofBasis(ref_trait,eb)

end

function scatter_free_and_dirichlet_values(f::ExtendedFESpace,fv,dv)

  cell_to_val = scatter_free_and_dirichlet_values(f.space,fv,dv)
  T = eltype(fv)
  void_to_val = Fill(T[],length(f.trian.void_to_oldcell))

  ExtendedVector(
    void_to_val,
    cell_to_val,
    f.trian.oldcell_to_cell,
    f.trian.void_to_oldcell,
    f.trian.cell_to_oldcell)
end

function gather_free_and_dirichlet_values!(fv,dv,f::ExtendedFESpace,cv)
  _cv = reindex(cv,f.trian.cell_to_oldcell)
  gather_free_and_dirichlet_values!(fv,dv,f.space,_cv)
end

# Delegated functions

function num_free_dofs(f::ExtendedFESpace)
  num_free_dofs(f.space)
end

function zero_free_values(f::ExtendedFESpace)
  zero_free_values(f.space)
end

function num_dirichlet_dofs(f::ExtendedFESpace)
  num_dirichlet_dofs(f.space)
end

function zero_dirichlet_values(f::ExtendedFESpace)
  zero_dirichlet_values(f.space)
end

function num_dirichlet_tags(f::ExtendedFESpace)
  num_dirichlet_tags(f.space)
end

function get_dirichlet_dof_tag(f::ExtendedFESpace)
  get_dirichlet_dof_tag(f.space)
end

function TrialFESpace(f::ExtendedFESpace)
  U = TrialFESpace(f.space)
  ExtendedFESpace(U,f.trian)
end


text:
# This list was obtained through AWK with Knitro 11.1 (and then pared down) with the following two commands:
# grep "#define" knitro.h | awk '{ printf("const %s = Cint(%s)\n",$2,$3) }'
# grep "#  define" knitro.h | awk '{ printf("const %s = Cint(%s)\n",$3,$4) }'
# For version up until 10.3, use following commands :
# grep "#define" knitro.h | awk '{ printf("const %s = Cint(%s)\n",$2,$3) }'
# grep "#  define" knitro.h | awk '{ printf("const %s = Cint(%s)\n",$3,$4) }'

const KNTRUE = Cint(1)
const KNFALSE = Cint(0)
const KN_INFINITY = Float64(0x1.fffffffffffffp+1023) # equal to DBL_MAX
const KN_PARAMTYPE_INTEGER = Cint(0)
const KN_PARAMTYPE_FLOAT = Cint(1)
const KN_PARAMTYPE_STRING = Cint(2)
const KN_OBJGOAL_MINIMIZE = Cint(0)
const KN_OBJGOAL_MAXIMIZE = Cint(1)
const KN_OBJTYPE_CONSTANT = Cint(-1)
const KN_OBJTYPE_GENERAL = Cint(0)
const KN_OBJTYPE_LINEAR = Cint(1)
const KN_OBJTYPE_QUADRATIC = Cint(2)
const KN_CONTYPE_CONSTANT = Cint(-1)
const KN_CONTYPE_GENERAL = Cint(0)
const KN_CONTYPE_LINEAR = Cint(1)
const KN_CONTYPE_QUADRATIC = Cint(2)
const KN_CONTYPE_CONIC = Cint(3)
const KN_RSDTYPE_CONSTANT = Cint(-1)
const KN_RSDTYPE_GENERAL = Cint(0)
const KN_RSDTYPE_LINEAR = Cint(1)
const KN_CCTYPE_VARVAR = Cint(0)
const KN_CCTYPE_VARCON = Cint(1)
const KN_CCTYPE_CONCON = Cint(2)
const KN_VARTYPE_CONTINUOUS = Cint(0)
const KN_VARTYPE_INTEGER = Cint(1)
const KN_VARTYPE_BINARY = Cint(2)
const KN_VAR_LINEAR = Cint(1)
const KN_OBJ_CONVEX = Cint(1)
const KN_OBJ_CONCAVE = Cint(2)
const KN_OBJ_CONTINUOUS = Cint(4)
const KN_OBJ_DIFFERENTIABLE = Cint(8)
const KN_OBJ_TWICE_DIFFERENTIABLE = Cint(16)
const KN_OBJ_NOISY = Cint(32)
const KN_OBJ_NONDETERMINISTIC = Cint(64)
const KN_CON_CONVEX = Cint(1)
const KN_CON_CONCAVE = Cint(2)
const KN_CON_CONTINUOUS = Cint(4)
const KN_CON_DIFFERENTIABLE = Cint(8)
const KN_CON_TWICE_DIFFERENTIABLE = Cint(16)
const KN_CON_NOISY = Cint(32)
const KN_CON_NONDETERMINISTIC = Cint(64)
const KN_DENSE = Cint(-1)
const KN_DENSE_ROWMAJOR = Cint(-2)
const KN_DENSE_COLMAJOR = Cint(-3)
const KN_RC_EVALFC = Cint(1)
const KN_RC_EVALGA = Cint(2)
const KN_RC_EVALH = Cint(3)
const KN_RC_EVALHV = Cint(7)
const KN_RC_EVALH_NO_F = Cint(8)
const KN_RC_EVALHV_NO_F = Cint(9)
const KN_RC_EVALR = Cint(10)
const KN_RC_EVALRJ = Cint(11)
const KN_RC_EVALFCGA = Cint(12)
const KN_RC_OPTIMAL_OR_SATISFACTORY = Cint(0)
const KN_RC_OPTIMAL = Cint(0)
const KN_RC_NEAR_OPT = Cint(-100)
const KN_RC_FEAS_XTOL = Cint(-101)
const KN_RC_FEAS_NO_IMPROVE = Cint(-102)
const KN_RC_FEAS_FTOL = Cint(-103)
const KN_RC_INFEASIBLE = Cint(-200)
const KN_RC_INFEAS_XTOL = Cint(-201)
const KN_RC_INFEAS_NO_IMPROVE = Cint(-202)
const KN_RC_INFEAS_MULTISTART = Cint(-203)
const KN_RC_INFEAS_CON_BOUNDS = Cint(-204)
const KN_RC_INFEAS_VAR_BOUNDS = Cint(-205)
const KN_RC_UNBOUNDED = Cint(-300)
if KNITRO_VERSION > v"12.0"
    const KN_RC_UNBOUNDED_OR_INFEAS = Cint(-301)
end
const KN_RC_ITER_LIMIT_FEAS = Cint(-400)
const KN_RC_TIME_LIMIT_FEAS = Cint(-401)
const KN_RC_FEVAL_LIMIT_FEAS = Cint(-402)
const KN_RC_MIP_EXH_FEAS = Cint(-403)
const KN_RC_MIP_TERM_FEAS = Cint(-404)
const KN_RC_MIP_SOLVE_LIMIT_FEAS = Cint(-405)
const KN_RC_MIP_NODE_LIMIT_FEAS = Cint(-406)
const KN_RC_ITER_LIMIT_INFEAS = Cint(-410)
const KN_RC_TIME_LIMIT_INFEAS = Cint(-411)
const KN_RC_FEVAL_LIMIT_INFEAS = Cint(-412)
const KN_RC_MIP_EXH_INFEAS = Cint(-413)
const KN_RC_MIP_SOLVE_LIMIT_INFEAS = Cint(-415)
const KN_RC_MIP_NODE_LIMIT_INFEAS = Cint(-416)
const KN_RC_CALLBACK_ERR = Cint(-500)
const KN_RC_LP_SOLVER_ERR = Cint(-501)
const KN_RC_EVAL_ERR = Cint(-502)
const KN_RC_OUT_OF_MEMORY = Cint(-503)
const KN_RC_USER_TERMINATION = Cint(-504)
const KN_RC_OPEN_FILE_ERR = Cint(-505)
const KN_RC_BAD_N_OR_F = Cint(-506)
const KN_RC_BAD_CONSTRAINT = Cint(-507)
const KN_RC_BAD_JACOBIAN = Cint(-508)
const KN_RC_BAD_HESSIAN = Cint(-509)
const KN_RC_BAD_CON_INDEX = Cint(-510)
const KN_RC_BAD_JAC_INDEX = Cint(-511)
const KN_RC_BAD_HESS_INDEX = Cint(-512)
const KN_RC_BAD_CON_BOUNDS = Cint(-513)
const KN_RC_BAD_VAR_BOUNDS = Cint(-514)
const KN_RC_ILLEGAL_CALL = Cint(-515)
const KN_RC_BAD_KCPTR = Cint(-516)
const KN_RC_NULL_POINTER = Cint(-517)
const KN_RC_BAD_INIT_VALUE = Cint(-518)
const KN_RC_LICENSE_ERROR = Cint(-520)
const KN_RC_BAD_PARAMINPUT = Cint(-521)
const KN_RC_LINEAR_SOLVER_ERR = Cint(-522)
const KN_RC_DERIV_CHECK_FAILED = Cint(-523)
const KN_RC_DERIV_CHECK_TERMINATE = Cint(-524)
const KN_RC_OVERFLOW_ERR = Cint(-525)
const KN_RC_BAD_SIZE = Cint(-526)
const KN_RC_BAD_VARIABLE = Cint(-527)
const KN_RC_BAD_VAR_INDEX = Cint(-528)
const KN_RC_BAD_OBJECTIVE = Cint(-529)
const KN_RC_BAD_OBJ_INDEX = Cint(-530)
const KN_RC_BAD_RESIDUAL = Cint(-531)
const KN_RC_BAD_RSD_INDEX = Cint(-532)
const KN_RC_INTERNAL_ERROR = Cint(-600)
const KN_PARAM_NEWPOINT = Cint(1001)
const KN_PARAM_HONORBNDS = Cint(1002)
const KN_PARAM_ALGORITHM = Cint(1003)
const KN_PARAM_ALG = Cint(1003)
const KN_PARAM_BAR_MURULE = Cint(1004)
const KN_PARAM_BAR_FEASIBLE = Cint(1006)
const KN_PARAM_GRADOPT = Cint(1007)
const KN_PARAM_HESSOPT = Cint(1008)
const KN_PARAM_BAR_INITPT = Cint(1009)
const KN_PARAM_ACT_LPSOLVER = Cint(1012)
const KN_PARAM_CG_MAXIT = Cint(1013)
const KN_PARAM_MAXIT = Cint(1014)
const KN_PARAM_OUTLEV = Cint(1015)
const KN_PARAM_OUTMODE = Cint(1016)
const KN_PARAM_SCALE = Cint(1017)
const KN_PARAM_SOC = Cint(1019)
const KN_PARAM_DELTA = Cint(1020)
const KN_PARAM_BAR_FEASMODETOL = Cint(1021)
const KN_PARAM_FEASTOL = Cint(1022)
const KN_PARAM_FEASTOLABS = Cint(1023)
const KN_PARAM_MAXTIMECPU = Cint(1024)
const KN_PARAM_BAR_INITMU = Cint(1025)
const KN_PARAM_OBJRANGE = Cint(1026)
const KN_PARAM_OPTTOL = Cint(1027)
const KN_PARAM_OPTTOLABS = Cint(1028)
const KN_PARAM_LINSOLVER_PIVOTTOL = Cint(1029)
const KN_PARAM_XTOL = Cint(1030)
const KN_PARAM_DEBUG = Cint(1031)
const KN_PARAM_MULTISTART = Cint(1033)
const KN_PARAM_MSENABLE = Cint(1033)
const KN_PARAM_MSMAXSOLVES = Cint(1034)
const KN_PARAM_MSMAXBNDRANGE = Cint(1035)
const KN_PARAM_MSMAXTIMECPU = Cint(1036)
const KN_PARAM_MSMAXTIMEREAL = Cint(1037)
const KN_PARAM_LMSIZE = Cint(1038)
const KN_PARAM_BAR_MAXCROSSIT = Cint(1039)
const KN_PARAM_MAXTIMEREAL = Cint(1040)
const KN_PARAM_CG_PRECOND = Cint(1041)
const KN_PARAM_BLASOPTION = Cint(1042)
const KN_PARAM_BAR_MAXREFACTOR = Cint(1043)
const KN_PARAM_LINESEARCH_MAXTRIALS = Cint(1044)
const KN_PARAM_BLASOPTIONLIB = Cint(1045)
const KN_PARAM_OUTAPPEND = Cint(1046)
const KN_PARAM_OUTDIR = Cint(1047)
const KN_PARAM_CPLEXLIB = Cint(1048)
const KN_PARAM_BAR_PENRULE = Cint(1049)
const KN_PARAM_BAR_PENCONS = Cint(1050)
const KN_PARAM_MSNUMTOSAVE = Cint(1051)
const KN_PARAM_MSSAVETOL = Cint(1052)
const KN_PARAM_PRESOLVEDEBUG = Cint(1053)
const KN_PARAM_MSTERMINATE = Cint(1054)
const KN_PARAM_MSSTARTPTRANGE = Cint(1055)
const KN_PARAM_INFEASTOL = Cint(1056)
const KN_PARAM_LINSOLVER = Cint(1057)
const KN_PARAM_BAR_DIRECTINTERVAL = Cint(1058)
const KN_PARAM_PRESOLVE = Cint(1059)
const KN_PARAM_PRESOLVE_TOL = Cint(1060)
const KN_PARAM_BAR_SWITCHRULE = Cint(1061)
const KN_PARAM_HESSIAN_NO_F = Cint(1062)
const KN_PARAM_MA_TERMINATE = Cint(1063)
const KN_PARAM_MA_MAXTIMECPU = Cint(1064)
const KN_PARAM_MA_MAXTIMEREAL = Cint(1065)
const KN_PARAM_MSSEED = Cint(1066)
const KN_PARAM_MA_OUTSUB = Cint(1067)
const KN_PARAM_MS_OUTSUB = Cint(1068)
const KN_PARAM_XPRESSLIB = Cint(1069)
const KN_PARAM_TUNER = Cint(1070)
const KN_PARAM_TUNER_OPTIONSFILE = Cint(1071)
const KN_PARAM_TUNER_MAXTIMECPU = Cint(1072)
const KN_PARAM_TUNER_MAXTIMEREAL = Cint(1073)
const KN_PARAM_TUNER_OUTSUB = Cint(1074)
const KN_PARAM_TUNER_TERMINATE = Cint(1075)
const KN_PARAM_LINSOLVER_OOC = Cint(1076)
const KN_PARAM_BAR_RELAXCONS = Cint(1077)
const KN_PARAM_MSDETERMINISTIC = Cint(1078)
const KN_PARAM_BAR_REFINEMENT = Cint(1079)
const KN_PARAM_DERIVCHECK = Cint(1080)
const KN_PARAM_DERIVCHECK_TYPE = Cint(1081)
const KN_PARAM_DERIVCHECK_TOL = Cint(1082)
const KN_PARAM_LINSOLVER_INEXACT = Cint(1083)
const KN_PARAM_LINSOLVER_INEXACTTOL = Cint(1084)
const KN_PARAM_MAXFEVALS = Cint(1085)
const KN_PARAM_FSTOPVAL = Cint(1086)
const KN_PARAM_DATACHECK = Cint(1087)
const KN_PARAM_DERIVCHECK_TERMINATE = Cint(1088)
const KN_PARAM_BAR_WATCHDOG = Cint(1089)
const KN_PARAM_FTOL = Cint(1090)
const KN_PARAM_FTOL_ITERS = Cint(1091)
const KN_PARAM_ACT_QPALG = Cint(1092)
const KN_PARAM_BAR_INITPI_MPEC = Cint(1093)
const KN_PARAM_XTOL_ITERS = Cint(1094)
const KN_PARAM_LINESEARCH = Cint(1095)
const KN_PARAM_OUT_CSVINFO = Cint(1096)
const KN_PARAM_INITPENALTY = Cint(1097)
const KN_PARAM_ACT_LPFEASTOL = Cint(1098)
const KN_PARAM_CG_STOPTOL = Cint(1099)
const KN_PARAM_RESTARTS = Cint(1100)
const KN_PARAM_RESTARTS_MAXIT = Cint(1101)
const KN_PARAM_BAR_SLACKBOUNDPUSH = Cint(1102)
const KN_PARAM_CG_PMEM = Cint(1103)
const KN_PARAM_BAR_SWITCHOBJ = Cint(1104)
const KN_PARAM_OUTNAME = Cint(1105)
const KN_PARAM_OUT_CSVNAME = Cint(1106)
const KN_PARAM_ACT_PARAMETRIC = Cint(1107)
const KN_PARAM_ACT_LPDUMPMPS = Cint(1108)
const KN_PARAM_ACT_LPALG = Cint(1109)
const KN_PARAM_ACT_LPPRESOLVE = Cint(1110)
const KN_PARAM_ACT_LPPENALTY = Cint(1111)
const KN_PARAM_BNDRANGE = Cint(1112)
const KN_PARAM_BAR_CONIC_ENABLE = Cint(1113)
const KN_PARAM_CONVEX = Cint(1114)
const KN_PARAM_OUT_HINTS = Cint(1115)
const KN_PARAM_EVAL_FCGA = Cint(1116)
const KN_PARAM_MIP_METHOD = Cint(2001)
const KN_PARAM_MIP_BRANCHRULE = Cint(2002)
const KN_PARAM_MIP_SELECTRULE = Cint(2003)
const KN_PARAM_MIP_INTGAPABS = Cint(2004)
const KN_PARAM_MIP_INTGAPREL = Cint(2005)
const KN_PARAM_MIP_MAXTIMECPU = Cint(2006)
const KN_PARAM_MIP_MAXTIMEREAL = Cint(2007)
const KN_PARAM_MIP_MAXSOLVES = Cint(2008)
const KN_PARAM_MIP_INTEGERTOL = Cint(2009)
const KN_PARAM_MIP_OUTLEVEL = Cint(2010)
const KN_PARAM_MIP_OUTINTERVAL = Cint(2011)
const KN_PARAM_MIP_OUTSUB = Cint(2012)
const KN_PARAM_MIP_DEBUG = Cint(2013)
const KN_PARAM_MIP_IMPLICATNS = Cint(2014)
const KN_PARAM_MIP_GUB_BRANCH = Cint(2015)
const KN_PARAM_MIP_KNAPSACK = Cint(2016)
const KN_PARAM_MIP_ROUNDING = Cint(2017)
const KN_PARAM_MIP_ROOTALG = Cint(2018)
const KN_PARAM_MIP_LPALG = Cint(2019)
const KN_PARAM_MIP_TERMINATE = Cint(2020)
const KN_PARAM_MIP_MAXNODES = Cint(2021)
const KN_PARAM_MIP_HEURISTIC = Cint(2022)
const KN_PARAM_MIP_HEUR_MAXIT = Cint(2023)
const KN_PARAM_MIP_HEUR_MAXTIMECPU = Cint(2024)
const KN_PARAM_MIP_HEUR_MAXTIMEREAL = Cint(2025)
const KN_PARAM_MIP_PSEUDOINIT = Cint(2026)
const KN_PARAM_MIP_STRONG_MAXIT = Cint(2027)
const KN_PARAM_MIP_STRONG_CANDLIM = Cint(2028)
const KN_PARAM_MIP_STRONG_LEVEL = Cint(2029)
const KN_PARAM_MIP_INTVAR_STRATEGY = Cint(2030)
const KN_PARAM_MIP_RELAXABLE = Cint(2031)
const KN_PARAM_MIP_NODEALG = Cint(2032)
const KN_PARAM_MIP_HEUR_TERMINATE = Cint(2033)
const KN_PARAM_MIP_SELECTDIR = Cint(2034)
const KN_PARAM_PAR_NUMTHREADS = Cint(3001)
const KN_PARAM_PAR_CONCURRENT_EVALS = Cint(3002)
const KN_PARAM_PAR_BLASNUMTHREADS = Cint(3003)
const KN_PARAM_PAR_LSNUMTHREADS = Cint(3004)
const KN_PARAM_PAR_MSNUMTHREADS = Cint(3005)
const KN_NEWPOINT_NONE = Cint(0)
const KN_NEWPOINT_SAVEONE = Cint(1)
const KN_NEWPOINT_SAVEALL = Cint(2)
const KN_HONORBNDS_NO = Cint(0)
const KN_HONORBNDS_ALWAYS = Cint(1)
const KN_HONORBNDS_INITPT = Cint(2)
const KN_ALG_AUTOMATIC = Cint(0)
const KN_ALG_AUTO = Cint(0)
const KN_ALG_BAR_DIRECT = Cint(1)
const KN_ALG_BAR_CG = Cint(2)
const KN_ALG_ACT_CG = Cint(3)
const KN_ALG_ACT_SQP = Cint(4)
const KN_ALG_MULTI = Cint(5)
const KN_BAR_MURULE_AUTOMATIC = Cint(0)
const KN_BAR_MURULE_AUTO = Cint(0)
const KN_BAR_MURULE_MONOTONE = Cint(1)
const KN_BAR_MURULE_ADAPTIVE = Cint(2)
const KN_BAR_MURULE_PROBING = Cint(3)
const KN_BAR_MURULE_DAMPMPC = Cint(4)
const KN_BAR_MURULE_FULLMPC = Cint(5)
const KN_BAR_MURULE_QUALITY = Cint(6)
const KN_BAR_FEASIBLE_NO = Cint(0)
const KN_BAR_FEASIBLE_STAY = Cint(1)
const KN_BAR_FEASIBLE_GET = Cint(2)
const KN_BAR_FEASIBLE_GET_STAY = Cint(3)
const KN_GRADOPT_EXACT = Cint(1)
const KN_GRADOPT_FORWARD = Cint(2)
const KN_GRADOPT_CENTRAL = Cint(3)
const KN_GRADOPT_USER_FORWARD = Cint(4)
const KN_GRADOPT_USER_CENTRAL = Cint(5)
const KN_HESSOPT_AUTO = Cint(0)
const KN_HESSOPT_EXACT = Cint(1)
const KN_HESSOPT_BFGS = Cint(2)
const KN_HESSOPT_SR1 = Cint(3)
const KN_HESSOPT_PRODUCT_FINDIFF = Cint(4)
const KN_HESSOPT_PRODUCT = Cint(5)
const KN_HESSOPT_LBFGS = Cint(6)
const KN_HESSOPT_GAUSS_NEWTON = Cint(7)
const KN_BAR_INITPT_AUTO = Cint(0)
const KN_BAR_INITPT_CONVEX = Cint(1)
const KN_BAR_INITPT_NEARBND = Cint(2)
const KN_BAR_INITPT_CENTRAL = Cint(3)
const KN_ACT_LPSOLVER_INTERNAL = Cint(1)
const KN_ACT_LPSOLVER_CPLEX = Cint(2)
const KN_ACT_LPSOLVER_XPRESS = Cint(3)
const KN_OUTLEV_NONE = Cint(0)
const KN_OUTLEV_SUMMARY = Cint(1)
const KN_OUTLEV_ITER_10 = Cint(2)
const KN_OUTLEV_ITER = Cint(3)
const KN_OUTLEV_ITER_VERBOSE = Cint(4)
const KN_OUTLEV_ITER_X = Cint(5)
const KN_OUTLEV_ALL = Cint(6)
const KN_OUTMODE_SCREEN = Cint(0)
const KN_OUTMODE_FILE = Cint(1)
const KN_OUTMODE_BOTH = Cint(2)
const KN_SCALE_NEVER = Cint(0)
const KN_SCALE_NO = Cint(0)
const KN_SCALE_USER_INTERNAL = Cint(1)
const KN_SCALE_USER_NONE = Cint(2)
const KN_SCALE_INTERNAL = Cint(3)
const KN_SOC_NO = Cint(0)
const KN_SOC_MAYBE = Cint(1)
const KN_SOC_YES = Cint(2)
const KN_DEBUG_NONE = Cint(0)
const KN_DEBUG_PROBLEM = Cint(1)
const KN_DEBUG_EXECUTION = Cint(2)
const KN_MULTISTART_NO = Cint(0)
const KN_MULTISTART_YES = Cint(1)
const KN_CG_PRECOND_NONE = Cint(0)
const KN_CG_PRECOND_CHOL = Cint(1)
const KN_BLASOPTION_KNITRO = Cint(0)
const KN_BLASOPTION_INTEL = Cint(1)
const KN_BLASOPTION_DYNAMIC = Cint(2)
const KN_OUTAPPEND_NO = Cint(0)
const KN_OUTAPPEND_YES = Cint(1)
const KN_BAR_PENRULE_AUTO = Cint(0)
const KN_BAR_PENRULE_SINGLE = Cint(1)
const KN_BAR_PENRULE_FLEX = Cint(2)
const KN_BAR_PENCONS_AUTO = Cint(0)
const KN_BAR_PENCONS_NONE = Cint(1)
const KN_BAR_PENCONS_ALL = Cint(2)
const KN_BAR_PENCONS_EQUALITIES = Cint(3)
const KN_BAR_PENCONS_INFEAS = Cint(4)
const KN_PRESOLVEDBG_NONE = Cint(0)
const KN_PRESOLVEDBG_BASIC = Cint(1)
const KN_PRESOLVEDBG_VERBOSE = Cint(2)
const KN_MSTERMINATE_MAXSOLVES = Cint(0)
const KN_MSTERMINATE_OPTIMAL = Cint(1)
const KN_MSTERMINATE_FEASIBLE = Cint(2)
const KN_MSTERMINATE_ANY = Cint(3)
const KN_LINSOLVER_AUTO = Cint(0)
const KN_LINSOLVER_INTERNAL = Cint(1)
const KN_LINSOLVER_HYBRID = Cint(2)
const KN_LINSOLVER_DENSEQR = Cint(3)
const KN_LINSOLVER_MA27 = Cint(4)
const KN_LINSOLVER_MA57 = Cint(5)
const KN_LINSOLVER_MKLPARDISO = Cint(6)
const KN_LINSOLVER_MA97 = Cint(7)
const KN_LINSOLVER_MA86 = Cint(8)
const KN_PRESOLVE_NONE = Cint(0)
const KN_PRESOLVE_BASIC = Cint(1)
const KN_PRESOLVE_ADVANCED = Cint(2)
if KNITRO_VERSION >= v"12.0"
    const KN_PRESOLVE_NO = Cint(0)
    const KN_PRESOLVE_YES = Cint(1)
    const KN_PRESOLVE_ADVANCED = Cint(2)
end
const KN_BAR_SWITCHRULE_AUTO = Cint(0)
const KN_BAR_SWITCHRULE_NEVER = Cint(1)
const KN_BAR_SWITCHRULE_MODERATE = Cint(2)
const KN_BAR_SWITCHRULE_AGGRESSIVE = Cint(3)
const KN_HESSIAN_NO_F_FORBID = Cint(0)
const KN_HESSIAN_NO_F_ALLOW = Cint(1)
const KN_MA_TERMINATE_ALL = Cint(0)
const KN_MA_TERMINATE_OPTIMAL = Cint(1)
const KN_MA_TERMINATE_FEASIBLE = Cint(2)
const KN_MA_TERMINATE_ANY = Cint(3)
const KN_MA_OUTSUB_NONE = Cint(0)
const KN_MA_OUTSUB_YES = Cint(1)
const KN_MS_OUTSUB_NONE = Cint(0)
const KN_MS_OUTSUB_YES = Cint(1)
const KN_TUNER_OFF = Cint(0)
const KN_TUNER_ON = Cint(1)
const KN_TUNER_OUTSUB_NONE = Cint(0)
const KN_TUNER_OUTSUB_SUMMARY = Cint(1)
const KN_TUNER_OUTSUB_ALL = Cint(2)
const KN_TUNER_TERMINATE_ALL = Cint(0)
const KN_TUNER_TERMINATE_OPTIMAL = Cint(1)
const KN_TUNER_TERMINATE_FEASIBLE = Cint(2)
const KN_TUNER_TERMINATE_ANY = Cint(3)
const KN_LINSOLVER_OOC_NO = Cint(0)
const KN_LINSOLVER_OOC_MAYBE = Cint(1)
const KN_LINSOLVER_OOC_YES = Cint(2)
const KN_BAR_RELAXCONS_NONE = Cint(0)
const KN_BAR_RELAXCONS_EQS = Cint(1)
const KN_BAR_RELAXCONS_INEQS = Cint(2)
const KN_BAR_RELAXCONS_ALL = Cint(3)
const KN_MSDETERMINISTIC_NO = Cint(0)
const KN_MSDETERMINISTIC_YES = Cint(1)
const KN_BAR_REFINEMENT_NO = Cint(0)
const KN_BAR_REFINEMENT_YES = Cint(1)
const KN_DERIVCHECK_NONE = Cint(0)
const KN_DERIVCHECK_FIRST = Cint(1)
const KN_DERIVCHECK_SECOND = Cint(2)
const KN_DERIVCHECK_ALL = Cint(3)
const KN_DERIVCHECK_FORWARD = Cint(1)
const KN_DERIVCHECK_CENTRAL = Cint(2)
const KN_LINSOLVER_INEXACT_NO = Cint(0)
const KN_LINSOLVER_INEXACT_YES = Cint(1)
const KN_DATACHECK_NO = Cint(0)
const KN_DATACHECK_YES = Cint(1)
const KN_DERIVCHECK_STOPERROR = Cint(1)
const KN_DERIVCHECK_STOPALWAYS = Cint(2)
const KN_BAR_WATCHDOG_NO = Cint(0)
const KN_BAR_WATCHDOG_YES = Cint(1)
const KN_ACT_QPALG_AUTO = Cint(0)
const KN_ACT_QPALG_BAR_DIRECT = Cint(1)
const KN_ACT_QPALG_BAR_CG = Cint(2)
const KN_ACT_QPALG_ACT_CG = Cint(3)
const KN_LINESEARCH_AUTO = Cint(0)
const KN_LINESEARCH_BACKTRACK = Cint(1)
const KN_LINESEARCH_INTERPOLATE = Cint(2)
const KN_OUT_CSVINFO_NO = Cint(0)
const KN_OUT_CSVINFO_YES = Cint(1)
const KN_BAR_SWITCHOBJ_NONE = Cint(0)
const KN_BAR_SWITCHOBJ_SCALARPROX = Cint(1)
const KN_BAR_SWITCHOBJ_DIAGPROX = Cint(2)
const KN_ACT_PARAMETRIC_NO = Cint(0)
const KN_ACT_PARAMETRIC_MAYBE = Cint(1)
const KN_ACT_PARAMETRIC_YES = Cint(2)
const KN_ACT_LPDUMPMPS_NO = Cint(0)
const KN_ACT_LPDUMPMPS_YES = Cint(1)
const KN_ACT_LPALG_DEFAULT = Cint(0)
const KN_ACT_LPALG_PRIMAL = Cint(1)
const KN_ACT_LPALG_DUAL = Cint(2)
const KN_ACT_LPALG_BARRIER = Cint(3)
const KN_ACT_LPPRESOLVE_OFF = Cint(0)
const KN_ACT_LPPRESOLVE_ON = Cint(1)
const KN_ACT_LPPENALTY_ALL = Cint(1)
const KN_ACT_LPPENALTY_NONLINEAR = Cint(2)
const KN_ACT_LPPENALTY_DYNAMIC = Cint(3)
const KN_BAR_CONIC_ENABLE_NONE = Cint(0)
const KN_BAR_CONIC_ENABLE_SOC = Cint(1)
const KN_CONVEX_NO = Cint(0)
const KN_CONVEX_YES = Cint(1)
const KN_OUT_HINTS_NO = Cint(0)
const KN_OUT_HINTS_YES = Cint(1)
const KN_EVAL_FCGA_NO = Cint(0)
const KN_EVAL_FCGA_YES = Cint(1)
const KN_MIP_METHOD_AUTO = Cint(0)
const KN_MIP_METHOD_BB = Cint(1)
const KN_MIP_METHOD_HQG = Cint(2)
const KN_MIP_METHOD_MISQP = Cint(3)
const KN_MIP_BRANCH_AUTO = Cint(0)
const KN_MIP_BRANCH_MOSTFRAC = Cint(1)
const KN_MIP_BRANCH_PSEUDOCOST = Cint(2)
const KN_MIP_BRANCH_STRONG = Cint(3)
const KN_MIP_SEL_AUTO = Cint(0)
const KN_MIP_SEL_DEPTHFIRST = Cint(1)
const KN_MIP_SEL_BESTBOUND = Cint(2)
const KN_MIP_SEL_COMBO_1 = Cint(3)
const KN_MIP_OUTLEVEL_NONE = Cint(0)
const KN_MIP_OUTLEVEL_ITERS = Cint(1)
const KN_MIP_OUTLEVEL_ITERSTIME = Cint(2)
const KN_MIP_OUTLEVEL_ROOT = Cint(3)
const KN_MIP_OUTSUB_NONE = Cint(0)
const KN_MIP_OUTSUB_YES = Cint(1)
const KN_MIP_OUTSUB_YESPROB = Cint(2)
const KN_MIP_DEBUG_NONE = Cint(0)
const KN_MIP_DEBUG_ALL = Cint(1)
const KN_MIP_IMPLICATNS_NO = Cint(0)
const KN_MIP_IMPLICATNS_YES = Cint(1)
const KN_MIP_GUB_BRANCH_NO = Cint(0)
const KN_MIP_GUB_BRANCH_YES = Cint(1)
const KN_MIP_KNAPSACK_NO = Cint(0)
const KN_MIP_KNAPSACK_NONE = Cint(0)
const KN_MIP_KNAPSACK_INEQ = Cint(1)
const KN_MIP_KNAPSACK_INEQ_EQ = Cint(2)
const KN_MIP_ROUND_AUTO = Cint(0)
const KN_MIP_ROUND_NONE = Cint(1)
const KN_MIP_ROUND_HEURISTIC = Cint(2)
const KN_MIP_ROUND_NLP_SOME = Cint(3)
const KN_MIP_ROUND_NLP_ALWAYS = Cint(4)
const KN_MIP_ROOTALG_AUTO = Cint(0)
const KN_MIP_ROOTALG_BAR_DIRECT = Cint(1)
const KN_MIP_ROOTALG_BAR_CG = Cint(2)
const KN_MIP_ROOTALG_ACT_CG = Cint(3)
const KN_MIP_ROOTALG_ACT_SQP = Cint(4)
const KN_MIP_ROOTALG_MULTI = Cint(5)
const KN_MIP_LPALG_AUTO = Cint(0)
const KN_MIP_LPALG_BAR_DIRECT = Cint(1)
const KN_MIP_LPALG_BAR_CG = Cint(2)
const KN_MIP_LPALG_ACT_CG = Cint(3)
const KN_MIP_TERMINATE_OPTIMAL = Cint(0)
const KN_MIP_TERMINATE_FEASIBLE = Cint(1)
const KN_MIP_HEURISTIC_AUTO = Cint(0)
const KN_MIP_HEURISTIC_NONE = Cint(1)
const KN_MIP_HEURISTIC_FEASPUMP = Cint(2)
const KN_MIP_HEURISTIC_MPEC = Cint(3)
const KN_MIP_PSEUDOINIT_AUTO = Cint(0)
const KN_MIP_PSEUDOINIT_AVE = Cint(1)
const KN_MIP_PSEUDOINIT_STRONG = Cint(2)
const KN_MIP_INTVAR_STRATEGY_NONE = Cint(0)
const KN_MIP_INTVAR_STRATEGY_RELAX = Cint(1)
const KN_MIP_INTVAR_STRATEGY_MPEC = Cint(2)
const KN_MIP_RELAXABLE_NONE = Cint(0)
const KN_MIP_RELAXABLE_ALL = Cint(1)
const KN_MIP_NODEALG_AUTO = Cint(0)
const KN_MIP_NODEALG_BAR_DIRECT = Cint(1)
const KN_MIP_NODEALG_BAR_CG = Cint(2)
const KN_MIP_NODEALG_ACT_CG = Cint(3)
const KN_MIP_NODEALG_ACT_SQP = Cint(4)
const KN_MIP_NODEALG_MULTI = Cint(5)
const KN_MIP_HEUR_TERMINATE_FEASIBLE = Cint(1)
const KN_MIP_HEUR_TERMINATE_LIMIT = Cint(2)
const KN_MIP_SELECTDIR_DOWN = Cint(0)
const KN_MIP_SELECTDIR_UP = Cint(1)
const KN_PAR_CONCURRENT_EVALS_NO = Cint(0)
const KN_PAR_CONCURRENT_EVALS_YES = Cint(1)
const KN_PAR_MSNUMTHREADS_AUTO = Cint(0)
if KNITRO_VERSION >= v"12.0"
    # Cuts parameters
    const KN_PARAM_MIP_ZEROHALF = Cint(2036)
    const KN_MIP_ZEROHALF_NONE = Cint(0)
    const KN_MIP_ZEROHALF_ROOT = Cint(1)
    const KN_MIP_ZEROHALF_TREE = Cint(2)
    const KN_MIP_ZEROHALF_ALL = Cint(3)
    const KN_PARAM_MIP_MIR = Cint(2037)
    const KN_MIP_MIR_NONE = Cint(0)
    const KN_MIP_MIR_TREE = Cint(1)
    const KN_PARAM_MIP_CLIQUE = Cint(2038)
    const KN_MIP_CLIQUE_NONE = Cint(0)
    const KN_MIP_CLIQUE_ROOT = Cint(1)
    const KN_MIP_CLIQUE_TREE = Cint(2)
    const KN_MIP_CLIQUE_ALL = Cint(3)
    # SOCP parameters
    const KN_PARAM_PAR_CONICNUMTHREADS = Cint(3006)
    # Presolve parameters
    const KN_PARAM_PRESOLVE_PASSES = Cint(1121)
    const KN_PARAM_PRESOLVE_LEVEL = Cint(1122)
    const KN_PRESOLVE_LEVEL_AUTO = Cint(-1)
    const KN_PRESOLVE_LEVEL_1 = Cint(1)
    const KN_PRESOLVE_LEVEL_2 = Cint(2)
end
if KNITRO_VERSION >= v"12.1"
    const KN_PARAM_FINDIFF_RELSTEPSIZE = Cint(1123)
    const KN_PARAM_INFEASTOL_ITERS = Cint(1124)
    # New options for MIP MIR cuts
    const KN_MIP_MIR_AUTO = Cint(-1)
    const KN_MIP_MIR_NLP = Cint(2)
    # New options for presolve
    const KN_PARAM_PRESOLVEOP_TIGHTEN = Cint(1125)
    const KN_PRESOLVEOP_TIGHTEN_AUTO = Cint(-1)
    const KN_PRESOLVEOP_TIGHTEN_NONE = Cint(0)
    const KN_PRESOLVEOP_TIGHTEN_VARBND = Cint(1)
end
if KNITRO_VERSION >= v"12.2"
    const KN_LINESEARCH_WEAKWOLFE = Cint(3)
    # Tune factorization of KKT system
    const KN_PARAM_BAR_LINSYS = Cint(1126)
    const KN_BAR_LINSYS_AUTO = Cint(-1)
    const KN_BAR_LINSYS_FULL = Cint(0)
    const KN_BAR_LINSYS_COMPACT1 = Cint(1)
    const KN_BAR_LINSYS_COMPACT2 = Cint(2)
    # Add new ways to handle init point during presolve
    const KN_PARAM_PRESOLVE_INITPT = Cint(1127)
    const KN_PRESOLVE_INITPT_AUTO = Cint(-1)
    const KN_PRESOLVE_INITPT_NOSHIFT = Cint(0)
    const KN_PRESOLVE_INITPT_LINSHIFT = Cint(1)
    const KN_PRESOLVE_INITPT_ANYSHIFT = Cint(2)
    # QP penalty in active set algorithm
    const KN_PARAM_ACT_QPPENALTY = Cint(1128)
    const KN_ACT_QPPENALTY_AUTO = Cint(-1)
    const KN_ACT_QPPENALTY_NONE = Cint(0)
    const KN_ACT_QPPENALTY_ALL = Cint(1)
end
if KNITRO_VERSION >= v"12.3"
    const KN_PARAM_BAR_LINSYS_STORAGE = Cint(1129)
    const KN_BAR_LINSYS_STORAGE_AUTO = Cint(-1)
    const KN_BAR_LINSYS_STORAGE_LOWMEM = Cint(1)
    const KN_BAR_LINSYS_STORAGE_NORMAL = Cint(2)
    const KN_PARAM_LINSOLVER_MAXITREF = Cint(1130)
    const KN_PARAM_BFGS_SCALING = Cint(1131)
    const KN_BFGS_SCALING_DYNAMIC = Cint(0)
    const KN_BFGS_SCALING_INVHESS = Cint(1)
    const KN_BFGS_SCALING_HESS = Cint(2)
end


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui
using ReduxImGui.OnOffButtons

function basic_button(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.OnOffButton(store, s->get_state(s).buttons["basic_button"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end
    # this works like a button with on/off state
    if OnOffButtons.is_on(get_state(store).buttons["basic_button"])
        CImGui.SameLine()
        CImGui.Text("Thanks for clicking me!")
    end
end

function radio_button_group(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    if ReduxImGui.RadioButtonGroup(store, s->get_state(s).radio_button_group)
        @info "This triggers $(@__FILE__):$(@__LINE__)."
    end
end

function color_buttons(store::AbstractStore, get_state=Redux.get_state)
    for (i, cb) in enumerate(get_state(store).colorful_buttons)
        label = ReduxImGui.ColorButtons.get_label(cb)
        dispatch!(store, ColorButtons.SetButtonColorTo(label, CImGui.HSV((i-1)/7.0, 0.6, 0.6)))
        dispatch!(store, ColorButtons.SetHoveredColorTo(label, CImGui.HSV((i-1)/7.0, 0.7, 0.7)))
        dispatch!(store, ColorButtons.SetActiveColorTo(label, CImGui.HSV((i-1)/7.0, 0.8, 0.8)))
        if ReduxImGui.ColorButton(store, s->get_state(s).colorful_buttons[i])
            @info "This triggers $(@__FILE__):$(@__LINE__)."
        end
        CImGui.SameLine()
    end
    CImGui.NewLine()
end


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

mutable struct UsageName <: SwaggerModel
    value::Any # spec type: Union{ Nothing, String } # spec name: value
    localizedValue::Any # spec type: Union{ Nothing, String } # spec name: localizedValue

    function UsageName(;value=nothing, localizedValue=nothing)
        o = new()
        validate_property(UsageName, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(UsageName, Symbol("localizedValue"), localizedValue)
        setfield!(o, Symbol("localizedValue"), localizedValue)
        o
    end
end # type UsageName

const _property_map_UsageName = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("localizedValue")=>Symbol("localizedValue"))
const _property_types_UsageName = Dict{Symbol,String}(Symbol("value")=>"String", Symbol("localizedValue")=>"String")
Base.propertynames(::Type{ UsageName }) = collect(keys(_property_map_UsageName))
Swagger.property_type(::Type{ UsageName }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_UsageName[name]))}
Swagger.field_name(::Type{ UsageName }, property_name::Symbol) =  _property_map_UsageName[property_name]

function check_required(o::UsageName)
    true
end

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


text:
"""
    rescore(x::Vector{Real})
Rescore the vector of item response.

# Usage
```julia
mapcols(c -> rescore(c), U)
```
"""
function rescore(u)
    obs = sort(unique(skipmissing(u)))
    rescored = 0:1:length(obs)-1
    d = Dict(zip(obs, rescored))
    return [o === missing ? missing : d[o] for o in u]
end