text:
function logit∂𝑙(y::AbstractFloat, ŷ::AbstractFloat)
    res = ŷ - y
    isnan(res) ? zero(res) : res
end

function logit∂²𝑙(ŷ::T) where {T<:AbstractFloat}
    res = max(ŷ * (one(T) - ŷ), eps(T))
    isnan(res) ? zero(T) : res
end

function logitraw(p::T) where {T<:AbstractFloat}
    -log(one(T) / p - one(T))
end

function logitraw(p::T, posweight::T) where {T<:AbstractFloat}
    -log(one(T) / p - one(T)) + log(posweight)
end

function sigmoid(x::T) where {T<:AbstractFloat}
    one(T) / (one(T) + exp(-x))
end

function logloss(y::AbstractFloat, ŷ::AbstractFloat)
    ϵ = eps(ŷ)
    if ŷ < ϵ
        -y * log(ϵ) - (one(y) - y)  * log(one(ϵ) - ϵ)
    elseif one(ŷ) - ŷ < ϵ
        -y * log(one(ϵ) - ϵ) - (one(y) - y)  * log(ϵ)
    else
        -y * log(ŷ) - (one(y) - y) * log(one(ŷ) - ŷ)
    end
end

function xgblogit(label::AbstractCovariate{S}, factors::Vector{<:AbstractFactor};
                  trainselector::AbstractBoolVariate = BoolVariate("", BitArray{1}(undef, 0)),
                  validselector::AbstractBoolVariate = BoolVariate("", BitArray{1}(undef, 0)),
                  μ::Real = 0.5, subsample::Real = 1.0,
                  η::Real = 0.3, λ::Real = 1.0, γ::Real = 0.0, maxdepth::Integer = 6, nrounds::Integer = 2, ordstumps::Bool = false, optsplit::Bool = false, pruning::Bool = false,
                  minchildweight::Real = 1.0, caching::Bool = true, filecaching::Bool = false, slicelength::Integer = 0, usefloat64::Bool = false, leafwise::Bool = false, maxleaves::Integer = 255,
                  singlethread::Bool = false) where {S<:AbstractFloat}

    T = usefloat64 ? Float64 : Float32
    factors = filter((f -> getname(f) != getname(label)), factors)
    factors = caching ? map(cache, factors) : (filecaching ? map(filecache, factors) : factors)
    label = caching ? cache(label) : (filecaching ? filecache(label) : label)
    slicelength = slicelength <= 0 ? length(label) : slicelength
    validselector = caching ? (validselector |> cache) : (filecaching ? filecache(validselector) : validselector)
    trainselector = length(trainselector) == 0 ? Vector{Bool}() : convert(Vector{Bool}, trainselector)

    λ = T(λ)
    γ = T(γ)
    η = T(η)
    minchildweight = T(minchildweight)
    μ = T(μ)
    subsample = T(subsample)
    f0 = Vector{T}(undef, length(label))
    fill!(f0, T(logitraw(μ)))
    ∂𝑙 = zeros(T, length(label))
    ∂²𝑙 = zeros(T, length(label))
    label = convert(Vector{S}, label)

    fm, trees = fold((f0, Vector{XGTree}()), Seq(1:nrounds)) do x, m
        fm, trees = x

        if length(trainselector) == 0
            @inbounds for i in 1:length(label)
                yhat = sigmoid(fm[i])
                y = label[i]
                ∂𝑙[i] = logit∂𝑙(y, yhat)
                ∂²𝑙[i] = logit∂²𝑙(yhat)
            end
        else
            @inbounds for i in 1:length(label)
                if trainselector[i]
                    yhat = sigmoid(fm[i])
                    y = label[i]
                    ∂𝑙[i] = logit∂𝑙(y, yhat)
                    ∂²𝑙[i] = logit∂²𝑙(yhat)
                end
            end
        end

        tree, predraw = growtree(factors, Covariate(∂𝑙), Covariate(∂²𝑙), maxdepth, λ, γ, leafwise, maxleaves, minchildweight, ordstumps, optsplit, pruning, slicelength, singlethread)

        fm .= muladd.(η, predraw, fm)
        push!(trees, tree)
        (fm, trees)
    end
    fm .= sigmoid.(fm)
    XGModel{T}(trees, λ, γ, η, minchildweight, maxdepth, fm)
end

function cvxgblogit(label::AbstractCovariate, factors::Vector{<:AbstractFactor}, nfolds::Integer;
                    aucmetric::Bool = true, loglossmetric::Bool = true, trainmetric::Bool = false, μ::Real = 0.5,
                    η::Real = 0.3, λ::Real = 1.0, γ::Real = 0.0, maxdepth::Integer = 6, nrounds::Integer = 2, ordstumps::Bool = false, optsplit::Bool = false, pruning::Bool = true,
                    minchildweight::Real = 1.0, caching::Bool = true, slicelength::Integer = 0, usefloat64::Bool = false,
                    singlethread::Bool = false)

    cvfolds = getnfolds(nfolds, false, length(label))
    trainaucfold = Vector{Float64}(undef, nfolds)
    trainloglossfold = Vector{Float64}(undef, nfolds)
    testaucfold = Vector{Float64}(undef, nfolds)
    testloglossfold = Vector{Float64}(undef, nfolds)
    for i in 1:nfolds
        trainselector = cvfolds .!= UInt8(i)
        testselector = cvfolds .== UInt8(i)
        model = xgblogit(label, factors; selector = BoolVariate("", trainselector), η = η, λ = λ, γ = γ, μ = μ, maxdepth = maxdepth,
                         nrounds = nrounds, ordstumps = ordstumps, optsplit = optsplit, pruning = pruning, minchildweight = minchildweight,
                         caching = caching, slicelength = slicelength, usefloat64 = usefloat64, singlethread = singlethread)
        if aucmetric
            testaucfold[i] = getauc(model.pred, label; selector = testselector)
            if trainmetric
                trainaucfold[i] = getauc(model.pred, label; selector = trainselector)
            end
        end
        if loglossmetric
            testloglossfold[i] = getlogloss(model.pred, label; selector = testselector)
            if trainmetric
                trainloglossfold[i] = getlogloss(model.pred, label; selector = trainselector)
            end
        end
    end
    res = Dict{String, Float64}()
    if aucmetric
        if trainmetric
            res["train_auc_mean"] = mean(trainaucfold)
            res["train_auc_std"] = std(trainaucfold)
        end
        res["test_auc_mean"] = mean(testaucfold)
        res["test_auc_std"] = std(testaucfold)
    end
    if loglossmetric
        if trainmetric
            res["train_logloss_mean"] = mean(trainloglossfold)
            res["train_logloss_std"] = std(trainloglossfold)
        end
        res["test_logloss_mean"] = mean(testloglossfold)
        res["test_logloss_std"] = std(testloglossfold)
    end
    res
end

function predict(model::XGModel{T}, dataframe::AbstractDataFrame; μ::Real = 0.5, posweight::Real = 1.0, nrounds::Integer = 0) where {T<:AbstractFloat}
    trees = nrounds == 0 ? model.trees : model.trees[1:nrounds]
    μ = T(μ)
    posweight = T(posweight)
    η = model.η
    f0 = Vector{T}(undef, length(dataframe))
    if posweight == one(T)
        fill!(f0, T(logitraw(μ)))
    else
        fill!(f0, T(logitraw(μ, posweight)))
    end

    factormap = Dict{AbstractFactor, Tuple{AbstractFactor, Dict{Int64, Int64}, Set{Int64}, Int64}}()
    for tree in trees
        for layer in tree.layers
            for node in layer.nodes
                if isa(node, SplitNode) && !(node.factor in keys(factormap))
                    mappedfactor = map(node.factor, dataframe)
                    levelmap = getlevelmap(node.factor, mappedfactor)
                    newind = getnewindices(node.factor, mappedfactor)
                    levelcount = length(getlevels(mappedfactor))
                    factormap[node.factor] = (mappedfactor, levelmap, newind, levelcount)
                end
            end
        end
    end

     mappedfactors = map(cache, collect(map((x -> x[1]), values(factormap))))
     for (i, f) in enumerate(keys(factormap))
         _, levelmap, newind, levelcount = factormap[f]
         factormap[f] = mappedfactors[i], levelmap, newind, levelcount
     end

    mappedtrees = [XGTree{T}(map((layer -> TreeLayer{T}([map(n, dataframe, factormap) for n in layer.nodes])), tree.layers), tree.λ, tree.γ, tree.min∂²𝑙, tree.maxdepth, tree.leafwise, tree.maxleaves, tree.slicelength, tree.singlethread) for tree in trees]

    for tree in mappedtrees
        predraw = predict(tree, dataframe)
        f0 .= muladd.(η, predraw, f0)
    end
    f0 .= sigmoid.(f0)
    f0
end

function getauc(pred::Vector{T}, label::AbstractCovariate{S}; slicelength::Integer = SLICELENGTH) where {T <: AbstractFloat} where {S <: AbstractFloat}
    labelslices = slice(label, 1, length(label), slicelength)

    len = length(pred)
    uniqcount = Dict{T, Int64}()
    labelagg = Dict{T, S}()

    fold(0, labelslices) do offset, slice
        for i in 1:length(slice)
            v = pred[offset + i]
            uniqcount[v] = get(uniqcount, v, 0) + 1
            labelagg[v] = get(labelagg, v, zero(S)) + slice[i]
        end
        offset + length(slice)
    end

    uniqpred = collect(keys(uniqcount))
    sort!(uniqpred; rev = true)
    ucount = map((v -> uniqcount[v]), uniqpred)

    sumlabel = sum(values(labelagg))
    cumlabel = cumsum(map((p -> labelagg[p]), uniqpred))
    cumcount = cumsum(ucount)
    map!((x -> convert(S, x)), cumcount, cumcount)
    tpr = (cumlabel ./ sumlabel)
    fpr = ((cumcount .- cumlabel) ./ (len .- sumlabel))
    len = length(tpr)
    fpr1 = view(fpr, 1:(len-1))
    fpr2 = view(fpr, 2:len)
    tpr1 = view(tpr, 1:(len - 1))
    tpr2 = view(tpr, 2:len)
    area0 = fpr[1] * tpr[1]
    0.5 * (sum((tpr1 .+ tpr2) .* (fpr2 .- fpr1)) + area0)
end

