text:
module FaultDetectionTools

using LinearAlgebra
using MatrixEquations
using MatrixPencils
using DescriptorSystems
using Polynomials
using Random

abstract type AbstractFDDObject end

import LinearAlgebra: BlasFloat, BlasReal, BlasComplex
import DescriptorSystems: DescriptorStateSpace, chess, rcond, gbalmr, gminreal, gpole
import Combinatorics: combinations
import Base: +, -, *

export FDIModel, fdimodset, FDFilter, FDFilterIF, FDIFilter, FDIFilterIF, fdIFeval, gbalmr
export fditspec, fdisspec, fdiscond, fdif2ngap, fdimmperf
export fdigenspec, fdichkspec
export efdsyn, efdisyn, efdbasesel, afdsyn, afdredsyn, afdisyn, afdbasesel, emmsyn, emmbasesel, ammsyn, ammbasesel
export fdhinfminus, fdhinfmax, binmat2dec, dec2binmat, fditspec_, fdisspec_, fdiscond_

export MDMModel, MDModel, mdmodset, MDFilter, MDFilterIF, mdIFeval
export mddist, mddist2c
export mdsspec, mdperf, mdmatch
export emdsyn, emdbasesel

const VRS = Union{Vector{Int}, UnitRange{Int}, Int}

include("types/FDIObjects.jl")
include("FDIanalysis.jl")
include("FDIperformance.jl")
include("FDIsynthesis.jl")
include("FDIutils.jl")
include("types/MDObjects.jl")
include("MDanalysis.jl")
include("MDperformance.jl")
include("MDsynthesis.jl")
end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
const PLANETOID_URL = "https://github.com/kimiyoung/planetoid/raw/master/data"
const PLANETOID_DATASETS = [:citeseer, :cora, :pubmed]
const EXTS = ["allx", "ally", "graph", "test.index", "tx", "ty", "x", "y"]
const DATAURLS = [joinpath(PLANETOID_URL, "ind.$(d).$(ext)") for d in PLANETOID_DATASETS, ext in EXTS]

planetoid_init() = register(DataDep(
    "Planetoid",
    """
    The citation network datasets "Cora", "CiteSeer", "PubMed" from
    "Revisiting Semi-Supervised Learning with Graph Embeddings"
    <https://arxiv.org/abs/1603.08861> paper.
    Nodes represent documents and edges represent citation links.
    """,
    reshape(DATAURLS, :),
    "f52b3d47f5993912d7509b51e8090b6807228c4ba8c7d906f946868005c61c18";
    post_fetch_method=preprocess_planetoid,
))

function preprocess_planetoid(local_path)
    for dataset in PLANETOID_DATASETS
        graph_file = @datadep_str "Planetoid/ind.$(dataset).graph"
        trainX_file = @datadep_str "Planetoid/ind.$(dataset).x"
        trainy_file = @datadep_str "Planetoid/ind.$(dataset).y"
        testX_file = @datadep_str "Planetoid/ind.$(dataset).tx"
        testy_file = @datadep_str "Planetoid/ind.$(dataset).ty"
        allX_file = @datadep_str "Planetoid/ind.$(dataset).allx"
        ally_file = @datadep_str "Planetoid/ind.$(dataset).ally"

        train_X = read_data(trainX_file)
        train_y = read_data(trainy_file)
        test_X = read_data(testX_file)
        test_y = read_data(testy_file)
        all_X = read_data(allX_file)
        all_y = read_data(ally_file)
        graph = read_graph(graph_file)

        num_V = length(graph)
        sg = to_simplegraph(graph, num_V)
        num_E = ne(sg)
        feat_dim = size(all_X, 2)
        label_dim = size(all_y, 2)
        train_X, train_y = sparse(train_X'), sparse(train_y')
        test_X, test_y = sparse(test_X'), sparse(test_y')
        all_X, all_y = sparse(all_X'), sparse(all_y')
        raw = Dict(:graph=>graph, :train_X=>train_X, :train_y=>train_y,
                   :test_X=>test_X, :test_y=>test_y, :all_X=>all_X, :all_y=>all_y)
        meta = (graph=(num_V=num_V, num_E=num_E),
                train=(features_dim=(feat_dim, size(train_X, 2)), labels_dim=(label_dim, size(train_y, 2))),
                test=(features_dim=(feat_dim, size(test_X, 2)), labels_dim=(label_dim, size(test_y, 2))),
                all=(features_dim=(feat_dim, size(all_X, 2)), labels_dim=(label_dim, size(all_y, 2)))
                )

        graphfile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).graph.jld2")
        trainfile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).train.jld2")
        testfile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).test.jld2")
        allfile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).all.jld2")
        rawfile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).raw.jld2")
        metadatafile = replace(graph_file, "ind.$(dataset).graph"=>"$(dataset).metadata.jld2")

        @save graphfile sg
        @save trainfile train_X train_y
        @save testfile test_X test_y
        @save allfile all_X all_y
        @save rawfile raw
        @save metadatafile meta
    end
end

function read_data(filename)
    py"""
    import pickle
    from scipy.sparse import csr_matrix

    with open($filename,"rb") as f:
        u = pickle._Unpickler(f)
        u.encoding = "latin1"
        data = u.load()

    if type(data) is csr_matrix:
        data = data.toarray()
    """
    return SparseMatrixCSC(Array(py"data"))
end

read_index(filename) = map(x -> parse(Int64, x), readlines(filename))

function read_graph(filename)
    py"""
    import pickle

    with open($filename,"rb") as f:
        u = pickle._Unpickler(f)
        u.encoding = "latin1"
        data = u.load()
    """
    return Dict(py"data")
end

struct Planetoid <: Dataset
    dataset::Symbol

    function Planetoid(ds::Symbol)
        @assert ds in PLANETOID_DATASETS "`dataset` should be one of citeseer, cora, pubmed."
        new(ds)
    end
end

function graphdata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).graph.jld2"
    @load file sg
    sg
end

function traindata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).train.jld2"
    @load file train_X train_y
    train_X, train_y
end

function testdata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).test.jld2"
    @load file test_X test_y
    test_X, test_y
end

function alldata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).all.jld2"
    @load file all_X all_y
    all_X, all_y
end

function rawdata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).raw.jld2"
    @load file raw
    raw
end

function metadata(pla::Planetoid)
    file = @datadep_str "Planetoid/$(pla.dataset).metadata.jld2"
    @load file meta
    meta
end

ext: jl


text:
module MLJTutorial

import IJulia
import Pluto

const DATASETS = [:horse, :house, :small]
const ROOT = joinpath(@__DIR__, "..")
const NOTEBOOKS = joinpath(ROOT, "notebooks")
const FIRST_PLUTO_NOTEBOOK =
    joinpath(NOTEBOOKS, "01_data_representation", "notebook.pluto.jl")

export jupyter, pluto, go, remove, load_house, load_horse, load_small

jupyter() = begin
    IJulia.notebook(dir=NOTEBOOKS)
end
const go = jupyter

pluto() = begin
    Pluto.run(notebook=FIRST_PLUTO_NOTEBOOK)
end

end # module


text:
@info "Loading test packages..."

using BioMotifInference, BioBackgroundModels, BioSequences, Distributions, Distributed, Random, Serialization, Test
import StatsFuns: logsumexp, logaddexp
import BioMotifInference:estimate_dirichlet_prior_on_wm, assemble_source_priors, init_logPWM_sources, wm_shift, permute_source_weights, get_length_params, permute_source_length, get_pwm_info, get_erosion_idxs, erode_source, init_mix_matrix, mixvec_decorrelate, mix_matrix_decorrelate, most_dissimilar, most_similar, revcomp_pwm, score_sources_ds!, score_sources_ss!, weave_scores_ss!, weave_scores_ds!, IPM_likelihood, consolidate_check, consolidate_srcs, pwm_distance, permute_source, permute_mix, perm_src_fit_mix, fit_mix, random_decorrelate, reinit_src, erode_model, reinit_src, shuffle_sources, accumulate_mix, distance_merge, similarity_merge, info_fill, converge_ensemble!, reset_ensemble!, Permute_Tuner, PRIOR_WT, TUNING_MEMORY, CONVERGENCE_MEMORY, tune_weights!, update_weights!, clamp_pvec!
import Distances: euclidean

@info "Beginning tests..."
using Random
Random.seed!(786)
O=1000;S=50

include("pwm_unit_tests.jl")
include("mix_matrix_unit_tests.jl")
include("likelihood_unit_tests.jl")
include("consolidate_unit_tests.jl")
include("permute_func_tests.jl")
include("permute_tuner_tests.jl")
include("ensemble_tests.jl")
@info "Tests complete!"


text:
devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticTransformer),
    :TT => DeviceModel(TapTransformer, StaticTransformer),
)
services = Dict{Symbol, ServiceModel}()

@testset "Operation Model kwargs with CopperPlatePowerModel base" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    c_sys5_re = build_system("c_sys5_re")
    c_sys14 = build_system("c_sys14")

    @test_throws ArgumentError OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        bad_kwarg = 10,
    )
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    moi_tests(op_problem, true, 120, 0, 120, 120, 24, false)
    op_problem =
        OperationsProblem(TestOpProblem, template, c_sys14; optimizer = OSQP_optimizer)
    moi_tests(op_problem, false, 120, 0, 120, 120, 24, false)
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        use_forecast_data = false,
        optimizer = GLPK_optimizer,
    )
    moi_tests(op_problem, false, 5, 0, 5, 5, 1, false)
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        use_forecast_data = false,
        use_parameters = false,
        optimizer = GLPK_optimizer,
    )
    moi_tests(op_problem, false, 5, 0, 5, 5, 1, false)

    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        optimizer = GLPK_optimizer,
        balance_slack_variables = true,
    )
    moi_tests(op_problem, false, 168, 0, 120, 120, 24, false)

end

@testset "Test optimization debugging functions" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    MOIU.attach_optimizer(op_problem.psi_container.JuMPmodel)
    constraint_indices = get_all_constraint_index(op_problem)
    for (key, index, moi_index) in constraint_indices
        val1 = get_con_index(op_problem, moi_index)
        val2 = op_problem.psi_container.constraints[key].data[index]
        @test val1 == val2
    end
    @test isnothing(get_con_index(op_problem, length(constraint_indices) + 1))

    var_indices = get_all_var_index(op_problem)
    for (key, index, moi_index) in var_indices
        val1 = get_var_index(op_problem, moi_index)
        val2 = op_problem.psi_container.variables[key].data[index]
        @test val1 == val2
    end
    @test isnothing(get_var_index(op_problem, length(var_indices) + 1))
end

@testset "Test passing custom JuMP model" begin
    my_model = JuMP.Model()
    my_model.ext[:PSI_Testing] = 1
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5,
        my_model;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    @test haskey(op_problem.psi_container.JuMPmodel.ext, :PSI_Testing)
    @test (:params in keys(op_problem.psi_container.JuMPmodel.ext)) == true
end

