text:
module GLPKInterfaceBase

using Compat
using Compat.SparseArrays
using Compat.LinearAlgebra
import GLPK

import MathProgBase
const MPB = MathProgBase

export GLPKMathProgModel

abstract type GLPKMathProgModel <: MPB.AbstractLinearQuadraticModel end

function MPB.loadproblem!(lpm::GLPKMathProgModel, filename::AbstractString)
    if endswith(filename, ".mps") || endswith(filename, ".mps.gz")
       GLPK.read_mps(lpm.inner, GLPK.MPS_FILE, filename)
    elseif endswith(filename, ".lp") || endswith(filename, ".lp.gz")
       GLPK.read_lp(lpm.inner, filename)
    elseif endswith(filename, ".prob") || endswith(filename, ".prob.gz")
       GLPK.read_prob(lpm.inner, filename)
    else
       error("unrecognized input format extension in $filename")
    end
end

nonnull(x) = (x != nothing && !isempty(x))

function MPB.loadproblem!(lpm::GLPKMathProgModel, A::AbstractMatrix, collb, colub, obj, rowlb, rowub, sense)
    lp = lpm.inner

    m, n = size(A)

    function checksize(x, l, str)
        if nonnull(x) && length(x) != l
            error("size of $str is incompatible with size of A")
        end
    end

    checksize(collb, n, "collb")
    checksize(colub, n, "colub")
    checksize(rowlb, m, "rowlb")
    checksize(rowub, m, "rowub")
    checksize(obj, n, "obj")

    (ia, ja, ar) = findnz(sparse(A))

    GLPK.erase_prob(lp)

    function getbounds(lb, ub, i)
        if nonnull(lb) && nonnull(ub) && lb[i] != -Inf && ub[i] != Inf
            if lb[i] != ub[i]
                return lb[i], ub[i], GLPK.DB
            else
                return lb[i], ub[i], GLPK.FX
            end
        elseif nonnull(lb) && lb[i] != -Inf
            return lb[i], Inf, GLPK.LO
        elseif nonnull(ub) && ub[i] != Inf
            return -Inf, ub[i], GLPK.UP
        else
            return -Inf, Inf, GLPK.FR
        end
    end

    m > 0 && GLPK.add_rows(lp, m)
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for r = 1:m
        #println("  r=$r b=$(b[r])")
        l, u, t = getbounds(rowlb, rowub, r)
        GLPK.set_row_bnds(lp, r, t, l, u)
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)

    n > 0 && GLPK.add_cols(lp, n)
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for c = 1:n
        #println("  r=$r b=$(b[r])")
        l, u, t = getbounds(collb, colub, c)
        GLPK.set_col_bnds(lp, c, t, l, u)
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)

    if nonnull(obj)
        for c = 1:n
            GLPK.set_obj_coef(lp, c, obj[c])
        end
    else
        for c = 1:n
            GLPK.set_obj_coef(lp, c, 0)
        end
    end

    m > 0 && n > 0 && GLPK.load_matrix(lp, ia, ja, ar)
    MPB.setsense!(lpm, sense)

    return lpm
end

#writeproblem(m, filename::AbstractString)

function MPB.getvarLB(lpm::GLPKMathProgModel)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    lb = Array{Float64}(undef, n)
    for c = 1:n
        l = GLPK.get_col_lb(lp, c)
        if l <= -floatmax(Float64)
            l = -Inf
        end
        lb[c] = l
    end
    return lb
end

function MPB.setvarLB!(lpm::GLPKMathProgModel, collb)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    if nonnull(collb) && length(collb) != n
        error("invalid size of collb")
    end
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for c = 1:n
        u = GLPK.get_col_ub(lp, c)
        if u >= floatmax(Float64)
            u = Inf
        end
        if nonnull(collb) && collb[c] != -Inf
            l = collb[c]
            if u < Inf
                if l != u
                    GLPK.set_col_bnds(lp, c, GLPK.DB, l, u)
                else
                    GLPK.set_col_bnds(lp, c, GLPK.FX, l, u)
                end
            else
                GLPK.set_col_bnds(lp, c, GLPK.LO, l, 0.0)
            end
        else
            if u < Inf
                GLPK.set_col_bnds(lp, c, GLPK.UP, 0.0, u)
            else
                GLPK.set_col_bnds(lp, c, GLPK.FR, 0.0, 0.0)
            end
        end
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)
end

function MPB.getvarUB(lpm::GLPKMathProgModel)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    ub = Array{Float64}(undef, n)
    for c = 1:n
        u = GLPK.get_col_ub(lp, c)
        if u >= floatmax(Float64)
            u = Inf
        end
        ub[c] = u
    end
    return ub
end

function MPB.setvarUB!(lpm::GLPKMathProgModel, colub)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    if nonnull(colub) && length(colub) != n
        error("invalid size of colub")
    end
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for c = 1:n
        l = GLPK.get_col_lb(lp, c)
        if l <= -floatmax(Float64)
            l = -Inf
        end
        if nonnull(colub) && colub[c] != Inf
            u = colub[c]
            if l > -Inf
                if l != u
                    GLPK.set_col_bnds(lp, c, GLPK.DB, l, u)
                else
                    GLPK.set_col_bnds(lp, c, GLPK.FX, l, u)
                end
            else
                GLPK.set_col_bnds(lp, c, GLPK.UP, 0.0, u)
            end
        else
            if l > -Inf
                GLPK.set_col_bnds(lp, c, GLPK.LO, l, 0.0)
            else
                GLPK.set_col_bnds(lp, c, GLPK.FR, 0.0, 0.0)
            end
        end
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)
end

function MPB.getconstrLB(lpm::GLPKMathProgModel)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)
    lb = Array{Float64}(undef, m)
    for r = 1:m
        l = GLPK.get_row_lb(lp, r)
        if l <= -floatmax(Float64)
            l = -Inf
        end
        lb[r] = l
    end
    return lb
end

function MPB.setconstrLB!(lpm::GLPKMathProgModel, rowlb)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)
    if nonnull(rowlb) && length(rowlb) != m
        error("invalid size of rowlb")
    end
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for r = 1:m
        u = GLPK.get_row_ub(lp, r)
        if u >= floatmax(Float64)
            u = Inf
        end
        if nonnull(rowlb) && rowlb[r] != -Inf
            l = rowlb[r]
            if u < Inf
                if l != u
                    GLPK.set_row_bnds(lp, r, GLPK.DB, l, u)
                else
                    GLPK.set_row_bnds(lp, r, GLPK.FX, l, u)
                end
            else
                GLPK.set_row_bnds(lp, r, GLPK.LO, l, 0.0)
            end
        else
            if u < Inf
                GLPK.set_row_bnds(lp, r, GLPK.UP, 0.0, u)
            else
                GLPK.set_row_bnds(lp, r, GLPK.FR, 0.0, 0.0)
            end
        end
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)
end

function MPB.getconstrUB(lpm::GLPKMathProgModel)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)
    ub = zeros(m)
    for r = 1:m
        u = GLPK.get_row_ub(lp, r)
        if u >= floatmax(Float64)
            u = Inf
        end
        ub[r] = u
    end
    return ub
end

function MPB.setconstrUB!(lpm::GLPKMathProgModel, rowub)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)
    if nonnull(rowub) && length(rowub) != m
        error("invalid size of rowub")
    end
    prev_preemptive_check = GLPK.jl_get_preemptive_check()
    GLPK.jl_set_preemptive_check(false)
    for r = 1:m
        l = GLPK.get_row_lb(lp, r)
        if l <= -floatmax(Float64)
            l = -Inf
        end
        if nonnull(rowub) && rowub[r] != Inf
            u = rowub[r]
            if l > -Inf
                if l != u
                    GLPK.set_row_bnds(lp, r, GLPK.DB, l, u)
                else
                    GLPK.set_row_bnds(lp, r, GLPK.FX, l, u)
                end
            else
                GLPK.set_row_bnds(lp, r, GLPK.UP, 0.0, u)
            end
        else
            if l > -Inf
                GLPK.set_row_bnds(lp, r, GLPK.LO, l, 0.0)
            else
                GLPK.set_row_bnds(lp, r, GLPK.FR, 0.0, 0.0)
            end
        end
    end
    GLPK.jl_set_preemptive_check(prev_preemptive_check)
end

function MPB.getconstrmatrix(lpm::GLPKMathProgModel)
    lp = lpm.inner
    m = GLPK.get_num_rows(lp)
    n = GLPK.get_num_cols(lp)
    colwise = [GLPK.get_mat_col(lp,i) for i in 1:n]
    nnz = 0
    for i in 1:n
        nnz += length(colwise[i][1])
    end
    colptr = Array{Int}(undef, n+1)
    rowval = Array{Int}(undef, nnz)
    nzval  = Array{Float64}(undef, nnz)
    cur_nnz = 1
    for i in 1:n
        colptr[i] = cur_nnz
        ind,vals = colwise[i]
        p = sortperm(ind) # indices must be sorted
        rowval[cur_nnz:(cur_nnz+length(ind)-1)] = ind[p]
        nzval[cur_nnz:(cur_nnz+length(ind)-1)] = vals[p]
        cur_nnz += length(ind)
    end
    colptr[n+1] = cur_nnz
    return SparseMatrixCSC(m,n,colptr,rowval,nzval)
end

function MPB.getobj(lpm::GLPKMathProgModel)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    return [GLPK.get_obj_coef(lp, i) for i in 1:n]
end

function MPB.setobj!(lpm::GLPKMathProgModel, obj)
    lp = lpm.inner
    n = GLPK.get_num_cols(lp)
    if nonnull(obj) && length(obj) != n
        error("invalid size of obj")
    end
    for c = 1:n
        if nonnull(obj)
            GLPK.set_obj_coef(lp, c, obj[c])
        else
            GLPK.set_obj_coef(lp, c, 0.0)
        end
    end
end

function MPB.addvar!(lpm::GLPKMathProgModel, rowidx::Vector, rowcoef::Vector, collb::Real, colub::Real, objcoef::Real)
    if length(rowidx) != length(rowcoef)
        error("rowidx and rowcoef have different legths")
    end
    lp = lpm.inner
    GLPK.add_cols(lp, 1)
    n = GLPK.get_num_cols(lp)
    GLPK.set_mat_col(lp, n, rowidx, rowcoef)
    if collb > -Inf && colub < Inf
        if collb != colub
            bt = GLPK.DB
        else
            bt = GLPK.FX
        end
    elseif collb > -Inf
        bt = GLPK.LO
    elseif colub < Inf
        bt = GLPK.UP
    else
        bt = GLPK.FR
    end
    GLPK.set_col_bnds(lp, n, bt, collb, colub)
    GLPK.set_obj_coef(lp, n, objcoef)
    return
end

function MPB.delvars!(lpm::GLPKMathProgModel, idx::Vector)
    GLPK.std_basis(lpm.inner)
    GLPK.del_cols(lpm.inner, length(idx), idx)
end

function MPB.addconstr!(lpm::GLPKMathProgModel, colidx::Vector, colcoef::Vector, rowlb::Real, rowub::Real)
    if length(colidx) != length(colcoef)
        error("colidx and colcoef have different legths")
    end
    lp = lpm.inner
    GLPK.add_rows(lp, 1)
    m = GLPK.get_num_rows(lp)
    GLPK.set_mat_row(lp, m, colidx, colcoef)
    if rowlb > -Inf && rowub < Inf
        if rowlb != rowub
            bt = GLPK.DB
        else
            bt = GLPK.FX
        end
    elseif rowlb > -Inf
        bt = GLPK.LO
    elseif rowub < Inf
        bt = GLPK.UP
    else
        bt = GLPK.FR
    end
    GLPK.set_row_bnds(lp, m, bt, rowlb, rowub)
    return
end

function MPB.delconstrs!(lpm::GLPKMathProgModel, idx::Vector)
    GLPK.std_basis(lpm.inner)
    GLPK.del_rows(lpm.inner, length(idx), idx)
end

function MPB.setsense!(lpm::GLPKMathProgModel, sense)
    lp = lpm.inner
    if sense == :Min
        GLPK.set_obj_dir(lp, GLPK.MIN)
    elseif sense == :Max
        GLPK.set_obj_dir(lp, GLPK.MAX)
    else
        error("Unrecognized objective sense $sense")
    end
end

function MPB.getsense(lpm::GLPKMathProgModel)
    lp = lpm.inner
    s = GLPK.get_obj_dir(lp)
    if s == GLPK.MIN
        return :Min
    elseif s == GLPK.MAX
        return :Max
    else
        error("Internal library error")
    end
end

MPB.numvar(lpm::GLPKMathProgModel) = GLPK.get_num_cols(lpm.inner)
MPB.numconstr(lpm::GLPKMathProgModel) = GLPK.get_num_rows(lpm.inner)

MPB.getrawsolver(lpm::GLPKMathProgModel) = lpm.inner

end


text:
using Test
import DiffOpt
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
import DelimitedFiles
import Ipopt
import HiGHS
import SCS

const VAF = MOI.VectorAffineFunction{Float64}
_vaf(c::Vector{Float64}) = VAF(MOI.ScalarAffineTerm{Float64}[], c)

@testset "MOI Unit" begin
    function test_runtests()
        model = DiffOpt.diff_optimizer(HiGHS.Optimizer)
        # `Variable.ZerosBridge` makes dual needed by some tests fail.
        MOI.Bridges.remove_bridge(model.optimizer.optimizer, MOI.Bridges.Variable.ZerosBridge{Float64})
        MOI.set(model, MOI.Silent(), true)
        config = MOI.Test.Config(atol = 1e-7)
        MOI.Test.runtests(model, config),
        return
    end
    test_runtests()
end