function getauc(pred::Vector{T}, label::AbstractCovariate{S}, trainselector::AbstractBoolVariate, validselector::AbstractBoolVariate; slicelength::Integer = SLICELENGTH) where {T <: AbstractFloat} where {S <: AbstractFloat}
    labelslices = slice(label, 1, length(label), slicelength)
    sellen = length(trainselector)

    trainlentrue = begin
        fold(0, slice(trainselector, 1, sellen, slicelength)) do acc, slice
            res = acc
            for v in slice
                if v
                    res += 1
                end
            end
            res
        end
    end

    validlen = length(validselector)
    validlentrue = validlen == 0 ? 0 : begin
        fold(0, slice(validselector, 1, sellen, slicelength)) do acc, slice
            res = acc
            for v in slice
                if v
                    res += 1
                end
            end
            res
        end
    end

    uniqcountin = Dict{T, Int64}()
    labelaggin = Dict{T, S}()
    uniqcountout = Dict{T, Int64}()
    labelaggout = Dict{T, S}()

    if validlen > 0
        zipslices = zip(labelslices, slice(trainselector, 1, sellen, slicelength), slice(validselector, 1, sellen, slicelength))
        fold(0, zipslices) do offset, slice
            labelslice, trainselslice, validselslice = slice
            for i in 1:length(labelslice)
                if trainselslice[i]
                    v = pred[offset + i]
                    uniqcountin[v] = get(uniqcountin, v, 0) + 1
                    labelaggin[v] = get(labelaggin, v, zero(S)) + labelslice[i]
                elseif validselslice[i]
                    v = pred[offset + i]
                    uniqcountout[v] = get(uniqcountout, v, 0) + 1
                    labelaggout[v] = get(labelaggout, v, zero(S)) + labelslice[i]
                end
            end
            offset + length(labelslice)
        end
    else
        zipslices = zip(labelslices, slice(trainselector, 1, sellen, slicelength))
        fold(0, zipslices) do offset, slice
            labelslice, trainselslice = slice
            for i in 1:length(labelslice)
                if trainselslice[i]
                    v = pred[offset + i]
                    uniqcountin[v] = get(uniqcountin, v, 0) + 1
                    labelaggin[v] = get(labelaggin, v, zero(S)) + labelslice[i]
                end
            end
            offset + length(labelslice)
        end
    end

    uniqpredin = collect(keys(uniqcountin))
    sort!(uniqpredin; rev = true)
    ucountin = map((v -> uniqcountin[v]), uniqpredin)

    if validlen > 0
        uniqpredout = collect(keys(uniqcountout))
        sort!(uniqpredout; rev = true)
        ucountout = map((v -> uniqcountout[v]), uniqpredout)
    end

    aucin = begin
        sumlabel = sum(values(labelaggin))
        cumlabel = cumsum(map((p -> labelaggin[p]), uniqpredin))
        cumcount = cumsum(ucountin)
        map!((x -> convert(S, x)), cumcount, cumcount)
        tpr = (cumlabel ./ sumlabel)
        fpr = ((cumcount .- cumlabel) ./ (trainlentrue .- sumlabel))
        len = length(tpr)
        fpr1 = view(fpr, 1:(len-1))
        fpr2 = view(fpr, 2:len)
        tpr1 = view(tpr, 1:(len - 1))
        tpr2 = view(tpr, 2:len)
        area0 = fpr[1] * tpr[1]
        0.5 * (sum((tpr1 .+ tpr2) .* (fpr2 .- fpr1)) + area0)
    end

    aucout = validlen == 0 ? NaN64 : begin
        sumlabel = sum(values(labelaggout))
        cumlabel = cumsum(map((p -> labelaggout[p]), uniqpredout))
        cumcount = cumsum(ucountout)
        map!((x -> convert(S, x)), cumcount, cumcount)
        tpr = (cumlabel ./ sumlabel)
        fpr = ((cumcount .- cumlabel) ./ (validlentrue .- sumlabel))
        len = length(tpr)
        fpr1 = view(fpr, 1:(len-1))
        fpr2 = view(fpr, 2:len)
        tpr1 = view(tpr, 1:(len - 1))
        tpr2 = view(tpr, 2:len)
        area0 = fpr[1] * tpr[1]
        0.5 * (sum((tpr1 .+ tpr2) .* (fpr2 .- fpr1)) + area0)
    end
    aucin, aucout
end

function getlogloss(pred::Vector{T}, label::AbstractCovariate{S}, trainselector::AbstractBoolVariate, validselector::AbstractBoolVariate; slicelength::Integer = SLICELENGTH) where {T <: AbstractFloat} where {S <: AbstractFloat}
    fromobs = 1
    toobs = length(pred)
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    predslices = slice(pred, fromobs, toobs, slicelength)
    labelslices = slice(label, fromobs, toobs, slicelength)
    trainslices = slice(trainselector, fromobs, toobs, slicelength)
    validslices = slice(validselector, fromobs, toobs, slicelength)
    zipslices = zip(predslices, labelslices, trainslices, validslices)
    trainlosssum, validlosssum, traincount, validcount = fold((0.0, 0.0, 0, 0), zipslices) do acc, zipslice
         trainlosssum, validlosssum, traincount, validcount = acc
         predslice, labelslice, trainslice, validslice = zipslice
         for i in 1:length(predslice)
             if trainslice[i]
                trainlosssum += logloss(labelslice[i], predslice[i])
                traincount += 1
             end
             if validslice[i]
                validlosssum += logloss(labelslice[i], predslice[i])
                validcount += 1
             end
         end
         trainlosssum, validlosssum, traincount, validcount
    end
    trainlosssum / traincount, validlosssum / validcount

    # label = convert(Vector{S}, label)
    # label::Vector{S} = length(selector) == 0 ? label : label[selector]
    # pred::Vector{T} = length(selector) == 0 ? pred : pred[selector]
    # mean(logloss.(label, pred))
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Revise
using SQP
using PowerModels, JuMP, Ipopt
using filterSQP
using CPLEX

PowerModels.silence()

include("acwr.jl")
include("init_opf.jl")

function build_opf(pm::PowerModels.AbstractPowerModel)
    PowerModels.variable_bus_voltage(pm)
    PowerModels.variable_gen_power(pm)
    PowerModels.variable_branch_power(pm)
    PowerModels.variable_dcline_power(pm)

    PowerModels.objective_min_fuel_and_flow_cost(pm)

    PowerModels.constraint_model_voltage(pm)

    for i in PowerModels.ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in PowerModels.ids(pm, :bus)
        PowerModels.constraint_power_balance(pm, i)
    end

    for i in PowerModels.ids(pm, :branch)
        PowerModels.constraint_ohms_yt_from(pm, i)
        PowerModels.constraint_ohms_yt_to(pm, i)

        # constraint_voltage_angle_difference(pm, i)

        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end

    for i in PowerModels.ids(pm, :dcline)
        PowerModels.constraint_dcline_power_losses(pm, i)
    end
end

build_acp(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), ACPPowerModel, build_opf)
build_acr(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), ACRPowerModel, build_opf)
build_iv(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), IVRPowerModel, PowerModels.build_opf_iv)
build_dcp(data_file::String) = instantiate_model(PowerModels.parse_file(data_file), DCPPowerModel, PowerModels.build_opf_iv)

##
function run_sqp_opf(data_file::String, max_iter::Int = 100)
    pm = build_acr(data_file)
    # init_vars(pm)
    # pm2 = build_acp(data_file)
    # JuMP.@objective(pm2.model, Min, 0)
    # init_vars_from_ipopt(pm, pm2)

    # choose an internal QP solver
    qp_solver = optimizer_with_attributes(
        Ipopt.Optimizer,
        "print_level" => 0,
        "warm_start_init_point" => "yes",
        "linear_solver" => "ma57",
        # "ma57_pre_alloc" => 5.0,
        # CPLEX.Optimizer,
        # "CPX_PARAM_SCRIND" => 1,
        # "CPX_PARAM_THREADS" => 1,
        # "CPXPARAM_OptimalityTarget" => 2, # 1: convex, 2: local, 3: global
        # "CPXPARAM_Barrier_ConvergeTol" => 1.0e-4,
    )

    result = optimize_model!(pm, optimizer = optimizer_with_attributes(
        SQP.Optimizer,
        "algorithm" => "SQP-TR",
        "external_optimizer" => qp_solver,
        "tol_infeas" => 1.e-6,
        "tol_residual" => 1.e-4,
        "max_iter" => max_iter,
        "use_soc" => true,
    ))
    return pm, result
end

run_sqp_opf("../data/case9.m", 50);

##
function run_ipopt!(data_file::String)
    pm = build_acp(data_file)
    init_vars(pm)
    # pm2 = build_acp(data_file)
    # JuMP.@objective(pm2.model, Min, 0)
    # init_vars_from_ipopt(pm, pm2)
    solver = optimizer_with_attributes(
        Ipopt.Optimizer,
        "warm_start_init_point" => "yes",
        "linear_solver" => "ma57",
    )
    optimize_model!(pm, optimizer = solver)
    return
end

run_ipopt!("../data/case2869pegase.m");

# ##
# function run_filter_sqp!(data_file::String)
#     pm = build_acp(data_file)
#     init_vars(pm)
#     # pm2 = build_acp(data_file)
#     # JuMP.@objective(pm2.model, Min, 0)
#     # init_vars_from_ipopt(pm, pm2)
#     solver = optimizer_with_attributes(
#         filterSQP.Optimizer,
#         "iprint" => 1,
#     )
#     optimize_model!(pm, optimizer = solver)
#     return
# end

# run_filter_sqp!("../data/case1354pegase.m");

ext: jl


text:
function dotplot(; x::AbstractVector = Real[], y::AbstractVector = Real[], group::AbstractVector = Real[], sorted::Bool = true)

    v = VegaVisualization(width = 800, height = 400)

    #Quick defaults
    add_data!(v, y = y, x = x, group = group)
    default_axes!(v)
    default_scales!(v)

    if group != Real[]
      legend!(v)
    end

    #Get unique table name
    table = v.data[1].name

    #Switch data source
    v.scales[1].domain.data = v.scales[2].domain.data = "aggregate"
    v.scales[1].domain.field = "average_x"

    #Switch y scale for categorical, padding and data source
    v.scales[2]._type = "ordinal"
    v.scales[2].padding = 1

    #Add gridlines
    ylab!(v, grid = true, layer = "back")

    #Calculate average
    push!(v.data, VegaData(name = "aggregate", source = table,
                          transform = [VegaTransform(Dict{Any, Any}("type" => "aggregate", "groupby"=> ["y", "group"], "summarize"=> Dict{Any, Any}("x" => ["average"])))
                                      ]
                )
       )

    if sorted
        v.scales[2].domain.sort = Dict{Any, Any}("field" => "average_x", "op" => "mean")
        v.scales[2].reverse = true
    end

    v.marks = [VegaMark(_type = "symbol",
                        from = VegaMarkFrom(data = "aggregate"),
                        properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "average_x"),
                                                                                y = VegaValueRef(scale = "y", field = "y"),
                                                                                size = VegaValueRef(value = 50),
                                                                                shape = VegaValueRef(value = "circle"),
                                                                                stroke = VegaValueRef(field = "group", scale = "group"),
                                                                                strokeWidth = VegaValueRef(value = 2)
                                                                                )
                                                       )
      )]

    return v
end

max_stars_repo_name: jpthor/vega


text:
abstract type AbstractStateSpaceProblem <: DiffEqBase.DEProblem end
abstract type AbstractPerturbationProblem <: AbstractStateSpaceProblem end

using DiffEqBase: get_concrete_tspan, get_concrete_u0, get_concrete_p, promote_u0, promote_tspan,
                  isconcreteu0
# Perturbation problesm don't have f, g
# In discrete time, tspan should not have a sensitivity so the concretization is less obvious
function DiffEqBase.get_concrete_problem(prob::AbstractPerturbationProblem, isadapt; kwargs...)
    p = DiffEqBase.get_concrete_p(prob, kwargs)
    tspan = prob.tspan #get_concrete_tspan(prob, isadapt, kwargs, p)
    u0 = get_concrete_u0(prob, isadapt, tspan[1], kwargs)
    u0_promote = promote_u0(u0, p, tspan[1])

    if isconcreteu0(prob, tspan[1], kwargs) &&
       typeof(u0_promote) === typeof(prob.u0) &&
       p === prob.p &&
       return prob
    else
        return remake(prob; u0 = u0_promote, p = p)
    end
end

SciMLBase.isinplace(prob::AbstractPerturbationProblem) = false  # necessary for the get_concrete_u0 overloads