@testset "Operation Model Constructors with Parameters" begin

    networks = [
        CopperPlatePowerModel,
        StandardPTDFModel,
        DCPPowerModel,
        NFAPowerModel,
        ACPPowerModel,
        ACRPowerModel,
        ACTPowerModel,
        DCPLLPowerModel,
        LPACCPowerModel,
        SOCWRPowerModel,
        QCRMPowerModel,
        QCLSPowerModel,
    ]

    thermal_gens = [
        ThermalStandardUnitCommitment,
        ThermalDispatch,
        ThermalRampLimited,
        ThermalDispatchNoMin,
    ]

    c_sys5 = build_system("c_sys5")
    c_sys5_re = build_system("c_sys5_re")
    c_sys5_bat = build_system("c_sys5_bat")
    c_sys5_pwl_ed = build_system("c_sys5_pwl_ed")
    systems = [c_sys5, c_sys5_re, c_sys5_bat, c_sys5_pwl_ed]
    for net in networks, thermal in thermal_gens, system in systems, p in [true, false]
        @testset "Operation Model $(net) - $(thermal) - $(system)" begin
            devices = Dict{Symbol, DeviceModel}(
                :Generators => DeviceModel(ThermalStandard, thermal),
                :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
            )
            branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(Line, StaticLine))
            template = OperationsProblemTemplate(net, devices, branches, services)
            op_problem = OperationsProblem(
                TestOpProblem,
                template,
                system;
                use_parameters = p,
                PTDF = build_PTDF5(),
                export_pwl_vars = true,
            )
            @test :nodal_balance_active in keys(op_problem.psi_container.expressions)
            @test (:params in keys(op_problem.psi_container.JuMPmodel.ext)) == p
        end
    end

    @testset "Operation Model CopperPlatePowerModel - ThermalDispatchNoMin - c_sys5_pwl_ed_nonconvex" begin
        c_sys5_pwl_ed_nonconvex = build_system("c_sys5_pwl_ed_nonconvex")
        devices = Dict{Symbol, DeviceModel}(
            :Generators => DeviceModel(ThermalStandard, ThermalDispatchNoMin),
            :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
        )
        branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(Line, StaticLine))
        template =
            OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
        @test_throws IS.InvalidValue OperationsProblem(
            TestOpProblem,
            template,
            c_sys5_pwl_ed_nonconvex;
            use_parameters = true,
            PTDF = build_PTDF5(),
            export_pwl_vars = true,
        )
    end

    @testset "Operations template constructors" begin
        c_sys5 = build_system("c_sys5")
        op_problem_ed = PSI.EconomicDispatchProblem(c_sys5)
        op_problem_uc = PSI.UnitCommitmentProblem(c_sys5)
        moi_tests(op_problem_uc, false, 480, 0, 240, 120, 144, true)
        moi_tests(op_problem_ed, false, 120, 0, 168, 120, 24, false)
        ED = PSI.run_economic_dispatch(c_sys5; optimizer = fast_lp_optimizer)
        UC = PSI.run_unit_commitment(c_sys5; optimizer = fast_lp_optimizer)
        @test ED.optimizer_log[:primal_status] == MOI.FEASIBLE_POINT
        @test UC.optimizer_log[:primal_status] == MOI.FEASIBLE_POINT
    end
end


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Eighteen, Snailfish."

#= This day really took me some time, I had to rewrite the code twice, because
 = I was inattentive to one of the conditions of the problem (the order of the
 = operations).  Alas, Julia’s typesystem didn’t helped me that much either.
 =#

include("Utils.jl")

using .Utils

"Calculate the magnitude of the `snailnum`."
function calc_magn(snailnum :: Vector{Complex{Int}})
    magn = 0
    depth = 0
    factors = Array{Int}([])

    for lex in snailnum[1:end - 1]
        if isreal(lex)
            magn += Int(lex)*prod(factors[1:depth])
            factors[depth] = if factors[depth] == 3 2 else 3 end
        end

        depth += imag(lex)

        if depth > length(factors)
            push!(factors, 3)
        elseif imag(lex) == -1
            factors[depth] = if factors[depth] == 3 2 else 3 end
        end
    end

    return magn
end

"""Explode a pair in `snailnum`, when applicable, return whether the vector was
mutated."""
function explode!(snailnum :: Vector{Complex{Int}})
    depth = 0

    i = 0
    while (i += 1) <= length(snailnum)
        lex = snailnum[i]

        if lex == im && depth >= 4
            (a, b) = snailnum[i + 1], snailnum[i + 2]

            j = i - 1
            while j >= 1
                if isreal(snailnum[j])
                    snailnum[j] += a
                    break
                end

                j -= 1
            end

            j = i + 4
            while j <= length(snailnum)
                if isreal(snailnum[j])
                    snailnum[j] += b
                    break
                end

                j += 1
            end

            splice!(snailnum, i:i + 3, 0)

            return true
        end

        depth += imag(lex)
    end

    return false
end

"""Split a number in `snailnum`, when applicable, return whether the vector was
mutated."""
function split!(snailnum :: Vector{Complex{Int}})
    i = 0
    while (i += 1) <= length(snailnum)
        lex = snailnum[i]

        if isreal(lex) && Int(lex) >= 10
            (m, r) = divrem(Int(lex), 2)
            splice!(snailnum, i, [im, m, m + r, -im])

            return true
        end
    end

    return false
end

"Reduce the `snailnum` (mutating!), return whether the vector was changed."
reduce!(snailnum :: Vector{Complex{Int}}) = explode!(snailnum) || split!(snailnum)

"Sum and reduce all `snailnums`, in order, and return the final magnitude."
function solve1(snailnums :: Vector{Vector{Complex{Int}}})
    res = first(snailnums)

    for addend in snailnums[2:end]
        res = vcat(im, res, addend, -im)

        while reduce!(res) end
    end

    return calc_magn(res)
end

"Find the maximum magnitude of any sum of two distinct `snailnums` in the set."
function solve2(snailnums :: Vector{Vector{Complex{Int}}})
    max_magn = 0

    for i in eachindex(snailnums)
        for j in eachindex(snailnums)
            if i == j
                continue
            end

            addition = solve1([snailnums[i], snailnums[j]])

            max_magn = max(max_magn, addition)
        end
    end

    return max_magn
end

"Parse `snailnum` string into a numerical format."
function parsenum(snailnum :: String)
    parsed = Array{Complex{Int}}([])

    i = 0
    while (i += 1) <= length(snailnum)
        lex = snailnum[i]

        #= Use the imaginary unit for the pair opening, and its additive
         = inverse for the closing.
         =#
        if lex == '['
            push!(parsed, im)
        elseif lex == ']'
            push!(parsed, -im)
        elseif isnumeric(lex)
            numstr = Array{Char}([lex])

            while isnumeric(snailnum[i + 1])
                push!(numstr, snailnum[i += 1])
            end

            push!(parsed, array2num(parse.(Int, numstr)))
        end
    end

    return parsed
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))
    real_data = parsenum.(data)

    println(solve1(real_data))
    println(solve2(real_data))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
using FrankWolfe
using LinearAlgebra
using Test

n = Int(1e3)
k = n

f(x) = dot(x, x)
function grad!(storage, x)
    @. storage = 2 * x
end

# pick feasible region
lmo = FrankWolfe.ProbabilitySimplexOracle{Rational{BigInt}}(1); # radius needs to be integer or rational

# compute some initial vertex
x0 = FrankWolfe.compute_extreme_point(lmo, zeros(n));

x, v, primal, dual_gap0, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    print_iter=k / 10,
    verbose=false,
    emphasis=FrankWolfe.blas,
)

xmem, vmem, primal, dual_gap, trajectory = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    copy(x0),
    max_iteration=k,
    line_search=FrankWolfe.Agnostic(),
    print_iter=k / 10,
    verbose=true,
    emphasis=FrankWolfe.memory,
)

@time xstep, _ = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.RationalShortstep(),
    L=2,
    print_iter=k / 10,
    verbose=true,
)

@test eltype(x) == eltype(xmem) == eltype(xstep) == Rational{BigInt}
@test xmem == x
@test abs(f(xstep) - f(x)) <= 1e-3
@test vmem == v


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SHTns"
version = v"3.5.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://bitbucket.org/nschaeff/shtns/downloads/shtns-$(version).tar.gz", "77f8a33b94df8786d2ce9b95cbfbe548f00443625b310b8c64012b22c8a7394f")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/shtns/
export CFLAGS="-fPIC"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
if [[ "${target}" == *-apple-* ]]; then
    # For some reasons, `configure` insists on setting CC2=gcc, also on macOS
    sed -i -e 's/gcc/cc/' -e 's/ -fno-tree-loop-distribute-patterns//' Makefile
fi
make -j${nproc}
make install
mkdir -p ${libdir}
cc -shared -o "${libdir}/libshtns.${dlext}" *.o -lfftw3
rm "${prefix}/lib/libshtns.a"
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
struct TSLSModel <: LinearModel
    y::Vector # outcome variable
    X::Matrix # Second stage variables (endog + exog)
    Z::Matrix # First stage variables (instrument + exog)
    Q_Z::Matrix
    R_Z::Matrix
    vcov::vcov
    X_names::Vector
    N::Int64
    K::Int64
    K_endog::Int64
    K_exog::Int64
    K_instrument::Int64
    function TSLSModel(y::Vector,
                       X_endog::Matrix,
                       X_exog::Matrix,
                       instruments::Matrix;
                       vcov::vcov = vcovIID(),
                       X_names = nothing)
        N = length(y)
        K_endog = size(X_endog, 2)
        K_exog = size(X_exog, 2)
        K_instrument = size(instruments, 2)
        # If only one column of 0s we drop the exog vars
        if K_exog == 1 & all(X_exog[:, 1] .== 0.0)
            Z = instruments
            X = X_endog
        else
            Z = hcat(instruments, X_exog)
            X = hcat(X_endog, X_exog)
        end
        Q_Z, R_Z = qr(Z)
        K = size(X, 2)
        if isnothing(X_names)
            X_names = vcat(
                "x_endog_" .* string.(1:size(X_endog, 2)),
                "x_exog_" .* string.(1:size(X_exog, 2))
            )
        end
        new(y, X, Z, Q_Z, R_Z, vcov, X_names, N, K, K_endog, K_exog, K_instrument)
    end
end

# Adding DataFrames functionality
function TSLSModel(y::Vector,
                    X_endog::DataFrame,
                    X_exog::DataFrame,
                    instruments::DataFrame;
                    vcov::vcov = vcovIID())
    X_names = vcat(names(X_endog), names(X_exog))
    X_endog = Matrix(X_endog)
    X_exog = Matrix(X_exog)
    instruments = Matrix(instruments)
    return TSLSModel(y, X_endog, X_exog, instruments, vcov = vcov, X_names = X_names)
end
struct FittedTSLSModel <: LinearModelFit
    y::Vector # outcome variable
    X::Matrix # Second stage variables (endog + exog)
    Z::Matrix # First stage variables (instrument + exog)
    Q::Matrix
    R::Matrix
    vcov::vcov
    X_names::Vector
    N::Int64
    K::Int64
    K_endog::Int64
    K_exog::Int64
    K_instrument::Int64
    modelfit::FitOutput
    P_Z::Matrix
end

struct FittedJIVEModel <:LinearModelFit
    y::Vector # outcome variable
    X::Matrix # Second stage variables (endog + exog)
    Z::Matrix # First stage variables (instrument + exog)
    Q::Matrix
    R::Matrix
    vcov::vcov
    X_names::Vector
    N::Int64
    K::Int64
    K_endog::Int64
    K_exog::Int64
    K_instrument::Int64
    modelfit::FitOutput
    P_Z::Matrix
    X_jive::Matrix
end