@testset "Testing forward on trivial QP" begin
    # using example on https://osqp.org/docs/examples/setup-and-solve.html
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [
         1.0 1.0
         1.0 0.0
         0.0 1.0
        -1.0 -1.0
        -1.0 0.0
        0.0 -1.0
    ]
    h = [1, 0.7, 0.7, -1, 0, 0]
    qp_test_with_solutions(
        Ipopt.Optimizer;
        Q = Q,
        q = q,
        G = G,
        h = h,
        dzb = ones(2),
        dQf = [1 -1; -1 1.0],
        dqf = [1, -1.0],
        dGf = ones(6, 2),
        dhf = ones(6),
        # Expected solutions
        z = [0.3, 0.7],
    )
end

@testset "Differentiating trivial QP 1" begin
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [1.0 1.0]
    h = [-1.0]

    model = DiffOpt.diff_optimizer(Ipopt.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)

    qp_test_with_solutions(
        Ipopt.Optimizer;
        Q = Q,
        q = q,
        G = G,
        h = h,
        dzb = ones(2),
        dQf = -ones(2, 2),
        dqf = ones(2),
        dGf = ones(1, 2),
        dhf = -ones(1),
        # Expected solutions
        z = [-0.25; -0.75],
        dhb = ones(1),
    )
end

# @testset "Differentiating a non-convex QP" begin
#     Q = [0.0 0.0; 1.0 2.0]
#     q = [1.0; 1.0]
#     G = [1.0 1.0;]
#     h = [-1.0;]

#     model = DiffOpt.diff_optimizer(Ipopt.Optimizer)
#     x = MOI.add_variables(model, 2)

#     # define objective
#     quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
#     for i in 1:2
#         for j in i:2 # indexes (i,j), (j,i) will be mirrored. specify only one kind
#             push!(
#                 quad_terms,
#                 MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j]),
#             )
#         end
#     end

#     objective_function = MOI.ScalarQuadraticFunction(
#                             quad_terms,
#                             MOI.ScalarAffineTerm.(q, x),
#                             0.0,
#                         )
#     MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
#     MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

#     # add constraint
#     MOI.add_constraint(
#         model,
#         MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[1, :], x), 0.0),
#         MOI.LessThan(h[1]),
#     )

#     @test_throws ErrorException MOI.optimize!(model) # should break
# end

@testset "Differentiating QP with inequality and equality constraints" begin
    # refered from: https://www.mathworks.com/help/optim/ug/quadprog.html#d120e113424
    # Find equivalent qpth program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_4_py.ipynb

    Q = [
         1.0 -1.0 1.0;
        -1.0  2.0 -2.0;
         1.0 -2.0 4.0
    ]
    q = [2.0, -3.0, 1.0]
    G = [0.0 0.0 1.0
         0.0 1.0 0.0
         1.0 0.0 0.0
         0.0 0.0 -1.0
         0.0 -1.0 0.0
         -1.0 0.0 0.0
    ]
    h = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0,]
    A = [1.0 1.0 1.0;]
    b = [0.5]

    qp_test_with_solutions(
        Ipopt.Optimizer;
        Q = Q,
        q = q,
        G = G,
        h = h,
        A = A,
        b = b,
        dzb = ones(3),
        dQf = ones(3, 3),
        dqf = ones(3),
        dGf = ones(6, 3),
        dhf = ones(6),
        dAf = ones(1, 3),
        dbf = ones(1),
        # Expected solutions
        z = [0.0, 0.5, 0.0],
        dQb = zeros(3, 3),
        dqb = zeros(3),
        dGb = zeros(6, 3),
        dhb = zeros(6),
        dAb = [0.0 -0.5 0.0],
        dbb = [1.0],
    )
end

# refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contquadratic.jl#L3
# Find equivalent CVXPYLayers and QPTH code here:
#               https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_1_py.ipynb
@testset "Differentiating MOI examples 1" begin
    # homogeneous quadratic objective
    # Min x^2 + xy + y^2 + yz + z^2
    # st  x + 2y + 3z >= 4 (c1)
    #     x +  y      >= 1 (c2)
    #     x, y, z \in R

    Q = [
        2.0 1.0 0.0
        1.0 2.0 1.0
        0.0 1.0 2.0
    ]
    q = zeros(3)
    G = [
        -1.0 -2.0 -3.0
        -1.0 -1.0 0.0
    ]
    h = [-4.0, -1.0]

    dQ = [-0.12244895  0.01530609 -0.11224488
           0.01530609  0.09183674  0.07653058
          -0.11224488  0.07653058 -0.06122449]
    dq = [-0.2142857,  0.21428567, -0.07142857]

    dG = [0.05102692 0.30612244 0.25510856
          0.06120519 0.36734693 0.30610315]
    dh = [-0.35714284; -0.4285714]

    qp_test_with_solutions(
        Ipopt.Optimizer;
        Q = Q,
        q = q,
        G = G,
        h = h,
        dzb = ones(3),
        dQf = dQ,
        dqf = dq,
        dGf = dG,
        dhf = dh,
        # Expected solutions
        dQb = dQ,
        dqb = dq,
        dGb = dG,
        dhb = dh,
    )
end

# refered from https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contquadratic.jl#L3
# Find equivalent CVXPYLayers and QPTH code here:
#               https://github.com/AKS1996/jump-gsoc-2020/blob/master/DiffOpt_tests_2_py.ipynb
@testset "Differentiating MOI examples 2" begin
    # non-homogeneous quadratic objective
    #    minimize 2 x^2 + y^2 + xy + x + y
    #       s.t.  x, y >= 0
    #             x + y = 1

    Q = [4 1.0
         1 2]
    q = [1, 1.0]
    G = [-1 0.0
         0 -1]
    h = [0, 0.0]
    A = [1 1.0]
    b = [1.0]

    dQ = [-0.05   -0.05
          -0.05    0.15]
    dq = [-0.2, 0.2]
    dG = zeros(2, 2)
    dh = zeros(2)
    dA = [0.375 -1.075]
    db = [0.7]

    qp_test_with_solutions(
        Ipopt.Optimizer;
        Q = Q,
        q = q,
        G = G,
        h = h,
        A = A,
        b = b,
        dzb = [1.3, 0.5],
        dQf = dQ,
        dqf = dq,
        dGf = dG,
        dhf = dh,
        dAf = dA,
        dbf = db,
        # Expected solutions
        dQb = dQ,
        dqb = dq,
        dGb = dG,
        dhb = dh,
        dAb = dA,
        dbb = db,
        z = [0.25, 0.75],
        dzf = [1.4875, -0.075],
        ∇zf = [-1.28125, 3.25625],
        ∇zb = [-0.2, 0.2],
        λ = zeros(2),
        dλf = zeros(2),
        dλb = zeros(2),
        ∇λb = [0.8, -0.8/3],
        ν = [-2.75],
        dνb = zeros(1),
        ∇νb = [-0.7],
    )
end

@testset "Differentiating non trivial convex QP MOI" begin
    nz = 10
    nineq_le = 25
    neq = 10

    # read matrices from files
    names = ["P", "q", "G", "h", "A", "b"]
    matrices = []

    for name in names
        push!(matrices, DelimitedFiles.readdlm(joinpath(dirname(dirname(pathof(DiffOpt))), "test", "data", name * ".txt"), ' ', Float64, '\n'))
    end

    Q, q, G, h, A, b = matrices
    q = vec(q)
    h = vec(h)
    b = vec(b)

    # read gradients from files
    names = ["dP", "dq", "dG", "dh", "dA", "db"]
    grads_actual = []

    for name in names
        push!(grads_actual, DelimitedFiles.readdlm(joinpath(@__DIR__, "data", name * ".txt"), ' ', Float64, '\n'))
    end

    dqb = vec(grads_actual[2])
    dhb = vec(grads_actual[4])
    dbb = vec(grads_actual[6])

    qp_test(
        Ipopt.Optimizer, true, true, true;
        Q = Q,
        q = q,
        G = G,
        h = h,
        A = A,
        b = b,
        dzb = ones(nz),
        dQf = ones(nz, nz),
        dqf = ones(nz),
        dGf = ones(length(h), nz),
        dhf = ones(length(h)),
        dAf = ones(length(b), nz),
        dbf = ones(length(b)),
        # Expected solutions
        dqb = dqb,
        dhb = dhb,
        dbb = dbb,
        atol = 1e-3, # The values in `data` seems to have low accuracy
        rtol = 1e-3, # The values in `data` seems to have low accuracy
    )
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # Issue #40 from Gurobi.jl
    # min  x
    # s.t. x >= 0
    #      x >= 3
    nz = 1
    qp_test_with_solutions(
        HiGHS.Optimizer;
        q = ones(nz),
        G = -ones(2, nz),
        h = [0.0, -3.0],
        dzb = ones(nz),
        dqf = ones(nz),
        # Expected solutions
        dGb = [0.0, 3.0],
        dhb = [0.0, -1.0],
    )
end

@testset "Differentiating a simple LP with GreaterThan constraint" begin
    # this is canonically same as above test
    # min  x
    # s.t. x >= 3
    nz = 1
    qp_test_with_solutions(
        Ipopt.Optimizer;
        q = ones(nz),
        G = -ones(1, nz),
        h = [-3.0],
        dzb = ones(nz),
        dqf = ones(nz),
        # Expected solutions
        dGb = [3.0],
        dhb = [-1.0],
    )
end

@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # refered from - https://en.wikipedia.org/wiki/Simplex_algorithm#Example

    # max 2x + 3y + 4z
    # s.t. 3x+2y+z <= 10
    #      2x+5y+3z <= 15
    #      x,y,z >= 0
    nz = 3
    qp_test_with_solutions(
        SCS.Optimizer;
        q = [-2.0, -3.0, -4.0],
        G = [
            3.0  2.0  1.0
            2.0  5.0  3.0
            -1.0  0.0  0.0
            0.0 -1.0  0.0
            0.0  0.0 -1.0
        ],
        h = [10.0, 15.0, 0.0, 0.0, 0.0],
        dzb = ones(nz),
        dqf = ones(nz),
        # Expected solutions
        dqb = zeros(nz),
        dGb = [
            0.0 0.0 0.0
            0.0 0.0 -5/3
            0.0 0.0 5/3
            0.0 0.0 -10/3
            0.0 0.0 0.0
        ],
        dhb = [0.0, 1/3, -1/3, 2/3, 0.0],
    )
end

@testset "Differentiating LP with variable bounds" begin
    # max 2x + 3y + 4z
    # s.t. 3x+2y+z <= 10
    #      2x+5y+3z <= 15
    #      x ≤ 3
    #      0 ≤ y ≤ 2
    #      z ≥ 2
    #      x,y,z >= 0
    # variant of previous test with same solution
    nz = 3
    qp_test_with_solutions(
        HiGHS.Optimizer;
        q = [-2.0, -3.0, -4.0],
        G = [
            3.0  2.0  1.0
            2.0  5.0  3.0
           -1.0  0.0  0.0
            0.0 -1.0  0.0
            0.0  0.0 -1.0
            1.0  0.0  0.0
            0.0  1.0  0.0
            0.0  0.0  1.0
        ],
        h = [10.0, 15.0, 0.0, 0.0, 0.0, 3.0, 2.0, 6.0],
        dzb = ones(nz),
        dqf = ones(nz),
        # Expected solutions
        dqb = zeros(nz),
        dGb = [
            0.0 0.0 0.0
            0.0 0.0 -5/3
            0.0 0.0 5/3
            0.0 0.0 -10/3
            0.0 0.0 0.0
            0.0 0.0 0.0
            0.0 0.0 0.0
            0.0 0.0 0.0
        ],
        dhb = [0.0, 1/3, -1/3, 2/3, 0.0, 0.0, 0.0, 0.0],
    )
end

@testset "Differentiating LP with variable bounds 2" begin
    nz = 3
    qp_test_with_solutions(
        HiGHS.Optimizer;
        q = [-2.0, -3.0, -4.0],
        G = [
            3.0  2.0  1.0
            2.0  5.0  3.0
            0.0 -1.0  0.0
            0.0  0.0 -1.0
        ],
        h = [10.0, 15.0, 0.0, 0.0],
        fix_indices = [1],
        fix_values = [0.0],
        dzb = ones(nz),
        dqf = ones(nz),
        # Expected solutions
        dqb = zeros(nz),
        dGb = [0.0 0.0 0.0
               0.0 0.0 -5/3
               0.0 0.0 -10/3
               0.0 0.0 0.0],
        dhb = [0.0, 1/3, 2/3, 0.0],
        dAb = [0.0 0.0 -5/3],
        dbb = [1/3],
    )
end

@testset "Differentiating LP with SAF, SV with LE, GE constraints" begin
    """
        max  2x + 3y + 4z
        s.t. 3x+2y+z <= 10
             2x+5y+3z <= 15
             x ≤ 3
             0 ≤ y ≤ 2
             6 ≥ z ≥ -1
             x, y, z >= 0
        variant of previous test with same solution
    """
    nz = 3
    qp_test_with_solutions(
        HiGHS.Optimizer;
        q = [-2.0, -3.0, -4.0],
        G = [
            3.0  2.0  1.0
            2.0  5.0  3.0
           -1.0  0.0  0.0
            0.0 -1.0  0.0
            0.0  0.0 -1.0
        ],
        h = [10.0, 15.0, 0.0, 0.0, 0.0], #5
        ub_indices = [1, 2, 3],
        ub_values = [3.0, 2.0, 6.0],
        lb_indices = [1],
        lb_values = [-1.0],
        dzb = ones(nz),
        dqb = zeros(nz),
        dGb = [0.0 0.0 0.0;
               0.0 0.0 -5/3;
               0.0 0.0 5/3;
               0.0 0.0 -10/3;
               0.0 0.0 0.0
               0.0 0.0 0.0
               0.0 0.0 0.0
               0.0 0.0 0.0
               0.0 0.0 0.0],
        dhb = [0.0, 1/3, -1/3, 2/3, 0.0, 0.0, 0.0, 0.0, 0.0],
    )