# the {iip} isn't relevant here at this point, but if we remove it then there are failures in the "remake" call above
# when using the Ensemble unit tests
struct LinearStateSpaceProblem{uType,uPriorType,tType,P,NP,F,AType,BType,CType,RType,ObsType,K} <:
       AbstractPerturbationProblem
    f::F # HACK: used only for standard interfaces/syms/etc., not used in calculations
    A::AType
    B::BType
    C::CType
    observables_noise::RType
    observables::ObsType
    u0::uType
    u0_prior::uPriorType
    tspan::tType
    p::P
    noise::NP
    kwargs::K
    @add_kwonly function LinearStateSpaceProblem{iip}(A, B, u0, tspan, p = NullParameters();
                                                      u0_prior = nothing, C = nothing,
                                                      observables_noise = nothing,
                                                      observables = nothing, noise = nothing,
                                                      syms = nothing,
                                                      f = ODEFunction{false}(((u, p, t) -> error("not implemented"));
                                                                             syms = syms),
                                                      kwargs...) where {iip}
        _tspan = promote_tspan(tspan)
        # _observables = promote_vv(observables)
        _observables = observables

        # Require integer distances between time periods for now.  Later could check with dt != 1
        @assert round(_tspan[2] - _tspan[1]) - (_tspan[2] - _tspan[1]) ≈ 0.0

        return new{typeof(u0),typeof(u0_prior),typeof(_tspan),typeof(p),typeof(noise),typeof(f),
                   typeof(A),typeof(B),typeof(C),typeof(observables_noise),typeof(_observables),
                   typeof(kwargs)}(f, A, B, C, observables_noise, _observables, u0, u0_prior,
                                   _tspan, p, noise, kwargs)
    end
end
# just forwards to a iip = false case
LinearStateSpaceProblem(args...; kwargs...) = LinearStateSpaceProblem{false}(args...; kwargs...)

# """
# u_f(t+1) = A_1 u_f(t) .+ B * noise(t+1)
# u(t+1) = A_0 + A_1 u(t) + quad(A_2, u_f(t)) .+ B noise(t+1)
# z(t) = C_0 + C_1 u(t) + quad(C_2, u_f(t))
# z_tilde(t) = z(t) + v(t+1)
# """
struct QuadraticStateSpaceProblem{uType,uPriorType,tType,P,NP,F,A0Type,A1Type,A2Type,BType,C0Type,
                                  C1Type,C2Type,RType,ObsType,K} <: AbstractPerturbationProblem
    f::F # HACK: used only for standard interfaces/syms/etc., not used in calculations
    A_0::A0Type
    A_1::A1Type
    A_2::A2Type
    B::BType
    C_0::C0Type
    C_1::C1Type
    C_2::C2Type
    observables_noise::RType
    observables::ObsType
    u0::uType
    u0_prior::uPriorType
    tspan::tType
    p::P
    noise::NP
    kwargs::K
    @add_kwonly function QuadraticStateSpaceProblem{iip}(A_0, A_1, A_2, B, u0, tspan,
                                                         p = NullParameters(); u0_prior = nothing,
                                                         C_0 = nothing, C_1 = nothing,
                                                         C_2 = nothing, observables_noise = nothing,
                                                         observables = nothing, noise = nothing,
                                                         syms = nothing,
                                                         f = ODEFunction{false}(((u, p, t) -> error("not implemented"));
                                                                                syms = syms),
                                                         kwargs...) where {iip}
        _tspan = promote_tspan(tspan)
        # _observables = promote_vv(observables)
        _observables = observables

        # Require integer distances between time periods for now.  Later could check with dt != 1
        @assert round(_tspan[2] - _tspan[1]) - (_tspan[2] - _tspan[1]) ≈ 0.0

        return new{typeof(u0),typeof(u0_prior),typeof(_tspan),typeof(p),typeof(noise),typeof(f),
                   typeof(A_0),typeof(A_1),typeof(A_2),typeof(B),typeof(C_0),typeof(C_1),
                   typeof(C_2),typeof(observables_noise),typeof(_observables),typeof(kwargs)}(f,
                                                                                A_0,
                                                                                A_1,
                                                                                A_2,
                                                                                B,
                                                                                C_0,
                                                                                C_1,
                                                                                C_2,
                                                                                observables_noise,
                                                                                _observables,
                                                                                u0,
                                                                                u0_prior,
                                                                                _tspan,
                                                                                p,
                                                                                noise,
                                                                                kwargs)
    end
end
# just forwards to a iip = false case
QuadraticStateSpaceProblem(args...; kwargs...) = QuadraticStateSpaceProblem{false}(args...;
                                                                                kwargs...)


text:
function test_add_reaction(model)
    value1 = 1.31592;
	value2 = 21.2434;
    tmp = deepcopy(model)

    # testing the first one
    test_model_1 = add_reaction(tmp,"test1", "pyr_e <=> pyr_c");
    @test_approx_eq value1 round(fba(test_model_1).obj,5)
    @test test_model_1 != tmp
    @test test_model_1 != model

    @test_approx_eq length(tmp.c) length(model.c)
    @test_approx_eq length(tmp.lb) length(model.lb)
    @test_approx_eq length(tmp.ub) length(model.ub)
    @test_approx_eq length(tmp.rxns) length(model.rxns)
    @test_approx_eq length(tmp.mets) length(model.mets)

    @test_approx_eq length(test_model_1.c) length(model.c) + 1
    @test_approx_eq length(test_model_1.lb) length(model.lb) + 1
    @test_approx_eq length(test_model_1.ub) length(model.ub) + 1
    @test_approx_eq length(test_model_1.rxns) length(model.rxns) + 1
    @test_approx_eq length(test_model_1.mets) length(model.mets)
    @test_approx_eq test_model_1.S.m tmp.S.m
    @test test_model_1.S != tmp.S
    @test_approx_eq test_model_1.S.n tmp.S.n + 1
    test_model_2 = add_reaction(tmp, "test1", "pyr_e <=> test_c");
    test_model_2 = add_reaction!(test_model_2, "test2", "test_e <=> test_c");
    test_model_2 = add_reaction!(test_model_2, "test3", "test_e <=> ");
    @test_approx_eq_eps(value2, fba(test_model_2).obj, 1e-4)
    @test test_model_2 != tmp

    @test_approx_eq length(tmp.c) length(model.c)
    @test_approx_eq length(tmp.lb) length(model.lb)
    @test_approx_eq length(tmp.ub) length(model.ub)
    @test_approx_eq length(tmp.rxns) length(model.rxns)
    @test_approx_eq length(tmp.mets) length(model.mets)
    @test test_model_2 != model
    @test_approx_eq length(test_model_2.c) length(model.c) + 3
    @test_approx_eq length(test_model_2.lb) length(model.lb) + 3
    @test_approx_eq length(test_model_2.ub) length(model.ub) + 3
    @test_approx_eq length(test_model_2.rxns) length(model.rxns) + 3
    @test_approx_eq length(test_model_2.mets) length(model.mets) + 2
    @test_approx_eq test_model_2.S.m tmp.S.m + 2
    @test test_model_2.S != tmp.S
    @test_approx_eq test_model_2.S.n tmp.S.n + 3
end

function test_change_objective(model)
	value = 20.0
	safe_model = deepcopy(model)

	safe_model.c = zeros(Float64, length(safe_model.c))
	safe_model.c[13] = 1
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 13, 1.0)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 13; objective = 1.0)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 13; objective = 1)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, model.rxns[13]; objective = 1)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, model.rxns[13], 1.0)).obj
	@test_approx_eq safe_model.c model.c

	safe_model.c = zeros(Float64, length(safe_model.c))
	safe_model.c[18] = 1
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 18, 1.0)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 18; objective = 1.0)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, 18; objective = 1)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, model.rxns[18]; objective = 1)).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, model.rxns[18], 1.0)).obj
	@test_approx_eq safe_model.c change_objective(model, model.rxns[18], 1).c
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, model.rxns[18], 1).c)

	safe_model.c = zeros(Float64, length(safe_model.c))
	safe_model.c[5] = 1
	safe_model.c[18] = 1
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, [5,18])).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, [5,18],[1,1])).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, [5,18],[1.0,1.0])).obj
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, [5,18]).c)
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, [5,18], [1,1]).c)
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, [5,18], [1.0,1.0]).c)

	safe_model.c = zeros(Float64, length(safe_model.c))
	safe_model.c[5] = 0.5
	safe_model.c[18] = 1
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, [5,18], [0.5, 1.0])).obj
	@test_approx_eq fba(safe_model).obj fba(change_objective(model, [5.0,18.0]; objectives = [0.5, 1.0])).obj
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, [5,18], [0.5, 1.0]).c)
	@test_approx_eq sum(safe_model.c) sum(change_objective(model, [5,18], [0.5, 1.0]).c)
end

function test_change_reaction_bounds(model)
	play_model = deepcopy(model)
	safe_model = deepcopy(model)
	safe_model.lb[13] = 0.5
	safe_model.ub[13] = 0.5
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5, "both")).obj
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "both").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "both").ub[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "l").ub[13] play_model.ub[13]
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "l").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "u").ub[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, "u").lb[13] play_model.lb[13]
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5, "u")).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5)).obj
	safe_model.ub[13] = 0.9
	@test abs(fba(safe_model).obj - fba(change_reaction_bounds(model, model.rxns[13], 0.5, "l")).obj) < 1e-6

	safe_model = deepcopy(model)
	safe_model.lb[13] = 0.5
	safe_model.ub[13] = 0.7
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "both")).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5, 0.7)).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, model.rxns[13], 0.5, 0.7, "u")).obj
	safe_model.ub[13] = 0.9
	@test abs(fba(safe_model).obj - fba(change_reaction_bounds(model, model.rxns[13], 0.5, 0.9, "l")).obj) < 1e-6

	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "both").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "both").ub[13] 0.7
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "l").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "l").ub[13] play_model.ub[13]
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "u").lb[13] play_model.lb[13]
	@test_approx_eq change_reaction_bounds(model, model.rxns[13], 0.5, 0.7,  "u").ub[13] 0.7

	safe_model = deepcopy(model)
	safe_model.lb[13] = 0.5
	safe_model.ub[13] = 0.7
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.7,  "both")).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.7,  "u")).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.7)).obj
	safe_model.lb[13] = 0.7
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.7, "both")).obj
	safe_model.ub[13] = 0.9
	@test abs(fba(safe_model).obj - fba(change_reaction_bounds(model, 13, 0.9, "u")).obj) < 1e-6
	@test abs(fba(safe_model).obj - fba(change_reaction_bounds(model, 13, 0.1, "l")).obj) < 1e-6
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "both").lb[13] 0.7
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "both").ub[13] 0.7
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "l").lb[13] 0.7
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "l").ub[13] play_model.ub[13]
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "u").lb[13] play_model.lb[13]
	@test_approx_eq change_reaction_bounds(model, 13, 0.7,  "u").ub[13] 0.7

	safe_model = deepcopy(model)
	safe_model.lb[13] = 0.5
	safe_model.ub[13] = 0.7
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.5, 0.7,  "both")).obj
	@test_approx_eq fba(safe_model).obj fba(change_reaction_bounds(model, 13, 0.5, 0.7,  "u")).obj
	safe_model.ub[13] = 0.9
	@test abs(fba(safe_model).obj - fba(change_reaction_bounds(model, 13, 0.5, 0.7,  "l")).obj) < 1e-6
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "both").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "both").ub[13] 0.7
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "l").lb[13] 0.5
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "l").ub[13] play_model.ub[13]
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "u").lb[13] play_model.lb[13]
	@test_approx_eq change_reaction_bounds(model, 13, 0.5, 0.7,  "u").ub[13] 0.7