function fit!(model::TSLSModel)
    y = model.y
    X = model.X
    Z = model.Z
    R = model.R_Z
    ZZ_inv = inv(cholesky(R' * R))
    P_Z = Z * ZZ_inv * Z'
    β = inv(X' * P_Z * X) * X' * P_Z * y
    resid = y - X*β
    se_β, pval, σ_sq, vcov_matrix = inference(resid,
                                              β,
                                              P_Z,
                                              model, model.vcov)
    return FitOutput(β, se_β, pval, resid, σ_sq, vcov_matrix), P_Z
end

function fit(model::TSLSModel)
    FittedTSLSModel(
        model.y,
        model.X,
        model.Z,
        model.Q_Z,
        model.R_Z,
        model.vcov,
        model.X_names,
        model.N,
        model.K,
        model.K_endog,
        model.K_exog,
        model.K_instrument,
        fit!(model)...
    )
end

function AndersonRubinCI!(β_AR,
                          fit_obj::Jeeves.FittedTSLSModel;
                          additional_controls = nothing,
                          signif_level = 0.05)
    K_instrument = fit_obj.K_instrument
    K_endog = fit_obj.K_endog
    X_endog = fit_obj.X[:, 1:K_endog]
    X_exog_only = fit_obj.Z[:, (K_instrument + 1):end]
    X_exog = fit_obj.Z
    AR_resid = fit_obj.y - X_endog*β_AR
    if !isnothing(additional_controls)
        X_exog = hcat(fit_obj.Z, additional_controls)
        X_exog_only = hcat(X_exog_only, additional_controls)
    end
    AR_model = OLSModel(AR_resid, X_exog, vcov = fit_obj.vcov)
    restricted_AR_model = OLSModel(AR_resid, X_exog_only, vcov = fit_obj.vcov)
    AR_fit = fit(AR_model)
    restricted_AR_fit = fit(restricted_AR_model)
    AR_F_vals = Ftest(AR_fit, restricted_AR_fit, signif_level = signif_level)
    return [β_AR..., AR_F_vals...]
end

function AndersonRubinCI(β_grid, fit_obj::FittedTSLSModel; additional_controls = nothing)
    β_grid = collect(β_grid)
    K_endog = fit_obj.K_endog
    K_β =  length(β_grid[1])
    K_β == K_endog || error("Grid dimension doesn't match K_endog")
    AR_results = Vector(undef, length(β_grid))
    for i in 1:length(β_grid)
        AR_results[i] = AndersonRubinCI!(
            [β_grid[i]...],
            fit_obj,
            additional_controls = additional_controls)
    end
    return AR_results
end

function jive!(model::TSLSModel)
    y = model.y
    X = model.X
    Z = model.Z
    R = model.R_Z
    ZZ_inv = inv(cholesky(R' * R))
    P_Z = Z * ZZ_inv * Z'
    h_i = diag(P_Z)
    π = ZZ_inv * (Z' * X)
    X_jive = (Z * π - h_i .* X)./(1 .- h_i)
    XX_inv = inv(X_jive' * X)
    β = XX_inv * (X_jive' * y)
    resid = y - X*β

    se_β, pval, σ_sq, vcov_matrix = inference(model.N,
                                              model.K,
                                              resid,
                                              β,
                                              XX_inv,
                                              X_jive,
                                              model.vcov)
    return FitOutput(β, se_β, pval, resid, σ_sq, vcov_matrix), P_Z, X_jive
end

function jive(model::TSLSModel)
    FittedJIVEModel(
        model.y,
        model.X,
        model.Z,
        model.Q_Z,
        model.R_Z,
        model.vcov,
        model.X_names,
        model.N,
        model.K,
        model.K_endog,
        model.K_exog,
        model.K_instrument,
        jive!(model)...
    )
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.
module Data

import Unitful
using Unitful.DefaultSymbols:mm,°

include("HumanEye.jl")

end #module


text:
module JointPosteriors

using   LineSearches: BackTracking
using   Optim
using   ForwardDiff
using   SparseQuadratureGrids
using   ConstrainedParameters
using   LogDensities
using   Interpolations
using   Distributions

import  Base: show, quantile, Val
import  LogDensities: log_density
import  SparseQuadratureGrids: SlidingVecFun

export  JointPosterior,
        JointPosteriorRaw,
        fit,
        marginal,
        Data,
        parameter,
        CovarianceMatrix,
        PositiveVector,
        ProbabilityVector,
        RealVector,
        Simplex,
        Model,
        construct,
        log_density,
        log_density!,
        log_jacobian!,
        quad_form,
        inv_det,
        inv_root_det,
        root_det,
        log_root_det,
        trace_inverse,
        lpdf_InverseWishart,
        lpdf_normal,
        Normal,
        Adaptive,
        AdaptiveRaw,
        Smolyak,
        pdf,
        cdf,
        quantile

include("joint_posterior.jl")
include("interp.jl")
include("marginal_posterior.jl")

end # module


text:
@test_throws ArgumentError OrbitalIndex(-1, -2)
@test_throws ArgumentError OrbitalIndex(3, -1)
@test_throws ArgumentError OrbitalIndex(3, 5)

@testset ">> Operators" begin
    @test OrbitalIndex(0, 0) == OrbitalIndex(0, 0)
    @test OrbitalIndex(2, 1) == OrbitalIndex(2, 1)
    @test OrbitalIndex(2, 1) != OrbitalIndex(2, 0)
    @test OrbitalIndex(2, 1) ≠ OrbitalIndex(2, 0)

    @test OrbitalIndex(2, 1) > OrbitalIndex(2, 0)
    @test OrbitalIndex(2, 1) > OrbitalIndex(1, 0)
    @test OrbitalIndex(2, 1) < OrbitalIndex(3, 0)
end

@testset ">> Linear indexing for first numbers" begin
    index = 1
    for n in 0:100
        for l in 0:n
            orb = OrbitalIndex(n, l)
            @test convert(Integer, orb) == index
            @test convert(OrbitalIndex, convert(Integer, orb)) == orb
            index += 1
        end
    end
end

@testset ">> conversions for n=$n" for n in rand(0:1000000, 5)
    @testset ">> conversions for l=$l" for l in rand(0:n, 5)
        orb = OrbitalIndex(n, l)
        back  = convert(OrbitalIndex, convert(Integer, orb))
        @test typeof(back) === typeof(orb)
        @test back == orb
    end
end

@testset ">> String indexing" begin
    @test nl"2, 2" == OrbitalIndex(2, 2)
	@test nl"3, 2" == OrbitalIndex(3, 2)
end

@testset ">> Printing" begin
	buffer = IOBuffer()
	print(buffer, nl"12, 4")
	seekstart(buffer)
	@test readstring(buffer) == "n=12, l=4"
end

@testset ">> looping by 1" begin
    index = 3
    for u in nl"1,1":nl"5, 2"
        @test convert(Integer, u) == index
        index += 1
    end
    @test convert(Integer, nl"5, 3") == index
end

@testset ">> looping by (n, 0)" begin
    n = 1
    index = 0
    for u in nl"1,1":nl"2, 0":nl"5, 2"
        @test u.l == 1
        @test u.n == n
        index = convert(Integer, u)
        n += 2
    end
    @test index == convert(Integer, nl"5, 1")
end

@testset ">> looping by (n, l)" begin
    expected = [OrbitalIndex(1, 1), OrbitalIndex(4, 3), OrbitalIndex(7, 5)]
    actual = [u for u in nl"1,1":nl"3, 2":nl"9, 2"]
    @test expected == actual
end


text:
# Original file: https://github.com/JuliaGraphics/Luxor.jl/blob/master/src/play.jl
# Replaced whileloop with forloop in play macro to limit number of displayed frames to n

using MiniFB

function onclick(window, button, mod, isPressed)::Cvoid
    if Bool(isPressed)
        println("mouse clicked")
    end
    mousex = mfb_get_mouse_x(window)
    mousey = mfb_get_mouse_y(window)
    println("x: $mousex y: $mousey")
end

macro play(w, h, n, body)
    quote
        window = mfb_open_ex("Visualization", $(esc(w)), $(esc(h)), MiniFB.WF_RESIZABLE)

        mfb_set_mouse_button_callback(window, onclick)

        buffer = zeros(UInt32, $(esc(w)), $(esc(h)))
        for i = 1:$(esc(n))
            Drawing($(esc(w)), $(esc(h)), :image)
            $(esc(body))
            m = permutedims(image_as_matrix!(buffer), (2, 1))
            finish()
            state = mfb_update(window, m)
            if state != MiniFB.STATE_OK
                break
            end
        end
        mfb_close(window)
    end
end


text:
function sciml_train(loss, θ, opt, adtype::DiffEqBase.AbstractADType = GalacticOptim.AutoZygote(), args...; kwargs...)
  @warn("sciml_train has been deprecated in favor of GalacticOptim.jl (https://github.com/SciML/GalacticOptim.jl)")
  optf = GalacticOptim.OptimizationFunction((x, p) -> loss(x), adtype)
  optfunc = GalacticOptim.instantiate_function(optf, θ, adtype, nothing)
  optprob = GalacticOptim.OptimizationProblem(optfunc, θ; kwargs...)
  GalacticOptim.solve(optprob, opt, args...; kwargs...)
end


text:
using MATLAB

x = range(-10.0, stop=10.0, length=500)
current_color="r"
mat"plot($x, sin($x), 'color',$current_color)"


text:
@testset "ch03" begin
    Random.seed!(1234)
    rng = MersenneTwister(1234)

    function create_toy_data(func, sample_size, std, domain = [0.0, 1.0])
        x = collect(range(domain[1], stop = domain[2], length = sample_size))
        shuffle!(rng, x)
        noise = rand(Normal(0.0, std), sample_size)
        return x, func(x) + noise
    end

    function sinusoidal(x)
        return sin.(2 * pi * x)
    end

    @testset "ch03_1" begin
        x = collect(range(-1.0, stop = 1.0, length = 100))
        means = collect(range(-1.0, stop = 1.0, length = 12))
        feature_poly = PolynomialFeature(11)
        feature_gauss = GaussianFeature(means, 0.1)
        feature_sigmoid = SigmoidalFeature(means, 10.0)

        X_polynomial = transform(feature_poly, x)
        X_gaussian = transform(feature_gauss, x)
        X_sigmoidal = transform(feature_sigmoid, x)
    end

    @testset "ch03_2" begin
        x_train, y_train = create_toy_data(sinusoidal, 10, 0.4)
        x_test = collect(range(0, stop = 1.0, length = 100))
        y_test = sinusoidal(x_test)
        means = collect(range(0, stop = 1.0, length = 7))
        feature = GaussianFeature(means, 0.1)

        X_train = transform(feature, x_train)
        X_test = transform(feature, x_test)
        model = LinearRegressor([0], 0)
        fitting(model, X_train, y_train)

        y, var = predict(model, X_test, true)

        model = RidgeRegressor([0.0], 1e-2)
        fitting(model, X_train, y_train)
        y = predict(model, X_test)
    end

    @testset "ch03_3" begin
        means = collect(range(0, stop = 1.0, length = 24))
        feature = GaussianFeature(means, 0.1)

        for a in [1e2, 1.0, 1e-9]
            y_list = []
            for i = 1:100
                x_train, y_train = create_toy_data(sinusoidal, 25, 0.25)
                x_test = collect(range(0, stop = 1.0, length = 100))
                X_train = transform(feature, x_train)
                X_test = transform(feature, x_test)
                model = BayesianRegressor(a, 1.0, 24)
                fitting(model, X_train, y_train)
                y = predict(model, X_test)
            end
        end
    end

    @testset "ch03_4" begin
        function linear(x)
            return -0.3 .+ 0.5 .* x
        end

        x_train, y_train = create_toy_data(linear, 20, 0.1, [-1.0, 1.0])
        x = collect(range(-1.0, stop = 1.0, length = 100))
        y = collect(range(-1.0, stop = 1.0, length = 100))
        #grids = [[i, j] for i in x, j in y];

        model = BayesianRegressor(1.0, 100.0, 2)
        feature = PolynomialFeature(1)
        X_train = transform(feature, x_train)
        X = transform(feature, x)

        dist = MvNormal(model._w_mean, inv(model._w_precision))
        vals = [Distributions.pdf(dist, [i, j]) for j in y, i in x]

        Y_random = predictSampling(model, X, 6)

        for (index, values) in enumerate([[1, 1], [2, 2], [3, 3], [4, 20]])
            first = values[1]
            last = values[2]
            fitting(model, X_train[first:last, :], y_train[first:last])
            dist = MvNormal(model._w_mean, inv(model._w_precision))
            vals = [Distributions.pdf(dist, [i, j]) for j in y, i in x]

            Y_random = predictSampling(model, X, 6)
        end
    end

    @testset "ch03_5" begin
        x_train, y_train = create_toy_data(sinusoidal, 25, 0.25)
        x_test = collect(range(0, stop = 1.0, length = 100))
        y_test = sinusoidal(x_test)

        means = collect(range(0, stop = 1.0, length = 9))
        feature = GaussianFeature(means, 0.1)
        X_train = transform(feature, x_train)
        X_test = transform(feature, x_test)

        model = BayesianRegressor(1e-3, 2.0, 9)

        for (index, ranges) in
            enumerate([[1, 1], [2, 2], [3, 4], [5, 9], [10, 15], [16, 25]])
            first, last = ranges[1], ranges[2]
            fitting(model, X_train[first:last, :], y_train[first:last])
            y, y_std = predict(model, X_test, true)
        end
    end

    @testset "ch03_6" begin
        function cubic(x)
            return x .* (x .- 5.0) .* (x .+ 5.0)
        end

        x_train, y_train = create_toy_data(cubic, 30, 10, [-5.0, 5.0])
        x_test = collect(range(-5.0, stop = 5.0, length = 100))
        evidences = []
        models = []

        for i = 1:6
            feature = PolynomialFeature(i)
            X_train = transform(feature, x_train)
            model = EmpiricalBayesianRegressor(100.0, 100.0, i + 1)
            fitting(model, X_train, y_train, 100)
            push!(evidences, log_evidence(model, X_train, y_train))
            push!(models, model)
        end

        degree = findall(isequal(minimum(evidences)), evidences)
        regressor = models[degree[1]]
        feature = PolynomialFeature(degree[1])

        X_test = transform(feature, x_test)
        y, y_std = predict(regressor, X_test, true)
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
type RenderStates
    ptr::Ptr{Void}

    function RenderStates(ptr::Ptr{Void})
        self = new(ptr)
        finalizer(self, destroy)
        self
    end
end

function RenderStates(blendmode::BlendMode=blend_alpha, shader::Shader=Shader(), texture::Texture=Texture())
    RenderStates(ccall((:sjRenderStates_create, "libjuliasfml"), Ptr{Void}, (BlendMode, Ptr{Void}, Ptr{Void}), blendmode, shader.ptr, texture.ptr))
end

RenderStates(s::Shader) = RenderStates(blend_alpha, s, Texture())
RenderStates(t::Texture) = RenderStates(blend_alpha, Shader(), t)

function set_texture(states::RenderStates, texture::Texture)
    ccall((:sjRenderStates_setTexture, "libjuliasfml"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), states.ptr, texture.ptr)
end

function destroy(states::RenderStates)
    ccall((:sjRenderStates_destroy, "libjuliasfml"), Void, (Ptr{Void},), states.ptr)
end


text:
module FilterOptimizer

const debug = false
const fullblown_optim = true  # false: LsqFit, true: Optim

using DSP
using LinearAlgebra
using NPZ
using Optim
using LineSearches
using LsqFit
if !debug
    using Plots
end
using Polynomials
using Printf
using Random
using Statistics

function main_simple_plot(; ensemble_size=8, cutoff_frac=0.05, noDC=false, killOutputs=1:0)
    # cutoff_frac = 0.05  # [fraction of Nyquist, i.e. of fs/2]
    num_freqs = 100

    filter = make_initial_filter(ensemble_size, cutoff_frac)

    # K88  Filter given from Matlab is transposed.
    # filter = [ 0.101928620359513  -0.144692907668316   0.001818686042631   0.034079418317095   0.017585080356053   0.001315214737816  -0.003444416529060  -0.001436438714533
    #           -0.182345096610246   0.360756599962924  -0.147994984888741  -0.059129118196972   0.002868426988251   0.015712495777574   0.007502730016391  -0.002913799498315
    #            0.325318278686291  -0.644192151825324   0.366462705585208  -0.039188548203213  -0.002501377627620   0.008040320659969   0.004771188361656  -0.001218136494786
    #            0.214429240579235   0.020907800279936  -0.640407310078860   0.438171965206373  -0.001984312130525   0.000672744414457   0.000815968852531   0.000020116755318
    #            0.039388771853078   0.158515140964055   0.021611283755305  -0.627238094500448   0.444963982609944  -0.001482793266190  -0.000658657133152   0.000289435629589
    #           -0.028159300455806   0.079366422909548   0.158022482575611   0.012192601874105  -0.632146169468256   0.444503863315704  -0.000536386191742   0.000149049199835
    #           -0.023962585339180   0.005858919704412   0.078943901425177   0.150008766716844   0.008032778760745  -0.632505199692507   0.445310958035830   0.000012891795086
    #           -0.006277091297826  -0.015051750417134   0.005747354831077   0.076845015215617   0.148923596940809   0.007947473300587  -0.632293311502353   0.445417719102865];
    # filter = filter';

    if noDC
        filter = -mean(filter, dims=1) .+ filter;
    end

    filter[:, killOutputs] .= 0

    # path = "$(homedir())/Downloads/"
    path = "$(homedir())/Code/ButterflyNetwork/software/imaging/imaging/recon/doppler/wallfilter/data/"
    filepath = path * "optimized_filter_cutoff$(@sprintf("%.2f", cutoff_frac))_ensemble$(ensemble_size).npy"
    npzwrite(filepath, filter)

    input_vec, freq_vec = make_input(ensemble_size)
    output_vec = apply_wallfilter(input_vec, filter)
    plot_response(freq_vec, make_desired_response(cutoff_frac), make_response(input_vec), make_response(output_vec))

    FilterOptimizer.PlotWallFilter(FilterOptimizer.EvaluateWallFilter(filter)...)
    filter
end

function main_objective(; noDC=false)
    ensemble_size = 8
    cutoff_frac = 0.10  # [fraction of Nyquist, i.e. of fs/2]
    num_freqs = 100

    input_vec, freq_vec = make_input(ensemble_size, num_freqs)
    input_response = make_response(input_vec)

    desired_response = make_desired_response(cutoff_frac, num_freqs)
    @show desired_response

    initial_filter = make_initial_filter(ensemble_size, cutoff_frac)
    # objective = objective_autocorr
    objective = objective_consistency

    # warm up and debug
    x = initial_filter
    f, fout = EvaluateWallFilter(x, num_freqs)
    output_vec = apply_wallfilter(input_vec, x)
    output_response = make_response(output_vec)
    loss = objective(desired_response, output_response, fout)
    if !debug
        PlotWallFilter(f, fout)
    end

    function optim_iter(x)
        _, fout = EvaluateWallFilter(x, num_freqs)
        output_vec = apply_wallfilter(input_vec, x)
        output_response = make_response(output_vec)
        if fullblown_optim
            loss = objective(desired_response, output_response, fout)
            return loss
        else
            # c = vec(output_response[1] .* exp.(1im .* output_response[2]))
            # return vec([real.(c) ; imag.(c)])
            return [vec(real.(fout)) ; vec(imag.(fout))]
        end
    end

    optim_iter(initial_filter)

    if fullblown_optim
        # Continuous, multivariate, optimization
        # Simulated Annealing, Evolutionary algo
        # CG, etc.: Compute / Auto Grad of f?
        # SPSA: black box, doesnt need gradient of f. Perturbs all params at the same time.
        # FD: black box, doesnt need gradient of f. Perturbs one param at a time.
        result = optimize(optim_iter, copy(initial_filter),
            # LBFGS(linesearch=LineSearches.BackTracking()),
            BFGS(linesearch=LineSearches.BackTracking()),
            # ConjugateGradient(linesearch=LineSearches.BackTracking()),
            # NelderMead(),
            debug ? Optim.Options(show_trace = true, iterations = 1) : Optim.Options(show_trace = true),
            ; autodiff = :forward
            )
        @show result

        optimized_filter = Optim.minimizer(result)
    else
        # t: array of independent variable
        # p: array of model parameters
        model(t, p) = optim_iter(reshape(p, size(initial_filter)))
        tdata = vec(input_vec)
        c = vec(desired_response[1] .* exp.(1im .* desired_response[2]))
        c = repeat(c, ensemble_size)
        ydata = [vec(real.(c)) ; vec(imag.(c))]
        p0 = vec(initial_filter)
        fit = curve_fit(model, tdata, ydata, p0)
        optimized_filter = reshape(fit.param, size(initial_filter))
    end

    if noDC
        optimized_filter = mean(optimized_filter, dims=1) .- optimized_filter;
    end

    npzwrite("$(homedir())/Downloads/optimized_filter_cutoff$(cutoff_frac)_ensemble$(ensemble_size).npy", optimized_filter)

    output_vec = apply_wallfilter(input_vec, optimized_filter)
    if !debug
        PlotWallFilter(EvaluateWallFilter(optimized_filter, num_freqs)...)
        plot_response(freq_vec, desired_response, input_response, make_response(output_vec))
    end

    optimized_filter
end

function make_from_iir_to_matrix(ensemble_size, cutoff_frac)
    responsetype = Highpass(cutoff_frac)
    n = 2  # poles
    rp = 0.5  # [dB] ripple in the passband
    rs = 60  # [dB] stopband attentuation
    # designmethod = Chebyshev1(n, rp)  # ripple in passband (no ripple in stopband)
    designmethod = Chebyshev2(n, rs)  # ripple in stopband (no ripple in passband)
    # designmethod = Elliptic(n, rp, rs)  # n pole elliptic (Cauer) filter with rp dB ripple in the passband and rs dB attentuation in the stopband

    filter = digitalfilter(responsetype, designmethod)
    @show filter
    filter = convert(Biquad, filter)
    # DSP.Filters.Biquad(b0, b1, b2, a1, a2) form:
    # H(z) = (b0 +b1z^−1 +b2z^−2) / (1 +a1z^−1 +a2z^−2)

    display(plot(20log10.(abs.(freqz(filter))), title="20log10.(abs.(freqz(filter)))"))
    display(plot(abs.(phasez(filter)), title="abs.(phasez(filter))"))

    # Put in form: Edward S. Chornoboy Lincoln Lab MIT Technical Report 828, 31 December 1990
    # H(z) = (α0 +α1z^-1 +α2z^-2) / (1 + β1z^-1 + β2z^-2)
    # by association:
    α0 = filter.b0
    α1 = filter.b1
    α2 = filter.b2
    β1 = filter.a1
    β2 = filter.a2

    A = [α1-α0*β1
         α2-α0*β2]
    B = [-β1  -β2
           1    0]
    C = [  1
           0]

    F = zeros(ensemble_size, size(A, 1))
    for n in 1:ensemble_size
        F[n, :] = A' * B^(n-1)
    end

    G = zeros(ensemble_size, ensemble_size)
    for n in 1:ensemble_size
        G[n,n] = α0
    end
    for irow in 1:ensemble_size
        for icol in 1:irow-1
            G[irow,icol] = A' * B^(irow-icol-1) * C
        end
    end

    Pf = F * (F' * F)^-1 * F'

    filter = (I - Pf) * G
    filter'
end

function objective_autocorr(desired_response, output_response, fout)
    desired_mag, desired_pha = desired_response
    output_mag, output_pha = output_response
    e = [desired_mag .- output_mag ; desired_pha .- output_pha]
    norm(e)
end

function objective_consistency(desired_response, output_response, fout)
    # size(fout) = (num_freqs, ensemble_size_out)
    desired_mag, desired_pha = desired_response
    output_mag, output_pha = output_response
    # e = desired_mag .- output_mag  # OK
    # e = [desired_mag .- output_mag ; desired_pha .- output_pha]  # OK
    # e = vec(abs.(fout) .- desired_mag)  # OK
    # e = vec(std(diff(angle.(fout), dims=2), dims=2))  # NaN
    # e = vec(std(angle.(fout), dims=2))  # NaN
    # e = vec(diff(angle.(fout), dims=2))  # OK
    # e = [vec(abs.(fout) .- desired_mag) ; vec(diff(angle.(fout), dims=2))]  # OK but not great filter.
    # e = vec(angle.(fout))  # OK. Bad result with UNWRAP.
    # e = [vec(abs.(fout) .- desired_mag) ; vec(angle.(fout))]  # OK but not great filter.
    # e = vec(diff(diff(angle.(fout), dims=2), dims=2))  # OK
    # e = [vec(abs.(fout) .- desired_mag) ; vec(diff(diff(angle.(fout), dims=2), dims=2))]  # NaN

    # e = [desired_mag .- output_mag ; desired_pha .- output_pha ;
    #      vec(abs.(fout) .- desired_mag) ; vec(diff(unwrap(angle.(fout), dims=1), dims=2))]  # OK

    e = [vec(abs.(fout) .- desired_mag) ; vec(diff(angle.(fout), dims=2))]

    # e = [desired_mag .- output_mag ; vec(angle.(fout))]
    # e = [desired_mag .- output_mag ; vec(std(angle.(fout), dims=2))]  # NaN
    # e = [desired_mag .- output_mag ; vec(std(diff(angle.(fout), dims=2), dims=2))]  # NaN
    # e = abs.(fout) .- desired_mag
    # e = [desired_mag .- output_mag ; vec(diff(abs.(fout), dims=2))]
    # e = [desired_mag .- output_mag ; vec(diff(abs.(fout), dims=2))]  # NaN
    # e = [desired_mag .- output_mag ; vec(diff(abs.(fout), dims=2)) ; vec(diff(angle.(fout), dims=2))]  # NaN
    # e = [vec(desired_mag .- fout) ; vec(diff(abs.(fout), dims=2)) ; vec(diff(angle.(fout), dims=2))]
    # println(length(e))
    norm(e)
end

function make_initial_filter(ensemble_size, cutoff_frac)
    # #SHIFTED TRUNCATED FIR
    # responsetype = Highpass(cutoff_frac)
    # # ripple = 3.0  # [dB] in the passband
    # # n = 2  # poles
    # # designmethod = Chebyshev1(n, ripple)
    # window = hanning(5)
    # designmethod = FIRWindow(window; scale=true)
    # digitalfilter(responsetype, designmethod)

    # POLYNOMIAL REGRESSION FILTER
    # filter_order = 1
    # polynomial_regression_filter(ensemble_size, filter_order)

    # FIR 2zero FILTER
    # fir_2zero_wallfilter(ensemble_size)

    # RANDOM matrix
    Random.seed!(1);
    # need randn and not rand b/c need some negative coefficients in order to get a highpass filter
    # randn(ensemble_size, ensemble_size)
    # rand uniform [-1, +1)
    2*rand(ensemble_size, ensemble_size) .- 1

    # IIR
    make_from_iir_to_matrix(ensemble_size, cutoff_frac)
end

function make_input(ensemble_size, num_freqs=100)
    # generate input freqs
    time_vec = linspace_with_endpoint(0, ensemble_size-1, ensemble_size)
    freq_vec = linspace_with_endpoint(0, 1 / 2, num_freqs)
    outer = freq_vec .* time_vec'
    input_vec = exp.(1im .* outer .* 2π)
    (input_vec, freq_vec)
end

linspace_with_endpoint(start, last, length) = range(start, last, length=length)

linspace_without_endpoint(start, stop, length) = range(start, stop, length=length+2)[2:end-1]

function apply_wallfilter(input_vec, filter_mat)
    input_vec * complex(filter_mat)
end

function make_response(vec)
    mag = abs.(autocorrelation(vec, lag=0))
    pha = angle.(autocorrelation(vec, lag=1))
    (mag, pha)
end

function make_desired_response(cutoff_frac, num_freqs=100)
    n_stops = Int(round(num_freqs * cutoff_frac))
    n_pass = num_freqs - n_stops
    mag = [repeat([0], n_stops) ; repeat([1], n_pass)]
    phase = collect(linspace_with_endpoint(0, π, num_freqs))
    (mag, phase)
end

"""
    Autocorrelation over 1st dimension for the given non-negative lag of 0, 1, 2, ...
    Arguments:
        data -- 2D matrix of complex samples where the first dimension is of
                ensemble size and the second dimension is of 1D reshaped pixels size.
        lag  -- Number of 1-way "lag" samples to offset in autocorrelation.
"""
function autocorrelation(data; lag=0)
    @assert (lag >= 0) "autocorrelation lag must be positive (currently $lag)"
    ens = size(data)[2]
    @assert (ens >= 2 * lag) "autocorrelation requires an ensemble of at least size $(lag*2) (currently $ens)"
    mean(data[:, 1+lag:end] .* conj(circshift(data, (0, lag))[:, 1+lag:end]), dims=2)
end

"""
    Generate the wall filter matrix that can reject polynomials up to
    order filter_order i.e. f[k]= sum over i from 0 to filter_order (a_i*k^i)
    Arguments:
    ensemble_size    -- Number of acquisitions in an ensemble
    filter_order     -- Maximum order of polynomial we want to reject.
"""
function polynomial_regression_filter(filter_length, filter_order)
    # Python:
    #  xx = np.linspace(-1, 1, num=filter_length)
    #  A = xx[:, np.newaxis] ** np.arange(filter_order + 1)
    xx = linspace_with_endpoint(-1, 1, filter_length)
    A = xx .^ (0:filter_order)'
    # Find the column basis of the matrix
    Uonly, = svd(A; full=true)
    # Construct a matrix where each column is a basis of the perp. space of Range(A).
    Uc = Uonly[:, (1 + filter_order+1):end]
    HiP = Uc * Uc'
end

function fir_2zero_wallfilter(ensemble_size)
    """Generate a 2 zero wall filter:  [-1 2 -1]/4
        Arguments:
            ensemble_size    -- Number of acquisitions in an ensemble
    """
    wallfilter = zeros(ensemble_size, ensemble_size)
    for k in 1:ensemble_size - 3
        wallfilter[k:k+2, k] = [-1, 2, -1] / 4.0
    end
    wallfilter
end

function plot_response(freq_vec, desired_response, input_response, output_response)
    desired_mag, desired_pha = desired_response
    input_mag, input_pha = input_response
    output_mag, output_pha = output_response

    # plot frequency response (should attenuate lower freqs)
    # fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
    p1 = plot(2freq_vec, [20log10.(desired_mag.+1e-6) 20log10.(input_mag) 20log10.(output_mag)],
        label=["Desired" "Input Ensemble" "WF'ed Ensemble"],
        title="Autocorrelation (frequency) response of wall filter",
        ylabel="Autocor0 abs [log a.u.]",
        xlims=(0, 1),
        ylims=(-100, Inf),
        lw=2,
    )

    # plot phase response
    p2 = plot(2freq_vec, [desired_pha input_pha output_pha],
        label=["Desired" "Input Ensemble" "WF'ed Ensemble"],
        xlabel="Frequency (0 to Nyquist) [fraction of fs/2]",
        ylabel="Autocor1 phase [radians]",
        xlims=(0, 1),
        ylims=(-π, π),
        lw=2,
    )

    # A Plot is only displayed when returned (a semicolon will suppress the return),
    # or if explicitly displayed with display(plt), gui(), or by adding show = true to your plot command.
    display(plot(p1, p2, size=(800, 500), layout = (2, 1), legend = true))
end

"""
    Func from Karl
"""
function EvaluateWallFilter(wf, num_freqs=100, PhaseOffset=0, PhaseMultiplier=1)
    # Matlab, old: Assume out (Nx1) = WF (NxM) x in (Mx1)
    # M = ensemble size, N = number of filters(output samples)
    # Assume in (1xM) x WF (MxN) = out (1xN)
    M, N = size(wf)

    # Assume Real
    f = linspace_with_endpoint(0, 1, num_freqs)
    xf = exp.(1im * π * f)
    P = length(f)

    # Fourier Output
    fout = zeros(Complex, P, N)
    for n=1:N
        coefs = wf[:,n]
        # Matlab coefs in descending powers, Julia in ascending powers!
        # fcoefs = polyval(coefs[:], xf[:])
        p = Polynomial(reverse(coefs[:]))
        fcoefs = p.(xf[:])
        fcoefs = fcoefs .* xf[:] .^(n-M+PhaseOffset)*PhaseMultiplier  # .* exp.(1im*(n-M/2).*xf[:])
        fout[:,n] = fcoefs[:]
    end

    f, fout
end

function PlotWallFilter(f, fout)
    # plot mag response
    p1 = plot(f, 20log10.(abs.(fout)),
        title="Magnitude Consistency Across Wall Filter Polyphases",
        ylabel="Abs [log a.u.]",
        xlims=(0, 1),
        ylims=(-80, Inf),
        lw=2,
    )
    # plot phase response
    uwpha = reverse(unwrap(reverse(angle.(fout .+ eps()), dims=1), dims=1), dims=1)
    p2 = plot(f, uwpha,
        title="Phase Consistency Across Wall Filter Polyphases",
        xlabel="Frequency (0 to Nyquist) [fraction of fs/2]",
        ylabel="Phase [radians]",
        xlims=(0, 1),
        ylims=(-π*1.1, π*1.1),
        lw=2,
    )
    # plot group delay (d phase / d w)
    p3 = plot(f[1:end-1], diff(uwpha, dims=1) ./ diff(f, dims=1) ./ π,
        title="Group Delay Consistency Across Wall Filter Polyphases",
        xlabel="Frequency (0 to Nyquist) [fraction of fs/2]",
        ylabel="Group delay [samples]",
        xlims=(0, 1),
        ylims=(-2, 2),
        lw=2,
    )
    display(plot(p1, p2, p3, size=(800, 800), layout = (3, 1), legend = true))
end

end # module

# wf = FilterOptimizer.main_simple_plot(noDC=true, killOutputs=1:0);
# wf = FilterOptimizer.main_simple_plot(ensemble_size=8, cutoff_frac=0.05, noDC=true, killOutputs=[1, 2, 7, 8]);

# wf = FilterOptimizer.main_simple_plot(ensemble_size=8, cutoff_frac=0.02, noDC=true, killOutputs=[1, 2, 3]);  # 0.05 in Bfly
wf = FilterOptimizer.main_simple_plot(ensemble_size=8, cutoff_frac=0.03, noDC=true, killOutputs=[1, 2, 3]);
# wf = FilterOptimizer.main_simple_plot(ensemble_size=8, cutoff_frac=0.04, noDC=true, killOutputs=[1, 2, 3]);
# wf = FilterOptimizer.main_simple_plot(ensemble_size=10, cutoff_frac=0.04, noDC=true, killOutputs=[1, 2, 3]);

# wf = FilterOptimizer.main_simple_plot(ensemble_size=8, cutoff_frac=0.03, noDC=false, killOutputs=[1, 2, 3]);

# wf = FilterOptimizer.main_simple_plot(ensemble_size=10, cutoff_frac=0.05, noDC=true, killOutputs=[1, 2, 8, 9, 10]);
# wf = FilterOptimizer.main_simple_plot(ensemble_size=10, cutoff_frac=0.05, noDC=true, killOutputs=[1, 2, 8, 9, 10]);
# wf = FilterOptimizer.main_simple_plot(ensemble_size=12, cutoff_frac=0.50, noDC=true, killOutputs=[1, 2, 10, 11, 12]);

# wf = FilterOptimizer.main_objective(noDC=true);

show(stdout, "text/plain", wf)


text:
# This interface was partially derived from the OSQP MathOptInterface written by T. Koolen and B. Stellato
# https://github.com/oxfordcontrol/OSQP.jl/blob/master/src/MathOptInterfaceOSQP.jl
# certain utility function are taken from the SCS  MathOptInterface:
# https://github.com/JuliaOpt/SCS.jl
using MathOptInterface

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const SparseTriplets{Tv} = Tuple{Vector{<:Integer}, Vector{<:Integer}, Vector{Tv}}

const SingleVariable = MOI.SingleVariable
const Affine = MOI.ScalarAffineFunction{<: AbstractFloat}
const Quadratic = MOI.ScalarQuadraticFunction{<: AbstractFloat}
const VectorOfVariables = MOI.VectorOfVariables
const VectorAffine = MOI.VectorAffineFunction{<: AbstractFloat}

const Interval = MOI.Interval{<: AbstractFloat}
const LessThan = MOI.LessThan{<: AbstractFloat}
const GreaterThan = MOI.GreaterThan{<: AbstractFloat}
const EqualTo = MOI.EqualTo{<: AbstractFloat}
const IntervalConvertible = Union{LessThan, GreaterThan, EqualTo, Interval}

const Zeros = MOI.Zeros
#Nonnegatives = MOI.Nonnegatives
const Nonpositives = MOI.Nonpositives
const SOC = MOI.SecondOrderCone
#const PSDS = MOI.PositiveSemidefiniteConeSquare
const PSDT = MOI.PositiveSemidefiniteConeTriangle
const PSD = Union{MOI.PositiveSemidefiniteConeSquare,MOI.PositiveSemidefiniteConeTriangle}
const SupportedVectorSets = Union{Zeros, MOI.Nonnegatives, Nonpositives, SOC, PSDT, MOI.ExponentialCone, MOI.DualExponentialCone, MOI.PowerCone, MOI.DualPowerCone}

#export sortSets, assign_constraint_row_ranges!, processconstraints, constraint_rows, processobjective, processlinearterms!, symmetrize!, processconstraints!, constant, processconstant!, processlinearpart!, processconstraintset!
export Optimizer

##############################
# MAIN INTERFACE OBJECTS AND FUNCTIONS
##############################

mutable struct Optimizer{T} <: MOI.AbstractOptimizer
    inner::COSMO.Model{T}
    hasresults::Bool
    results::COSMO.Result{T}
    is_empty::Bool
    sense::MOI.OptimizationSense
    objconstant::T
    set_constant::Vector{T}
    constr_constant::Vector{T}
    rowranges::Dict{Int, UnitRange{Int}}
    idxmap::MOIU.IndexMap

    function Optimizer{T}(; user_settings...) where {T <: AbstractFloat}
        inner = COSMO.Model{T}()
        set_MOI_default_settings!(inner.settings)
        length(user_settings) > 0 && set_user_settings!(inner, user_settings)
        hasresults = false
        results = COSMO.Result{T}()
        is_empty = true
        sense = MOI.MIN_SENSE
        objconstant = zero(T)
        set_constant = T[]
        constr_constant = T[]
        rowranges = Dict{Int, UnitRange{Int}}()
        idxmap = MOIU.IndexMap()
        new(inner, hasresults, results, is_empty, sense, objconstant, set_constant, constr_constant, rowranges, idxmap)
    end
end
Optimizer(args...; kwargs...) = Optimizer{DefaultFloat}(args...; kwargs...)

function printIdxmap(idxmap::MOIU.IndexMap)
    println(">>Variable Map with $(length(idxmap.varmap)) entries:")
    dkeys = collect(keys(idxmap.varmap))
    dvalues = collect(values(idxmap.varmap))
    for i=1:length(dkeys)
        println("i=$(i): $(dkeys[i].value) => $(dvalues[i].value)")
    end

     println(">>Constraint Map with $(length(idxmap.conmap)) entries:")
    dkeys = collect(keys(idxmap.conmap))
    dvalues = collect(values(idxmap.conmap))
    for i=1:length(dkeys)
        println("i=$(i): $(dkeys[i].value) => $(dvalues[i].value)")
    end
end

function Base.show(io::IO, obj::Optimizer)
    if obj.is_empty
        print(io,"Empty COSMO - Optimizer")
    elseif obj.hasresults
        print(io,"COSMO - Optimizer\n- Has results: $(obj.hasresults)\n- Objective constant: $(obj.objconstant)\n- Problem status: $(obj.results.status)\n- Optimal objective: $(round(obj.results.obj_val, digits = 3))\n- Sense: $(obj.sense)\n- Iterations: $(obj.results.iter)\n- Solve time: $(round.(obj.results.times.solver_time*1000,digits=2))ms")
    else
        print(io,"COSMO - Optimizer\n- Has results: $(obj.hasresults)\n- Objective constant: $(obj.objconstant)\n- Sense: $(obj.sense)")
    end
end

function set_user_settings!(inner::COSMO.Model, user_settings)
    for (k, v) in user_settings
        !in(k, fieldnames(typeof(inner.settings))) && error("The user setting $(k) is not defined.")
        setfield!(inner.settings, k, v)
    end
end

# Function to set MOI specific default values from COSMO that are different from plain COSMO
# the default behaviour for MOI solver is verbose printing turned on
function set_MOI_default_settings!(settings)
    settings.verbose = true
end

# MG: function to reset otimizer
function MOI.empty!(optimizer::Optimizer{T}) where {T <: AbstractFloat}
    COSMO.empty_model!(optimizer.inner)
    optimizer.hasresults = false
    optimizer.results = COSMO.Result{T}()
    optimizer.is_empty = true
    optimizer.sense = MOI.MIN_SENSE # model parameter, so needs to be reset
    optimizer.objconstant = zero(T)
    optimizer.set_constant = T[]
    optimizer.constr_constant = T[]   # the 5 in (3 * x1 + 2 * x2 + 5 ≤ 10 )
    optimizer.idxmap = MOIU.IndexMap()
    optimizer.rowranges = Dict{Int, UnitRange{Int}}()
    optimizer
end

# MG: check if optimizer object is empty
MOI.is_empty(optimizer::Optimizer) = optimizer.is_empty

##############################
# MODEL --> SOLVER FORMAT FUNCTIONS
##############################

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; copy_names = false)
    copy_names && error("Copying names is not supported.")
    MOI.empty!(dest)
    idxmap = MOIU.IndexMap(dest, src)
    assign_constraint_row_ranges!(dest.rowranges, idxmap, src)
    processconstraints!(dest, src, idxmap, dest.rowranges)
    pre_allocate_variables!(dest.inner)

    dest.is_empty = false
    dest.idxmap = idxmap
    # pass attributes, e.g. objective, warm strarting vars, etc.
    dest.sense = MOI.get(src, MOI.ObjectiveSense())
    COSMO.pass_attributes!(dest, src, idxmap)

    # if no cost function is provided, allocate P, q with correct sizes
    dest.sense == MOI.FEASIBILITY_SENSE && COSMO.allocate_cost_variables!(dest)
    return idxmap
end
function MOI.optimize!(optimizer::Optimizer)
    optimizer.results = COSMO.optimize!(optimizer.inner)
    optimizer.hasresults = true
    nothing
end

function MOIU.IndexMap(dest::Optimizer, src::MOI.ModelLike)
    idxmap = MOIU.IndexMap()
    # map model variable indices to solver variable indices
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    for i in eachindex(vis_src)
        idxmap[vis_src[i]] = MOI.VariableIndex(i)
    end
    # map model constraint indices to solver constraint indices. For now this is important since solver expects following
    # order: {0}-variables, R+-variables, SOCP cones, psd cones
    # LOCs = MOI.get(src, MOI.ListOfConstraints())
    # sort!(LOCs, by=x-> sort_sets(x[2]))
    i = 0
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        MOI.supports_constraint(dest, F, S) || throw(MOI.UnsupportedConstraint{F, S})
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        for ci in cis_src
            i += 1
            idxmap[ci] = CI{F, S}(i)
        end
    end
    idxmap
end

# This is important for set merging
sort_sets(s::Union{Type{<: MOI.EqualTo}, Type{<: MOI.Zeros}}) = 1
sort_sets(s::Union{Type{ <: LessThan}, Type{ <: GreaterThan}, Type{ <: MOI.Nonnegatives}, Type{ <: MOI.Nonpositives}}) = 2
sort_sets(s::Type{MOI.Interval}) = 3
sort_sets(s::Type{<: MOI.AbstractSet}) = 4

# this has to be in the right order
# returns a Dictionary that maps a constraint index to a range
function assign_constraint_row_ranges!(rowranges::Dict{Int, UnitRange{Int}}, idxmap::MOIU.IndexMap, src::MOI.ModelLike)
    startrow = 1
    LOCs = MOI.get(src, MOI.ListOfConstraints())
    sort!(LOCs, by=x-> sort_sets(x[2]))
    for (F, S) in LOCs
        # Returns Array of constraint indexes that match F,S, each constraint index is just a type-safe wrapper for Int64
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        for ci_src in cis_src
            set = MOI.get(src, MOI.ConstraintSet(), ci_src)
            ci_dest = idxmap[ci_src]
            endrow = startrow + MOI.dimension(set) - 1
            rowranges[ci_dest.value] = startrow : endrow
            startrow = endrow + 1
        end
    end
end

function constraint_rows(rowranges::Dict{Int, UnitRange{Int}}, ci::CI{<:Any, <:MOI.AbstractScalarSet})
    rowrange = rowranges[ci.value]
    length(rowrange) == 1 || error()
    first(rowrange)
end
constraint_rows(rowranges::Dict{Int, UnitRange{Int}}, ci::CI{<:Any, <:MOI.AbstractVectorSet}) = rowranges[ci.value]

# this is handled explicitly in `copy_to` before the objective function is processed
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveSense, sense) = nothing

function allocate_cost_variables!(optimizer::Optimizer)
    m, n = size(optimizer.inner.p.A)
    optimizer.inner.p.P = spzeros(n, n)
    optimizer.inner.p.q = zeros(n)
    return nothing
end

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<:Union{MOI.SingleVariable, Affine, Quadratic}}) = true
function MOIU.load(dest::Optimizer{T}, ::MOI.ObjectiveFunction, f::MOI.SingleVariable) where {T <: AbstractFloat}
    idxmap = dest.idxmap
    n = MOI.get(dest, MOI.NumberOfVariables())
    dest.inner.p.q = zeros(T, n)
    dest.inner.p.q[idxmap[f.variable].value] = one(T)
    dest.inner.p.P = spzeros(T, n, n)
    dest.objconstant = zero(T)
    dest.objconstant = apply_sense!(dest.sense, dest.inner.p.P, dest.inner.p.q, dest.objconstant)
    return nothing
end

function MOIU.load(dest::Optimizer{T}, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction{T}) where {T <: AbstractFloat}
    idxmap = dest.idxmap
    n = MOI.get(dest, MOI.NumberOfVariables())
    dest.inner.p.P = spzeros(T, n, n)

    dest.inner.p.q = zeros(T, n)
    processlinearterms!(dest.inner.p.q, f.terms)# idxmap)

    dest.objconstant = f.constant
    dest.objconstant = apply_sense!(dest.sense, dest.inner.p.P, dest.inner.p.q, dest.objconstant)
    return nothing
end

function MOIU.load(dest::Optimizer{T}, ::MOI.ObjectiveFunction, f::MOI.ScalarQuadraticFunction{T}) where {T <: AbstractFloat}
    idxmap = dest.idxmap
    n = MOI.get(dest, MOI.NumberOfVariables())

    I = [Int(idxmap[term.variable_index_1].value) for term in f.quadratic_terms]
    J = [Int(idxmap[term.variable_index_2].value) for term in f.quadratic_terms]
    V = [term.coefficient for term in f.quadratic_terms]
    symmetrize!(I, J, V)
    dest.inner.p.P = sparse(I, J, V, n, n)

    dest.inner.p.q = zeros(T, n)
    processlinearterms!(dest.inner.p.q, f.affine_terms)#, idxmap)
    dest.objconstant = f.constant

    dest.objconstant = apply_sense!(dest.sense, dest.inner.p.P, dest.inner.p.q, dest.objconstant)
    return nothing
end

function MOIU.load(dest::Optimizer, ::MOI.ObjectiveFunction, f)
    throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction{typeof(f)}()))
end

function apply_sense!(sense::MOI.OptimizationSense, P::AbstractMatrix{T}, q::AbstractVector{T}, c::T) where {T <: AbstractFloat}
    if sense == MOI.MAX_SENSE
        LinearAlgebra.rmul!(P, -one(T))
        LinearAlgebra.rmul!(q, -one(T))
        c = -c
    end
    return c
end

function processlinearterms!(q, terms::Vector{<:MOI.ScalarAffineTerm}, idxmapfun::Function = identity)
    for term in terms
        var = term.variable_index
        coeff = term.coefficient
        q[var.value] += coeff
    end
end

function symmetrize!(I::Vector{Ti}, J::Vector{Ti}, V::Vector{Tv}) where {Tv <: AbstractFloat, Ti <: Integer}
    n = length(V)
    (length(I) == length(J) == n) || error()
    for i = 1 : n
        if I[i] != J[i]
            push!(I, J[i])
            push!(J, I[i])
            push!(V, V[i])
        end
    end
end

# The following utility functions are from https://github.com/JuliaOpt/SCS.jl

# Scale coefficients depending on rows index
# rows: List of row indices
# coef: List of corresponding coefficients
# minus: if true, multiply the result by -1
# d: dimension of set
# rev: if true, we unscale instead (e.g. divide by √2 instead of multiply for PSD cone)
_scalecoef(rows, coef::Union{Vector{T}, T}, minus, ::Type{<:MOI.AbstractSet}, rev) where {T <: AbstractFloat} = minus ? -coef : coef
_scalecoef(rows, coef::Union{Vector{T}, T}, minus, ::Union{Type{<:MOI.LessThan}, Type{<:MOI.Nonpositives}, Type{<:MOI.EqualTo}}, rev) where {T <: AbstractFloat} = minus ? coef : -coef
function _scalecoef(rows, coef::Union{Vector{T}, T}, minus, ::Type{MOI.PositiveSemidefiniteConeTriangle}, rev) where {T <: AbstractFloat}
    scaling = minus ? -one(T) : one(T)
    scaling2 = rev ? scaling / sqrt(T(2)) : scaling * sqrt(T(2))
    output = copy(coef)
    idx = 0
    for i in 1:length(output)
        # See https://en.wikipedia.org/wiki/Triangular_number#Triangular_roots_and_tests_for_triangular_numbers
        is_diagonal_index = isinteger(sqrt(8 * rows[i] + 1))
        if is_diagonal_index
            output[i] *= scaling
        else
            output[i] *= scaling2
        end
    end
    output
end
# Unscale the coefficients in `coef` with respective rows in `rows` for a set `s` and multiply by `-1` if `minus` is `true`.
scalecoef(rows, coef, minus, s) = _scalecoef(rows, coef, minus, typeof(s), false)
# Unscale the coefficients in `coef` with respective rows in `rows` for a set of type `S`
unscalecoef(rows, coef, S::Type{<:MOI.AbstractSet}) = _scalecoef(rows, coef, false, S, true)

# This helper function is to provide scale- and unscalecoef with the nomimal rows in the case of at PSDTriangle constraint
# this is because for _scalecoef the actual row in A doesn't matter, what matters is the row in the upper triangle matrix
nom_rows(rows, s::Type{MOI.PositiveSemidefiniteConeTriangle}) = 1:length(rows)
nom_rows(rows, s::Type{<:MOI.AbstractSet}) = rows

output_index(t::MOI.VectorAffineTerm) = t.output_index
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
get_var_index(t::MOI.ScalarAffineTerm) = t.variable_index
get_var_index(t::MOI.VectorAffineTerm) = get_var_index(t.scalar_term)
coefficient(t::MOI.ScalarAffineTerm) = t.coefficient
coefficient(t::MOI.VectorAffineTerm) = coefficient(t.scalar_term)

function processconstraints!(optimizer::Optimizer{T}, src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}}) where {T <: AbstractFloat}

    m = mapreduce(length, +, values(rowranges), init=0)

    b = zeros(T, m)
    constant = zeros(T, m)
    I = Int[]
    J = Int[]
    V = T[]
    COSMOconvexSets = Array{COSMO.AbstractConvexSet{T}}(undef, 0)
    set_constant = zeros(T, m)
    # loop over constraints and modify A, l, u and constants
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        processconstraints!((I, J, V), b, COSMOconvexSets, constant, set_constant, src, idxmap, rowranges, F, S)
    end
    optimizer.set_constant = set_constant
    # subtract constant from right hand side
    n = MOI.get(src, MOI.NumberOfVariables())
    A = sparse(I, J, V, m, n)

    # in a second step: merge some of the sets
    convex_sets = COSMO.merge_sets(COSMOconvexSets)

    # store constraint matrices in inner model
    model = optimizer.inner
    model.p.A = A
    model.p.b = b
    model.p.C = CompositeConvexSet{T}(deepcopy(convex_sets))
    model.p.model_size = [size(A, 1); size(A,2 )]
    optimizer.constr_constant = constant
    return nothing
end

function processconstraints!(triplets::SparseTriplets{T}, b::Vector{T}, COSMOconvexSets::Array{COSMO.AbstractConvexSet{T}}, constant::Vector{T}, set_constant::Vector{T},
        src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}},
        F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet}) where {T <: AbstractFloat}
    cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
    # loop over all constraints of same (F,S)
    for ci in cis_src
        s = MOI.get(src, MOI.ConstraintSet(), ci)
        f = MOI.get(src, MOI.ConstraintFunction(), ci)
        rows = constraint_rows(rowranges, idxmap[ci])
        processConstant!(b, rows, f, s)
        constant[rows] = b[rows]
        if typeof(s) <: MOI.AbstractScalarSet && !(typeof(s) <: MOI.Interval)
            set_constant[rows] =  MOI.constant(s)
        end
        processConstraint!(triplets, f, rows, idxmap, s)
        processSet!(b, rows, COSMOconvexSets, s)
    end
    nothing