end

# TODO: split file here
# above is QP Back
# below is conic forw

function test_simple_socp(eq_vec::Bool)
    # referred from _soc2test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L1355
    # find equivalent diffcp python program here: https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_socp_1_py.ipynb

    # Problem SOC2
    # min  x
    # s.t. y ≥ 1/√2
    #      x² + y² ≤ 1
    # in conic form:
    # min  x
    # s.t.  -1/√2 + y ∈ R₊
    #        1 - t ∈ {0}
    #      (t,x,y) ∈ SOC₃

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x,y,t = MOI.add_variables(model, 3)

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), 1.0x)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    if eq_vec
        ceq  = MOI.add_constraint(model, MOIU.vectorize([-1.0t + 1.0]), MOI.Zeros(1))
    else
        ceq  = MOI.add_constraint(model, -1.0t, MOI.EqualTo(-1.0))
    end
    cnon = MOI.add_constraint(model, MOIU.vectorize([1.0y - 1/√2]), MOI.Nonnegatives(1))
    csoc = MOI.add_constraint(model, MOIU.vectorize([1.0t, 1.0x, 1.0y]), MOI.SecondOrderCone(3))

    MOI.optimize!(model)

    if eq_vec
        MOI.set(model, DiffOpt.ForwardConstraintFunction(), ceq, MOIU.vectorize([1.0 * x]))
    else
        MOI.set(model, DiffOpt.ForwardConstraintFunction(), ceq, 1.0 * x)
    end
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), cnon, MOIU.vectorize([1.0 * y]))
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), csoc, MOIU.operate(vcat, Float64, 1.0 * t, 0.0, 0.0))

    DiffOpt.forward_differentiate!(model)

    # these matrices are benchmarked with the output generated by diffcp
    # refer the python file mentioned above to get equivalent python source code
    @test model.diff.model.x ≈ [-1/√2; 1/√2; 1.0] atol=ATOL rtol=RTOL
    if eq_vec
        @test model.diff.model.s ≈ [0.0, 0.0, 1.0, -1/√2, 1/√2] atol=ATOL rtol=RTOL
        @test model.diff.model.y ≈ [√2, 1.0, √2, 1.0, -1.0] atol=ATOL rtol=RTOL
    else
        @test model.diff.model.s ≈ [0.0, 1.0, -1/√2, 1/√2, 0.0] atol=ATOL rtol=RTOL
        @test model.diff.model.y ≈ [1.0, √2, 1.0, -1.0, √2] atol=ATOL rtol=RTOL
    end

    dx = [1.12132144; 1/√2; 1/√2]
    for (i, vi) in enumerate([x, y, t])
        @test dx[i] ≈ MOI.get(model, DiffOpt.ForwardVariablePrimal(), vi) atol=ATOL rtol=RTOL
    end
    # @test dx ≈ [1.12132144; 1/√2; 1/√2] atol=ATOL rtol=RTOL
    # @test ds ≈ [0.0; 0.0; -2.92893438e-01;  1.12132144e+00; 7.07106999e-01]  atol=ATOL rtol=RTOL
    # @test dy ≈ [2.4142175; 5.00000557; 3.8284315; √2; -4.00000495] atol=ATOL rtol=RTOL
end

@testset "Differentiating simple SOCP" begin
    for eq_vec in [false, true]
        test_simple_socp(eq_vec)
    end
end

# refered from _psd0test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L3919
# find equivalent diffcp program here: https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_1_py.ipynb

# min X[1,1] + X[2,2]    max y
#     X[2,1] = 1         [0   y/2     [ 1  0
#                         y/2 0    <=   0  1]
#     X >= 0              y free
# Optimal solution:
#
#     ⎛ 1   1 ⎞
# X = ⎜       ⎟           y = 2
#     ⎝ 1   1 ⎠
function simple_psd(solver)
    model = DiffOpt.diff_optimizer(solver)
    MOI.set(model, MOI.Silent(), true)
    X = MOI.add_variables(model, 3)
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(
        model,
        VAF(vov),
        MOI.PositiveSemidefiniteConeTriangle(2)
    )

    c  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, X[2]))],
            [-1.0]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            1.0 * X[1] + 1.0 * X[end])
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    x = MOI.get(model, MOI.VariablePrimal(), X)

    cone_types = unique([S for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())])
    conic_form = DiffOpt.ConicProgram.Form{Float64}()
    cones = conic_form.constraints.sets
    DiffOpt.set_set_types(cones, cone_types)
    index_map = MOI.copy_to(conic_form, model)

    # s = DiffOpt.map_rows((ci, r) -> MOI.get(model.optimizer, MOI.ConstraintPrimal(), ci), model.optimizer, cones, index_map, DiffOpt.Flattened{Float64}())
    # y = DiffOpt.map_rows((ci, r) -> MOI.get(model.optimizer, MOI.ConstraintDual(), ci), model.optimizer, cones, index_map, DiffOpt.Flattened{Float64}())

    @test x ≈ ones(3) atol=ATOL rtol=RTOL
    # @test s ≈ [0.0; ones(3)] atol=ATOL rtol=RTOL
    # @test y ≈ [2.0, 1.0, -1.0, 1.0]  atol=ATOL rtol=RTOL

    # test1: changing the constant in `c`, i.e. changing value of X[2]
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, _vaf([1.0]))

    DiffOpt.forward_differentiate!(model)

    dx = -ones(3)
    for (i, vi) in enumerate(X)
        @test dx[i] ≈ MOI.get(model, DiffOpt.ForwardVariablePrimal(), vi) atol=ATOL rtol=RTOL
    end

    # @test dx ≈ -ones(3) atol=ATOL rtol=RTOL  # will change the value of other 2 variables
    # @test ds[2:4] ≈ -ones(3)  atol=ATOL rtol=RTOL  # will affect PSD constraint too

    # test2: changing X[1], X[3] but keeping the objective (their sum) same
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, MOIU.zero_with_output_dimension(MOI.VectorAffineFunction{Float64}, 1))
    MOI.set(model, DiffOpt.ForwardObjectiveFunction(), -1.0X[1] + 1.0X[3])

    DiffOpt.forward_differentiate!(model)

    # @test dx ≈ [1.0, 0.0, -1.0] atol=ATOL rtol=RTOL  # note: no effect on X[2]
    dx = [1.0, 0.0, -1.0]
    for (i, vi) in enumerate(X)
        @test dx[i] ≈ MOI.get(model, DiffOpt.ForwardVariablePrimal(), vi) atol=ATOL rtol=RTOL
    end
end

@testset "Differentiating simple PSD program" begin
    simple_psd(SCS.Optimizer)
end

@testset "Differentiating conic with PSD and SOC constraints" begin
    # similar to _psd1test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L4054
    # find equivalent diffcp example here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_2_py.ipynb

    #     | 2 1 0 |
    # min | 1 2 1 | . X + x1
    #     | 0 1 2 |
    #
    #
    # s.t. | 1 0 0 |
    #      | 0 1 0 | . X + x1 = 1
    #      | 0 0 1 |
    #
    #      | 1 1 1 |
    #      | 1 1 1 | . X + x2 + x3 = 1/2
    #      | 1 1 1 |
    #
    #      (x1,x2,x3) in C^3_q
    #      X in C_psd
    #
    # The dual is
    # max y1 + y2/2
    #
    # s.t. | y1+y2    y2    y2 |
    #      |    y2 y1+y2    y2 | in C_psd
    #      |    y2    y2 y1+y2 |
    #
    #      (1-y1, -y2, -y2) in C^3_q

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    y2 = 1 - ε*δ
    y1 = 1 - √2*y2
    obj = y1 + y2/2
    k = -2*δ/ε
    x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    α = √(3-2obj-4x2)/2
    β = k*α

    X = MOI.add_variables(model, 6)
    x = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(model,
        VAF(vov),
        MOI.PositiveSemidefiniteConeTriangle(3))
    cx = MOI.add_constraint(model,
        VAF(MOI.VectorOfVariables(x)),
        MOI.SecondOrderCone(3))

    c1 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 1., 1., 1.], [X[1], X[3], X[end], x[1]])),
            [-1.0]
        ),
        MOI.Zeros(1)
    )
    c2 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]])),
            [-0.5]
        ),
        MOI.Zeros(1),
    )

    # this is a useless constraint - refer the tests below
    # even if we comment this, it won't affect the optimal values
    c_extra = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.(ones(3), x)),
            [100.0]
        ),
        MOI.Nonnegatives(1)
    )

    objXidx = [1:3; 5:6]
    objXcoefs = 2*ones(5)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    _x = MOI.get(model, MOI.VariablePrimal(), x)
    _X = MOI.get(model, MOI.VariablePrimal(), X)

    cone_types = unique([S for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())])
    conic_form = DiffOpt.ConicProgram.Form{Float64}()
    cones = conic_form.constraints.sets
    DiffOpt.set_set_types(cones, cone_types)
    index_map = MOI.copy_to(conic_form, model)

    # s = DiffOpt.map_rows((ci, r) -> MOI.get(model.optimizer, MOI.ConstraintPrimal(), ci), model.optimizer, cones, index_map, DiffOpt.Flattened{Float64}())
    # y = DiffOpt.map_rows((ci, r) -> MOI.get(model.optimizer, MOI.ConstraintDual(), ci), model.optimizer, cones, index_map, DiffOpt.Flattened{Float64}())

    @test _X ≈ [0.21725121; -0.25996907;  0.31108582;  0.21725009; -0.25996907;  0.21725121
        ] atol=ATOL rtol=RTOL
    @test _x ≈ [0.2544097;   0.17989425;  0.17989425
        ] atol=ATOL rtol=RTOL

    # @test x ≈ [ 0.21725121; -0.25996907;  0.31108582;  0.21725009; -0.25996907;  0.21725121;
    #             0.2544097;   0.17989425;  0.17989425] atol=ATOL rtol=RTOL
    # @test s ≈ [
    #     0.0, 0.0, 100.614,
    #     0.254408, 0.179894, 0.179894,
    #     0.217251, -0.25997, 0.31109, 0.217251, -0.25997, 0.217251,
    # ] atol=ATOL rtol=RTOL
    # TODO: it should be 100, not 100.614, its surely a residual error
    # Joaquim: it is not an error it is sum(_x)
    # there seemss to be an issue with this test (note the low precision)
    #
    # @test y ≈ [
    #     0.544758, 0.321905, 0.0,
    #     0.455242, -0.321905, -0.321905,
    #     1.13334, 0.678095, 1.13334, -0.321905, 0.678095, 1.13334,
    # ]  atol=ATOL rtol=RTOL

    # test c_extra
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c_extra, _vaf([1.0]))

    DiffOpt.forward_differentiate!(model)

    # a small change in the constant in c_extra should not affect any other variable or constraint other than c_extra itself
    for (i, vi) in enumerate(X)
        @test 0.0 ≈ MOI.get(model,
            DiffOpt.ForwardVariablePrimal(), vi) atol=1e-2 rtol=RTOL
    end
    for (i, vi) in enumerate(x)
        @test 0.0 ≈ MOI.get(model,
            DiffOpt.ForwardVariablePrimal(), vi) atol=1e-2 rtol=RTOL
    end
    # @test dx ≈ zeros(9) atol=1e-2
    # @test dy ≈ zeros(12) atol=0.012
    # @test [ds[1:2]; ds[4:end]] ≈ zeros(11) atol=1e-2
    # @test ds[3] ≈ 1.0 atol=1e-2   # except c_extra itself
end