end
function test_fba(model)
	value = 0.87392

    fbaA = fba(model).obj
	@test_approx_eq_eps value fbaA 1e-5
	#@test value == round(fba(model,"max","gurobi").obj,5)
	@test_approx_eq_eps value fba(model,"max").obj 1e-5
end
function test_find_blocked_reactions(model::Model)
	index = [26,27,29,34,45,47,52,63];
	values = find_blocked_reactions(model);
	@test values == model.rxns[index];
end
function test_find_deadend_metabolites(model)
	ans = [30;32;37;49;60]

	@test ans == find_deadend_metabolites(model)
end
function test_find_exchange_reactions(model::Model)
	index = [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39];
	values = find_exchange_reactions(model);

	@test index == values
end

function test_find_reactions_from_gene(model::Model)
	index = [2,14,58,69,70]
	values = find_reactions_from_gene(model, "s0001")

	@test index == values
end
function test_find_reactions_from_metabolite(model::Model)
	index = [4,5,11,12,13,15,16,18,41,46,52,53,54,58,62,76,79,81]
	value = find_reactions_from_metabolite(model, "h2o_c")

	@test index == value
end
function test_fva(model::Model)

    # this is a necessary fix
    # the ecoli_core.mat model
    # which the /CBM/test/test_vars/fva_test.mat
    # values test values are calculated
    # from has these two values swapped
    # relative to the e_coli_core.json model
    # so we must do that aswell
    # in addition, the ecoli_core.mat has
    # the upper limit for ATPM at 1000 while
    # e_coli_core.json has it at 8.39
    tmp = deepcopy(model)
    tmp.S.rowval[111] = 42
    tmp.S.rowval[112] = 44
    tmp.ub[11] = 1000

    fva_test = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/fva_test.mat")

    minF0 = fva_test["minF0"]
    maxF0 = fva_test["maxF0"]

    minF90 = fva_test["minF90"]
    maxF90 = fva_test["maxF90"]

    minF100 = fva_test["minF100"]
    maxF100 = fva_test["maxF100"]

    minF, maxF = fva(tmp, 0)
    @test_approx_eq minF minF0
    @test_approx_eq maxF maxF0

    minF, maxF = fva(tmp, 0.9)
    @test_approx_eq_eps minF minF90 1e-3
    @test_approx_eq_eps maxF maxF90 1e-3

    minF, maxF = fva(tmp, 1)
    @test_approx_eq_eps minF minF100 1e-3
    @test_approx_eq_eps maxF maxF100 1e-3
end
function test_gene_deletion(model)
    single_gene = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/sngl_gene.mat")

    matsol = single_gene["sol"]
    genes = single_gene["genelist"]

    sol = gene_deletion(model,1).g_f

    @test sum(matsol) - sum(collect(values(sol))) < 1e-6

    solution_order = []

    for gene in genes
        push!(solution_order, sol[[gene]])
    end

    @test !any(abs(matsol - solution_order) .> 1e-3)

    double_matrix = zeros(Float64, length(genes), length(genes))

    sol = gene_deletion(model,2).g_f

    for i in 1:length(genes)
        for j in i+1:length(genes)
            if haskey(sol, [genes[i], genes[j]])
                double_matrix[i,j] = sol[[genes[i], genes[j]]]
            elseif haskey(sol, [genes[j], genes[i]])
                double_matrix[i,j] = sol[[genes[j], genes[i]]]
            else
                warn("something wrong")
                break
            end
        end
    end

    double_matrix += double_matrix'

    for i in 1:length(genes)
        double_matrix[i,i] = sol[[genes[i]]]
    end

    double_gene = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/dbl_gene.mat")["ss"]

    @test !any(abs(double_gene - double_matrix) .> 1e-3)
end
function test_remove_reaction(model)
    eps = 1e-2

    test = deepcopy(model)
    tmp = deepcopy(model)

    tmp.lb[12] = 0
    tmp.ub[12] = 0
    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, model.rxns[12])).obj, eps)
    remove_reaction!(test, model.rxns[12])
    @test_approx_eq_eps(fba(tmp).obj, fba(test).obj, eps)
    @test_approx_eq length(model.rxns) length(test.rxns) + 1
    @test_approx_eq length(model.lb) length(test.lb) + 1
    @test_approx_eq length(model.ub) length(test.ub) + 1
    @test_approx_eq length(model.c) length(test.c) + 1
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 1
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 1
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 1

    tmp.lb[8] = 0
    tmp.ub[8] = 0

    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, model.rxns[8])).obj, eps)

    remove_reaction!(test, model.rxns[8])
    @test_approx_eq length(model.rxns) length(test.rxns) + 2
    @test_approx_eq length(model.lb) length(test.lb) + 2
    @test_approx_eq length(model.ub) length(test.ub) + 2
    @test_approx_eq length(model.c) length(test.c) + 2
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 2
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 2
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 2

    test = deepcopy(model)

    tmp = deepcopy(model)
    tmp.lb[12] = 0
    tmp.ub[12] = 0

    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, 12)).obj, eps)

    remove_reaction!(test, 12)
    @test_approx_eq length(model.rxns) length(test.rxns) + 1
    @test_approx_eq length(model.lb) length(test.lb) + 1
    @test_approx_eq length(model.ub) length(test.ub) + 1
    @test_approx_eq length(model.c) length(test.c) + 1
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 1
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 1
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 1

    tmp.lb[8] = 0
    tmp.ub[8] = 0
    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, 8)).obj, eps)
    remove_reaction!(test, 8)
    @test_approx_eq_eps(fba(tmp).obj, fba(test).obj, eps)
    @test_approx_eq length(model.rxns) length(test.rxns) + 2
    @test_approx_eq length(model.lb) length(test.lb) + 2
    @test_approx_eq length(model.ub) length(test.ub) + 2
    @test_approx_eq length(model.c) length(test.c) + 2
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 2
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 2
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 2

    test = deepcopy(model)
    tmp = deepcopy(model)
    tmp.lb[12] = 0
    tmp.ub[12] = 0
    tmp.lb[8] = 0
    tmp.ub[8] = 0
    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, [test.rxns[8], test.rxns[12]])).obj, eps)
    remove_reaction!(test, [test.rxns[8], test.rxns[12]])
    @test_approx_eq_eps(fba(tmp).obj, fba(test).obj, eps)
    @test_approx_eq length(model.rxns) length(test.rxns) + 2
    @test_approx_eq length(model.lb) length(test.lb) + 2
    @test_approx_eq length(model.ub) length(test.ub) + 2
    @test_approx_eq length(model.c) length(test.c) + 2
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 2
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 2
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 2

    test = deepcopy(model)
    tmp = deepcopy(model)
    tmp.lb[12] = 0
    tmp.ub[12] = 0
    tmp.lb[8] = 0
    tmp.ub[8] = 0
    @test_approx_eq_eps(fba(tmp).obj,fba(remove_reaction(test, [8,12])).obj, eps)
    remove_reaction!(test, [8,12])
    @test_approx_eq_eps(fba(tmp).obj, fba(test).obj, eps)
    @test_approx_eq length(model.rxns) length(test.rxns) + 2
    @test_approx_eq length(model.lb) length(test.lb) + 2
    @test_approx_eq length(model.ub) length(test.ub) + 2
    @test_approx_eq length(model.c) length(test.c) + 2
    @test_approx_eq length(model.rxn_rules) length(test.rxn_rules) + 2
    @test_approx_eq length(model.rxn_name) length(test.rxn_name) + 2
    @test_approx_eq length(model.rxn_subsystem) length(test.rxn_subsystem) + 2
end
function test_lp(lp, model)
    tmp = clone(model)
    tol = 1e-6
    @test_approx_eq_eps answer_lp(lp) fba(tmp).obj tol

    change_objective_coef(lp, 11)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 11)).obj tol
    change_objective_coef(lp, 28)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 28)).obj tol
    change_objective_coef(lp, 7)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 7)).obj tol
    change_objective_coef(lp, 13)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 13)).obj tol

    change_objective_coef(lp, 11, 0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 11, 0.5)).obj tol
    change_objective_coef(lp, 28, 0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 28, 0.5)).obj tol
    change_objective_coef(lp, 7, 0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 7, 0.5)).obj tol
    change_objective_coef(lp, 13, 0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 13, 0.5)).obj tol

    change_objective_coef(lp, [5,7])
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, [5,7])).obj tol
    @test get_objective(lp) == change_objective(tmp, [5,7]).c

    change_objective_coef(lp, [5,7],[1,1])
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, [5,7],[1,1])).obj tol
    @test get_objective(lp) == change_objective(tmp, [5,7],[1,1]).c

    change_objective_coef(lp, [5,7],[0.5,0.7])
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, [5,7],[0.5,0.7])).obj tol
    @test get_objective(lp) == change_objective(tmp, [5,7],[0.5,0.7]).c

    change_objective_coef(lp, [13,16],[0.5,0.7])
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, [13,16],[0.5,0.7])).obj tol
    @test get_objective(lp) == change_objective(tmp, [13,16],[0.5,0.7]).c

    change_objective_coef(lp, [7,5],[0.5,0.7])
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, [7,5],[0.5,0.7])).obj tol
    @test get_objective(lp) == change_objective(tmp, [7,5],[0.5,0.7]).c

    change_objective_coef(lp, 13)
    @test_approx_eq_eps answer_lp(lp) fba(change_objective(tmp, 13)).obj tol
    @test get_objective(lp) == change_objective(tmp, 13).c

    @test get_variable_bounds(lp) == (tmp.lb, tmp.ub)

    set_col_bounds(lp, 5, -500,500)
    @test get_variable_bounds(lp)[1] == (change_reaction_bounds(tmp, 5, -500, 500).lb)
    @test get_variable_bounds(lp)[2] == (change_reaction_bounds(tmp, 5, -500, 500).ub)

    set_col_bounds(lp, 5, 0,0)
    @test get_variable_bounds(lp)[1] == (change_reaction_bounds(tmp, 5, 0, 0).lb)
    @test get_variable_bounds(lp)[2] == (change_reaction_bounds(tmp, 5, 0, 0).ub)

    set_col_bounds(lp, 5, 0)
    @test get_variable_bounds(lp)[1] == (change_reaction_bounds(tmp, 5, 0).lb)
    @test get_variable_bounds(lp)[2] == (change_reaction_bounds(tmp, 5, 0).ub)

    set_col_bounds(lp, 5, tmp.lb[5], tmp.ub[5])
    @test_approx_eq_eps answer_lp(lp) fba(tmp).obj tol

    set_col_bounds(lp, 13, 0.5,0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_reaction_bounds(tmp, 13, 0.5, 0.5)).obj tol

    set_col_bounds(lp, 13, tmp.lb[13], tmp.ub[13])

    set_col_bounds(lp, 13, 0.5)
    @test_approx_eq_eps answer_lp(lp) fba(change_reaction_bounds(tmp, 13, 0.5)).obj tol

    set_col_bounds(lp, 13, tmp.lb[13], tmp.ub[13])
    @test_approx_eq_eps answer_lp(lp) fba(tmp).obj tol
end

function test_robustness_analysis(model)
    matsol = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/robustness.mat")
    single = matsol["single"]
    double = matsol["double"]

    sol = robustness_analysis(model, [8]).result
    @test_approx_eq sol single

    sol = robustness_analysis(model, [5,8]).result'
    @test_approx_eq sol double
end
function test_fast_cc()
    testfile = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/fastcctest.mat")
    A = testfile["testA"]

    tmp = load_matlab(Pkg.dir() * "/CBM/test/test_vars/consistRecon1.mat")

    @test A == fast_cc(tmp, 1e-3)