end

##############################
# PROCESS FUNCTIONS
##############################
constant(f::MOI.SingleVariable) = 0
constant(f::MOI.ScalarAffineFunction) = f.constant

# process constant for functions Union{Affine, SingleVariable}
function processConstant!(b::AbstractVector{T}, row::Int, f::Affine, s) where {T <: AbstractFloat}
    b[row] = scalecoef(row, constant(f), false, s)
    nothing
end

function processConstant!(b::AbstractVector{T}, rows::UnitRange{Int}, f::VectorOfVariables, s) where {T <: AbstractFloat}
    b[rows] .= zero(T)
    nothing
end

# process constant for functions VectorAffineFunction{Float64}
function processConstant!(b::AbstractVector{T}, rows::UnitRange{Int}, f::VectorAffine, s) where {T <: AbstractFloat}
        b[rows] = scalecoef(nom_rows(rows, typeof(s)), f.constants, false, s)
    nothing
end

# process function like f(x) = coeff*x_1
function processConstraint!(triplets::SparseTriplets{T}, f::MOI.ScalarAffineFunction, row::Int, idxmap::MOIU.IndexMap, s::MOI.AbstractSet) where {T <: AbstractFloat}
    (I, J, V) = triplets
    for term in f.terms
        push!(I, row)
        push!(J, idxmap[term.variable_index].value)
        push!(V, scalecoef(row, term.coefficient, true, s))
    end