@testset "Differentiating conic with PSD and POS constraints" begin
    # refer psdt2test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L4306
    # find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variables(model, 7)
    @test MOI.get(model, MOI.NumberOfVariables()) == 7

    η = 10.0

    c1  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(-1.0, x[1:6])),
            [η]
        ),
        MOI.Nonnegatives(1)
    )
    c2 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                1:6, MOI.ScalarAffineTerm.(1.0, x[1:6])), zeros(6)),
            MOI.Nonnegatives(6))
    α = 0.8
    δ = 0.9
    c3 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [fill(1, 7); fill(2, 5);     fill(3, 6)],
                MOI.ScalarAffineTerm.(
                    [ δ/2,       α,   δ, δ/4, δ/8,      0.0, -1.0,
                        -δ/(2*√2), -δ/4, 0,     -δ/(8*√2), 0.0,
                        δ/2,     δ-α,   0,      δ/8,      δ/4, -1.0],
                    [x[1:7];     x[1:3]; x[5:6]; x[1:3]; x[5:7]])),
                zeros(3)), MOI.PositiveSemidefiniteConeTriangle(2))
    c4 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1,
                MOI.ScalarAffineTerm.(0.0, [x[1:3]; x[5:6]])),
            [0.0]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x[7])], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(model)

    # dc = ones(7)
    MOI.set(model, DiffOpt.ForwardObjectiveFunction(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(7), x), 0.0))
    # db = ones(11)
    # dA = ones(11, 7)
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c1, MOIU.vectorize(ones(1, 7) * x + ones(1)))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c2, MOIU.vectorize(ones(6, 7) * x + ones(6)))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c3, MOIU.vectorize(ones(3, 7) * x + ones(3)))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c4, MOIU.vectorize(ones(1, 7) * x + ones(1)))

    DiffOpt.forward_differentiate!(model)

    @test model.diff.model.x ≈ [20/3., 0.0, 10/3., 0.0, 0.0, 0.0, 1.90192379] atol=ATOL rtol=RTOL
    @test model.diff.model.s ≈ [0.0, 0.0, 20/3.0,  0.0,  10/3.0,  0.0,  0.0,  0.0, 4.09807621, -2.12132,  1.09807621] atol=ATOL rtol=RTOL
    @test model.diff.model.y ≈ [0.0, 0.19019238, 0., 0.12597667, 0., 0.14264428, 0.14264428, 0.01274047, 0.21132487, 0.408248, 0.78867513] atol=ATOL rtol=RTOL

    atol = 0.3
    rtol = 0.01

    # compare these with https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb
    # results are not exactly as: 1. there is some residual error   2. diffcp results are SCS specific, hence scaled
    dx = [-39.6066, 10.8953, -14.9189, 10.9054, 10.883, 10.9118, -21.7508]
    for (i, vi) in enumerate(x)
        @test dx[i] ≈ MOI.get(model,
            DiffOpt.ForwardVariablePrimal(), vi) atol=atol rtol=rtol
    end
    # @test dy ≈ [0.0, -3.56905, 0.0, -0.380035, 0.0, -0.41398, -0.385321, -0.00743119, -0.644986, -0.550542, -2.36765] atol=atol rtol=rtol
    # @test ds ≈ [0.0, 0.0, -50.4973, 0.0, -25.8066, 0.0, 0.0, 0.0, -7.96528, -1.62968, -2.18925] atol=atol rtol=rtol

    # TODO: future example, how to differentiate wrt a specific constraint/variable, refer QPLib article for more
    dA = zeros(11, 7)
    dA[3:8, 1:6] = Matrix{Float64}(LinearAlgebra.I, 6, 6)  # differentiating only wrt POS constraint c2
    db = zeros(11)
    dc = zeros(7)

    # db = zeros(11)
    # dA = zeros(11, 7)
    # dA[3:8, 1:6] = Matrix{Float64}(LinearAlgebra.I, 6, 6)  # differentiating only wrt POS constraint c2
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c1, MOIU.zero_with_output_dimension(VAF, 1))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c2, MOIU.vectorize(ones(6) .* x[1:6]))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c3, MOIU.zero_with_output_dimension(VAF, 3))
    MOI.set(model,
        DiffOpt.ForwardConstraintFunction(), c4, MOIU.zero_with_output_dimension(VAF, 1))

    DiffOpt.forward_differentiate!(model)

    # for (i, vi) in enumerate(X)
    #     @test 0.0 ≈ MOI.get(model,
    #         DiffOpt.ForwardVariablePrimal(), vi) atol=1e-2 rtol=RTOL
    # end

    # TODO add a test here, probably on duals

    # # note that there's no change in the PSD slack values or dual optimas
    # @test dy ≈ [0.0, 0.0, 0.0, 0.125978, 0.0, 0.142644, 0.142641, 0.0127401, 0.0, 0.0, 0.0] atol=atol rtol=RTOL
    # @test ds ≈ [0.0, 0.0, -6.66672, 0.0, -3.33336, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] atol=atol rtol=RTOL
end

@testset "Differentiating a simple PSD" begin
    # refer _psd3test, https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L4484
    # find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_0_py.ipynb

    # min x
    # s.t. [x 1 1]
    #      [1 x 1] ⪰ 0
    #      [1 1 x]

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variable(model)

    func = MOIU.operate(vcat, Float64, x, one(Float64), x, one(Float64), one(Float64), x)

    # do not confuse this constraint with the matrix `c` in the conic form (of the matrices A, b, c)
    c = MOI.add_constraint(model, func, MOI.PositiveSemidefiniteConeTriangle(3))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    MOI.optimize!(model)

    # SCS/Mosek specific
    # @test s' ≈ [1.         1.41421356 1.41421356 1.         1.41421356 1.        ] atol=ATOL rtol=RTOL
    # @test y' ≈ [ 0.33333333 -0.23570226 -0.23570226  0.33333333 -0.23570226  0.33333333]  atol=ATOL rtol=RTOL

    # dA = zeros(6, 1)
    # db = ones(6)
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, _vaf(ones(6)))
    # dc = zeros(1)

    DiffOpt.forward_differentiate!(model)

    @test model.diff.model.x ≈ [1.0] atol=10ATOL rtol=10RTOL
    @test model.diff.model.s ≈ ones(6) atol=ATOL rtol=RTOL
    @test model.diff.model.y ≈ [1/3, -1/6, 1/3, -1/6, -1/6, 1/3]  atol=ATOL rtol=RTOL

    @test -0.5 ≈ MOI.get(model,
        DiffOpt.ForwardVariablePrimal(), x) atol=1e-2 rtol=RTOL

    # @test dx ≈ [-0.5] atol=ATOL rtol=RTOL
    # @test dy ≈ zeros(6) atol=ATOL rtol=RTOL
    # @test ds ≈ [0.5, 1.0, 0.5, 1.0, 1.0, 0.5] atol=ATOL rtol=RTOL

    # test 2
    dA = zeros(6, 1)
    db = zeros(6)
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, MOIU.zero_with_output_dimension(VAF, 6))
    MOI.set(model, DiffOpt.ForwardObjectiveFunction(), 1.0 * x)

    DiffOpt.forward_differentiate!(model)

    @test 0.0 ≈ MOI.get(model,
        DiffOpt.ForwardVariablePrimal(), x) atol=1e-2 rtol=RTOL

    # @test dx ≈ zeros(1) atol=ATOL rtol=RTOL
    # @test dy ≈ [0.333333, -0.333333, 0.333333, -0.333333, -0.333333, 0.333333] atol=ATOL rtol=RTOL
    # @test ds ≈ zeros(6) atol=ATOL rtol=RTOL
end

@testset "Verifying cache after differentiating a QP" begin
    Q = [
        4.0 1.0
        1.0 2.0
    ]
    q = [1.0, 1.0]
    G = [1.0 1.0]
    h = [-1.0]

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 2)

    # define objective
    quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
    for i in 1:2
        for j in i:2 # indexes (i,j), (j,i) will be mirrored. specify only one kind
            push!(
                quad_terms,
                MOI.ScalarQuadraticTerm(Q[i,j], x[i], x[j])
            )
        end
    end
    objective_function = MOI.ScalarQuadraticFunction(
        quad_terms,
        MOI.ScalarAffineTerm.(q, x),
        0.0,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)

    # add constraint
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[1, :], x), 0.0),
        MOI.LessThan(h[1]),
    )
    MOI.optimize!(model)

    x_sol = MOI.get(model, MOI.VariablePrimal(), x)
    @test x_sol ≈ [-0.25, -0.75] atol=ATOL rtol=RTOL

    MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, ones(2))

    @test model.diff === nothing
    DiffOpt.reverse_differentiate!(model)

    grad_wrt_h = MOI.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), c))
    @test grad_wrt_h ≈ -1.0 atol=2ATOL rtol=RTOL

    # adding two variables invalidates the cache
    y = MOI.add_variables(model, 2)
    MOI.delete(model, y)

    @test model.diff === nothing
    MOI.optimize!(model)

    MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, ones(2))
    DiffOpt.reverse_differentiate!(model)

    grad_wrt_h = MOI.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), c))
    @test grad_wrt_h ≈ -1.0 atol=2ATOL rtol=RTOL

    # adding single variable invalidates the cache
    y0 = MOI.add_variable(model)
    @test model.diff === nothing
    MOI.add_constraint(model, MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, y0)], 0.0), MOI.EqualTo(42.0))

    MOI.optimize!(model)
    @test model.diff === nothing

    MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, ones(2))
    DiffOpt.reverse_differentiate!(model)
    grad_wrt_h = MOI.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), c))
    @test grad_wrt_h ≈ -1.0 atol=5e-3 rtol=RTOL
    @test model.diff.model.gradient_cache isa DiffOpt.QuadraticProgram.Cache

    # adding constraint invalidates the cache
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], x), 0.0),
        MOI.LessThan(0.0),
    )
    @test model.diff === nothing
    MOI.optimize!(model)

    MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, ones(2))
    DiffOpt.reverse_differentiate!(model)
    grad_wrt_h = MOI.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), c))
    @test grad_wrt_h ≈ -1.0 atol=5e-3 rtol=RTOL
    # second constraint inactive
    grad_wrt_h = MOI.constant(MOI.get(model, DiffOpt.ReverseConstraintFunction(), c2))
    @test grad_wrt_h ≈ 0.0 atol=5e-3 rtol=RTOL
    @test model.diff.model.gradient_cache isa DiffOpt.QuadraticProgram.Cache
end

@testset "Verifying cache on a PSD" begin

    model = DiffOpt.diff_optimizer(SCS.Optimizer)
    MOI.set(model, MOI.Silent(), true)

    x = MOI.add_variable(model)

    func = MOIU.operate(vcat, Float64, x, 1.0, x, 1.0, 1.0, x)

    c = MOI.add_constraint(model, func, MOI.PositiveSemidefiniteConeTriangle(3))

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x]), 0.0)
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test model.diff === nothing
    MOI.optimize!(model)
    @test model.diff === nothing

    dA = zeros(6, 1)
    db = ones(6)
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, _vaf(ones(6)))
    dc = zeros(1)

    DiffOpt.forward_differentiate!(model)

    @test model.diff.model.x ≈ [1.0] atol=10ATOL rtol=10RTOL
    @test model.diff.model.s ≈ ones(6) atol=ATOL rtol=RTOL
    @test model.diff.model.y ≈ [1/3,  -1/6,  1/3,  -1/6,  -1/6,  1/3]  atol=ATOL rtol=RTOL

    @test -0.5 ≈ MOI.get(model,
    DiffOpt.ForwardVariablePrimal(), x) atol=1e-2 rtol=RTOL

    @test model.diff.model.gradient_cache isa DiffOpt.ConicProgram.Cache

    DiffOpt.forward_differentiate!(model)

    @test -0.5 ≈ MOI.get(model,
        DiffOpt.ForwardVariablePrimal(), x) atol=1e-2 rtol=RTOL

    # test 2
    MOI.set(model, DiffOpt.ForwardConstraintFunction(), c, _vaf(zeros(6)))
    MOI.set(model, DiffOpt.ForwardObjectiveFunction(), 1.0 * x)

    DiffOpt.forward_differentiate!(model)

    @test 0.0 ≈ MOI.get(model,
        DiffOpt.ForwardVariablePrimal(), x) atol=1e-2 rtol=RTOL

end

ext: jl


text:
module LinearExpressions

export Symbolic, LinExpr, AbstractVariable, BasicVariable, TypedVariable,
    RealVariable, getcoeff, setcoeff!, differentiate

import Base: hash, show, print, showcompact, convert, promote_rule, zero, one,
             conj, (==), (!=), (+), (-), (*), (.*), transpose

# Borrowed from Calculus.jl
abstract type Symbolic end
abstract type AbstractVariable <: Symbolic end
struct BasicVariable <: AbstractVariable
  sym::Symbol
end
show(io::IO, x::BasicVariable) = print(io, x.sym)
(==)(x::BasicVariable, y::BasicVariable) = x.sym == y.sym

# Symbolic variables with assosiated type (experimental)
struct TypedVariable{T} <: AbstractVariable
  sym::Symbol
end
show(io::IO, x::TypedVariable) = print(io, x.sym)
(==)(x::TypedVariable{T}, y::TypedVariable{T}) where {T} = x.sym == y.sym
const RealVariable = TypedVariable{Real}

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

const ScalarCoeff = Real
const ArrayCoeff{T <: Real} = AbstractArray{T}
const Coeff = Union{ScalarCoeff,ArrayCoeff}

mutable struct LinExpr{Tc<:Coeff,Tv<:AbstractVariable}
  constt::Tc
  coeffs::Dict{Tv,Tc}

  LinExpr{Tc,Tv}(c::Tc, vc::Dict{Tv,Tc}) where {Tc,Tv} = new(c, filter(vc -> vc.second != zero(vc.second), vc))
  LinExpr{Tc,Tv}(c::Tc) where {Tc,Tv} = new(c, Dict{Tv,Tc}())
  LinExpr{Tc,Tv}() where {Tc,Tv} = new(zero(Tc), Dict{Tv,Tc}())
end