end
function test_fast_core()
    testfile = open_mat_file(Pkg.dir() * "/CBM/test/test_vars/fastcoretest.mat")
    C = testfile["testC"]
    A = testfile["testA"]

    tmp = load_matlab(Pkg.dir() * "/CBM/test/test_vars/consistRecon1.mat")

    @test A == fast_core(C, tmp, 1e-3)
end
function test_solvers(model)
    supported_solvers = ["GLPK", "CPLEX", "Gurobi", "Clp"]
    available_solvers = map(x -> any(x .== readdir(Pkg.dir())), supported_solvers)

    if available_solvers[1]
        test_lp(setup_glpk(model), model)
    end
    if available_solvers[2]
        test_lp(setup_cplex(model), model)
    end
    if available_solvers[3]
        test_lp(setup_gurobi(model), model)
    end
    if available_solvers[4]
        test_lp(setup_clp(model), model)
    end
end

max_stars_repo_name: isebarn/CBM

max_issues_repo_name: isebarn/CBM

max_forks_repo_name: isebarn/CBM


text:
__precompile__()

module BayesianSumProductNetworks

    import Base.eltype
    import Base.length
    import Base.rand
    import Base.copy

    # Required packages
    using Reexport
    @reexport using Distributions
    using AxisArrays
    using SpecialFunctions
    using StatsFuns
    import StatsFuns.normlogpdf

    # Required base packages.
    using Statistics
    using Random
    using SparseArrays

    # Exports
    export NormalInverseGamma
    export DirichletSufficientStats, NormalInvGammaSufficientStats
    export GammaSufficienStats, BetaSufficientStats
    export AbstractSufficientStats
    export postparams, update!
    export RegionGraphNode, PartitionGraphNode
    export FactorizedDistributionGraphNode, AbstractRegionGraphNode
    export FactorizedMixtureGraphNode, InfiniteSumNode
    export haschildren, topoligicalorder, id
    export scope, randomscopes!, nscope, hasscope, setscope!, updatescope!
    export haschildren, children
    export logpdf!, logmllh!, logpdf

    ######################
    ## Type Definitions ##
    ######################

    abstract type AbstractSufficientStats end

    """
        Gamma prior sufficient stats.
    """
    mutable struct GammaSufficienStats <: AbstractSufficientStats
        n::Int
        sn::Float64
    end

    GammaSufficienStats() = GammaSufficienStats(0, 0.0)

    function update!(s::GammaSufficienStats, x::AbstractVector)
        if isempty(x)
            s.n = 0
            s.sn = 0.0
        else
            @assert all(x .>= 0) "Not all positive, $x"
            s.n = length(x)
            s.sn = sum(x)
        end
        return s
    end

    @inline postparams(d::Gamma, T::GammaSufficienStats) = (d.α+T.sn, (1/d.θ)+T.n)

    """
        Beta prior sufficient stats.
    """
    mutable struct BetaSufficientStats <: AbstractSufficientStats
        n::Int
        count::Int
    end

    function update!(s::BetaSufficientStats, x::T) where {T<:AbstractVector{<:Real}}
        if isempty(x)
            s.n = 0
            s.count = 0
        else
            @assert all((y) -> isone(y) | iszero(y), x)
            s.n = length(x)
            s.count = sum(x)
        end
        return s
    end
    @inline BetaSufficientStats() = BetaSufficientStats(0,0)
    @inline postparams(d::Beta, T::BetaSufficientStats) = (T.count+d.α,(T.n-T.count)+d.β)

    """
        Dirichlet prior sufficient stats.
    """
    struct DirichletSufficientStats <: AbstractSufficientStats
        counts::Vector{Int}
    end
    function update!(s::DirichletSufficientStats, x::AbstractVector)
        if isempty(x)
            fill!(s.counts, 0)
        else
            s.counts[:] = map(k -> sum(x .== k), 1:length(s.counts))
        end
        return s
    end
    @inline DirichletSufficientStats(n::Int) = DirichletSufficientStats(zeros(Int, n))
    @inline postparams(d::Dirichlet, T::DirichletSufficientStats) = T.counts.+d.alpha

    ## Gaussian with Normal Inverse Gamma prior ##
    struct NormalInverseGamma <: Distribution{Multivariate,Continuous}
        μ::Float64
        ν::Float64
        a::Float64
        b::Float64
    end

    length(d::NormalInverseGamma) = 2

    function logpdf(d::NormalInverseGamma, θ...)
        (μ, σ²) = θ
        l = _invgammalogpdf(d.a, d.b, σ²)
        l += normlogpdf(d.μ, sqrt(σ²) / d.ν, μ)
        return l
    end
    @inline _invgammalogpdf(a::Real, b::Real, x::Real) = a*log(b)-lgamma(a)-(a+1)*log(x)-b/x

    """
        Normal-Inverse-Gamma prior sufficient stats.
    """
    mutable struct NormalInvGammaSufficientStats <: AbstractSufficientStats
        n::Int
        s::Float64
        sdiff::Float64
    end

    NormalInvGammaSufficientStats() = NormalInvGammaSufficientStats(0,0.0,0.0)

    function update!(s::NormalInvGammaSufficientStats, x::AbstractVector)
        if isempty(x)
            s.n = 0
            s.s = 0
            s.sdiff = 0
        else
            s.n = length(x)
            s.s = sum(x)
            s.sdiff = sum((x .- mean(x)).^2)
        end
        s
    end

    function rand(d::NormalInverseGamma)
        σ² = rand(InverseGamma(d.a, d.b))
        μ = rand(Normal(d.μ, sqrt(σ² / d.ν )))
        return [μ, σ²]
    end

    @inline params(d::NormalInverseGamma) = (d.μ, d.ν, d.a, d.b)
    function postparams(d::NormalInverseGamma, T::NormalInvGammaSufficientStats)

        if T.n == 0
            return (d.μ, d.ν, d.a, d.b)
        end

        μ = (d.ν*d.μ + T.s) / (d.ν + T.n)
        ν = d.ν + T.n
        a = d.a + (T.n / 2)
        b = d.b
        b += T.sdiff/2 + ((T.n * d.ν) / (d.ν + T.n)) * (((T.s / T.n) - d.μ)^2 / 2)

        return (μ, ν, a, b)
    end

    # include implementations
    include("mathutil.jl")
    include("distUtils.jl")
    include("templateStructure.jl")
    include("regiongraph.jl")
    include("nonparametrics.jl")
    include("ancestralsampler.jl")
    include("initialization.jl")
    include("structureinference.jl")

end # module

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
using Optim

## initial values
x0 = -5.0*ones(1)

## Our objective function
objfun(x) = x[1]^4 - 3.0*x[1]^3 + 2.0

## Find the optimum!
result = optimize(objfun, x0, LBFGS())
println("minimum value is: ",result.minimizer)
println("obj value at minimum is: ",result.minimum)

## Now do Particle Swarm (derivative-free method)
result = optimize(objfun, x0, ParticleSwarm())
println("minimum value is: ",result.minimizer)
println("obj value at minimum is: ",result.minimum)


text:
################################################################################

# Parasagittal/supra-sylvian electrodes
ParasagittalL = [
"Fp1",
"F3",
"C3",
"P3",
"O1",
]

ParasagittalR = [
"Fp2",
"F4",
"C4",
"P4",
"O2",
]

# Lateral/temporal electrodes
TemporalL = [
"F7",
"T3", # "T7",
"T5", # "P7",
]

TemporalR = [
"F8",
"T4", # "T8",
"T6", # "P8",
]

# Midline electrodes
Midline = [
"Fz",
"Cz",
"Pz",
]

# Earlobe electrodes
Earlobe = [
"A1",
"A2",
]

Bipolar = [
"FP1-F7",
"F7-T7",
"T7-P7",
"P7-O1",
"FP1-F3",
"F3-C3",
"C3-P3",
"P3-O1",
"FP2-F4",
"F4-C4",
"C4-P4",
"P4-O2",
"FP2-F8",
"F8-T8",
"T8-P8",
"P8-O2",
"FZ-CZ",
"CZ-PZ",
"P7-T7",
"T7-FT9",
"FT9-FT10",
"FT10-T8",
]

elecID = [ParasagittalL;
ParasagittalR;
TemporalL;
TemporalR;
Midline;
Earlobe;
Bipolar
]

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

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


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

using Markdown
using InteractiveUtils

# ╔═╡ 9353168e-bd57-446f-9fc3-ce063aa9c0f5
using Plots, PlutoTest

# ╔═╡ 4049e274-f948-468f-9e43-71855a426603
plotly()

# ╔═╡ 1b64f670-b017-11eb-140d-ebd5ac033872
md"""
# [Day 6: Probably a Fire Hazard](https://adventofcode.com/2015/day/6)
Because your neighbors keep defeating you in the holiday house decorating contest year after year, you've decided to deploy one million lights in a 1000x1000 grid.

Furthermore, because you've been especially nice this year, Santa has mailed you instructions on how to display the ideal lighting configuration.

Lights in your grid are numbered from 0 to 999 in each direction; the lights at each corner are at `0,0`, `0,999`, `999,999`, and `999,0`. The instructions include whether to `turn on`, `turn off`, or `toggle` various inclusive ranges given as coordinate pairs. Each coordinate pair represents opposite corners of a rectangle, inclusive; a coordinate pair like `0,0 through 2,2` therefore refers to 9 lights in a 3x3 square. The lights all start turned off.

To defeat your neighbors this year, all you have to do is set up your lights by doing the instructions Santa sent you in order.

For example:

-   `turn on 0,0 through 999,999` would turn on (or leave on) every light.
-   `toggle 0,0 through 999,0` would toggle the first line of 1000 lights, turning off the ones that were on, and turning on the ones that were off.
-   `turn off 499,499 through 500,500` would turn off (or leave off) the middle four lights.

**After following the instructions, how many lights are lit?**
"""