end

function processConstraint!(triplets::SparseTriplets{T}, f::MOI.VectorOfVariables, rows::UnitRange{Int}, idxmap::MOIU.IndexMap, s::MOI.AbstractSet) where {T <: AbstractFloat}
    (I, J, V) = triplets
    cols = zeros(length(rows))
    for (i, var) in enumerate(f.variables)
        cols[i] = idxmap[var].value
    end
    append!(I, rows)
    append!(J, cols)
    append!(V, scalecoef(nom_rows(rows, typeof(s)), ones(T, length(rows)), true, s))
    nothing
end

function processConstraint!(triplets::SparseTriplets{T}, f::MOI.VectorAffineFunction, rows::UnitRange{Int}, idxmap::MOIU.IndexMap, s::MOI.AbstractSet) where {T <: AbstractFloat}
    (I, J, V) = triplets

    vis_src = get_var_index.(f.terms)
    vis_dest = map(vi -> idxmap[vi], vis_src)

    A = sparse(output_index.(f.terms), map(x-> x.value, vis_dest), coefficient.(f.terms))
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    A_I, A_J, A_V = findnz(A)

    offset = first(rows) - 1

    append!(J, A_J)
    append!(V, scalecoef(A_I, A_V, true, s))
    append!(I, A_I .+ offset)