LinExpr{Tc,Tv}(vc::Dict{Tv,Tc}) where {Tc,Tv} = LinExpr{Tc,Tv}(zero(Tc), vc)
LinExpr(c::Tc, vc::Dict{Tv,Tc}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(c, vc)
LinExpr(vc::Dict{Tv,Tc}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(zero(Tc), vc)

function linexpr_show(io::IO, e::LinExpr{Tc,Tv}, compact::Bool) where {Tc <: ScalarCoeff, Tv <: AbstractVariable}
  start = true
  if e.constt != zero(e.constt) || isempty(e.coeffs)
    print(io, e.constt)
    start = false
  end
  for (s, c) in e.coeffs
    if start
      start = false
    else
      if c < zero(c)
        c = -c
        print(io, compact ? "-" : " - ")
      else
        print(io, compact ? "+" : " + ")
      end
    end
    if c == -one(c)
      print(io, "-")
    elseif c != one(c)
      print(io, c)
    end
    print(io, s)
  end
end
show(io::IO, e::LinExpr{Tc,Tv}) where {Tc <: ScalarCoeff, Tv <: AbstractVariable} = linexpr_show(io, e, false)
showcompact(io::IO, e::LinExpr{Tc,Tv}) where {Tc <: ScalarCoeff, Tv <: AbstractVariable} = linexpr_show(io, e, true)

function show(io::IO, e::LinExpr{Tc,Tv}) where {Tc <: ArrayCoeff, Tv <: AbstractVariable}
  start = true
  if e.constt != zero(e.constt) || isempty(e.coeffs)
    print(io, e.constt)
    start = false
  end
  for (s, c) in e.coeffs
    if start
      start = false
    else
      print(io, "\n+\n\n")
    end
    print(io, c)
    print(io, s)
    print(io, "\n")
  end
end

function ==(a::LinExpr, b::LinExpr)
  a.constt == b.constt && a.coeffs == b.coeffs
end

function getcoeff(e::LinExpr{Tc,Tv}, v::Tv) where {Tc <: Coeff, Tv <: AbstractVariable}
    get(e.coeffs, v, zero(e.constt))
end
function setcoeff!(e::LinExpr{Tc,Tv}, v::Tv, x::Tc) where {Tc <: Coeff, Tv <: AbstractVariable}
    e.coeffs[v] = x
end

one(e::LinExpr{Tc,Tv}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(one(e.constt))
zero(e::LinExpr{Tc,Tv}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(zero(e.constt))

one(::Type{LinExpr{Tc,Tv}}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(one(Tc))
zero(::Type{LinExpr{Tc,Tv}}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(zero(Tc))

conj(x::TypedVariable{T}) where {T <: Real} = x
conj(e::LinExpr{Tc,TypedVariable{T}}) where {Tc <: Real, T <: Real} = e

function convert(::Type{Tc}, e::LinExpr{Tc,Tv}) where {Tc <: Coeff, Tv <: AbstractVariable}
    if isempty(e.coeffs)
        e.constt
    else
        throw(InexactError(:convert, Tc, e))
    end
end

function convert(::Type{Tv}, e::LinExpr{Tc,Tv}) where {Tc <: ScalarCoeff, Tv <: AbstractVariable}
    if e.constt != zero(Tc) || length(e.coeffs) != 1
        throw(InexactError(:convert, Tv, e))
    else
        s = collect(keys(e.coeffs))[1]
        if e.coeffs[s] != one(Tc)
            throw(InexactError(:convert, Tv, e))
        else
            s
        end
    end
end

convert(::Type{LinExpr{Tc,Tv}}, e::LinExpr{Tc,Tv}) where {Tc <: Coeff, Tv <: AbstractVariable} = e
function convert(::Type{LinExpr{Tc2,Tv}}, e::LinExpr{Tc1,Tv}) where {Tc1 <: Coeff, Tc2 <: Coeff, Tv <: AbstractVariable}
    LinExpr{Tc2,Tv}(convert(Tc2, e.constt), Dict(v => convert(Tc2, c) for (v, c) in e.coeffs))
end

convert(::Type{LinExpr{Tc1,Tv}}, x::Tc2) where {Tc1 <: Coeff, Tv <: AbstractVariable, Tc2 <: Coeff} = LinExpr{Tc1,Tv}(convert(Tc1, x), Dict{Tv,Tc1}())
convert(::Type{LinExpr{Tc,Tv}}, x::Tv) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(zero(Tc), Dict(x => one(Tc)))

promote_rule(::Type{LinExpr{Tc1,Tv}}, ::Type{LinExpr{Tc2,Tv}}) where {Tc1 <: Coeff, Tc2 <: Coeff, Tv <: AbstractVariable} = LinExpr{promote_type(Tc1, Tc2),Tv}
promote_rule(::Type{LinExpr{Tc1,Tv}}, ::Type{Tc2}) where {Tc1 <: Coeff, Tv <: AbstractVariable, Tc2 <: Coeff} = LinExpr{promote_type(Tc1, Tc2),Tv}
promote_rule(::Type{Tc}, ::Type{Tv}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}
promote_rule(::Type{LinExpr{Tc,Tv}}, ::Type{Tv}) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}

function promote_rule(::Type{LinExpr{Array{T1,N},Tv}}, ::Type{LinExpr{Array{T2,N},Tv}}) where {T1 <: ScalarCoeff, T2 <: ScalarCoeff, N, Tv <: AbstractVariable}
  LinExpr{Array{promote_type(T1, T2),N},Tv}
end
function promote_rule(::Type{LinExpr{Array{T1,N},Tv}}, ::Type{Array{T2,N}}) where {T1 <: ScalarCoeff, T2 <: ScalarCoeff, N, Tv <: AbstractVariable}
  LinExpr{Array{promote_type(T1, T2),N},Tv}
end

function *(c::Tc, e::LinExpr{Tc,Tv}) where {Tc <: Coeff, Tv <: AbstractVariable}
  constt = e.constt * c
  coeffs = empty(e.coeffs)
  for (s, coeff) in e.coeffs
    coeffs[s] = c * coeff
  end
  LinExpr{Tc,Tv}(constt, coeffs)
end
*(e::LinExpr{Tc,Tv}, c::Tc) where {Tc <: Coeff, Tv <: AbstractVariable} = c * e

function *(c::Tc2, e::LinExpr{Tc1,Tv}) where {Tc1 <: Coeff, Tv <: AbstractVariable, Tc2 <: Coeff}
    Tc_common = promote_type(Tc1, Tc2)
    convert(Tc_common, c) * convert(LinExpr{Tc_common,Tv}, e)
end
*(e::LinExpr{Tc1,Tv}, c::Tc2) where {Tc1 <: Coeff, Tv <: AbstractVariable, Tc2 <: Coeff} = c * e

*(c::Tc, s::Tv) where {Tc <: Coeff, Tv <: AbstractVariable} = LinExpr{Tc,Tv}(zero(c), Dict(s => c))
*(s::Tv, c::Tc) where {Tc <: Coeff, Tv <: AbstractVariable} = c * s

const VarOrLinExpr = Union{AbstractVariable,LinExpr}

# Base.broadcast(::typeof(*), a::ScalarCoeff, e::VarOrLinExpr) = a * e
# Base.broadcast(::typeof(*), e::VarOrLinExpr, a::ScalarCoeff) = e * a
#
Base.broadcast(::typeof(*), e::VarOrLinExpr, a::ArrayCoeff) = map(x -> e * x, a)
# Base.broadcast(::typeof(*), a::ArrayCoeff, e::VarOrLinExpr) = e .* a

function +(a::T, b::T) where {T <: LinExpr}
  constt = a.constt + b.constt
  coeffs = empty(a.coeffs)
  for (v, c) in a.coeffs
    coeffs[v] = c
  end
  for (v, c) in b.coeffs
    coeffs[v] = c + get(coeffs, v, zero(constt))
  end
  LinExpr(constt, coeffs)
end

+(a::LinExpr) = a

function -(a::LinExpr)
  b = LinExpr(-a.constt, empty(a.coeffs))
  for (v, c) in a.coeffs
    b.coeffs[v] = -c
  end
  b
end

-(a::T, b::T) where {T <: LinExpr} = a + (-b)

-(x::AbstractVariable) = -1 * x
+(x::AbstractVariable) = 1 * x

+(x::VarOrLinExpr, y::VarOrLinExpr) = +(promote(+x, y)...)
-(x::VarOrLinExpr, y::VarOrLinExpr) = x + (-y)

+(x::VarOrLinExpr, y::Coeff) = +(promote(x, y)...)
+(x::Coeff, y::VarOrLinExpr) = y + x

-(x::VarOrLinExpr, y::Coeff) = -(promote(x, y)...)
-(x::Coeff, y::VarOrLinExpr) = x + (-y)

transpose(x::RealVariable) = x
transpose(x::LinExpr{<:Real, <:RealVariable}) = x

==(x::VarOrLinExpr, y::Coeff) = ==(promote(x, y)...)
!=(x::VarOrLinExpr, y::Coeff) = !(==(promote(x, y)...))

==(::T1, ::T2) where {T1 <: AbstractVariable, T2 <: AbstractVariable} = false  # safe for different variable types
==(x::VarOrLinExpr, y::VarOrLinExpr) = ==(promote(x, y)...)
!=(x::VarOrLinExpr, y::VarOrLinExpr) = !(x==y)

differentiate(s::AbstractVariable, d::AbstractVariable) = s == d ? one(1) : zero(0)
differentiate(s::TypedVariable{T}, d::TypedVariable{T}) where {T} = s == d ? one(T) : zero(T)
differentiate(c::Coeff, d::AbstractVariable) = zero(c)
differentiate(e::LinExpr{Tc,Tv}, d::Tv) where {Tc, Tv} = get(e.coeffs, d, zero(e.constt))

end # module


text:
module MyProject

using ForwardDiff
export my_calc, my_deri
greet() = print("Hello World!")

include("calc_file.jl")

end


text:
abstract type AbstractMarkovChain <: AbstractMeasure end

abstract type AbstractMarkovChainPrior <: AbstractMeasure end

function Dists.fit(mctype::Type{<:AbstractMarkovChain}, args...; kwargs...)
    return fit_mle(mctype, args..., kwargs...)
end

function Dists.fit_mle(mctype::Type{<:AbstractMarkovChain}, args...; kwargs...)
    ss = suffstats(mctype, args...; kwargs...)
    return fit_mle(mctype, ss)
end


text:
module hyporheicBiogeochemistry

using DifferentialEquations
using Turing
using DiffEqOperators
using QuadGK
using CSV
using DataFrames
using ModelingToolkit
using Distributions
using MCMCChains

include("residenceTimeDistributions.jl")
include("advect.jl")
include("advectReact.jl")
include("prepData.jl")
include("reactionKinetics.jl")
include("mesocosmModels.jl")
include("statisticalModels.jl")
include("streamModels.jl")
include("simulateMesocosm.jl")
include("generateSyntheticMesocosmData.jl")

export qCalc_powerLaw
export prepData
export firstOrder
export firstOrderProductionAndConsumption
export E_powerLaw
export fit_sam_model
export simulateMesocosm
export generateSyntheticMesocosmData
export @advect
export @advectReact
export @model
export @sam_model
export @stream_model

end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
using PackageEvaluator
println(evalPkg("JuMP", false))


text:
#=
  Bayesian Data Analysis of Wine dataset

  I'm continue to testingMvNormal models to see how well it does on
  various datasets.

  From http://archive.ics.uci.edu/ml/datasets/Wine
  """
    Original Owners:

    Forina, M. et al, PARVUS -
    An Extendible Package for Data Exploration, Classification and Correlation.
    Institute of Pharmaceutical and Food Analysis and Technologies, Via Brigata Salerno,
    16147 Genoa, Italy.

    Donor:

    Stefan Aeberhard, email: stefan '@' coral.cs.jcu.edu.au

    Data Set Information:

    These data are the results of a chemical analysis of wines grown in the same region in Italy
    but derived from three different cultivars. The analysis determined the quantities of 13
    constituents found in each of the three types of wines.

    I think that the initial data set had around 30 variables, but for some reason I only have the
    13 dimensional version. I had a list of what the 30 or so variables were,
    but a.) I lost it, and b.), I would not know which 13 variables are
    included in the set.

    The attributes are (dontated by Riccardo Leardi, riclea '@' anchem.unige.it )
    1) Alcohol
    2) Malic acid
    3) Ash
    4) Alcalinity of ash
    5) Magnesium
    6) Total phenols
    7) Flavanoids
    8) Nonflavanoid phenols
    9) Proanthocyanins
    10)Color intensity
    11)Hue
    12)OD280/OD315 of diluted wines
    13)Proline

    In a classification context, this is a well posed problem with "well behaved" class structures.
    A good data set for first testing of a new classifier, but not very challenging.
  """

  The dataset is in the file wine.csv. The first entry is the class: 1, 2, 3 (Green,Red,Yellow)

  Summary:
  Num correct: 171.0
  Num incorrect: 7.0
  Accuracy: 0.9606741573033708
  Confusion matrix:
  3×3 Matrix{Float64}:
   56.0   3.0   0.0
    1.0  69.0   1.0
    0.0   2.0  46.0

  Bad classfied:
  7-element Vector{Any}:
    5 => (1 => 2)
    22 => (1 => 2)
    44 => (1 => 2)
    62 => (2 => 3)
    122 => (2 => 1)
    131 => (3 => 2)
    142 => (3 => 2)
  num bad: 7

  Note that this use the full dataset in training and testing. Comparing to
  other algorithms (e.g. Weka's Logistic, etc they most result in 100% correct predictions. )

=#
using Turing, StatsPlots, Distributions
using DataFrames, StatsBase
using CSV
include("jl_utils.jl")

#
# Read dataset and preprocessing
#
df = CSV.read("wine.csv",DataFrame,header=false)

display(describe(df))
num_rows,num_cols = size(df)
data = Matrix(df)
x = data[:,2:num_cols]
y = data[:,1]

categories = unique(y) |> sort
println("Mean value for each species: ", categories)
display(categories .|> s -> (s,mean(x[y.==s,:], dims=1)))
println("\n")

# Standardize the dataset (better and faster)
# x = StatsBase.standardize(UnitRangeTransform,x) # This give very bad result
means = mean(x,dims=1)
stds = std(x,dims=1)
maxes = [maximum(x[:,i]) for i in 1:size(x)[2]]
x = (x .- means) ./ stds
# x = (x .- means) ./ maxes'
# x = x ./ maxes'

println("mean(x): ", mean(x,dims=1))
println("std(x): ", std(x,dims=1))

@model function wine_model(x,y=missing, ::Type{T} = Float64) where {T}
    n,p = size(x)

    if y === missing
      y = Vector{T}(undef, n)
    end

    a ~ Normal(0,1)
    params ~ filldist(Normal(0,1),p)
    sigma ~ Uniform(0,1)

    mu = a .+ x*params
    y ~ MvNormal(mu, sigma)

    return y

end

model = wine_model(x,y)

# println("Prior:")
# chns = sample(model, Prior(), 10_000)
# display(chns)

println("\nModel:")
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
# chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, HMC(0.1,6), 1_000)
chns = sample(model, NUTS(), 1_000)

display(chns)

# display(plot(chns))
println("\nPredictions:")

# mpred = iris_model(x[50:51,:],missing)
mpred = wine_model(x,missing)
pred = predict(mpred,chns)
display(pred)

cats = Dict(1=>1, 2=>2, 3=>3)
confusion, bads = confusion_matrix(y,pred,cats,true)
println("Confusion matrix:")
display(confusion)
println("\nBad classfied:")
display(bads)
println("num bad: $(length(bads))")
# println("\nbads: credible_intervals")
# for b in bads
#   println("b:",b)
#   ix = b[1]
#   println("x[$ix]: $(x[ix,:])  true cat (y[$ix]): $(y[ix])  (pred mean: ", mean(pred["y[$ix]"]))
#   # credible_interval(pred,"y[$(b[1])]",0.3)
# end

# println("\nMean value for each category: ", 1:3)
# display(1:3 .|> s -> (s,mean(x[y.==s,:], dims=1)))
# println("\n")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#!me TODO:
#===
#
#  Trait trick is probably needed to decomplect solve my type issues
f(a, storage::Contiguous) = ...
f(a, storage::Strided) = ...
f(a) = f(a, getstorage(a))

...

getstorage(::MyType) = Contiguous()
#
===#
#==
# iterators and path Reverse need to be sorted out
==#

export
EmptyNode,
Node,
Branch,
Leaves,
BinaryTree,
Traverse,
TraverseWhen,
DepthFirst,
BreadthFirst,
BinarySearch,
Path,
ReversePath,
empty_node,
isempty_node,
left,
right,
data,
key,
value,
insert,
search,
tree,
leaves,
isleaf,
rotate_left,
rotate_right,
path,
reverse_path,
replace_node,
delete,
update,
fullness,
minheight

using Base.Iterators
using Base.Iterators:Reverse
using Base.Iterators:Drop

∂ = (f,x)-> (y...)->f(x,y...)

abstract type BinaryTree{T} end
const EmptyNode = Nothing
const empty_node = nothing

const Branch{T} = Union{BinaryTree{T},EmptyNode}

abstract type WrappedData{T} end

left(tree::EmptyNode) = nothing
right(tree::EmptyNode) = nothing
key(tree::EmptyNode) = nothing
value(tree::EmptyNode) = nothing
data(tree::EmptyNode) = nothing
keyval(tree::EmptyNode) = nothing

struct Node{T} <: BinaryTree{T}
  data::T
  left::Branch{T}
  right::Branch{T}
end

Node(data::T) where T = Node(data,empty_node,empty_node)
Node{T}(data::T) where T = Node(data)

left(tree::Node{T}) where T = tree.left
right(tree::Node{T}) where T = tree.right
data(tree::Node{T}) where T = tree.data
key(tree::Node{T}) where T = key(tree.data)
value(tree::Node{T}) where T = value(tree.data)
keyval(tree::Node{T}) where T = tree.data
leaf(data::T) where T = Node(data)

key(d::Any) = d
value(d::Any) = d

#==
internal_nodedata(::Type{T}) where T = one(T)
internal_nodedata(::Type{T}) where {T<:AbstractArray} = T()
internal_nodedata(::Type{W}) where {T,W<:WrappedData{T}} = internal_nodedata(T)
==#

#  Traits  #

# a full tree always has either 0 or 2 children per branch
abstract type Fullness  end
struct FullTree <: Fullness end
struct NotFullTree <: Fullness end

fullness(::Type{BT}) where {T, BT<:BinaryTree{T}} = NotFullTree
fullness(::Type{BT}) where {T, W<:WrappedData{T}, BT<:BinaryTree{W}} = fullness(BinaryTree{T})

# AbstractTrees interface
import AbstractTrees.children
import AbstractTrees.printnode
import AbstractTrees.Leaves

children( n::Nothing ) = ()
children( n::BinaryTree{T} ) where T = isleaf(n) ? () : (left(n), right(n))

printnode(io::IO, n::BinaryTree{T}) where T = show(data(n))
printnode(io::IO, n::Nothing ) = nothing

struct DepthFirst{T}
  tree::BinaryTree{T}
end

rest(a::Vector{T}) where T = Vector{T}(a[2:end])
rest(a::NTuple{T}) where T = a[2:end]

isempty_node(n::Nothing) = true
isempty_node(n) = false

isleaf(a::BinaryTree{T}) where T = (isempty_node(left(a)) && isempty_node(right(a)))

const PathElement{BT} = Tuple{BT,Symbol}
const Path{BT} = Vector{PathElement{BT}}
const ReversePath{BT} = Union{Reverse{Path{BT}},Drop{Reverse{Path{BT}}}}  # awkward.

#==== General Traversal ===#
# Transducers should give better performance than Channel

function traverse(node::BinaryTree{T}, c::Channel, pre, post, inorder) where T
  if pre != nothing
    push!(c, pre(node))
  end

  traverse(left(node), c, pre, post, inorder)
  if inorder != nothing
    push!(c, inorder(node))
  end

  traverse(right(node), c, pre, post, inorder)
  if post != nothing
    push!(c, post(node))
  end
end

function traverse(node::Nothing, c::Channel, f...)
end

Traverse(tree::BinaryTree{T}, pre, post = nothing, inorder = nothing) where T = Channel() do c
  traverse(tree, c, pre, post, inorder)
end

function traversewhen(when::Function, node::BinaryTree{T}, c::Channel, pre, post, inorder) where T
  if !when(node)
    return
  end

  if pre != nothing
    push!(c, pre(node))
  end

  traversewhen(when, left(node), c, pre, post, inorder)
  if inorder != nothing
    push!(c, inorder(node))
  end

  traversewhen(when, right(node), c, pre, post, inorder)
  if post != nothing
    push!(c, post(node))
  end
end

function traversewhen(when::Function, node::Nothing, c::Channel, f...)
end

TraverseWhen(when, tree::BinaryTree{T}, pre, post = nothing, inorder = nothing) where T = Channel() do c
  traversewhen(when, tree, c, pre, post, inorder)
end

#==== DepthFirst Traversal ===#
function Base.iterate(iter::DepthFirst{T}, children::Vector{BinaryTree{T}}) where T
  if length(children) == 0
    return nothing
  end
  child = first(children)
  (child, vcat(Vector{BinaryTree{T}}(filter(x->x!=nothing,[left(child), right(child)])), rest(children)))
end

function Base.iterate(iter::DepthFirst{T}) where T
  (iter.tree,Vector{BinaryTree{T}}(filter(x->x!=nothing,[left(iter.tree), right(iter.tree)])))
end

Base.length(iter::DepthFirst{T}) where T = reduce( (acc,_)->acc+1, iter; init=0)

Base.iterate(tree::BinaryTree{T}) where T = iterate(DepthFirst(tree))
Base.iterate(iter::BinaryTree{T}, children::Vector{BinaryTree{T}}) where T = iterate(DepthFirst(iter), children)
Base.length(iter::BinaryTree{T}) where T = length(DepthFirst{T}(iter))

#===== BreadthFirst Traversal ===#
struct BreadthFirst{T}
  tree::BinaryTree{T}
end

function Base.iterate(iter::BreadthFirst{T}, children::Vector{BinaryTree{T}}) where T
  if length(children) == 0
    return nothing
  end
  child = last(children)
  (child, vcat(Vector{BinaryTree{T}}(filter(x->x!=nothing,[right(child), left(child)])), children[1:end-1]))
end

function Base.iterate(iter::BreadthFirst{T}) where T
 (iter.tree,Vector{BinaryTree{T}}(filter(x->x!=nothing,[right(iter.tree), left(iter.tree)])))
end

Base.length(iter::BreadthFirst{T}) where T = length(DepthFirst(iter.tree))

#===== BinarySearch Traversal ===#
struct BinarySearch{T}
  tree::BinaryTree{T}
  comparison::Function
end

function pathfinder( k::T ) where T
  st->begin
    if key(k) == key(st)
      :found
    elseif key(k) < key(st)
      :left
    else
      :right
    end
  end
end

BinarySearch(tree::BinaryTree{T}, k::T) where {T} =
  BinarySearch{T}( tree, pathfinder(k) )

function Base.iterate(iter::BinarySearch{T}, subtree::Branch{T}) where T
  if subtree == empty_node
    return nothing
  end

  if iter.comparison(subtree) == :found
    ((subtree,:found),nothing)
  elseif iter.comparison(subtree) == :left
    ((subtree,:left), left(subtree))
  else
    ((subtree,:right), right(subtree))
  end
end

function Base.iterate(iter::BinarySearch{T}) where T
  iterate(iter, iter.tree)
end

Base.length(iter::BinarySearch{T}) where T = reduce( (acc,_)->acc+1, iter; init=0)

#height(t::BinaryTree{T}) where T = reduce(max,length.(∂(path,t).(key.(leaves(t)))))
height(::EmptyNode, h::Int) = h
height(t::Branch{T}, h::Int) where T = max(height(left(t),h+1), height(right(t), h+1))
height(t::BinaryTree{T}) where T = height(t, 0)

#===== Leaf Traversal ===#
Leaves( a::BinaryTree{T}, descend::Function = x->true ) where {T} = Channel(ctype=BinaryTree{T}) do c
  nodestack = Vector{BinaryTree{T}}()

  if descend(a)
    push!(nodestack,a)
  end

  while( !isempty(nodestack) )
    ancestor = pop!(nodestack)

    if( isleaf(ancestor) )
      push!(c, ancestor)
    else
      childL = left(ancestor)
      childR = right(ancestor)

      if childR != empty_node && descend(childR)
        push!(nodestack,childR)
      end
      if childL != empty_node && descend(childL)
        push!(nodestack,childL)
      end
    end
  end
end

leaves( a::BinaryTree{T}, descend::Function = x->true ) where T = Leaves( a, descend )

"""
  search nodes until ftest passes. returns the subtree rooted there
"""
search( ftest::Function, tree::BinaryTree{T} ) where T = Iterators.filter(ftest,tree) |> first

"""
  search nodes until we find the node with given value. returns the subtree rooted there
"""
search( k::T, tree::BinaryTree{T} ) where T = search(x->key(x)==k, tree)

"""
  search nodes until ftest passes, return the path to that node
"""
path( k::T, tree::BinaryTree{T} ) where T = path(BinarySearch(tree, k))
path( tree::BinaryTree{T}, k::T ) where T = path(k,tree)
path( bs::BinarySearch{T} ) where T =  (bs|>∂(collect,PathElement{BinaryTree{T}}))

"""
  reverse the path from top down to bottom up
"""
reverse( p::Path{BT} ) where BT = Reverse(p)
reverse( bs::BinarySearch{T} ) where T = Reverse(path(bs))

reverse_path = reverse∘path

"""
  returns a function to compare one node to another. for binary search. i.e. return :left if first < second
"""
branch_compare( d::S ) where {S} = st-> d < key(st) ? :left : :right

"""
  return a new tree with new node inserted on branch where data < parent.data decends to left
"""
insert( ::Type{NotFullTree}, t::N, d::T ) where {T, N<:BinaryTree{T}} =
  rbuild(reverse_path(BinarySearch(t,
                                   branch_compare( key(d) ))),
         N(d,empty_node,empty_node))

function insert( ::Type{FullTree}, t::N, d::T ) where {T, N<:BinaryTree{T}}
  rpath = reverse_path(BinarySearch(t,
                                    branch_compare( key(d) )))

  # no entities are stored on branches, so create new Branch with existing Leaf right and new Leaf left
  cleaf = (take(rpath,1) |> first)[1]
  rpath = drop(rpath,1)

  twig = rebuild_leaf( FullTree, :both, cleaf, d )
  rbuild( rpath, twig )
end

insert( t::N, d::T ) where {T, N<:BinaryTree{T}} = insert( fullness(N), t, d )
#insert( t::N, d::T ) where {T, N<:Node{T}} = insert( fullness(N), t, d )
#insert( t::N, d::T ) where {T, N<:BinaryTree{T}} = insert( fullness(N), t, d )
"""
  rebuild a tree given a path.  path will be created, all other branches preserved
"""
build( path::Path{BT} ) where BT =  rbuild( Reverse(path) )

rebuild_leaf( ::Type{NotFullTree}, direction::Symbol, parent::N, d::T ) where {N,T} =
  rebuild_node( Val(direction), empty_node, empty_node, parent, d )

function rebuild_leaf(::Type{FullTree},
                      direction::Symbol,
                      fatherbrother::NODE,
                      d::T) where {T,NODE<:BinaryTree{T}}
#  newleaf = rebuild_node( Val(direction), empty_node, empty_node, fatherbrother, d )
  newleaf = NODE(d)
  twigd = rebuild_parent_data( Val(:both), keyval(newleaf), keyval(fatherbrother), nothing )
  # well now, this is incestuous. the last fatherbrother is ignored usually.
  #!me maybe that last fatherbrother should be nothing?  then change rebuild_node to accept Branch there
  #!me might result in ambiguous dispatch for wrapped type?
  twig = rebuild_node(Val(:both), newleaf, fatherbrother, fatherbrother, twigd)
end

#WIP
rebuild_parent_data(dir,#::Union{Val{:both}, Val{:left}, Val{:right}, Val{:either}},
                    lcd::Union{T,Nothing}, rcd::Union{T,Nothing},
                    parent_d::T ) where T = parent_d

#==
rebuild_parent_data(dir, #::Union{Val{:both}, Val{:left}, Val{:right}, Val{:either}},
                    lcd::Union{T,Nothing}, rcd::Union{T,Nothing},
                    _::Nothing ) where {T<:Union{AbstractString,Number,AbstractArray}} = internal_nodedata(T)
==#

# interval tree, called when BinaryTree{Vector{T}} has fullness set to FullTree
function rebuild_parent_data(Any,
  lcd::Vector{T},
  rcd::Vector{T},
  parent::Union{Nothing, Vector{T}}) where T
  Vector{T}(sort(vcat(lcd, rcd))[[1,end]])
end

rebuild_node(::Union{Val{:both}, Val{:left}, Val{:right}},
             lc::Branch{T}, rc::Branch{T},
             n::Node{T}, d::T ) where T =
  Node(d, lc, rc)

function rebuild_node( direction::Symbol, lc::Branch{T}, rc::Branch{T}, n::BinaryTree{T} ) where T
  d = rebuild_parent_data( Val(direction), keyval(lc), keyval(rc), keyval(n) )
  rebuild_node( Val(direction), lc, rc, n, d )
end

function rebuild_node(direction::Symbol, child::Branch{T}, n::BinaryTree{T} ) where T
  if direction == :left
    rebuild_node( direction, child, right(n), n )
  else
    rebuild_node( direction, left(n), child, n )
  end
end

update(n::Branch{T}) where T = n
#WIP

"""
  rebuild a tree given a reversed (bottom up) path.  path will be created, all other branches preserved
  to implement a new type of BinaryTree, add a new method with this signature for rebuilding
"""
function rbuild( rpath::ReversePath{BT}, leaf::Branch{T} ) where {T, BT<:BinaryTree{T}}
  reduce((st,(n, direction)) -> rebuild_node( direction, st, n ),
         rpath ;
         init=update(leaf))
end

rbuild( rpath::ReversePath{BT} ) where BT = rbuild(drop(rpath,1), (first∘first)(rpath))

"""
  replace the deepest node in dest path with source subtree
"""
function replace_node( dest::ReversePath{BT}, source::Branch{T} ) where {T, BT<:BinaryTree{T}}
  rpath = drop(dest, 1)
  rbuild(rpath,source)
end

function delete( rpath::ReversePath{BT} ) where BT
  node,dir = first(rpath)
  dir == :found || return first(last(rpath))
  lc = left(node)
  rc = right(node)

  if lc == empty_node && rc == empty_node
    replace_node( rpath, empty_node )
  elseif lc == empty_node
    replace_node(rpath, rc)
  elseif rc == empty_node
    replace_node(rpath, lc)
  else  # both children present
    # choose left or right?  right for now ( right then take leftmost )
    # the new root of this subtree will be leftmost leaf of right branch
    goleft = x->:left
    successor_path = reverse_path(BinarySearch(rc,goleft))
    successor = (first∘first)(successor_path)
    #rtree = replace_node(successor_path, empty_node)
    rtree = replace_node(successor_path, right(successor))
    #!me NG
    #replace_node(rpath, Node(data(predecessor),lc,rtree))
    #newsource = rebuild_node(:left, lc,
    newelder = rebuild_node(:both, lc, rtree, successor)
    replace_node(rpath, newelder)
  end

end

delete( t::BinaryTree{T}, key::T ) where T = delete(reverse_path(key,t))
delete( key::T, t::BinaryTree{T} ) where T = delete(t, key)
delete( t::BinaryTree{W}, d::T ) where {T,W<:WrappedData{T}} = delete(t, W(d))

tree( d::T ) where T = Node(d,empty_node,empty_node)

tree( dc::Vector{T} ) where T = reduce(insert,rest(dc); init = tree(first(dc)))

#==
function rotate_right( parent::Node{T} ) where T
  child = left(parent)
  Node(data(child), left(child), Node(data(parent), right(child), right(parent)))
end

function rotate_left( parent::Node{T} ) where T
  child = right(parent)
  Node(data(child), Node(data(parent), left(parent), left(child)), right(child))
end
==#

function rotate_right( parent::N ) where {T, N<:BinaryTree{T}}
  child = left(parent)
  rebuild_node(:right, left(child), rebuild_node(:left, right(child), right(parent), parent), child)
  #Node(data(child), left(child), Node(data(parent), right(child), right(parent)))
end

function rotate_left( parent::N ) where {T, N<:BinaryTree{T}}
  child = right(parent)
  rebuild_node(:left, rebuild_node(:right, left(parent), left(child), parent), right(child), child)
  #Node(data(child), Node(data(parent), left(parent), left(child)), right(child))
end

update( t::BinaryTree{T}, data::T, newdata::T ) where T = insert(delete(t, key(data)), newdata)

""" minimum height of a binary tree with n leaves """
minheight(n) = log(2,n)+1


text:
# Functions that find special events, like e.g. lyrics, or tracknames
# are contained here

export trackname, addtrackname!, textevent, findtextevents
export tracknames

const NOTRACKNAME = "No track name found"

"""
    trackname(track::MIDI.MIDITrack)

Return the name of the given `track` as a string,
by finding the "track name" `MetaEvent`.

If no such event exists, `"No track name found"` is returned.
"""
function trackname(track::MIDI.MIDITrack)

    pos = findtrackname(track)
    if pos == 0
        return NOTRACKNAME
    # check if there really is a name
    elseif length(track.events[pos].data) == 0
        return NOTRACKNAME
    else
        event = track.events[pos]
        # extract the name (string(Char()) takes care of ASCII encoding)
        trackname = string(Char(event.data[1]))
        for c in event.data[2:end]
            trackname *= string(Char(c))
        end
        return trackname
    end
end

"`tracknames(m::MIDIFile) = trackname.(m.tracks)`"
tracknames(m::MIDIFile) = trackname.(m.tracks)

"""
    addtrackname!(track::MIDI.MIDITrack, name::String)

Add a name to the given `track` by attaching the
"track name" `MetaEvent` to the start of the `track`.
"""
function addtrackname!(track::MIDI.MIDITrack, name::String)
    # construct fitting name event
    data = UInt8[]
    for i = 1:length(name)
        push!(data, UInt8(name[i]))
    end
    meta = MetaEvent(0,0x03,data)

    # remove existing name
    prev = findtrackname(track)
    if prev != 0
        deleteat!(track.events, prev)
    end

    addevent!(track, 0, meta)
end

function findtrackname(track::MIDI.MIDITrack)
    position = 0
    for (i,event) in enumerate(track.events)
        if isa(event, MIDI.MetaEvent) && event.metatype == 0x03
            position = i
            break
        end
    end
    return position
end

"""
    textevent(eventtype, text)
Create an event using the string `text`.
The `eventtype` can be `:text, :lyric, :marker`, which will create the
appropriate type of `MetaEvent`.

The returned event can be added to a [`MIDITrack`](@ref) via either
[`addevent!`](@ref) or [`addevents!`](@ref) for multiple events.

*Notice* - Cubase can read the marker events and MuseScore can read the lyrics
events. We haven't seen any editor that can read the text events, so far.
"""
function textevent(eventtype, text)
    data = [UInt8(a) for a in text]
    event = MetaEvent(0, s_to_text[eventtype], data)
end

const s_to_text = Dict(
:text => TEXTEV, :lyric => LYRICEV, :marker => MARKEREV
)

"""
    findtextevents(eventtype, track)
Find all "text" events specifield by `eventtype` in the `track`.
The `eventtype` can be `:text, :lyric, :marker`, which will find the
appropriate `MetaEvent`s.

For convenience, this function does not return the events themselves.
Instead, it returns three vectors: the first is the strings of the events,
the second is the indices of the events in the `track` and the third is
the absolute position of the events (since start of `track`).

*Notice* - common music score editors like e.g. MuseScore, GuitarPro, etc., do
not export the lyrics and text information when exporting midi files.
"""
function findtextevents(eventtype, track)
    etype = s_to_text[eventtype]
    events = String[]
    idxs = Int[]
    for (i, event) ∈ enumerate(track.events)
        if typeof(event) == MetaEvent && event.metatype == etype
            # lol without copy the events are lost!
            push!(events, String(copy(event.data)))
            push!(idxs, i)
        end
    end
    length(idxs) == 0 && error("Found no events of such type")
    abspos = get_abs_pos(idxs, track)
    return events, idxs, abspos
end


text:
# https://github.com/jingzbu/InverseVIsTraffic/blob/master/08_develop_new_OD_demand_estimator_MA_uni_class_cdc16/06_demands_adjustment_MA.ipynb
# REMEMBER THAT THE FOLDER MUST BE EMPTY!!
#Importing parameters
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np
@pyimport json
@pyimport os
@pyimport pickle
@pyimport pandas as pd
@pyimport collections
@pyimport parameters_julia
@pyimport utils_julia

out_dir = parameters_julia.out_dir
files_ID = parameters_julia.files_ID
month_w = parameters_julia.month_w
year = parameters_julia.year
instances_1 = parameters_julia.instances_ID
deg_grid = parameters_julia.deg_grid
c_grid = parameters_julia.c_grid
lamb_grid = parameters_julia.lamb_grid
week_day_Apr_list = parameters_julia.week_day_list

include("Julia_files/initia_data.jl");
include("prepare_data.jl");
#include("Julia_files/inverseVI.jl");
include("Julia_files/demands_adjustment_gradi.jl");

function demandsDictFixed(demandsDict, flow_observ, link_vector, graph, ta_data, link_dic, day, gamma1,
    gamma2, out_dir, files_ID, month_w, instance, key_, free_flow_time, capacity, start_node, en_node, numZones, cnt)
    #day = 4  # day of April
    # observed flow vector
    xl = flow_observ[:, cnt]

    tapFlows = Dict()
    for i = 1:length(ta_data.start_node)
        key = link_vector[i]
        tapFlows[key] = xl[i]
    end
    #println(tapFlows)
    tapFlowVect = xl;

    # get observed flow vector (corresponding to ground truth demands and ground truth costs)
    tapFlowDicDict[0], tapFlowVecDict[0] = tapFlows, tapFlowVect;

    # get arcs data corresponding to ground truth demands and flows !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    network_data_file = files_ID * "_net_" * month_w * "_full_" * instance * ".txt"
    arcsDict[0] = observFlow(out_dir  * "data_traffic_assignment_uni-class/" * network_data_file, tapFlowDicDict[0]);

    coeffs_dict_Apr_weekend_ = readstring(out_dir * "coeffs_dict_"* month_w * "_" * instance * ".json")
    coeffs_dict_Apr_weekend_ = JSON.parse(coeffs_dict_Apr_weekend_)
    fcoeffs = coeffs_dict_Apr_weekend_[key_]

    demandsVecDict[1] = demandsDicToVec(demandsDict[1]);
    #println(demandsVecDict[1])

    objFunDict[1] = objF(graph, ta_data, link_dic, gamma1, gamma2, demandsVecDict[1], demandsVecDict[1], tapFlowVecDict, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones);

    # get initial flow vector (corresponding to initial demands)
    tapFlowDicDict[1], tapFlowVecDict[1] = tapMSA(graph, ta_data, link_dic, demandsDict[1], fcoeffs, free_flow_time, capacity, start_node, en_node, numZones)
    #tapMSA(graph, ta_data, link_dic, demandsDict[1], fcoeffs, free_flow_time, capacity, start_node, en_node,  numIter=500, tol=1e-6);

    demandsVecDict[0] = demandsDicToVec(demandsDict[0]);

    # get arcs data corresponding to initial demands and flows !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    network_data_file = files_ID * "_net_" * month_w * "_full_" * instance * ".txt"
    arcsDict[1] = observFlow(out_dir  * "data_traffic_assignment_uni-class/" * network_data_file, tapFlowDicDict[1]);

    linkCostDicDict[1] = tapFlowVecToLinkCostDict(tapFlowVecDict[1], fcoeffs, free_flow_time, capacity);

    linkCostDicDict[1]["0"], link_length_list[1]

    jacobiSpiessDict[1] = Compute_Jacobian.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                                  link_list_js, [linkCostDicDict[1]["$(i)"] for i=0:numLinks-1]);

    # maximum number of iterations
    N = 100;

    # Armijo rule parameters
    rho = 2;
    M = 10;

    # search direction parameter
    epsilon_1 = 0;

    # stop criterion parameter
    epsilon_2 = 1e-20;
    obj = 9999999999999

    las = 0
    for l = 1:N

        jacobDict[l] = jacobiSpiessDict[l]

        descDirecDict[l] = descDirec(gamma1, gamma2, demandsVecDict[l], demandsVecDict[1],  tapFlowVecDict[l],
            tapFlowVecDict[0], jacobDict[l], numODpairs, numLinks);

        demandsVecDict[l] = demandsDicToVec(demandsDict[l]);

        searchDirecDict[l] = searchDirec(demandsVecDict[l], descDirecDict[l], epsilon_1);

        thetaMaxDict[l] = thetaMax(demandsVecDict[l], searchDirecDict[l]);

        demandsVecDict[l+1] = similar(demandsVecDict[0]);

        demandsVecDict[l+1], objFunDict[l+1] = armijo(gamma1, gamma2, objFunDict[l], demandsVecDict[l], demandsVecDict[0], tapFlowVecDict, fcoeffs,
        	searchDirecDict[l], thetaMaxDict[l], rho, M, graph, ta_data, link_dic, free_flow_time, capacity, start_node, en_node, numZones);

        demandsDict[l+1] = demandsVecToDic(demandsVecDict[l+1]);

        tapFlowDicDict[l+1], tapFlowVecDict[l+1] = tapMSA(graph, ta_data, link_dic, demandsDict[l+1], fcoeffs, free_flow_time, capacity, start_node, en_node, numZones);

        arcsDict[l+1] = observFlow(out_dir  * "data_traffic_assignment_uni-class/" * network_data_file, tapFlowDicDict[l+1]);

        linkCostDicDict[l+1] = tapFlowVecToLinkCostDict(tapFlowVecDict[l+1], fcoeffs, free_flow_time, capacity);

        jacobiSpiessDict[l+1] = Compute_Jacobian.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                                  link_list_js, [linkCostDicDict[l+1]["$(i)"] for i=0:numLinks-1]);

        demandsDiffDict[l+1] = norm(demandsVecDict[l+1] - demandsVecDict[0]) / norm(demandsVecDict[0]);

        obj =  objFunDict[1] - objFunDict[l]

        las = l
        # stopping criterion
        if (objFunDict[l] - objFunDict[l+1]) / objFunDict[1] < epsilon_2
            break
        end

        println("iteration $(l) finished...")

    end

    # normalize objective function value
    for a = 1:(length(objFunDict))
        norObjFunDict[a] = objFunDict[a] / objFunDict[1];
    end

    outfile = open(out_dir * "demandsDict/demandsVecDict$(day)_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, demandsVecDict)

    close(outfile)

    outfile = open(out_dir * "demandsDict/demandsDict$(day)_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, demandsDict)

    close(outfile)

    outfile = open(out_dir * "demandsDict/tapFlowDicDict$(day)_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, tapFlowDicDict)

    close(outfile)

    outfile = open(out_dir * "demandsDict/tapFlowVecDict$(day)_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, tapFlowVecDict)

    close(outfile)

    outfile = open(out_dir * "demandsDict/jacobi$(day)_" * month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, jacobiSpiessDict[las])

    close(outfile)

    demandsDict[length(demandsDict)-1]

    demandsDict_ = Dict{}()

    for key in keys(demandsDict[length(demandsDict)-1])
        demandsDict_[key] = demandsDict[length(demandsDict)-1][key]
    end

    outfile = open(out_dir * "demandsDict/demandsDictFixed$(day)_"* month_w * "_" * instance * ".json", "w")

    JSON.print(outfile, demandsDict_)

    close(outfile)

    return obj

end

function socialObj(linkFlowVec, free_flow_time, polyDeg, fcoeffs, capacity, numLinks)
    objVal = sum([sum([free_flow_time[a] * fcoeffs[i] * linkFlowVec[a]^i / capacity[a]^(i-1) for i=1:polyDeg])
        for a = 1:numLinks])
    return objVal
end

#key_ = "(6, 1.5, 0.1, 1)"
using PyCall
function biblev(instance)
	#instance = "AM"
	instance1 = instance
	open(out_dir * "instance_comm.txt", "w") do f
	    write(f, instance1)
	end

	unshift!(PyVector(pyimport("sys")["path"]), "");
	@pyimport parameters_julia
	@pyimport GLS_julia
	@pyimport Compute_Jacobian

	numNodes = Compute_Jacobian.numNodes;
	numLinks = Compute_Jacobian.numLinks;
	numODpairs = Compute_Jacobian.numODpairs;
	numZones = Compute_Jacobian.numZones
	od_pairs = Compute_Jacobian.od_pairs;
	link_list_js = Compute_Jacobian.link_list_js;
	link_length_list = Compute_Jacobian.link_length_list;

	flow_observ, link_vector = GLS_julia.GLS_juliaf()

	include("extract_data.jl");
	include("Julia_files/tap_MSA.jl");
	include("Julia_files/demands_adjustment_gradi.jl");

	cnt = 0
	PoA_dict = Dict();
	tapSocialFlowDicDict = Dict();
	tapSocialFlowVecDict = Dict();
	user_sol_dict = Dict();
	social_sol_dict = Dict();
	obj_dict = Dict();

	for day in week_day_Apr_list
		cnt = cnt + 1
	#day = 9
		demandsDict = extract_demandDict(out_dir, files_ID, month_w ,day, instance1, numZones)
		numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data_Apr_PM, start_node, en_node = extract_dataf(out_dir, files_ID, month_w ,day, instance1, demandsDict)

		# preparing a graph
		graph = create_graph(start_node, en_node);
		link_dic = sparse(start_node, en_node, 1:numLinks);

		demandsVecDict[0] = demandsDicToVec(demandsDict[0]);
		demandsDiffDict[1] = norm(demandsDicToVec(demandsDict[1]) - demandsDicToVec(demandsDict[0]))/
		                     norm(demandsDicToVec(demandsDict[0]));
		gamma1 = 0
		gamma2 = 100

		if isdir(out_dir * "demandsDict") == false
		    mkdir(out_dir * "demandsDict")
		end

	    ta_data = load_ta_network_(out_dir, files_ID, month_w, day, instance1);

	    key_ = readstring(out_dir * "cross_validation_best_key/cross_validation_best_key_" * month_w * "_" * string(day) * "_" * instance * ".json")
	    key_ = JSON.parse(key_)
	    #key_ = "(7, 0.5, 1000.0, 1)"
	    obj_dict[day] = demandsDictFixed(demandsDict, flow_observ,link_vector, graph, ta_data, link_dic, day, gamma1, gamma2,
	                        out_dir, files_ID, month_w, instance, key_, free_flow_time, capacity, start_node, en_node, numZones, cnt)

		#end

		coeffs_dict_ = readstring(out_dir * "coeffs_dict_" * month_w * "_" * instance1 *".json")
		coeffs_dict_ = JSON.parse(coeffs_dict_)
		fcoeffs = coeffs_dict_[key_]
		polyDeg = length(fcoeffs)

		#for day in week_day_Apr_list

		demandsDict = readstring(out_dir * "demandsDict/demandsDictFixed$(day)_" * month_w * "_" * instance1 * ".json");
		demandsDict = JSON.parse(demandsDict)

		demandsDict_ = Dict()
		for key in keys(demandsDict)
		    key_2 = (parse(Int, split(split(key, ',')[1], '(')[2]), parse(Int, split(split(key, ',')[2], ')')[1]))
		    demandsDict_[key_2] = demandsDict[key]
		end

		#     demandsDict_

		     tapFlowDicDictP = Dict()
		     tapFlowVecDictP = Dict()
		     #tapFlowDicDict[day], tapFlowVecDict[day] = tapMSA(graph, ta_data, link_dic, demandsDict_, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones)

		#     tapFlowVecDict[day]

		    tapSocialFlowDicDict[day], tapSocialFlowVecDict[day] = tapMSASocial(demandsDict_, fcoeffs, ta_data, graph, link_dic, start_node, en_node, free_flow_time, capacity, numLinks, numZones);
	#demands, fcoeffs, graph, ta_data, link_dic, start_node, en_node, free_flow_time, capacity, numLinks, numIter=1000, tol=1e-6
		#     tapSocialFlowVecDict[day]

		#     flow_observ[:, day]

	        tapFlowVecDictP = readstring(out_dir * "demandsDict/tapFlowVecDict$(day)_" * month_w * "_" * instance * ".json")
	        tapFlowVecDictP = JSON.parse(tapFlowVecDictP)
	        k_ = keys(tapFlowVecDictP)
	        max_k = maximum(map(x->parse(Int,x),k_))
	        #println(max_k)
		    # PoA_dict[day] = socialObj(tapFlowVecDict[day]) / socialObj(tapSocialFlowVecDict[day])
	        user_sol_dict[day] = socialObj(flow_observ[:, cnt], free_flow_time, polyDeg, fcoeffs, capacity, numLinks) ;
	        #user_sol_dict[day] = socialObj(tapFlowVecDictP[string(max_k)], free_flow_time, polyDeg, fcoeffs, capacity, numLinks) ;
	        social_sol_dict[day] = socialObj(tapSocialFlowVecDict[day], free_flow_time, polyDeg, fcoeffs, capacity, numLinks);
		    PoA_dict[day] = user_sol_dict[day] / social_sol_dict[day];
		#end
		println("day $(day) finished...")
	    println("PoA for day $(day) and instance " * instance1 * " is " * string(PoA_dict[day]))
	end

	outfile =  open(out_dir * "PoA_dict_" * month_w * "_" * instance1 * ".json", "w")
	JSON.print(outfile, PoA_dict)
	close(outfile)

	outfile = open(out_dir * "tapSocialFlowVecDict_" * month_w * "_" * instance1 * ".json", "w")
	JSON.print(outfile, tapSocialFlowVecDict)
	close(outfile)

	outfile = open(out_dir * "cong_" * month_w * "_" * instance1 * ".json", "w")
	JSON.print(outfile, user_sol_dict)
	close(outfile)

	outfile = open(out_dir * "obj_dict" * month_w * "_" * instance1 * ".json", "w")
	JSON.print(outfile, obj_dict)
	close(outfile)
end

y = @parallel vcat for instan = ["AM", "MD", "PM" , "NT"]# function biblev(instance)
	biblev(instan);
end;

#=

#for day = 1: length(week_day_Apr_list)

    demandsDict = readstring(out_dir * "demandsDict/demandsDictFixed$(day)_" * month_w * "_" * instance * ".json");
    demandsDict = JSON.parse(demandsDict)

    demandsDict_ = Dict()
    for key in keys(demandsDict)
        key_ = (parse(Int, split(split(key, ',')[1], '(')[2]), parse(Int, split(split(key, ',')[2], ')')[1]))
        demandsDict_[key_] = demandsDict[key]
    end

#     demandsDict_

#     tapFlowDicDict = Dict()
#     tapFlowVecDict = Dict()
#     tapFlowDicDict[day], tapFlowVecDict[day] = tapMSA(demandsDict_, fcoeffs);

#     tapFlowVecDict[day]

    tapSocialFlowDicDict[day], tapSocialFlowVecDict[day] = tapMSASocial(demandsDict_, fcoeffs);

#     tapSocialFlowVecDict[day]

#     flow_observ[:, day]

    # PoA_dict[day] = socialObj(tapFlowVecDict[day]) / socialObj(tapSocialFlowVecDict[day])

    PoA_dict[day] = socialObj(flow_observ[:, day]) / socialObj(tapSocialFlowVecDict[day]);
    println(socialObj(flow_observ[:, day]))
#end

#PoA_dict

max_stars_repo_path: 05. TAP POA 2.jl

max_issues_repo_path: 05. TAP POA 2.jl

max_forks_repo_path: 05. TAP POA 2.jl


text:
# Suit of functions to create LCModels of Catalecticant matrices.
# For example, for ternary quartics:
#
# Σ = Catalecticant(2, 3)
#
# To compute the ML-degree of Σ
#
# ml_degree_witness(Σ; dual=true)
#

using LinearCovarianceModels
import HomotopyContinuation
const HC = HomotopyContinuation

function nvectors_sumingupm(n, m)
    ntuples = Iterators.product(ntuple(_ -> 0:m, n)...)
    return filter(I -> sum(I) == m, collect(ntuples))
end

function Catalecticantindices(r, s)
    rowindeces = collect.(nvectors_sumingupm(s, r))
    m = binomial(r + s - 1, s - 1)
    M = reduce(hcat, fill(rowindeces, m))
    return M + permutedims(M)
end

function Catalecticantmatrix(r, s)
    return [HC.Variable(:θ, I...) for I in Catalecticantindices(r, s)]
end

Catalecticant = LCModel ∘ Catalecticantmatrix

max_stars_repo_name: LauraBMo/InvertingCats

max_issues_repo_name: LauraBMo/InvertingCats

max_forks_repo_name: LauraBMo/InvertingCats


text:
module SearchLightCacheAdapter

using Genie, Logger, App, SearchLight, Util

type StorageCache <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  ### fields
  id::Nullable{SearchLight.DbId}
  name::String
  val::String
  expires_at::Int

  ### constructor
  StorageCache(;
    id = Nullable{SearchLight.DbId}(),
    name = "",
    val = "",
    expires_at = time_to_unixtimestamp()
  ) = new("storage_caches", "id", id, name, val, expires_at)
end

"""
    to_cache(key::Union{String,Symbol}, content::Any; dir = "") :: Void

Persists `content` onto the DB under the `key` key.
"""
function to_cache(key::Union{String,Symbol}, content::Any; dir = "") :: Void
  try
    io = IOBuffer()
    serialize(io, content)

    SearchLight.update_by_or_create!!(StorageCache(name = key, val = base64encode(io.data), expires_at = time_to_unixtimestamp()), :name)
  catch ex
    Logger.log("Error when serializing cache in $(@__FILE__):$(@__LINE__)", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    rethrow(ex)
  end

  nothing
end

"""
    from_cache(key::Union{String,Symbol}, expiration::Int) :: Nullable

Retrieves from cache the object stored under the `key` key if the `expiration` delta (in seconds) is in the future.
"""
function from_cache(key::Union{String,Symbol}, expiration::Int; dir = "") :: Nullable
  try
    cache_info = SearchLight.find_one_by!!(StorageCache, :name, key)

    App.config.log_cache && Logger.log("Found cache for $key", :info)

    cache_info.expires_at + expiration < time_to_unixtimestamp() && return Nullable()

    io = IOBuffer()
    io.data = base64decode(cache_info.val)
    io.size = size(io.data)[1]
    seekstart(io)
    cache = deserialize(io)

    Nullable(cache)
  catch ex
    Logger.log("Can't read cache", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    Nullable()
  end
end

"""
    purge(key::Union{String,Symbol}) :: Void

Removes the cache data stored under the `key` key.
"""
function purge(key::Union{String,Symbol}; dir = "") :: Void
  cache_info = SearchLight.find_one_by!!(StorageCache, :name, key)
  SearchLight.delete(cache_info)

  nothing
end

"""
    purge_all() :: Void

Removes all cached data.
"""
function purge_all(; dir = "") :: Void
  SearchLight.delete_all(StorageCache, truncate = true, reset_sequence = true)

  nothing
end

"""
    cache_path(key::Union{String,Symbol}; dir = "") :: String

Computes the path to a cache `key` based on current cache settings.
"""
function cache_path(key::Union{String,Symbol}; dir = "") :: String
  error("Not applicable!")
end

### Migrations module

module Migrations

module CreateTableStorageCaches

import Migration: create_table, column, column_id, add_index, drop_table

using Configuration, App

"""

"""
function up()
  create_table(App.config.cache_table) do
    [
    column_id()
    column(:name, :string, "UNIQUE", limit = 64)
    column(:val, :text)
    column(:expires_at, :integer)
    ]
  end

  add_index(App.config.cache_table, :expires_at)
end

function down()
  drop_table(App.config.cache_table)
end

end # module CreateTableStorageCaches

end # module Migrations

end