# ╔═╡ c87dc276-6a3a-4903-955e-6ee8063b0f53
puzzle_input = "turn on 887,9 through 959,629
turn on 454,398 through 844,448
turn off 539,243 through 559,965
turn off 370,819 through 676,868
turn off 145,40 through 370,997
turn off 301,3 through 808,453
turn on 351,678 through 951,908
toggle 720,196 through 897,994
toggle 831,394 through 904,860
toggle 753,664 through 970,926
turn off 150,300 through 213,740
turn on 141,242 through 932,871
toggle 294,259 through 474,326
toggle 678,333 through 752,957
toggle 393,804 through 510,976
turn off 6,964 through 411,976
turn off 33,572 through 978,590
turn on 579,693 through 650,978
turn on 150,20 through 652,719
turn off 782,143 through 808,802
turn off 240,377 through 761,468
turn off 899,828 through 958,967
turn on 613,565 through 952,659
turn on 295,36 through 964,978
toggle 846,296 through 969,528
turn off 211,254 through 529,491
turn off 231,594 through 406,794
turn off 169,791 through 758,942
turn on 955,440 through 980,477
toggle 944,498 through 995,928
turn on 519,391 through 605,718
toggle 521,303 through 617,366
turn off 524,349 through 694,791
toggle 391,87 through 499,792
toggle 562,527 through 668,935
turn off 68,358 through 857,453
toggle 815,811 through 889,828
turn off 666,61 through 768,87
turn on 27,501 through 921,952
turn on 953,102 through 983,471
turn on 277,552 through 451,723
turn off 64,253 through 655,960
turn on 47,485 through 734,977
turn off 59,119 through 699,734
toggle 407,898 through 493,955
toggle 912,966 through 949,991
turn on 479,990 through 895,990
toggle 390,589 through 869,766
toggle 593,903 through 926,943
toggle 358,439 through 870,528
turn off 649,410 through 652,875
turn on 629,834 through 712,895
toggle 254,555 through 770,901
toggle 641,832 through 947,850
turn on 268,448 through 743,777
turn off 512,123 through 625,874
turn off 498,262 through 930,811
turn off 835,158 through 886,242
toggle 546,310 through 607,773
turn on 501,505 through 896,909
turn off 666,796 through 817,924
toggle 987,789 through 993,809
toggle 745,8 through 860,693
toggle 181,983 through 731,988
turn on 826,174 through 924,883
turn on 239,228 through 843,993
turn on 205,613 through 891,667
toggle 867,873 through 984,896
turn on 628,251 through 677,681
toggle 276,956 through 631,964
turn on 78,358 through 974,713
turn on 521,360 through 773,597
turn off 963,52 through 979,502
turn on 117,151 through 934,622
toggle 237,91 through 528,164
turn on 944,269 through 975,453
toggle 979,460 through 988,964
turn off 440,254 through 681,507
toggle 347,100 through 896,785
turn off 329,592 through 369,985
turn on 931,960 through 979,985
toggle 703,3 through 776,36
toggle 798,120 through 908,550
turn off 186,605 through 914,709
turn off 921,725 through 979,956
toggle 167,34 through 735,249
turn on 726,781 through 987,936
toggle 720,336 through 847,756
turn on 171,630 through 656,769
turn off 417,276 through 751,500
toggle 559,485 through 584,534
turn on 568,629 through 690,873
toggle 248,712 through 277,988
toggle 345,594 through 812,723
turn off 800,108 through 834,618
turn off 967,439 through 986,869
turn on 842,209 through 955,529
turn on 132,653 through 357,696
turn on 817,38 through 973,662
turn off 569,816 through 721,861
turn on 568,429 through 945,724
turn on 77,458 through 844,685
turn off 138,78 through 498,851
turn on 136,21 through 252,986
turn off 2,460 through 863,472
turn on 172,81 through 839,332
turn on 123,216 through 703,384
turn off 879,644 through 944,887
toggle 227,491 through 504,793
toggle 580,418 through 741,479
toggle 65,276 through 414,299
toggle 482,486 through 838,931
turn off 557,768 through 950,927
turn off 615,617 through 955,864
turn on 859,886 through 923,919
turn on 391,330 through 499,971
toggle 521,835 through 613,847
turn on 822,787 through 989,847
turn on 192,142 through 357,846
turn off 564,945 through 985,945
turn off 479,361 through 703,799
toggle 56,481 through 489,978
turn off 632,991 through 774,998
toggle 723,526 through 945,792
turn on 344,149 through 441,640
toggle 568,927 through 624,952
turn on 621,784 through 970,788
toggle 665,783 through 795,981
toggle 386,610 through 817,730
toggle 440,399 through 734,417
toggle 939,201 through 978,803
turn off 395,883 through 554,929
turn on 340,309 through 637,561
turn off 875,147 through 946,481
turn off 945,837 through 957,922
turn off 429,982 through 691,991
toggle 227,137 through 439,822
toggle 4,848 through 7,932
turn off 545,146 through 756,943
turn on 763,863 through 937,994
turn on 232,94 through 404,502
turn off 742,254 through 930,512
turn on 91,931 through 101,942
toggle 585,106 through 651,425
turn on 506,700 through 567,960
turn off 548,44 through 718,352
turn off 194,827 through 673,859
turn off 6,645 through 509,764
turn off 13,230 through 821,361
turn on 734,629 through 919,631
toggle 788,552 through 957,972
toggle 244,747 through 849,773
turn off 162,553 through 276,887
turn off 569,577 through 587,604
turn off 799,482 through 854,956
turn on 744,535 through 909,802
toggle 330,641 through 396,986
turn off 927,458 through 966,564
toggle 984,486 through 986,913
toggle 519,682 through 632,708
turn on 984,977 through 989,986
toggle 766,423 through 934,495
turn on 17,509 through 947,718
turn on 413,783 through 631,903
turn on 482,370 through 493,688
turn on 433,859 through 628,938
turn off 769,549 through 945,810
turn on 178,853 through 539,941
turn off 203,251 through 692,433
turn off 525,638 through 955,794
turn on 169,70 through 764,939
toggle 59,352 through 896,404
toggle 143,245 through 707,320
turn off 103,35 through 160,949
toggle 496,24 through 669,507
turn off 581,847 through 847,903
turn on 689,153 through 733,562
turn on 821,487 through 839,699
turn on 837,627 through 978,723
toggle 96,748 through 973,753
toggle 99,818 through 609,995
turn on 731,193 through 756,509
turn off 622,55 through 813,365
turn on 456,490 through 576,548
turn on 48,421 through 163,674
turn off 853,861 through 924,964
turn off 59,963 through 556,987
turn on 458,710 through 688,847
toggle 12,484 through 878,562
turn off 241,964 through 799,983
turn off 434,299 through 845,772
toggle 896,725 through 956,847
turn on 740,289 through 784,345
turn off 395,840 through 822,845
turn on 955,224 through 996,953
turn off 710,186 through 957,722
turn off 485,949 through 869,985
turn on 848,209 through 975,376
toggle 221,241 through 906,384
turn on 588,49 through 927,496
turn on 273,332 through 735,725
turn on 505,962 through 895,962
toggle 820,112 through 923,143
turn on 919,792 through 978,982
toggle 489,461 through 910,737
turn off 202,642 through 638,940
turn off 708,953 through 970,960
toggle 437,291 through 546,381
turn on 409,358 through 837,479
turn off 756,279 through 870,943
turn off 154,657 through 375,703
turn off 524,622 through 995,779
toggle 514,221 through 651,850
toggle 808,464 through 886,646
toggle 483,537 through 739,840
toggle 654,769 through 831,825
turn off 326,37 through 631,69
turn off 590,570 through 926,656
turn off 881,913 through 911,998
turn on 996,102 through 998,616
turn off 677,503 through 828,563
turn on 860,251 through 877,441
turn off 964,100 through 982,377
toggle 888,403 through 961,597
turn off 632,240 through 938,968
toggle 731,176 through 932,413
turn on 5,498 through 203,835
turn on 819,352 through 929,855
toggle 393,813 through 832,816
toggle 725,689 through 967,888
turn on 968,950 through 969,983
turn off 152,628 through 582,896
turn off 165,844 through 459,935
turn off 882,741 through 974,786
turn off 283,179 through 731,899
toggle 197,366 through 682,445
turn on 106,309 through 120,813
toggle 950,387 through 967,782
turn off 274,603 through 383,759
turn off 155,665 through 284,787
toggle 551,871 through 860,962
turn off 30,826 through 598,892
toggle 76,552 through 977,888
turn on 938,180 through 994,997
toggle 62,381 through 993,656
toggle 625,861 through 921,941
turn on 685,311 through 872,521
turn on 124,934 through 530,962
turn on 606,379 through 961,867
turn off 792,735 through 946,783
turn on 417,480 through 860,598
toggle 178,91 through 481,887
turn off 23,935 through 833,962
toggle 317,14 through 793,425
turn on 986,89 through 999,613
turn off 359,201 through 560,554
turn off 729,494 through 942,626
turn on 204,143 through 876,610
toggle 474,97 through 636,542
turn off 902,924 through 976,973
turn off 389,442 through 824,638
turn off 622,863 through 798,863
turn on 840,622 through 978,920
toggle 567,374 through 925,439
turn off 643,319 through 935,662
toggle 185,42 through 294,810
turn on 47,124 through 598,880
toggle 828,303 through 979,770
turn off 174,272 through 280,311
turn off 540,50 through 880,212
turn on 141,994 through 221,998
turn on 476,695 through 483,901
turn on 960,216 through 972,502
toggle 752,335 through 957,733
turn off 419,713 through 537,998
toggle 772,846 through 994,888
turn on 881,159 through 902,312
turn off 537,651 through 641,816
toggle 561,947 through 638,965
turn on 368,458 through 437,612
turn on 290,149 through 705,919
turn on 711,918 through 974,945
toggle 916,242 through 926,786
toggle 522,272 through 773,314
turn on 432,897 through 440,954
turn off 132,169 through 775,380
toggle 52,205 through 693,747
toggle 926,309 through 976,669
turn off 838,342 through 938,444
turn on 144,431 through 260,951
toggle 780,318 through 975,495
turn off 185,412 through 796,541
turn on 879,548 through 892,860
turn on 294,132 through 460,338
turn on 823,500 through 899,529
turn off 225,603 through 483,920
toggle 717,493 through 930,875
toggle 534,948 through 599,968
turn on 522,730 through 968,950
turn off 102,229 through 674,529"

# ╔═╡ 462f3b05-97a4-4bc7-b817-0b74be141531
dims = (1000, 1000)

# ╔═╡ bdb2415c-bef7-4d31-81cf-4175679b8814
parse_input(input) = split(input, "\n")

# ╔═╡ cfeb0f48-8879-4c32-882a-6cfc5fbe5e78
grid = zeros(dims)

# ╔═╡ b64e92eb-2505-4203-a512-6292235ccea5
apply1!(grid::BitMatrix, command::AbstractString) = begin
	parts = split(command, " ")
	parse_value(str::AbstractString) = split(str, ",") |>
			Base.Fix1(map, s->parse(Int, s)) |>
			Base.Fix1(map, s->s+1)
	not(x::Bool) = !x

	if parts[1] == "turn"
		from = parse_value(parts[3])
		to = parse_value(parts[5])

		if parts[2] == "on"
			grid[from[1]:to[1], from[2]:to[2]] .= 1
		elseif parts[2] == "off"
			grid[from[1]:to[1], from[2]:to[2]] .= 0
		else
			throw(ErrorException("invalid syntax"))
		end
	elseif parts[1] == "toggle"
		from = parse_value(parts[2])
		to = parse_value(parts[4])
		grid[from[1]:to[1], from[2]:to[2]] = not.(grid[from[1]:to[1], from[2]:to[2]])
	else
		throw(ErrorException("invalid syntax"))
	end

	grid
end

# ╔═╡ 238fc1ab-c208-4f54-9d48-552fba424020
puzzle_commands = parse_input(puzzle_input)

# ╔═╡ fa569771-87e4-4534-83d3-22da3fe425f7
@assert sum(apply1!(falses(dims), "turn on 0,0 through 999,999")) == 1000*1000

# ╔═╡ 2d924161-974a-4209-bf72-5632dcc3e5b5
@assert sum(apply1!(falses(dims), "toggle 0,0 through 999,0")) == 1000

# ╔═╡ e3e5aa80-54bf-4d7d-912a-00d1d7061ea7
@assert sum(apply1!(trues(dims), "toggle 0,0 through 999,0")) == 1000*1000-1000

# ╔═╡ ea5f89aa-012b-4d23-8516-3d26269e27c5
@assert sum(apply1!(trues(dims), "turn off 499,499 through 500,500")) == 1000*1000-4

# ╔═╡ 41e0a047-f7a4-4e7c-a9b8-6a3891e52465
function mass_apply1(commands::Vector)
	grid = falses(dims)
	for command in commands
		apply1!(grid, command)
	end
	grid
end

# ╔═╡ 175908f4-2ef8-44a9-81c0-fb21e4bc10b7
grid1 = mass_apply1(split(puzzle_input, "\n"))