end

##############################
# PROCESS SETS
##############################

# process the following sets Union{LessThan, GreaterThan, EqualTo}
function processSet!(b::Vector{T}, row::Int, cs, s::LessThan) where {T <: AbstractFloat}
    b[row] += s.upper
    push!(cs, COSMO.Nonnegatives{T}(1))
    nothing
end
function processSet!(b::Vector{T}, row::Int, cs, s::GreaterThan) where {T <: AbstractFloat}
    b[row] -= s.lower
    push!(cs, COSMO.Nonnegatives{T}(1))
    nothing
end
function processSet!(b::Vector{T}, row::Int, cs, s::EqualTo) where {T <: AbstractFloat}
    b[row] += s.value
    push!(cs, COSMO.ZeroSet{T}(1))
    nothing
end

function processSet!(b::Vector{T}, row::Int, cs, s::MOI.Interval) where {T <: AbstractFloat}
    push!(cs, COSMO.Box{T}([s.lower], [s.upper]))
    nothing
end

# process the following sets Zeros
function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::Zeros) where {T <: AbstractFloat}
    push!(cs, COSMO.ZeroSet{T}(length(rows)))
    nothing
end

# process the following sets Union{Zeros, Nonnegatives, Nonpositives}
function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.Nonnegatives) where {T <: AbstractFloat}
    push!(cs, COSMO.Nonnegatives{T}(length(rows)))
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.Nonpositives) where {T <: AbstractFloat}
    push!(cs, COSMO.Nonnegatives{T}(length(rows)))
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::SOC) where {T <: AbstractFloat}
    push!(cs, COSMO.SecondOrderCone{T}(length(rows)))
    nothing