# ╔═╡ 184bb620-2d32-4b98-a02d-f5b64ae88eb6
heatmap(grid1)

# ╔═╡ 9840d540-607e-47d3-874c-9df7ec814d51
part1 = Int(sum(grid1))

# ╔═╡ 08f7886a-9cef-4b23-9885-8dab7066640a
@test part1 == 377891

# ╔═╡ 8ada2165-6667-4850-8131-e56b1fde7abe
md"Your puzzle answer was `377891`."

# ╔═╡ b3b7099c-b2c2-4ff4-b5c2-dfa6ec3cda87
md"""
# Part Two

You just finish implementing your winning light pattern when you realize you mistranslated Santa's message from Ancient Nordic Elvish.

The light grid you bought actually has individual brightness controls; each light can have a brightness of zero or more. The lights all start at zero.

The phrase `turn on` actually means that you should increase the brightness of those lights by `1`.

The phrase `turn off` actually means that you should decrease the brightness of those lights by `1`, to a minimum of zero.

The phrase `toggle` actually means that you should increase the brightness of those lights by `2`.

**What is the total brightness of all lights combined after following Santa's instructions?**

For example:

-   `turn on 0,0 through 0,0` would increase the total brightness by `1`.
-   `toggle 0,0 through 999,999` would increase the total brightness by `2000000`.

"""

# ╔═╡ eea139c9-3d08-499a-a9fb-df4f20bd129a
apply2!(grid::Matrix, command::AbstractString) = begin
	parts = split(command, " ")
	parse_value(str::AbstractString) = split(str, ",") |>
			Base.Fix1(map, s->parse(Int, s)) |>
			Base.Fix1(map, s->s+1)
	reduce_to_zero(n::Float64) = n > 0 ? n - 1 : 0

	if parts[1] == "turn"
		from = parse_value(parts[3])
		to = parse_value(parts[5])

		if parts[2] == "on"
			grid[from[1]:to[1], from[2]:to[2]] .+= 1
		elseif parts[2] == "off"
			grid[from[1]:to[1], from[2]:to[2]] = reduce_to_zero.(grid[from[1]:to[1], from[2]:to[2]])
		else
			throw(ErrorException("invalid syntax"))
		end
	elseif parts[1] == "toggle"
		from = parse_value(parts[2])
		to = parse_value(parts[4])
		grid[from[1]:to[1], from[2]:to[2]] .+= 2
	else
		throw(ErrorException("invalid syntax"))
	end

	grid
end

# ╔═╡ ee1b5584-8e75-4ebf-9514-0e5985b9869f
@assert sum(apply2!(zeros(dims), "turn on 0,0 through 0,0")) == 1

# ╔═╡ 13f0e9af-d1dd-40ec-99a8-c56b91ab6ce8
@assert sum(apply2!(zeros(dims), "toggle 0,0 through 999,999")) == 1000*1000*2

# ╔═╡ 7b0a6f85-0637-4ec0-8d17-3984bd9d8969
function mass_apply2(commands::Vector)
	grid = zeros(dims)
	for command in commands
		apply2!(grid, command)
	end
	grid
end

# ╔═╡ 4d5a5013-fd59-4fea-8f73-396fba062fc5
grid2 = mass_apply2(split(puzzle_input, "\n"))

# ╔═╡ fbae4e2b-2ea5-414a-b3fd-01cd5630bfda
heatmap(grid2)

# ╔═╡ f245c971-fa2d-482e-a355-67a51963c0cb
part2 = Int(sum(grid2))

# ╔═╡ 6e303fd4-5539-4ec6-9bb8-72f366904f1f
@test part2 == 14110788

# ╔═╡ 0f036ef3-4901-42e4-bf0c-0580215b20cb
md"Your puzzle answer was `14110788`."

# ╔═╡ Cell order:
# ╠═9353168e-bd57-446f-9fc3-ce063aa9c0f5
# ╠═4049e274-f948-468f-9e43-71855a426603
# ╟─1b64f670-b017-11eb-140d-ebd5ac033872
# ╟─c87dc276-6a3a-4903-955e-6ee8063b0f53
# ╠═462f3b05-97a4-4bc7-b817-0b74be141531
# ╠═bdb2415c-bef7-4d31-81cf-4175679b8814
# ╠═cfeb0f48-8879-4c32-882a-6cfc5fbe5e78
# ╠═b64e92eb-2505-4203-a512-6292235ccea5
# ╠═238fc1ab-c208-4f54-9d48-552fba424020
# ╠═fa569771-87e4-4534-83d3-22da3fe425f7
# ╠═2d924161-974a-4209-bf72-5632dcc3e5b5
# ╠═e3e5aa80-54bf-4d7d-912a-00d1d7061ea7
# ╠═ea5f89aa-012b-4d23-8516-3d26269e27c5
# ╠═41e0a047-f7a4-4e7c-a9b8-6a3891e52465
# ╠═175908f4-2ef8-44a9-81c0-fb21e4bc10b7
# ╠═184bb620-2d32-4b98-a02d-f5b64ae88eb6
# ╠═9840d540-607e-47d3-874c-9df7ec814d51
# ╠═08f7886a-9cef-4b23-9885-8dab7066640a
# ╟─8ada2165-6667-4850-8131-e56b1fde7abe
# ╟─b3b7099c-b2c2-4ff4-b5c2-dfa6ec3cda87
# ╠═eea139c9-3d08-499a-a9fb-df4f20bd129a
# ╠═ee1b5584-8e75-4ebf-9514-0e5985b9869f
# ╠═13f0e9af-d1dd-40ec-99a8-c56b91ab6ce8
# ╠═7b0a6f85-0637-4ec0-8d17-3984bd9d8969
# ╠═4d5a5013-fd59-4fea-8f73-396fba062fc5
# ╠═fbae4e2b-2ea5-414a-b3fd-01cd5630bfda
# ╠═f245c971-fa2d-482e-a355-67a51963c0cb
# ╠═6e303fd4-5539-4ec6-9bb8-72f366904f1f
# ╟─0f036ef3-4901-42e4-bf0c-0580215b20cb


text:
# include all my files, and setup the environment -
#include("Include.jl")

# main -
function main()

    # generate the problem dictionary -
    problem_dictionary = generate_problem_dictionary()

    # pass the problem dictionary to the GLPK wrapper -
    return calculate_optimal_flux_distribution(problem_dictionary)
end

# execute main -
results = main()


text:
# myfunctions_pids.jl

using Distributed

const PATH_DB = "./test3_invoicing.txt"
const PATH_DB_PAID = "./test3_invoicing_paid.txt"
const PATH_CSV = "./bank.csv"
const PATH_JOURNAL = "./test3_journal.txt"
const PATH_LEDGER = "./test3_ledger.txt"