end

# function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeSquare)
#     push!(cs, COSMO.PsdCone{Float64}(length(rows)))
#     nothing
# end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeTriangle) where {T <: AbstractFloat}
    push!(cs, COSMO.PsdConeTriangle{T}(length(rows)))
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.ExponentialCone) where {T <: AbstractFloat}
    push!(cs, COSMO.ExponentialCone{T}())
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.DualExponentialCone) where {T <: AbstractFloat}
    push!(cs, COSMO.DualExponentialCone{T}())
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.PowerCone) where {T <: AbstractFloat}
    push!(cs, COSMO.PowerCone{T}(T(s.exponent)))
    nothing
end

function processSet!(b::Vector{T}, rows::UnitRange{Int}, cs, s::MOI.DualPowerCone) where {T <: AbstractFloat}
    push!(cs, COSMO.DualPowerCone{T}(T(s.exponent)))
    nothing
end

# to reduce function calls combine, individual ZeroSets, Nonnegatives and Box constraints
# This assumes that C is ordered by set type
function merge_sets(C::Array{COSMO.AbstractConvexSet{T}}) where {T <: AbstractFloat}

    z_ind = findall(set -> typeof(set) <: COSMO.ZeroSet, C)
    nn_ind = findall(set -> typeof(set) <: COSMO.Nonnegatives, C)
    box_ind = findall(set -> typeof(set) <: COSMO.Box, C)
    other_ind = findall( x-> !in(x, union(z_ind, nn_ind, box_ind)), collect(1:1:length(C)))

    # if none of these sets are present, do nothing
    length(other_ind) == length(C) && return C

    length(z_ind) > 0 ? (num_z = 1) : (num_z = 0)
    length(nn_ind) > 0 ? (num_nn = 1) : (num_nn = 0)
    length(box_ind) > 0 ? (num_box = 1) : (num_box = 0)

    num_merged_sets = length(other_ind) + num_z + num_nn + num_box
    merged_sets = Array{COSMO.AbstractConvexSet{T}}(undef, num_merged_sets)
    set_ind = 1
    length(z_ind) > 0 && (set_ind = COSMO.merge_set!(merged_sets, z_ind, C, set_ind, COSMO.ZeroSet{T}))
    length(nn_ind) > 0 && (set_ind = COSMO.merge_set!(merged_sets, nn_ind, C, set_ind, COSMO.Nonnegatives{T}))
    length(box_ind) > 0 && (set_ind = COSMO.merge_box!(merged_sets, box_ind, C, set_ind))

    for other in other_ind
        merged_sets[set_ind] = C[other]
        set_ind += 1
    end
    return merged_sets