# =================================
# task_0 - get orders from Sales
# =================================
function task_0(rx, pid)
    tx = Channel(32)
    @async while true
        if isready(tx)
            start = take!(tx)
            @info("task_0 (master): $(typeof(start)) - $start - $(start == "START")")
            if start == "START"
                @info("task_0 (master) will start the process")
                orders = @fetchfrom pid AppliSales.process()
                @info("task_0 (master) will put $(length(orders)) the orders on rx channel")
                put!(rx, orders)
                @info("task_0 (master) has put $(length(orders)) the orders on rx channel")
            end
        else
            @info("task_0 (master) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_0

# =================================
# task_1 - process the orders
# =================================
function task_1(rx, pid)
    tx = Channel(32)
    @async while true
        if isready(tx)
            orders = take!(tx)
            @info("task 1 (process the orders): $(typeof(orders))")
            if orders isa Array{AppliSales.Order, 1}
                @info("task_1 (process the orders) will process $(length(orders)) orders remotely")
                result = @fetchfrom pid AppliAR.process(orders; path=PATH_DB)
                @info("task_1 (process the orders) will put $(length(result)) journal entries on rx channel")
                put!(rx, result)
            end
        else
            @info("task_1 (process the orders) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_1

# =================================
# task_2 - process journal entries
# =================================
function task_2(rx, pid)
    tx = Channel(32)
    @async while true
        if isready(tx)
            entries = take!(tx)
            @info("task_2 (process journal entries): $(typeof(entries))")
            if entries isa Array{AppliGeneralLedger.JournalEntry,1}
                @info("task_2 (process journal entries) will process $(length(entries)) journal entries remotely")
                result = @fetchfrom pid AppliGeneralLedger.process(entries; path_journal=PATH_JOURNAL, path_ledger=PATH_LEDGER)
                @info("task_2 (process journal entries) saved the journal entries")
                #put!(tx, result)
            end
        else
            @info("task_2 (process journal entries) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task2_2

# =================================
# task_3 - process payments
# =================================
function task_3(rx, pid)
    tx = Channel(32)
    @async while true
        if isready(tx)
            stms = take!(tx)
            @info("task_3 (process payments): $(typeof(stms))")
            if stms isa Array{AppliAR.BankStatement,1}
                @info("task_3 (process payments) will match unpaid invoices with bank statements")
                result = @fetchfrom pid begin
                    unpaid_invoices = retrieve_unpaid_invoices(; path=PATH_DB)
                    try
                        AppliAR.process(unpaid_invoices, stms; path=PATH_DB_PAID)
                    catch e
                        @info(e)
                    end
                end
                @info("task_3 (process payments) will put $(length(result)) journal entries on rx channel")
                put!(rx, result)
            end
        else
            @info("task_3 (process payments) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_3

# =================================
# task dispatcher
# =================================
function dispatcher(p, q)
    rx = Channel(32)

    # instantiate tasks
    tx0 = task_0(rx, p) # get the orders
    tx1 = task_1(rx, q) # process the orders
    tx2 = task_2(rx, p) # process the journal entries
    tx3 = task_3(rx, q) # process the unpaid invoices

    # definition Holy traits pattern Dispatcher in domain.jl
    dispatch(x::T) where {T} = dispatch(Dispatcher(T), x)

    dispatch(::T0, x) = put!(tx0, x)
    dispatch(::T1, x) = put!(tx1, x)
    dispatch(::T2, x) = put!(tx2, x)
    dispatch(::T3, x) = put!(tx3, x)
    # dispatch(::T4, x) = put!(tx0, x)

    @async while true
        if isready(rx)
            value = take!(rx)
            @info("Dispatcher received $(typeof(value))")
            dispatch(value)
        else
            wait(rx)
        end
    end
    return rx
end # dispatcher


text:
# # Support vector machine
# ### Support Vector Machine (SVM)
# We are given two sets of points in ${\bf R}^n$, $\{x_1, \ldots, x_N\}$ and $\{y_1, \ldots, y_M\}$, and wish to find a function $f(x) = w^T x - b$ that linearly separates the points, i.e. $f(x_i) \geq 1$ for $i = 1, \ldots, N$ and $f(y_i) \leq -1$ for $i = 1, \ldots, M$. That is, the points are separated by two hyperplanes, $w^T x - b = 1$ and $w^T x - b = -1$.
#
# Perfect linear separation is not always possible, so we seek to minimize the amount that these inequalities are violated. The violation of point $x_i$ is $\text{max} \{1 + b - w^T x_i, 0\}$, and the violation of point $y_i$ is $\text{max} \{1 - b + w^T y_i, 0\}$. We tradeoff the error $\sum_{i=1}^N \text{max} \{1 + b - w^T x_i, 0\} + \sum_{i=1}^M \text{max} \{1 - b + w^T y_i, 0\}$ with the distance between the two hyperplanes, which we want to be large, via minimizing $\|w\|^2$.
#
# We can write this problem as
#
# $$
# \begin{array}{ll}
#     \text{minimize}   & \|w\|^2 + C * (\sum_{i=1}^N \text{max} \{1 + b - w^T x_i, 0\} + \sum_{i=1}^M \text{max} \{1 - b + w^T y_i, 0\})
# \end{array},
# $$
#
# where $w \in {\bf R}^n$ and $b \in {\bf R}$ are our optimization variables.
#
# We can solve the problem as follows.

using Convex, SCS

#-

## Generate data.
n = 2; # dimensionality of data
C = 10; # inverse regularization parameter in the objective
N = 10; # number of positive examples
M = 10; # number of negative examples

using Distributions: MvNormal
## positive data points
pos_data = rand(MvNormal([1.0, 2.0], 1.0), N);
## negative data points
neg_data = rand(MvNormal([-1.0, 2.0], 1.0), M);

#-

function svm(pos_data, neg_data, solver=() -> SCS.Optimizer(verbose=0))
    ## Create variables for the separating hyperplane w'*x = b.
    w = Variable(n)
    b = Variable()
    ## Form the objective.
    obj = sumsquares(w) + C*sum(max(1+b-w'*pos_data, 0)) + C*sum(max(1-b+w'*neg_data, 0))
    ## Form and solve problem.
    problem = minimize(obj)
    solve!(problem, solver)
    return evaluate(w), evaluate(b)
end;

#-

w, b = svm(pos_data, neg_data);

#-

## Plot our results.
using Plots
# Generate the separating hyperplane
line_x = -2:0.1:2;
line_y = (-w[1] * line_x .+ b)/w[2];
# Plot the positive points, negative points, and separating hyperplane.
plot(pos_data[1,:], pos_data[2,:], st=:scatter, label="Positive points")
plot!(neg_data[1,:], neg_data[2,:], st=:scatter, label="Negative points")
plot!(line_x, line_y, label="Separating hyperplane")


text:
const RGAS = 8.31446261815324

unicode_subscript(a::Int) = join('₀'+d for d in reverse(digits(a)))

function molar_to_weight(value,x,mw)
    return value/(mw*x)*1000.0
end

# mixing rule:
#= =
an operation on a vector of properties that returns a single sumber.
is a weighted operation on p, by the weights x.
#

= =#
"""
    mixing_rule(op, x, p)

returns an efficient implementation of:
`sum(x[i]*x[j]*op(p[i],p[j]) for i = 1:n , j = 1:n)`
where `op(p[i],p[j]) == op(p[j],p[i])`

"""
function mixing_rule(op, x, p,A=nothing)
    if isnothing(A)
        N = length(x)
        @boundscheck checkbounds(p, N)
        @inbounds begin
            res1 = zero(eltype(x))
            for i = 1:N
                res1 += p[i] * x[i]^2
                for j = 1:i - 1
                    res1 += 2 * x[i] * x[j] * op(p[i], p[j])
                end
            end
        end
        return res1
    else
        return mixing_rulea(op,x,p,A)
    end
end
#
# Abstract mixing rule, generates a mixing rule,based on
# an operation, so mij = xi*xj*op(pi,pj)*Aij
# example: mixing_rule(geometric_mean_rule,x,Tc,1.-K)
#

function mixing_rulea(op, x, p, A)
    N = length(x)
    checkbounds(A, N, N)
    @boundscheck checkbounds(p, N)
    @inbounds begin
        res1 = zero(eltype(x))
        for i = 1:N
            res1 += p[i] * x[i]^2
            for j = 1:i - 1
                res1 += 2 * x[i] * x[j] * op(p[i], p[j]) * A[i, j]
            end
        end
    end
    return res1
end
#
# Abstract asymetric mixing rule
# Adds a Asymetric matrix A_asym, and a op_sim(xi,xj,Aij)
# the mayor example is the GERG2008 equation, where
# op_asym(xi,xj,Aij) = (xi+xj)/(Aij^2*xi + xj)
#
"""
    mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)

returns an efficient implementation of:
` sum(A[i,j] * x[i] * x[j] * op(p[i],p[j]) * op_asym(x[i],x[j],A_asym[i,j])) for i = 1:n , j = 1:n)`
where `op(p[i],p[j]) == op(p[j],p[i])` , op_asym doesn't follow this symmetry.

"""
function mixing_rule_asymetric(op, op_asym, x, p, A, A_asym)
    N = length(x)
    checkbounds(A, N, N)
    checkbounds(A_asym, N, N)
    @boundscheck checkbounds(p, N)
    @inbounds begin
        res1 = zero(eltype(x))

        for i = 1:N
            x[i] != 0 && begin
                res1 += p[i] * x[i]^2
                for j = 1:i - 1
                    res1 +=
                        2 *
                        x[i] *
                        x[j] *
                        op(p[i], p[j]) *
                        A[i, j] *
                        op_asym(x[i], x[j], A_asym[i, j])
                end
            end
        end
    end
    return res1
end

"""
    mixing_matrix!(A, op, p)

returns a matrix of size nxn , where A[i,j] = op(p[i],p[j]), modifying A inplace

"""
function mixing_matrix!(A, op, p)
    N = length(p)
    @boundscheck size(A) == (N, N)
    @inbounds begin
        res1 = zero(eltype(p))
        for i = 1:N
            A[i, i] = p[i]
            for j = 1:i - 1
                A[i, j] = op(p[i], p[j])
                A[j, i] = op(p[i], p[j])
            end
        end
    end
    return A

end

"""
    mixing_matrix!(op, p)

returns a matrix of size nxn , where A[i,j] = op(p[i],p[j])
"""
function mixing_matrix(op, p)
    N = length(p)
    A = Matrix{eltype(p)}(undef, N, N)
    return mixing_matrix!(A, op, p)
end

#a_in_b(t1, t2) = all(in(t2), t1)
#tuple_comparison(t1, t2) = ((length(t1) == length(t2)) && a_in_b(t1, t2) && a_in_b(t2, t1))

#returns f(x)/(∂f(x)/∂x), for Newton solvers, using ForwardDiff+StaticArrays+DiffResults for one single AD pass

"""
    autonewton(f,x)

returns f/(df/fx) evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate f and dfdx in one pass
"""
function autonewton(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.GradientResult(x_vec)
    _∂f =  ForwardDiff.gradient!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    return fx/∂f∂x
end

function f_dfdx(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.GradientResult(x_vec)
    _∂f =  ForwardDiff.gradient!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    return (fx,∂f∂x)
end
"""
    normalizefrac!!(x)

mutates x, so that sum(x) == 1. also removes negative values and replaces Inf with big numbers

"""
function normalizefrac!!(x)
    x0 = zero(eltype(x))
    function f(xi)
        if xi < x0
         return -x0
        elseif isinf(xi)
            return inv(eps(typeof(xi)))
        else
            return xi
        end
    end
    xn = one(x0)/sum(f,x)
    @! x .= f.(x) .* xn
    return x
end

function nan_to_num(xi)
    if x == typemax(xi)
        return prevfloat(xi)
    elseif isnan(xi)
        return zero(xi)
    else
        return xi
    end
end

function remove_minimums!!(x,lo=zero(eltype(x)))
    return @! x .= max.(x,lo)
end

function remove_maximums!!(x,hi=inv(eps(eltype(x))))
    return @! x .= min.(x,hi)
end

function cardan_reals(res)
    nan = real(zero(first(res))/zero(first(res)))
    a1, a2, a3 = res
    isreal1 = abs(imag(a1) < 4*eps(imag(a1)))
    isreal2 = abs(imag(a2) < 4*eps(imag(a2)))
    isreal3 = abs(imag(a3) < 4*eps(imag(a3)))
    r1 = real(a1)
    r2 = real(a2)
    r3 = real(a3)
    if isreal1 & isreal2 & isreal3
        return (r1,r2,r3)
    elseif !isreal1 & !isreal2
        return (r3,r3,r3)
    elseif  !isreal3 & !isreal2
        return (r1,r1,r1)
    elseif  !isreal3 & !isreal1
        return (r2,r2,r2)
    elseif !isreal1
        return (r3,r2,r2)
    elseif !isreal2
        return (r1,r3,r3)
    elseif !isreal3
        return (r1,r3,r3)
    else
        return (nan,nan,nan)
    end
end

function is_liquid(x::Symbol)::Bool
    x in (:liquid,:l,:L,:LIQUID)
end

function is_gas(x::Symbol)::Bool
    x in (:gas,:g,:G,:GAS,:vapor,:VAPOR,:V,:v)
end


text:
module decoder

using conv_code
using utils
using transitions
using code_structs
using cost_functions

export decode

function decode(p::Float64, code::ConvolutionalCode, in_sequence::Array{Array{Float64, 1}, 1})

	# O nível da raiz é 1
	current_leaves::Array{TreeNode} = [
		TreeNode(nothing, 0, 0, -1) for i in range(1; length=2^code.quant_mem)
	]
	current_leaves[1].level = 1

	next_leaves::Array{TreeNode} = [
		TreeNode(nothing, 0, 0, -1) for i in range(1; length=2^code.quant_mem)
	]

	# No índice level, calcular as folhas do nível level+1
	for level=1:length(in_sequence)
		# Sequência recebida
		curr_seq = in_sequence[level]

		# Calcular os filhos
		for curr_leaf_index=1:length(current_leaves)
			curr_leaf_node = current_leaves[curr_leaf_index]

			if curr_leaf_node.level < level
				continue
			end

			for transition=transition_list(code, curr_leaf_index)
				child = TreeNode(
					curr_leaf_node,
					level+1,
					euclidean_cost(p, curr_seq, curr_leaf_node, transition),
					transition.in_bit
				)
				# Se visto pela primeira vez, ou mais barato que o já
				# colocado na árvore, inserir `child` na árvore
				if next_leaves[transition.next_state].level < level+1 ||
					next_leaves[transition.next_state].cost > child.cost

					next_leaves[transition.next_state] = child

				end
			end
		end
		@swap!(current_leaves, next_leaves)
	end

	# Selecionar a folha mais barata do último nível
	best_leaf = current_leaves[1]
	for leaf=current_leaves
		if leaf.level != length(in_sequence)+1
			continue
		end
		if leaf.cost < best_leaf.cost
			best_leaf = leaf
		end
	end

	# Subir na árvore até a raiz
	node = best_leaf
	decoded_bits = Array{Int64}(undef, length(in_sequence))
	for i=length(in_sequence):-1:1
		decoded_bits[i] = node.info_bit
		node = node.parent
	end

	return decoded_bits
end

end


text:
# Automatically generated using Clang.jl

# const CSFML_SYSTEM_API = CSFML_API_IMPORT
const sfInputStreamReadFunc = Ptr{Cvoid}
const sfInputStreamSeekFunc = Ptr{Cvoid}
const sfInputStreamTellFunc = Ptr{Cvoid}
const sfInputStreamGetSizeFunc = Ptr{Cvoid}

struct sfInputStream
    read::sfInputStreamReadFunc
    seek::sfInputStreamSeekFunc
    tell::sfInputStreamTellFunc
    getSize::sfInputStreamGetSizeFunc
    userData::Ptr{Cvoid}
end

struct sfTime
    microseconds::sfInt64
end

const sfClock = Cvoid
const sfMutex = Cvoid
const sfThread = Cvoid

struct sfVector2i
    x::Cint
    y::Cint
end

struct sfVector2u
    x::UInt32
    y::UInt32
end

struct sfVector2f
    x::Cfloat
    y::Cfloat
end

struct sfVector3f
    x::Cfloat
    y::Cfloat
    z::Cfloat
end