end

function merge_set!(merged_sets::Array{COSMO.AbstractConvexSet{T}, 1}, ind::Array{Int64, 1}, C::Array{<: COSMO.AbstractConvexSet, 1}, set_ind::Int64, set_type::DataType) where {T <: AbstractFloat}
        if length(ind) > 1
            combined_dim = sum(x -> x.dim, C[ind])
        else
            combined_dim = C[ind[1]].dim
        end
        merged_sets[set_ind] = set_type(combined_dim)
        return set_ind + 1
end

function merge_box!(merged_sets::Array{COSMO.AbstractConvexSet{T}, 1}, ind::Array{Int64, 1}, C::Array{<: COSMO.AbstractConvexSet{T}, 1}, set_ind::Int64) where {T <: AbstractFloat}
        if length(ind) > 1
            combined_dim = sum(x -> x.dim, C[ind])
            l = zeros(T, combined_dim)
            u = zeros(T, combined_dim)
            row = 1
            for box in C[ind]
                l[row:row + box.dim - 1] = box.l
                u[row:row + box.dim - 1] = box.u
                row += box.dim
            end
            merged_sets[set_ind] = COSMO.Box{T}(l, u)
        else
            merged_sets[set_ind] = C[ind[1]]
        end

        return set_ind + 1
end

function pass_attributes!(dest::Optimizer{T}, src::MOI.ModelLike, idxmap::MOIU.IndexMap, pass_attr::Function=MOI.set) where {T <: AbstractFloat}

    # Copy model attributes, e.g. ObjectiveFunction
    attrs = MOI.get(src, MOI.ListOfModelAttributesSet())
    attrs = filter(attr -> !(attr isa MOI.Name), attrs)

    MOIU._pass_attributes(dest, src, idxmap, attrs, tuple(), tuple(), tuple(), MOIU.load)

    # Copy variable attributes, e.g. VariablePrimalStart
    var_attr = MOI.get(src, MOI.ListOfVariableAttributesSet())
    # We don't support MOI.VariableNames at this point
    var_attr = filter(attr -> !(attr isa MOI.VariableName), var_attr)
    if !isempty(var_attr)
        vis_src = MOI.get(src, MOI.ListOfVariableIndices())
        vis_dest = map(vi -> idxmap[vi], vis_src)
        MOIU._pass_attributes(dest, src, idxmap, var_attr, (VI,), (vis_src,), (vis_dest,), MOIU.load)
    end

    # Copy constraint attributes, e.g. ConstraintPrimalStart
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        attrs = MOI.get(src, MOI.ListOfConstraintAttributesSet{F, S}())
        # we don't support MOI.ConstraintName at this point
        attrs = filter(attr -> !(attr isa MOI.ConstraintName), attrs)
        if !isempty(attrs)
            cis_dest = map(ci -> idxmap[ci], cis_src)
            MOIU._pass_attributes(dest, src, idxmap, attrs, (CI{F, S},), (cis_src,), (cis_dest,), MOIU.load)
        end
    end
    return nothing
end

##############################
# GET / SET : OPTIMIZER ATTRIBUTES
##############################

## Supported Optimizer attributes, which  objective functions are supported:
# allow returning the inner model
MOI.get(optimizer::Optimizer, ::MOI.RawSolver) = optimizer.inner
# allow returning the number of results available
MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = optimizer.hasresults ? 1 : 0
MOI.get(optimizer::Optimizer, ::MOI.NumberOfVariables) = length(optimizer.inner.vars.x)
function MOI.get(optimizer::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(optimizer.inner.vars.x)]
end

# Get objective function
function MOI.get(optimizer::Optimizer, a::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, a)
    rawobj = optimizer.results.obj_val
    if optimizer.sense == MOI.MAX_SENSE
        return -rawobj - optimizer.objconstant
    else
        return rawobj + optimizer.objconstant
    end
end

MOI.get(optimizer::Optimizer, a::MOI.SolveTime) = optimizer.results.times.solver_time
MOI.get(optimizer::Optimizer, a::MOI.SolverName) = "COSMO"

# Get Termination Status
function MOI.get(optimizer::Optimizer, a::MOI.TerminationStatus)
    # Note that the :Dual_infeasible and :Primal_infeasible are mapped to MOI.Success
    # because COSMO can return a proof of infeasibility. For the same reason,
    # :Primal_infeasible_inaccurate is mapped to MOI.AlmostSuccess
    optimizer.hasresults || return MOI.OPTIMIZE_NOT_CALLED

    status = optimizer.results.status

    if status == :Dual_infeasible
        return MOI.DUAL_INFEASIBLE
    elseif status == :Primal_infeasible
        return MOI.INFEASIBLE
    elseif status == :Max_iter_reached
        return MOI.ITERATION_LIMIT
    elseif status == :Solved
        return MOI.OPTIMAL
    elseif status == :Time_limit_reached
        return MOI.TIME_LIMIT
    else
        return MOI.NUMERICAL_ERROR
    end
end

# Get Primal Status
function MOI.get(optimizer::Optimizer, a::MOI.PrimalStatus)
    if !(1 <= a.N <= MOI.get(optimizer, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end

    status = optimizer.results.status
    if status == :Primal_infeasible
        return MOI.INFEASIBLE_POINT
    elseif status == :Solved
        return MOI.FEASIBLE_POINT
    elseif status == :Max_iter_reached
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Dual_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    else
        return MOI.NO_SOLUTION
    end
end

# Get Dual Status
function MOI.get(optimizer::Optimizer, a::MOI.DualStatus)
    if !(1 <= a.N <= MOI.get(optimizer, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end

    status = optimizer.results.status
    if status == :Dual_infeasible
        return MOI.INFEASIBLE_POINT
    elseif status == :Primal_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == :Solved
        return MOI.FEASIBLE_POINT
    elseif status == :Max_iter_reached
        return MOI.NEARLY_FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return String(optimizer.results.status)
end

_unshift(optimizer::Optimizer, offset, value, s) = value
_unshift(optimizer::Optimizer, offset, value, s::Type{<:MOI.AbstractScalarSet}) = value + optimizer.set_constant[offset]
_shift(optimizer::Optimizer, offset, value, s) = value
_shift(optimizer::Optimizer, offset, value, s::Type{<:MOI.AbstractScalarSet}) = value - optimizer.set_constant[offset]

function MOI.get(optimizer::Optimizer, a::MOI.ConstraintPrimal, ci_src::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, a)
    # offset = constroffset(optimizer, ci)
    # rows = constrrows(optimizer, ci)
    rows = COSMO.constraint_rows(optimizer.rowranges, ci_src)
    c_primal = unscalecoef(nom_rows(rows, S), optimizer.results.s[rows], S)
    # (typeof(c_primal) <: AbstractArray && length(c_primal) == 1) && (c_primal = first(c_primal))
    return _unshift(optimizer, rows, c_primal, S)
end

function MOI.get(optimizer::Optimizer, a::MOI.ConstraintDual, ci_src::CI{<:MOI.AbstractFunction, S}) where S <: MOI.AbstractSet
    MOI.check_result_index_bounds(optimizer, a)
    rows = constraint_rows(optimizer.rowranges, ci_src)
    return unscalecoef(nom_rows(rows, S), optimizer.results.y[rows], S)
end

## Variable attributes:
function MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(optimizer, a)
    return optimizer.results.x[vi.value]
end

## Warm starting:
MOI.supports(::Optimizer, a::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}) = true
function MOIU.load(optimizer::Optimizer, a::MOI.VariablePrimalStart, vi::VI, value::Real)
    MOI.is_empty(optimizer) && throw(MOI.CannotSetAttribute(a))
    COSMO.warm_start_primal!(optimizer.inner, value, vi.value)
end
MOIU.load(optimizer::Optimizer, a::MOI.VariablePrimalStart, vi::VI, value::Nothing) = nothing

MOI.supports(::Optimizer, a::MOI.ConstraintPrimalStart, ::Type{<:MOI.ConstraintIndex}) = true
function MOIU.load(optimizer::Optimizer, a::MOI.ConstraintPrimalStart, ci::CI{<:MOI.AbstractFunction, S}, value) where S <: MOI.AbstractSet
    (value == nothing || isa(value, Array{Nothing, 1})) && return nothing
    MOI.is_empty(optimizer) && throw(MOI.CannotSetAttribute(a))
    rows = constraint_rows(optimizer.rowranges, ci)
    # this undoes the scaling and shifting that was used in get(MOI.ConstraintPrimal)
    # Off-diagonal entries of slack variable of a PSDTriangle constraint has to be scaled by sqrt(2)
    value = _shift(optimizer, rows, value, S)
    COSMO.warm_start_slack!(optimizer.inner, _scalecoef(nom_rows(rows, S), value, false, S, false), rows)
end

MOI.supports(::Optimizer, a::MOI.ConstraintDualStart, ::Type{<:MOI.ConstraintIndex}) = true
function MOIU.load(optimizer::Optimizer, a::MOI.ConstraintDualStart, ci::CI{<:MOI.AbstractFunction, S}, value) where S <: MOI.AbstractSet
    (value == nothing || isa(value, Array{Nothing, 1})) && return nothing
    MOI.is_empty(optimizer) && throw(MOI.CannotSetAttribute(a))
    rows = constraint_rows(optimizer.rowranges, ci)
    # this undoes the scaling that was used in get(MOI.ConstraintDual)
    # Off-diagonal entries of dual variable of a PSDTriangle constraint has to be scaled by sfqrt(2)
    COSMO.warm_start_dual!(optimizer.inner, _scalecoef(nom_rows(rows, S), value, false, S, false), rows)
    nothing
end

 MOI.supports(::Optimizer, ::MOI.RawParameter) = true
function MOI.set(optimizer::Optimizer, p::MOI.RawParameter, value)
    setfield!(optimizer.inner.settings, Symbol(p.name), value)
end

function MOI.get(optimizer::Optimizer, p::MOI.RawParameter)
    if in(Symbol(p.name), fieldnames(typeof(optimizer.inner.settings)))
        return getfield(optimizer.inner.settings, Symbol(p.name))
    end
    error("RawParameter with name $(p.name) is not set.")
end

# MOI.Silent == COSMO.Settings.verbose = false
MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.inner.settings.verbose = !value
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = !optimizer.inner.settings.verbose

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(optimizer::Optimizer{T}, ::MOI.TimeLimitSec, value::Real) where {T <: AbstractFloat}
    MOI.set(optimizer, MOI.RawParameter("time_limit"), T(value))
end

function MOI.set(optimizer::Optimizer{T}, attr::MOI.TimeLimitSec, ::Nothing) where {T <: AbstractFloat}
    MOI.set(optimizer, MOI.RawParameter("time_limit"), zero(T))
end
function MOI.get(optimizer::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(optimizer, MOI.RawParameter("time_limit"))
end

##############################
# SUPPORT OBJECTIVE AND SET FUNCTIONS
##############################

## Supported Constraints
MOI.supports_constraint(optimizer::Optimizer, ::Type{<:Affine}, ::Type{<:IntervalConvertible}) = true
MOI.supports_constraint(optimizer::Optimizer, ::Type{<:Union{VectorOfVariables, VectorAffine}}, ::Type{<:SupportedVectorSets}) = true


text:
include("BacksteppingPositionController.jl")
include("LQR.jl")
include("PID.jl")