text:
module Poisson_examples
using FinEtools

function Poisson_FE_example()
    println("""

    Heat conduction example described by Amuthan A. Ramabathiran
    http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia:
    Unit square, with known temperature distribution along the boundary,
    and uniform heat generation rate inside.  Mesh of regular linear TRIANGLES,
    in a grid of 1000 x 1000 edges (2M triangles, 1M degrees of freedom).
    Version: 05/29/2017
    """
    )
    t0 = time()

    A = 1.0 # dimension of the domain (length of the side of the square)
    thermal_conductivity =  [i==j ? one(FFlt) : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
    Q = -6.0; # internal heat generation rate
    tempf(x) = (1.0 .+ x[:,1].^2 .+ 2.0 .* x[:,2].^2);#the exact distribution of temperature
    N = 1000;# number of subdivisions along the sides of the square domain

    println("Mesh generation")
    @time fens,fes =T3block(A, A, N, N)

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))

    println("Searching nodes  for BC")
    @time l1 = selectnode(fens; box=[0. 0. 0. A], inflate = 1.0/N/100.0)
    @time l2 = selectnode(fens; box=[A A 0. A], inflate = 1.0/N/100.0)
    @time l3 = selectnode(fens; box=[0. A 0. 0.], inflate = 1.0/N/100.0)
    @time l4 = selectnode(fens; box=[0. A A A], inflate = 1.0/N/100.0)
    List = vcat(l1, l2, l3, l4)
    @time setebc!(Temp, List, true, 1, tempf(geom.values[List,:])[:])
    @time applyebc!(Temp)
    @time numberdofs!(Temp)

    t1 = time()

    material = MatHeatDiff(thermal_conductivity)

    femm = FEMMHeatDiff(IntegData(fes, TriRule(1)), material)

    println("Conductivity")
    @time K = conductivity(femm, geom, Temp)
    println("Nonzero EBC")
    @time F2 = nzebcloadsconductivity(femm, geom, Temp);
    println("Internal heat generation")
    # function getsource!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    #   forceout[1] = Q; #heat source
    # end
    # fi = ForceIntensity(FFlt, getsource!);# alternative  specification
    fi = ForceIntensity(FFlt[Q]);
    @time F1 = distribloads(FEMMBase(IntegData(fes, TriRule(1))), geom, Temp, fi, 3);

    println("Factorization")
    @time K = cholesky(K)
    println("Solution of the factorized system")
    @time U = K\(F1+F2)
    scattersysvec!(Temp,U[:])

    println("Total time elapsed = $(time() - t0) [s]")
    println("Solution time elapsed = $(time() - t1) [s]")

    Error= 0.0
    for k=1:size(fens.xyz,1)
        Error = Error + abs.(Temp.values[k,1] - tempf(reshape(fens.xyz[k,:], (1,2)))[1])
    end
    println("Error =$Error")

    # File =  "a.vtk"
    # MeshExportModule.vtkexportmesh (File, fes.conn, [geom.values Temp.values], MeshExportModule.T3; scalars=Temp.values, scalars_name ="Temperature")

    true
end

function Poisson_FE_example_algo()
    A= 1.0
    thermal_conductivity =  [i==j ? one(FFlt) : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
    magn = -6.0; #heat source
    truetempf(x)=1.0 .+ x[1].^2 .+ 2.0 .* x[2].^2;
    N=20;

    println("""
    Heat conduction example described by Amuthan A. Ramabathiran
    http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia:
    Unit square, with known temperature distribution along the boundary,
    and uniform heat generation rate inside.  Mesh of regular TRIANGLES,
    in a grid of $N x $N edges.
    This version uses the FinEtools algorithm module.
    """
    )
    t0 = time()

    fens,fes = T3block(A, A, N, N)

    # Define boundary conditions
    l1  = selectnode(fens; box=[0. 0. 0. A], inflate = 1.0/N/100.0)
    l2  = selectnode(fens; box=[A A 0. A], inflate = 1.0/N/100.0)
    l3  = selectnode(fens; box=[0. A 0. 0.], inflate = 1.0/N/100.0)
    l4  = selectnode(fens; box=[0. A A A], inflate = 1.0/N/100.0)

    essential1 = FDataDict("node_list"=>vcat(l1, l2, l3, l4),
    "temperature"=>truetempf);
    material = MatHeatDiff(thermal_conductivity)
    femm = FEMMHeatDiff(IntegData(fes, TriRule(1)), material)
    region1 = FDataDict("femm"=>femm, "Q"=>magn)
    # Make model data
    modeldata= FDataDict("fens"=> fens,
    "regions"=>[region1],
    "essential_bcs"=>[essential1]);

    # Call the solver
    modeldata = FinEtools.AlgoHeatDiffModule.steadystate(modeldata)

    println("Total time elapsed = ",time() - t0,"s")

    geom=modeldata["geom"]
    Temp=modeldata["temp"]
    femm=modeldata["regions"][1]["femm"]

    function errfh(loc,val)
        exact = truetempf(loc)
        return ((exact[1]-val[1])*exact)
    end
    E = integratefieldfunction(femm, geom, Temp, errfh, 0.0, m=3)
    println("Error=$E")

end # Poisson_FE_example_algo

function Poisson_FE_example_csys_1()
    println("""
    Heat conduction example described by Amuthan A. Ramabathiran
    http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia:
    Unit square, with known temperature distribution along the boundary,
    and uniform heat generation rate inside.  Mesh of regular linear TRIANGLES,
    in a grid of 1000 x 1000 edges (2M triangles, 1M degrees of freedom).
    The material response is defined in a local coordinate system.
    Version: 05/29/2017
    """
    )
    t0 = time()

    A = 1.0 # dimension of the domain (length of the side of the square)
    thermal_conductivity =  [i==j ? one(FFlt) : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
    Q = -6.0; # internal heat generation rate
    function getsource!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        forceout[1] = Q; #heat source
    end
    tempf(x) = (1.0 .+ x[:,1].^2 .+ 2.0 .* x[:,2].^2);#the exact distribution of temperature
    N = 1000;# number of subdivisions along the sides of the square domain
    Rm=[-0.9917568452513019 -0.12813414805267656
    -0.12813414805267656 0.9917568452513019]
    Rm=[-0.8020689950104449 -0.5972313850116512
    -0.5972313850116512 0.8020689950104447]

    println("Mesh generation")
    @time fens,fes =T3block(A, A, N, N)

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))

    println("Searching nodes  for BC")
    @time l1 = selectnode(fens; box=[0. 0. 0. A], inflate = 1.0/N/100.0)
    @time l2 = selectnode(fens; box=[A A 0. A], inflate = 1.0/N/100.0)
    @time l3 = selectnode(fens; box=[0. A 0. 0.], inflate = 1.0/N/100.0)
    @time l4 = selectnode(fens; box=[0. A A A], inflate = 1.0/N/100.0)
    List = vcat(l1, l2, l3, l4)
    @time setebc!(Temp, List, true, 1, tempf(geom.values[List,:])[:])
    @time applyebc!(Temp)
    @time numberdofs!(Temp)

    t1 = time()

    material = MatHeatDiff(thermal_conductivity)

    femm = FEMMHeatDiff(IntegData(fes, TriRule(1)), CSys(Rm), material)

    println("Conductivity")
    @time K = conductivity(femm, geom, Temp)
    println("Nonzero EBC")
    @time F2 = nzebcloadsconductivity(femm, geom, Temp);
    println("Internal heat generation")
    fi = ForceIntensity(FFlt[Q]);
    @time F1 = distribloads(femm, geom, Temp, fi, 3);

    println("Factorization")
    @time K = cholesky(K)
    println("Solution of the factorized system")
    @time U = K\(F1+F2)
    scattersysvec!(Temp,U[:])

    println("Total time elapsed = $(time() - t0) [s]")
    println("Solution time elapsed = $(time() - t1) [s]")

    Error= 0.0
    for k=1:size(fens.xyz,1)
        Error = Error + abs.(Temp.values[k,1] - tempf(reshape(fens.xyz[k,:], (1,2)))[1])
    end
    println("Error =$Error")

    # File =  "a.vtk"
    # MeshExportModule.vtkexportmesh (File, fes.conn, [geom.values Temp.values], MeshExportModule.T3; scalars=Temp.values, scalars_name ="Temperature")

    true

end # Poisson_FE_example_csys_1

function Poisson_FE_Q4_example()
    println("""

    Heat conduction example described by Amuthan A. Ramabathiran
    http://www.codeproject.com/Articles/579983/Finite-Element-programming-in-Julia:
    Unit square, with known temperature distribution along the boundary,
    and uniform heat generation rate inside.  Mesh of regular four-node QUADRILATERALS,
    in a grid of 1000 x 1000 edges (1M quads, 1M degrees of freedom).
    Version: 05/29/2017
    """
    )
    t0 = time()

    A = 1.0
    thermal_conductivity =  [i==j ? one(FFlt) : zero(FFlt) for i=1:2, j=1:2]; # conductivity matrix
    function getsource!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        forceout[1] = -6.0; #heat source
    end
    tempf(x) = (1.0 .+ x[:,1].^2 .+ 2.0 .* x[:,2].^2);#the exact distribution of temperature
    N = 1000;

    println("Mesh generation")
    @time fens,fes = Q4block(A, A, N, N)

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))

    println("Searching nodes  for BC")
    @time l1 = selectnode(fens; box=[0. 0. 0. A], inflate = 1.0/N/100.0)
    @time l2 = selectnode(fens; box=[A A 0. A], inflate = 1.0/N/100.0)
    @time l3 = selectnode(fens; box=[0. A 0. 0.], inflate = 1.0/N/100.0)
    @time l4 = selectnode(fens; box=[0. A A A], inflate = 1.0/N/100.0)
    List = vcat(l1, l2, l3, l4);
    setebc!(Temp, List, true, 1, tempf(geom.values[List,:])[:])
    applyebc!(Temp)

    numberdofs!(Temp)

    t1 = time()

    m = MatHeatDiff(thermal_conductivity)
    femm = FEMMHeatDiff(IntegData(fes, GaussRule(2, 2)), m)

    println("Conductivity")
    @time K=conductivity(femm, geom, Temp)
    #Profile.print()

    println("Nonzero EBC")
    @time F2 = nzebcloadsconductivity(femm, geom, Temp);
    println("Internal heat generation")
    fi = ForceIntensity(FFlt, 1, getsource!);
    @time F1 = distribloads(femm, geom, Temp, fi, 3);

    println("Factorization")
    @time K = cholesky(K)
    println("Solution of the factorized system")
    @time U=  K\(F1+F2)
    scattersysvec!(Temp, U[:])

    println("Total time elapsed = $(time() - t0) [s]")
    println("Solution time elapsed = $(time() - t1) [s]")

    # using MeshExportModule

    # File =  "a.vtk"
    # MeshExportModule.vtkexportmesh (File, fes.conn, [geom.values Temp.values], MeshExportModule.Q4; scalars=Temp.values, scalars_name ="Temperature")

    Error= 0.0
    for k=1:size(fens.xyz,1)
        Error = Error + abs.(Temp.values[k,1] - tempf(reshape(fens.xyz[k,:], (1,2)))[1])
    end
    println("Error =$Error")

    true
end # Poisson_FE_Q4_example

function allrun()
    println("#####################################################")
    println("# Poisson_FE_example ")
    Poisson_FE_example()
    println("#####################################################")
    println("# Poisson_FE_example_algo ")
    Poisson_FE_example_algo()
    println("#####################################################")
    println("# Poisson_FE_example_csys_1 ")
    Poisson_FE_example_csys_1()
    println("#####################################################")
    println("# Poisson_FE_Q4_example ")
    Poisson_FE_Q4_example()
end # function allrun

end # module Poisson_examples


text:
module CassetteV2

using Mixtape
import Mixtape: CompilationContext, allow, transform
using CodeInfoTools
using BenchmarkTools

# I had to do it.

module Target

function baz(y::Float64)
    return y + 40.0
end

function foo(x::Float64)
    return baz(x + 20.0)
end

end

# Controls when the transform on lowered code applies.
# I don't want to apply the recursive transform in type inference.
# So when stacklevel > 1, don't apply.
# I'll just dispatch to `Mixtape.call` at runtime.
mutable struct Mix <: CompilationContext
    stacklevel::Int
end
Mix() = Mix(1)

# Allow the transform on our Target module.
allow(ctx::Mix, m::Module, fn, args...) = m == Target

# Basically == a Cassette context.
struct Context
    d::Dict
    Context() = new(Dict())
end

# Our version of overdub.
overdub(::Context, f, args...) = f(args...)
function overdub(ctx::Context, ::typeof(+), args...)
    ret = (+)(args...)
    ctx.d[(Base.:(+), args...)] = ret
    return ret
end

function overdub(ctx::Context, fn::Function, args...)
    ret, inner = call(fn, args...; ctx = Mix())
    Base.merge!(ctx.d, inner.d)
    return ret
end

# The transform inserts state, then wraps calls in (overdub).
# Then, anytime there's a return value --
# create a tuple of (ret, state)
# and returns that.
# Consider monadic lifting f: R -> R => trans => R -> (R, state).
swap(r, e) = e
function swap(r, e::Expr)
    e.head == :call || return e
    return Expr(:call, overdub, r, e.args[1:end]...)
end

function transform(mix::Mix, src, sig)
    b = CodeInfoTools.Builder(src)
    mix.stacklevel == 1 || return src
    q = push!(b, Expr(:call, Context))
    rets = Any[]
    for (v, st) in b
        b[v] = swap(q, st)
        st isa Core.ReturnNode && push!(rets, v => st)
    end
    for (n, ret) in rets
        v = insert!(b, n, Expr(:call, Base.tuple, ret.val, q))
        b[n] = Core.ReturnNode(v)
    end
    mix.stacklevel += 1
    return CodeInfoTools.finish(b)
end

Mixtape.@load_abi()
ret = call(Target.foo, 5.0; ctx = Mix(1))
ret = call(Target.foo, 5.0; ctx = Mix(1), optlevel = 1)
ret = call(Target.foo, 5.0; ctx = Mix(1), optlevel = 2)
ret = call(Target.foo, 5.0; ctx = Mix(1), optlevel = 3)
display(ret)
@btime call(Target.foo, 5.0; ctx = Mix(1))

src = emit(Target.foo, Tuple{Float64};
           ctx = Mix(1), opt = true)
display(src)
display(src.ssavaluetypes)

end # module

ext: jl


text:
mutable struct cpop_data
    n # number of variables
    nb # number of binary variables
    m # number of constraints
    numeq # number of equality constraints
    x # set of variables
    pop # polynomial optimization problem
    gb # Grobner basis
    leadsupp # leader terms of the Grobner basis
    supp # support data
    coe # coefficient data
    basis # monomial bases
    ksupp # extending support at the k-th step
    sb # sizes of different blocks
    numb # numbers of different blocks
    blocks # block structure
    cl # numbers of blocks
    blocksize # sizes of blocks
    GramMat # Gram matrix
    solver # SDP solver
    tol # tolerance to certify global optimality
    flag # 0 if global optimality is certified; 1 otherwise
end

"""
    opt,sol,data = tssos_first(pop, x, d; nb=0, numeq=0, quotient=true, basis=[],
    reducebasis=false, TS="block", merge=false, md=3, solver="Mosek", QUIET=false, solve=true,
    MomentOne=false, Gram=false, solution=false, tol=1e-4)

Compute the first step of the TSSOS hierarchy for constrained polynomial optimization with
relaxation order `d`.
If `quotient=true`, then exploit the quotient ring structure defined by the equality constraints.
Return the optimum, the (near) optimal solution (if `solution=true`) and other auxiliary data.

# Arguments
- `pop`: the vector of the objective function, inequality constraints, and equality constraints.
- `x`: the set of variables.
- `d`: the relaxation order of the moment-SOS hierarchy.
- `nb`: the number of binary variables in `x`.
- `md`: the tunable parameter for merging blocks.
- `numeq`: the number of equality constraints.
"""
function tssos_first(pop, x, d; nb=0, numeq=0, quotient=true, basis=[], reducebasis=false,
    TS="block", merge=false, md=3, solver="Mosek", QUIET=false, solve=true, MomentOne=false, Gram=false,
    solution=false, tol=1e-4)
    println("***************************TSSOS***************************")
    println("TSSOS is launching...")
    n = length(x)
    if quotient == true
        cpop = copy(pop)
        gb = cpop[end-numeq+1:end]
        cpop = cpop[1:end-numeq]
        if numeq > 0 && QUIET == false
            println("Starting to compute the Gröbner basis...")
            println("This might take time. You can set quotient=false to close it.")
        end
        SemialgebraicSets.gröbnerbasis!(gb)
        cpop[1] = rem(cpop[1], gb)
        lead = leadingmonomial.(gb)
        llead = length(lead)
        leadsupp = zeros(UInt8, n, llead)
        for i = 1:llead, j = 1:n
            @inbounds leadsupp[j,i] = MultivariatePolynomials.degree(lead[i], x[j])
        end
    else
        cpop = pop
        gb = []
        leadsupp = []
    end
    m = length(cpop)-1
    dg = zeros(Int, m)
    coe = Vector{Vector{Float64}}(undef, m+1)
    supp = Vector{Array{UInt8,2}}(undef, m+1)
    for k = 1:m+1
        mon = monomials(cpop[k])
        coe[k] = coefficients(cpop[k])
        supp[k] = zeros(UInt8,n,length(mon))
        for i = 1:length(mon), j = 1:n
            @inbounds supp[k][j,i] = MultivariatePolynomials.degree(mon[i],x[j])
        end
    end
    isupp = supp[1]
    for i = 2:m+1
        dg[i-1] = maxdegree(pop[i])
        isupp = [isupp supp[i]]
    end
    if basis == []
        basis = Vector{Array{UInt8,2}}(undef, m+1)
        basis[1] = get_basis(n, d, nb=nb, lead=leadsupp)
        for k = 1:m
            basis[k+1] = get_basis(n, d-Int(ceil(dg[k]/2)), nb=nb, lead=leadsupp)
        end
    end
    tsupp = [isupp bin_add(basis[1],basis[1],nb)]
    tsupp = sortslices(tsupp, dims=2)
    tsupp = unique(tsupp, dims=2)
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,_ = get_cblocks(m, tsupp, supp[2:end], basis, nb=nb, TS=TS, QUIET=QUIET, merge=merge, md=md)
    if reducebasis == true && quotient == false
        gsupp = get_gsupp(n, m, supp, basis[2:end], blocks[2:end], cl[2:end], blocksize[2:end], nb=nb)
        psupp = [supp[1] zeros(UInt8,n)]
        psupp = [psupp gsupp]
        basis[1],flag = reducebasis!(psupp, basis[1], blocks[1], cl[1], blocksize[1], nb=nb)
        if flag == 1
            tsupp = [isupp bin_add(basis[1],basis[1],nb)]
            tsupp = sortslices(tsupp, dims=2)
            tsupp = unique(tsupp, dims=2)
            blocks,cl,blocksize,sb,numb,_ = get_cblocks(m, tsupp, supp[2:end], basis, nb=nb, TS=TS, QUIET=QUIET, merge=merge, md=md)
        end
    end
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp,moment,GramMat = blockcpop(n, m, supp, coe, basis, blocks, cl, blocksize, nb=nb, numeq=numeq, gb=gb, x=x, lead=leadsupp, solver=solver, QUIET=QUIET, solve=solve, solution=solution, MomentOne=MomentOne, Gram=Gram)
    data = cpop_data(n, nb, m, numeq, x, pop, gb, leadsupp, supp, coe, basis, ksupp, sb, numb, blocks, cl, blocksize, GramMat, solver, tol, 1)
    sol = nothing
    if solution == true
        sol,gap,data.flag = extract_solutions(moment, opt, pop, x, numeq=numeq, tol=tol)
        if data.flag == 1
            if gap > 0.5
                sol = randn(n)
            end
            sol,ub,gap = refine_sol(opt, sol, data, QUIET=true)
            if gap != nothing
                if gap < tol
                    data.flag = 0
                else
                    println("Found a local optimal solution giving an upper bound: $ub and a relative optimality gap: $gap.")
                end
            end
        end
    end
    return opt,sol,data
end

function tssos_higher!(data::cpop_data; TS="block", merge=false, md=3, QUIET=false, solve=true, MomentOne=false, Gram=false, solution=false)
    n = data.n
    nb = data.nb
    m = data.m
    numeq = data.numeq
    x = data.x
    pop = data.pop
    gb = data.gb
    leadsupp = data.leadsupp
    supp = data.supp
    coe = data.coe
    basis = data.basis
    ksupp = data.ksupp
    sb = data.sb
    numb = data.numb
    blocks = data.blocks
    cl = data.cl
    blocksize = data.blocksize
    solver = data.solver
    tol = data.tol
    ksupp = sortslices(ksupp, dims=2)
    ksupp = unique(ksupp, dims=2)
    if QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,status = get_cblocks(m, ksupp, supp[2:end], basis, blocks=blocks, cl=cl, blocksize=blocksize, sb=sb, numb=numb, nb=nb, TS=TS, QUIET=QUIET, merge=merge, md=md)
    end
    opt = nothing
    sol = nothing
    if status == 1
        if QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
        end
        opt,ksupp,moment,GramMat = blockcpop(n, m, supp, coe, basis, blocks, cl, blocksize, nb=nb, numeq=numeq, gb=gb, x=x, lead=leadsupp, solver=solver, QUIET=QUIET, solve=solve, solution=solution, MomentOne=MomentOne, Gram=Gram)
        if solution == true
            sol,gap,data.flag = extract_solutions(moment, opt, pop, x, numeq=numeq, tol=tol)
            if data.flag == 1
                if gap > 0.5
                    sol = randn(n)
                end
                sol,ub,gap = refine_sol(opt, sol, data, QUIET=true)
                if gap != nothing
                    if gap < tol
                        data.flag = 0
                    else
                        println("Found a local optimal solution giving an upper bound: $ub and a relative optimality gap: $gap.")
                    end
                end
            end
        end
        data.ksupp = ksupp
        data.sb = sb
        data.numb = numb
        data.blocks = blocks
        data.cl = cl
        data.blocksize = blocksize
        data.GramMat = GramMat
    end
    return opt,sol,data
end

function get_gsupp(n, m, supp, gbasis, blocks, cl, blocksize; nb=0)
    gsupp = zeros(UInt8, n, sum(size(supp[k+1],2)*Int(sum(Int.(blocksize[k]).^2+blocksize[k])/2) for k=1:m))
    l = 1
    for k = 1:m, i = 1:cl[k], j = 1:blocksize[k][i], r = j:blocksize[k][i], s = 1:size(supp[k+1],2)
        @inbounds bi = bin_add(gbasis[k][:,blocks[k][i][j]], gbasis[k][:,blocks[k][i][r]], nb)
        @inbounds bi = bin_add(bi, supp[k+1][:,s], nb)
        @inbounds gsupp[:,l] = bi
        l += 1
    end
    return gsupp
end

function reducebasis!(supp, basis, blocks, cl, blocksize; nb=0)
    esupp = supp[:, all.(iseven, eachcol(supp))]
    init,flag,check = 0,0,0
    while init==0 || check>0
        init,check = 1,0
        tsupp = esupp
        for i = 1:cl
            if blocksize[i] > 1
                for j = 1:blocksize[i], r = j+1:blocksize[i]
                    @inbounds bi = bin_add(basis[:,blocks[i][j]], basis[:,blocks[i][r]], nb)
                    tsupp = [tsupp bi]
                end
            end
        end
        tsupp = unique(tsupp, dims=2)
        tsupp = sortslices(tsupp, dims=2)
        ltsupp = size(tsupp, 2)
        for i = 1:cl
            lo = blocksize[i]
            indexb = [k for k=1:lo]
            j = 1
            while lo >= j
                bi = bin_add(basis[:,blocks[i][indexb[j]]], basis[:,blocks[i][indexb[j]]], nb)
                Locb = bfind(tsupp, ltsupp, bi)
                if Locb == 0
                   check,flag = 1,1
                   deleteat!(indexb, j)
                   lo -= 1
                else
                   j += 1
                end
            end
            blocks[i] = blocks[i][indexb]
            blocksize[i] = lo
        end
    end
    if flag == 1
       indexb = blocks[1]
       for i = 2:cl
           indexb = append!(indexb, blocks[i])
       end
       sort!(indexb)
       unique!(indexb)
       return basis[:,indexb],flag
    else
       return basis,flag
    end
end

function get_cgraph(tsupp::Array{UInt8, 2}, supp::Array{UInt8, 2}, basis::Array{UInt8, 2}; nb=0)
    lb = size(basis, 2)
    G = SimpleGraph(lb)
    ltsupp = size(tsupp, 2)
    for i = 1:lb, j = i+1:lb
        r = 1
        while r <= size(supp, 2)
            bi = bin_add(basis[:,i], basis[:,j], nb)
            bi = bin_add(bi, supp[:,r], nb)
            if bfind(tsupp, ltsupp, bi) != 0
               break
            else
                r += 1
            end
        end
        if r <= size(supp, 2)
           add_edge!(G, i, j)
        end
    end
    return G
end

function get_cblocks(m, tsupp, supp, basis; blocks=[], cl=[], blocksize=[], sb=[], numb=[], nb=0,
    TS="block", QUIET=true, merge=false, md=3)
    if isempty(blocks)
        blocks = Vector{Vector{Vector{UInt16}}}(undef, m+1)
        blocksize = Vector{Vector{UInt16}}(undef, m+1)
        cl = Vector{UInt16}(undef, m+1)
    end
    if TS == false
        for k = 1:m+1
            lb = ndims(basis[k])==1 ? length(basis[k]) : size(basis[k],2)
            blocks[k] = [[i for i=1:lb]]
            blocksize[k] = [lb]
            cl[k] = 1
        end
        status = 1
        nsb = Int.(blocksize[1])
        nnumb = [1]
    else
        G = get_graph(tsupp, basis[1], nb=nb)
        if TS == "block"
            blocks[1] = connected_components(G)
            blocksize[1] = length.(blocks[1])
            cl[1] = length(blocksize[1])
        else
            blocks[1],cl[1],blocksize[1] = chordal_cliques!(G, method=TS, minimize=false)
            if merge == true
                blocks[1],cl[1],blocksize[1] = clique_merge!(blocks[1], d=md, QUIET=true)
            end
        end
        nsb = sort(Int.(unique(blocksize[1])), rev=true)
        nnumb = [sum(blocksize[1].== i) for i in nsb]
        if isempty(sb) || nsb!=sb || nnumb!=numb
            status = 1
            if QUIET == false
                println("------------------------------------------------------")
                println("The sizes of PSD blocks:\n$nsb\n$nnumb")
                println("------------------------------------------------------")
            end
            for k = 1:m
                G = get_cgraph(tsupp, supp[k], basis[k+1], nb=nb)
                if TS == "block"
                    blocks[k+1] = connected_components(G)
                    blocksize[k+1] = length.(blocks[k+1])
                    cl[k+1] = length(blocksize[k+1])
                else
                    blocks[k+1],cl[k+1],blocksize[k+1] = chordal_cliques!(G, method=TS, minimize=false)
                    if merge == true
                        blocks[k+1],cl[k+1],blocksize[k+1] = clique_merge!(blocks[k+1], d=md, QUIET=true)
                    end
                end
            end
        else
            status = 0
            if QUIET == false
                println("No higher TSSOS hierarchy!")
            end
        end
    end
    return blocks,cl,blocksize,nsb,nnumb,status
end

function blockcpop(n, m, supp, coe, basis, blocks, cl, blocksize; nb=0, numeq=0, gb=[],
    x=[], lead=[], solver="Mosek", QUIET=true, solve=true, solution=false, MomentOne=false, Gram=false)
    ksupp = zeros(UInt8, n, Int(sum(Int.(blocksize[1]).^2+blocksize[1])/2))
    k = 1
    for i = 1:cl[1], j = 1:blocksize[1][i], r = j:blocksize[1][i]
        @inbounds bi = bin_add(basis[1][:,blocks[1][i][j]], basis[1][:,blocks[1][i][r]], nb)
        @inbounds ksupp[:,k] = bi
        k += 1
    end
    objv = nothing
    moment = nothing
    GramMat = nothing
    if solve == true
        tsupp = ksupp
        if m > 0
            gsupp = get_gsupp(n, m, supp, basis[2:end], blocks[2:end], cl[2:end], blocksize[2:end], nb=nb)
            tsupp = [tsupp gsupp]
        end
        if MomentOne == true || solution == true
            tsupp = [tsupp get_basis(n, 2, nb=nb)]
        end
        tsupp = unique(tsupp, dims=2)
        if !isempty(gb)
            nsupp = zeros(UInt8, n)
            llead = size(lead, 2)
            for i = 1:size(tsupp, 2)
                if divide(tsupp[:,i], lead, n, llead)
                    _,temp,_ = reminder(tsupp[:,i], x, gb, n)
                    nsupp = [nsupp temp]
                else
                    nsupp = [nsupp tsupp[:,i]]
                end
            end
            tsupp = nsupp
        end
        tsupp = sortslices(tsupp, dims=2)
        tsupp = unique(tsupp, dims=2)
        ltsupp = size(tsupp, 2)
        if QUIET == false
            println("Assembling the SDP...")
            println("There are $ltsupp affine constraints.")
        end
        if solver == "Mosek"
            model = Model(optimizer_with_attributes(Mosek.Optimizer))
        elseif solver == "COSMO"
            model = Model(optimizer_with_attributes(COSMO.Optimizer, "max_iter" => 10000))
        elseif solver == "SDPT3"
            model = Model(optimizer_with_attributes(SDPT3.Optimizer))
        else
            @error "The solver is currently not supported!"
            return nothing,nothing,nothing,nothing
        end
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        cons = [AffExpr(0) for i=1:ltsupp]
        pos = Vector{Union{VariableRef,Symmetric{VariableRef}}}(undef, cl[1])
        for i = 1:cl[1]
            if MomentOne == true || solution == true
                pos0 = @variable(model, [1:n+1, 1:n+1], PSD)
                for j = 1:n+1, k = j:n+1
                    @inbounds bi = bin_add(basis[1][:,j], basis[1][:,k], nb)
                    if !isempty(gb) && divide(bi, lead, n, llead)
                        bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                        for l = 1:bi_lm
                            Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                            if j == k
                               @inbounds add_to_expression!(cons[Locb], bi_coe[l], pos0[j,k])
                            else
                               @inbounds add_to_expression!(cons[Locb], 2*bi_coe[l], pos0[j,k])
                            end
                        end
                    else
                        Locb = bfind(tsupp,ltsupp,bi)
                        if j == k
                           @inbounds add_to_expression!(cons[Locb], pos0[j,k])
                        else
                           @inbounds add_to_expression!(cons[Locb], 2, pos0[j,k])
                        end
                    end
                end
            end
            bs = blocksize[1][i]
            if bs == 1
               @inbounds pos[i] = @variable(model, lower_bound=0)
               @inbounds bi = bin_add(basis[1][:,blocks[1][i][1]], basis[1][:,blocks[1][i][1]], nb)
               if !isempty(gb) && divide(bi, lead, n, llead)
                   bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                   for l = 1:bi_lm
                       Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                       @inbounds add_to_expression!(cons[Locb], bi_coe[l], pos[i])
                   end
               else
                   Locb = bfind(tsupp,ltsupp,bi)
                   @inbounds add_to_expression!(cons[Locb], pos[i])
               end
            else
               @inbounds pos[i] = @variable(model, [1:bs, 1:bs], PSD)
               for j = 1:bs, r = j:bs
                   @inbounds bi = bin_add(basis[1][:,blocks[1][i][j]], basis[1][:,blocks[1][i][r]], nb)
                   if !isempty(gb) && divide(bi, lead, n, llead)
                       bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                       for l = 1:bi_lm
                           Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                           if j == r
                              @inbounds add_to_expression!(cons[Locb], bi_coe[l], pos[i][j,r])
                           else
                              @inbounds add_to_expression!(cons[Locb], 2*bi_coe[l], pos[i][j,r])
                           end
                       end
                   else
                       Locb = bfind(tsupp, ltsupp, bi)
                       if j == r
                          @inbounds add_to_expression!(cons[Locb], pos[i][j,r])
                       else
                          @inbounds add_to_expression!(cons[Locb], 2, pos[i][j,r])
                       end
                   end
               end
            end
        end
        gpos = Vector{Vector{Union{VariableRef,Symmetric{VariableRef}}}}(undef, m)
        for k = 1:m
            gpos[k] = Vector{Union{VariableRef,Symmetric{VariableRef}}}(undef, cl[k+1])
            for i = 1:cl[k+1]
                bs = blocksize[k+1][i]
                if bs == 1
                    if !isempty(gb) || k <= m-numeq
                        gpos[k][i] = @variable(model, lower_bound=0)
                    else
                        gpos[k][i] = @variable(model)
                    end
                    for s = 1:size(supp[k+1],2)
                        @inbounds bi = bin_add(basis[k+1][:,blocks[k+1][i][1]], basis[k+1][:,blocks[k+1][i][1]], nb)
                        @inbounds bi = bin_add(bi, supp[k+1][:,s], nb)
                        if !isempty(gb) && divide(bi, lead, n, llead)
                            bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                            for l = 1:bi_lm
                                Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                                @inbounds add_to_expression!(cons[Locb], coe[k+1][s]*bi_coe[l], gpos[k][i])
                            end
                        else
                            Locb = bfind(tsupp, ltsupp, bi)
                            @inbounds add_to_expression!(cons[Locb], coe[k+1][s], gpos[k][i])
                        end
                    end
                else
                    if !isempty(gb) || k <= m-numeq
                       gpos[k][i] = @variable(model, [1:bs, 1:bs], PSD)
                    else
                       gpos[k][i] = @variable(model, [1:bs, 1:bs], Symmetric)
                    end
                    for j = 1:bs, r = j:bs, s = 1:size(supp[k+1],2)
                        @inbounds bi = bin_add(basis[k+1][:,blocks[k+1][i][j]], basis[k+1][:,blocks[k+1][i][r]], nb)
                        @inbounds bi = bin_add(bi,supp[k+1][:,s], nb)
                        if !isempty(gb) && divide(bi, lead, n, llead)
                            bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                            for l = 1:bi_lm
                                Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                                if j == r
                                   @inbounds add_to_expression!(cons[Locb], coe[k+1][s]*bi_coe[l], gpos[k][i][j,r])
                                else
                                   @inbounds add_to_expression!(cons[Locb], 2*coe[k+1][s]*bi_coe[l], gpos[k][i][j,r])
                                end
                            end
                        else
                            Locb = bfind(tsupp, ltsupp, bi)
                            if j == r
                               @inbounds add_to_expression!(cons[Locb], coe[k+1][s], gpos[k][i][j,r])
                            else
                               @inbounds add_to_expression!(cons[Locb], 2*coe[k+1][s], gpos[k][i][j,r])
                            end
                        end
                    end
                end
            end
        end
        bc = zeros(ltsupp)
        for i = 1:size(supp[1],2)
            Locb = bfind(tsupp, ltsupp, supp[1][:,i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing,nothing,nothing
            else
               bc[Locb] = coe[1][i]
           end
        end
        @variable(model, lower)
        cons[1] += lower
        @constraint(model, con[i=1:ltsupp], cons[i]==bc[i])
        @objective(model, Max, lower)
        end
        if QUIET==false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
        if Gram == true
            GramMat = Vector{Vector{Union{Float64,Matrix{Float64}}}}(undef, m+1)
            GramMat[1] = [value.(pos[i]) for i = 1:cl[1]]
            for k = 1:m
                GramMat[k+1] = [value.(gpos[k][i]) for i = 1:cl[k+1]]
            end
        end
        if solution == true
            dual_var = -dual.(con)
            moment = zeros(Float64, n+1, n+1)
            for j = 1:n+1, k = j:n+1
                bi = bin_add(basis[1][:,j], basis[1][:,k], nb)
                if !isempty(gb) && divide(bi, lead, n, llead)
                    bi_lm,bi_supp,bi_coe = reminder(bi, x, gb, n)
                    moment[j,k] = 0
                    for l = 1:bi_lm
                        Locb = bfind(tsupp, ltsupp, bi_supp[:,l])
                        moment[j,k] += bi_coe[l]*dual_var[Locb]
                    end
                else
                    Locb = bfind(tsupp,ltsupp,bi)
                    moment[j,k] = dual_var[Locb]
                end
            end
            moment = Symmetric(moment,:U)
        end
    end
    return objv,ksupp,moment,GramMat
end


text:
__precompile__(false)

@static if Sys.is_windows() using WinRPM end

"""
TODO
"""
function find_system_gcc()
    # On Windows, check to see if WinRPM is installed, and if so, see if gcc is installed
    @static if Sys.is_windows()
        try
            winrpmgcc = joinpath(WinRPM.installdir,"usr","$(Sys.ARCH)-w64-mingw32", "sys-root","mingw","bin","gcc.exe")
            if success(`$winrpmgcc --version`)
                return winrpmgcc
            end
        end
    end

    # See if `gcc` exists
    @static if is_unix()
        try
            if success(`gcc -v`)
                return "gcc"
            end
        end
    end

    error( "GCC not found on system: " * (Sys.is_windows() ? "GCC can be installed via `Pkg.add(\"WinRPM\"); WinRPM.install(\"gcc\")`" : "" ))
end

"""
TODO
"""
function gcc_compile(gcc,file,libname,env)
  gcc = find_system_gcc()
  C_cmd = setenv(`$gcc -fPIC -O3 -msse3 -xc -shared $file -Wl,-Bstatic -lopengl32 -o $(libname * "." * Libdl.dlext)`, env)
  Base.run(C_cmd)

  if !(!isempty(libname) && Libdl.dlopen_e(libname) != C_NULL) # HAS LIB?
    error("Cannot open lib $(libname)!")
  end
end

"""
TODO
"""
function write_c_file(libname)
  # write the C code inside a raw string
  C_code = raw"""

	__declspec(dllexport) unsigned long NvOptimusEnablement = 0x00000001; // NVIDIA
	//__declspec(dllexport) int AmdPowerXpressRequestHighPerformance = 1; // AMD Radeon

  int add(unsigned int a, unsigned int b){
      return a+b;
  }
  """

  cfile = libname * ".c"

  open(cfile, "w") do f
      print( f, C_code)
  end

  cfile
end

"""
TODO
"""
function compiler_setPaths(gcc,env_path)
  # set paths
  binary_path = dirname(gcc)
  inlcude_path = joinpath(abspath(binary_path,"../"),"include")
  lib_path = joinpath(abspath(binary_path,"../"),"lib")

  ENV2 = deepcopy(ENV)
  ENV2["CPATH"] = ""
  ENV2["LIBRARY_PATH"] = ""

  ENV2["PATH"] *= ";" * env_path
  ENV2["PATH"] *= ";" * binary_path

  ENV2["CPATH"] *= ";" * inlcude_path

  ENV2["LIBRARY_PATH"] *= ";" * env_path
  ENV2["LIBRARY_PATH"] *= ";" * binary_path
  ENV2["LIBRARY_PATH"] *= ";" * lib_path

  ENV2
end

### COMPILE WITH GCC ###
"""
TODO
"""
function compileWithGCC()
  # obtain the environment variable
  env_path = abspath(dirname(@__FILE__))

  # it MUST be declared as a `const` for `ccall` to work
  global const clibname = joinpath(env_path, "../", "build", "mylibByGCC")
  global const optimizer = joinpath(env_path,"../", "build", "JuliaOptimizer")

  # write the C code inside a raw string
  cfile = write_c_file(clibname)

  # find gcc compiler
  gcc = find_system_gcc()

  # set paths
  ENVS = compiler_setPaths(gcc,env_path)

  # the path to store the C library file
  const Clib = joinpath(env_path,clibname)

  # compile the C code into a shared library
  gcc_compile(gcc,cfile,Clib,ENVS)
end

### LINK LIB FUNCTIONS ###
"""
TODO
"""
function createLoop(index, array, func)
  cfunc = (x)->func(unsafe_pointer_to_objref(x))
  cLoopFunc = cfunction(cfunc, Void, Tuple{Ptr{Void},})
  ccall((:createLoop, optimizer),
  Ptr{Void},
  (Cuint,Ptr{Ptr{Void}},Cuint,Ptr{Void}),
  index,pointer(array),length(array),cLoopFunc)
end

"""
TODO
"""
function loopByIndex(index)
  ccall((:loopByIndex, optimizer),Void,(Cuint,),(index))
end

"""
TODO
"""
function loopByObject(pointer)
  ccall((:loopByObject, optimizer),Void,(Ptr{Void},),(pointer))
end

"""
TODO
"""
function prepareStaticLoop(x,a)
  f = cfunction(x, Void, Tuple{Ptr{Float32},})
  ccall((:prepare, optimizer),Void,(Ptr{Void},Ptr{Ptr{Void}},Cuint),f,pointer(a),length(a))
end

"""
TODO
"""
function staticloop()
  ccall((:loop, optimizer),Void,())
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
#=
analyze natvsext benchmark results
see natvsext/README.md
=#

using Printf
using CSV
using DataFrames

include(joinpath(@__DIR__, "../../examples/Examples.jl"))
import Main.Examples

bench_file = joinpath(@__DIR__, "raw", "bench.csv")
output_dir = mkpath(joinpath(@__DIR__, "analysis"))
tex_dir = mkpath(joinpath(output_dir, "tex"))
stats_dir = mkpath(joinpath(output_dir, "stats"))
csv_dir = mkpath(joinpath(output_dir, "csvs"))

print_table_solvers =
    true # add solver results to tex tables
    # false # just put formulation sizes in tex tables

# uncomment examples to process
examples_params = Dict(
    #= natural formulations paper =#
    "densityest" => (
        [:m, :twok], [2, 3],
        [:SEP,], [:nu_nat, :n_nat, :n_SEP]
        ),
    "doptimaldesign" => (
        [:logdet, :k], [5, 1],
        # [:EP,], [:n_EP, :q_nat, :q_EP]
        [:EP, :SEP],
        Symbol[]
        # [:nu_nat, :n_nat, :q_nat, :nu_EP, :n_EP, :q_EP, :nu_SEP, :n_SEP, :q_SEP]
        ),
    "matrixcompletion" => (
        [:m, :k], [1, 2],
        [:EP, :SEP],
        Symbol[]
        # [:nu_nat, :n_nat, :p_nat, :q_nat, :nu_EP, :n_EP, :q_EP, :nu_SEP,
        # :n_SEP, :q_SEP]
        ),
    "matrixregression" => (
        [:m, :k], [2, 1],
        [:SEP,], [:n_SEP, :q_nat]
        ),
    "polymin" => (
        [:m, :k], [1, 2],
        [:SEP,], [:nu_nat, :n_nat, :q_SEP]
        ),
    "portfolio" => (
        [:k], [1],
        [:SEP,], Symbol[]
        ),
    "shapeconregr" => (
        [:m, :twok], [1, 5],
        [:SEP,], [:nu_nat, :n_nat, :n_SEP, :q_nat]
        ),
    #= WSOS cones paper =#
    "polynorm" => (
        [:L1, :n, :d, :m], [5, 1, 3, 4],
        [:SEP,], Symbol[]
        ),
    #= spectral function cones paper =#
    "centralpolymat" => (
        [:func, :m, :k], [3, 1, 2],
        [:nat, :ext], Symbol[]
        ),
    "classicalquantum" => (
        [:d], [1],
        [:nat, :ext], Symbol[]
        ),
    "experimentdesign" => (
        [:func, :d], [2, 1],
        [:nat, :natlog, :ext], Symbol[]
        ),
    "nonparametricdistr" => (
        [:func, :d], [2, 1],
        [:nat, :natlog, :vecext], Symbol[]
        ),
    )

@info("analyzing examples: $(keys(examples_params))")

function post_process()
    all_df = CSV.read(bench_file, DataFrame)
    replace!(all_df.extender, extender_map...)
    all_df[!, :nu] = convert.(Union{Int, Missing}, all_df[!, :nu])

    for (ex_name, ex_params) in examples_params
        println()
        ex_df = filter(t -> (t.example == ex_name) && (ismissing(t.extender) ||
            t.extender in string.(ex_params[3])), all_df)
        if isempty(ex_df)
            @info("no data for $ex_name with params: $ex_params")
            continue
        end

        @info("starting $ex_name with params: $ex_params")
        (ex_df_wide, inst_solvers) = make_wide_csv(ex_df, ex_name, ex_params)
        make_table_tex(ex_name, ex_params, ex_df_wide, inst_solvers)
        make_plot_csv(ex_name, ex_params, ex_df_wide, inst_solvers)
        @info("finished $ex_name")
    end

    println()
    @info("finished all")
end

extender_map = Dict(
    "ExpPSD" => "EP",
    "SOCExpPSD" => "SEP",
    )

status_map = Dict(
    "SetupModelCaughtError" => "m",
    "SetupModelKilledTime" => "m",
    "SetupModelKilledMemory" => "m",
    "Optimal" => "co",
    "NearOptimal" => "ne",
    "TimeLimit" => "tl",
    "IterationLimit" => "il",
    "SolveCheckKilledTime" => "tl",
    "SolveCheckKilledMemory" => "rl",
    "SolveCheckCaughtError" => "rl", # rare, seems to be memory error
    "SlowProgress" => "sp",
    "NumericalFailure" => "er",
    "SkippedSolveCheck" => "sk",
    )

residual_tol_satisfied(a, tol = 1e-5) =
    (all(isfinite, a) && (maximum(abs, a) < tol))

relative_tol_satisfied(a::T, b::T, tol::T = 1e-4) where {T <: Real} =
    (abs(a - b) / (1 + max(abs(a), abs(b))) < tol)

get_name(x::Any) = string(x)
get_name(x::Main.Examples.SpectralExtender) = Main.Examples.nat_name(x)

function make_wide_csv(ex_df, ex_name, ex_params)
    @info("making wide csv for $ex_name")
    inst_keys = ex_params[1]
    str(x) = split(x[2:(end - 1)], ", ")

    # add columns
    inst_ext_name(inst_set, ext) = (ismissing(ext) ? inst_set : ext)
    transform!(ex_df, [:inst_set, :extender] => ((x, y) ->
        inst_ext_name.(x, y)) => :inst_ext)

    inst_solver_name(inst_ext, solver) = (inst_ext * "_" * solver)
    transform!(ex_df,
        [:inst_ext, :solver] => ((x, y) -> inst_solver_name.(x, y)) =>
        :inst_solver,
        [:x_viol, :y_viol, :z_viol, :rel_obj_diff] => ByRow((res...) ->
        residual_tol_satisfied(coalesce.(res, NaN))) => :converged,
        [:status, :script_status] => ByRow((x, y) ->
        (y == "Success" ? status_map[x] : status_map[y])) => :status)

    for (name, pos) in zip(inst_keys, ex_params[2])
        transform!(ex_df, :inst_data => ByRow(x ->
            get_name(eval(Meta.parse(x))[pos])) => name)
    end

    # get solver combinations and reorder
    s_temp = unique(ex_df[:, :inst_solver])
    exts = string.(ex_params[3])
    inst_solvers = filter(s -> !any(startswith(s, e) for e in exts), s_temp)
    for e in exts
        s_e = filter(s -> startswith(s, e * "_"), s_temp)
        append!(inst_solvers, s_e)
    end
    @assert length(inst_solvers) == length(s_temp)
    @info("instance set - solver combinations: $inst_solvers")

    # check objectives if solver claims optimality
    for group_df in groupby(ex_df, inst_keys)
        # check all pairs of verified converged results
        co_idxs = findall((group_df[:, :status] .== "co") .&
            group_df[:, :converged])
        (length(co_idxs) < 2) && continue

        first_optval = group_df[co_idxs[1], :primal_obj]
        other_optvals = group_df[co_idxs[2:end], :primal_obj]
        if !all(relative_tol_satisfied.(other_optvals, first_optval))
            dat = group_df[!, :inst_data][1]
            println("$ex_name $dat primal optimal objective values disagree:")
            show(group_df[!, [:inst_set, :inst_num, :extender, :solver,
                :status, :primal_obj, :dual_obj, :rel_obj_diff]],
                show_row_number=false, eltypes=false)
            println("\n")
        end
    end

    unstacked_dims = [
        unstack(ex_df, inst_keys, :inst_ext, v, renamecols = x ->
        Symbol(v, :_, x), allowduplicates=true)
        for v in [:nu, :n, :p, :q]
        ]
    unstacked_res = [
        unstack(ex_df, inst_keys, :inst_solver, v, renamecols = x ->
        Symbol(v, :_, x), allowduplicates=true)
        for v in [:status, :converged, :iters, :solve_time]
        ]
    ex_df_wide = outerjoin(unstacked_dims..., unstacked_res..., on = inst_keys)
    (:func in inst_keys) && sort!(ex_df_wide, :func)

    CSV.write(joinpath(stats_dir, ex_name * "_wide.csv"), ex_df_wide)

    return (ex_df_wide, inst_solvers)
end

process_entry(::Missing) = "\$\\ast\$"

process_entry(::Missing, ::Missing) = "sk" # no data so instance was skipped

function process_entry(x::Float64)
    isnan(x) && return "\$\\ast\$"
    @assert x > 0
    if x < 10
        return @sprintf("%.1f", x)
    else
        return @sprintf("%.0f.", x)
    end
end

process_entry(st::String, converged::Bool) =
    (converged ? "\\underline{$(st)}" : st)

process_entry(x) = string(x)

function process_inst_solver(row, inst_solver)
    sep = " & "
    row_str = sep * process_entry(row[Symbol(:status_, inst_solver)],
        row[Symbol(:converged_, inst_solver)])
    row_str *= sep * process_entry(row[Symbol(:iters_, inst_solver)])
    row_str *= sep * process_entry(row[Symbol(:solve_time_, inst_solver)])
    return row_str
end

function make_table_tex(ex_name, ex_params, ex_df_wide, inst_solvers)
    @info("making table tex for $ex_name")
    inst_keys = ex_params[1]
    print_sizes = ex_params[4]

    sep = " & "
    ex_tex = open(joinpath(tex_dir, ex_name * "_table.tex"), "w")
    header_str = prod(string(s) * sep for s in vcat(inst_keys, print_sizes))
    if print_table_solvers
        header_str *= prod(s * " & & & " for s in inst_solvers)
    end
    header_str = replace(header_str[1:(end - 2)], "_" => " ") * "\\\\"
    println(ex_tex, header_str)

    for row in eachrow(ex_df_wide)
        row_str = process_entry(row[1])
        for i in 2:length(inst_keys)
            row_str *= sep * process_entry(row[i])
        end
        for s in print_sizes
            row_str *= sep * process_entry(row[s])
        end
        if print_table_solvers
            for inst_solver in inst_solvers
                row_str *= process_inst_solver(row, inst_solver)
            end
        end
        row_str *= " \\\\"
        println(ex_tex, row_str)
    end

    close(ex_tex)

    return
end

function transform_plot_cols(ex_df_wide, inst_solver)
    old_cols = Symbol.([:converged_, :solve_time_], inst_solver)
    transform!(ex_df_wide, old_cols => ByRow((x, y) ->
        ((!ismissing(x) && x) ? y : missing)) => inst_solver)
end

function make_plot_csv(ex_name, ex_params, ex_df_wide, inst_solvers)
    @info("making plot csv for $ex_name")
    inst_keys = ex_params[1]
    num_params = length(inst_keys)

    for inst_solver in inst_solvers
        transform_plot_cols(ex_df_wide, inst_solver)
    end

    plot_file_start = joinpath(csv_dir, ex_name * "_plot")
    axis_name = last(inst_keys)
    if num_params == 1
        success_df = select(ex_df_wide, axis_name, inst_solvers...)
        CSV.write(plot_file_start * ".csv", success_df)
    else
        group_name = first(inst_keys)
        success_df = select(ex_df_wide, axis_name, group_name, inst_solvers...)
        for (group_id, group_df) in pairs(groupby(success_df, group_name))
            CSV.write(plot_file_start * "_$(group_id[1]).csv",
                select(group_df, Not(group_name)))
        end
    end

    return
end

post_process()
;


text:
import Base: IdentityUnitRange, @propagate_inbounds

# from https://github.com/JuliaArrays/CustomUnitRanges.jl/blob/master/src/ZeroRange.jl
"""
    ZeroRange(n)
Defines an `AbstractUnitRange` corresponding to the half-open interval
[0,n), equivalent to `0:n-1` but with the lower bound guaranteed to be
zero by Julia's type system.
"""
struct ZeroRange{T<:Integer} <: AbstractUnitRange{T}
    len::T
    ZeroRange{T}(len) where {T} = new{T}(max(zero(T), len))
end
ZeroRange(len::T) where {T<:Integer} = ZeroRange{T}(len)

Base.length(r::ZeroRange) = r.len

Base.length(r::ZeroRange{T}) where {T<:Union{Int,Int64}} = T(r.len)

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.length(r::ZeroRange{T}) where {T<:smallint} = Int(r.len)
end

Base.first(r::ZeroRange{T}) where {T} = zero(T)
Base.last(r::ZeroRange{T})  where {T} = r.len-one(T)

function Base.iterate(r::ZeroRange{T}) where T
    r.len <= 0 && return nothing
    x = oftype(one(T) + one(T), zero(T))
    return (x,x)
end
function Base.iterate(r::ZeroRange{T}, i) where T
    x = i + one(T)
    return (x >= oftype(x, r.len) ? nothing : (x,x))
end

@inline function Base.getindex(v::ZeroRange{T}, i::Integer) where T
    @boundscheck ((i > 0) & (i <= length(v))) || Base.throw_boundserror(v, i)
    convert(T, i-1)
end

@inline function Base.getindex(r::ZeroRange{R}, s::AbstractUnitRange{S}) where {R,S<:Integer}
    @boundscheck checkbounds(r, s)
    R(first(s)-1):R(last(s)-1)
end

Base.intersect(r::ZeroRange, s::ZeroRange) = ZeroRange(min(r.len,s.len))

Base.promote_rule(::Type{ZeroRange{T1}},::Type{ZeroRange{T2}}) where {T1,T2} = ZeroRange{promote_type(T1,T2)}
Base.convert(::Type{ZeroRange{T}}, r::ZeroRange{T}) where {T<:Real} = r
Base.convert(::Type{ZeroRange{T}}, r::ZeroRange) where {T<:Real} = ZeroRange{T}(r.len)
ZeroRange{T}(r::ZeroRange) where {T} = convert(ZeroRange{T}, r)

Base.show(io::IO, r::ZeroRange) = print(io, typeof(r).name, "(", r.len, ")")
Base.show(io::IO, ::MIME"text/plain", r::ZeroRange) = show(io, r)

"""
    ZVector{T} <: AbstractVector{T}
Simple 0-based Vector. For a more generic solution, see OffsetArrays.jl.
"""
struct ZVector{T} <: AbstractVector{T}
    x::Vector{T}
end
ZVector() = ZVector(Vector{Any}(undef, 0))
ZVector(::UndefInitializer, m::Integer) = ZVector(Vector{Any}(undef, Int(m)))

Base.parent(A::ZVector) = A.x
Base.size(A::ZVector) = size(parent(A))
Base.size(A::ZVector, d) = size(parent(A), d)

@inline Base.axes(A::ZVector) = (IdentityUnitRange(ZeroRange(length(parent(A)))),)
@inline Base.axes(A::ZVector, d) = d == 1 ? IdentityUnitRange(ZeroRange(length(parent(A)))) : IdentityUnitRange(0:0)

Base.IndexStyle(::Type{T}) where {T<:ZVector} = IndexLinear()
Base.eachindex(::IndexLinear, A::ZVector) = axes(A, 1)

Base.length(A::ZVector) = prod(size(A))

@inline @propagate_inbounds function Base.getindex(A::ZVector{T}, i::Int) where {T}
    @boundscheck checkbounds(A, i)
    @inbounds ret = parent(A)[i + 1]
    ret
end

@inline @propagate_inbounds function Base.getindex(A::ZVector{T}, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(parent(A), (I .+ 1)...)
    @inbounds ret = parent(A)[(I .+ 1)...]
    ret
end

@inline @propagate_inbounds function Base.setindex!(A::ZVector{T}, val, i::Int) where {T}
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[i + 1] = val
    val
end

@inline @propagate_inbounds function Base.setindex!(A::ZVector{T}, val, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(parent(A), (I .+ 1)...)
    @inbounds parent(A)[(I .+ 1)...] = val
    val
end

indexlength(r::AbstractRange) = length(r)
indexlength(i::Integer) = i

ZeroRangeIUR{T} = IdentityUnitRange{ZeroRange{T}}
Base.similar(A::ZVector) = similar(A, eltype(A), size(A))
Base.similar(A::ZVector, ::Type{T}) where {T} = similar(A, T, size(A))
Base.similar(A::ZVector, dims::Dims{1}) = similar(A, eltype(A), dims)
Base.similar(A::ZVector, ::Type{T}, dims::Dims{1}) where {T} = ZVector{T}(undef, first(dims))
Base.similar(A::AbstractVector, ::Type{T}, inds::Tuple{ZeroRangeIUR,Vararg{ZeroRangeIUR}}) where {T} = ZVector(similar(A, T, map(indexlength, inds)))
Base.similar(::Type{T}, shape::Tuple{ZeroRangeIUR,Vararg{ZeroRangeIUR}}) where {T<:AbstractArray} = ZVector(T(undef, map(indexlength, shape)))


text:
using HDF5
using BSON

abstract type DataManager end

# save data, overwriting existing
save(self::DataManager, path::String, data::Dict) = raise(DomainError("save() not defined!"))

# save data, adding to existing
save!(self::DataManager, path::String, data::Dict) = raise(DomainError("save!() not defined!"))

# load data
load(self::DataManager, path::String) = raise(DomainError("load() not defined!"))

# get manager filetype extension
extension() = raise(DomainError("extension() not defined!"))

# ===============
# --- H D F 5 ---
# ===============

struct HDF5Manager <: DataManager end

# Saving/Loading data
function _save(self::HDF5Manager, path::String, data::Dict, writeMode::String)
    h5open(path, writeMode) do f
      for (k,v) in data
          write(f, k, v)
      end
    end
end

extension(self::HDF5Manager) = ".h5"

save!(self::HDF5Manager, path::String, data::Dict) = _save(self, path, data, "cw")
save(self::HDF5Manager, path::String, data::Dict) = _save(self, path, data, "w")

function load(self::HDF5Manager, path::String)
    data = Dict()
    h5open(path) do f
        keys = names(f)
        for k in keys
            data[k] = read(f[k])
        end
    end
    return data
end

# ===============
# --- B S O N ---
# ===============

struct BSONManager <: DataManager end

extension(self::BSONManager) = ".bson"

function save(self::BSONManager, path::String, data::Dict)
    bson(path, data)
end

function save!(self::BSONManager, path::String, data::Dict)
    try
        priorData = BSON.load(path)
        newData = merge(data,BSON.load(path))
        bson(path, newData)
    catch
        bson(path, data)
    end
end

function load(self::BSONManager, path::String)
    return BSON.load(path)
end


text:
using Combinatorics
using StatsBase

using Base.Iterators

function njet(::Val{N}, ::typeof(sin), x₀) where {N}
    (s, c) = sincos(x₀)
    Jet(x₀, s, tuple(take(cycle((c, -s, -c, s)), N)...))
end

function njet(::Val{N}, ::typeof(cos), x₀) where {N}
    (s, c) = sincos(x₀)
    Jet(x₀, s, tuple(take(cycle((-s, -c, s, c)), N)...))
end

function njet(::Val{N}, ::typeof(exp), x₀) where {N}
    exped = exp(x₀)
    Jet(x₀, exped, tuple(take(repeated(exped), N)...))
end

jeval(j, x) = j(x)
for f in (sin, cos, exp)
    function (∂☆ₙ::∂☆{N})(fb::ZeroBundle{N, typeof(f)}, x::TaylorBundle{N}) where {N}
        njet(Val{N}(), primal(fb), primal(x))(x)
    end
    function (∂⃖ₙ::∂⃖{N})(∂☆ₘ::∂☆{M}, fb::ZeroBundle{M, typeof(f)}, x::TaylorBundle{M}) where {N, M}
        ∂⃖ₙ(jeval, njet(Val{N+M}(), primal(fb), primal(x)), x)
    end
end

function (::Diffractor.∂☆new{N})(B::ATB{N, Type{T}}, args::ATB{N}...) where {N, T<:SArray}
    error("Should have intercepted the constructor")
end


text:
# test for subst.jl which has operator over substitution
using Test

#include("subst.jl")

@testset "apply subst" begin
 @test apply([],:(P(x)),[]) == :(P(x))
 @test apply([:x],:x,[:a]) == :a
 @test apply([:x],:(x),[:a]) == :(a)
 @test apply([:x],:(P(x)),[:a]) == :(P(a))
 @test apply([:x,:y],:(P(x,f(y))),[:a,:b]) == :(P(a,f(b)))

 @test apply([:x,:y],:(P(x,y)),[:y,:x]) == :(P(y,x))
 @test apply([:x,:y],:(P(x,y)),[:a,:x]) == :(P(a,x))
 @test apply([:x,:y],:(P(x,y)),[:a,:b]) == :(P(a,b))
 @test apply([:x,:y,:z],:(P(x,y,z)),[:y,:z,:x]) == :(P(y,z,x))
 @test apply([:x,:y,:z],:(P(x,y,z)),[:(f(y,z)),:(f(y,x)),:(f(g(x),h(y)))]) == :(P(f(y,z), f(y,x), f(g(x),h(y))))
end

@testset "apply2 subst" begin
 @test apply2([], [], []) == []
 @test apply2([:x], [:x], [:(f(y))]) == [:(f(y))]
 @test apply2([:x], [:(g(x))], [:(f(y))]) == [:(g(f(y)))]
 @test apply2([:x,:y], [:y,:x], [:x,:y]) == [:y,:x]
 @test apply2([:x,:y], [:(f(y,x)),:x], [:x,:a]) == [:(f(a,x)),:x]
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
## Taken from version 0.4
cld{T<:Integer }(x::T, y::T) = div(x,y)+(!signbit(x$y)&(rem(x,y)!=0))

function makedata!(y::Vector{Float64}, x::Matrix{Float64})
    n = length(y)
    size(x)[1] == n || error("x and y have incompatible dimensions")
    randn!(y)
    for i in 1:n
        x[i] = 1
    end
    if length(x) > n
        for i in (n+1):length(x)
            x[i] = randn()
        end
    end
end

## Let's just do one for forecast unbiasedness first.
function oosstat!(f::Vector{Float64}, coef::Vector{Float64},
                  y::Vector{Float64}, x::Matrix{Float64})
    R = length(y) - length(f)
    for t = R:(length(y) - 1)
        coef[:] = x[1:t,:] \ y[1:t]
        f[t-R+1] = y[t+1] - (x[t+1,:] * coef)[1]
    end
    mean(f)
end

function bootmean!(ffull::Vector{Float64}, coef::Vector{Float64},
                   y::Vector{Float64}, x::Matrix{Float64})
    coef[:] = x \ y
    ffull[:] = y - x * coef
    mean(ffull)
end

## Let's just do one for absolute forecast error.
function oosstat2!(f::Vector{Float64}, coef::Vector{Float64},
                  y::Vector{Float64}, x::Matrix{Float64})
    R = length(y) - length(f)
    for t = R:(length(y) - 1)
        coef[:] = x[1:t,:] \ y[1:t]
        f[t-R+1] = abs(y[t+1] - (x[t+1,:] * coef)[1])
    end
    mean(f)
end

function bootmean2!(ffull::Vector{Float64}, coef::Vector{Float64},
                   y::Vector{Float64}, x::Matrix{Float64})
    coef[:] = x \ y
    ffull[:] = abs(y - x * coef)
    mean(ffull)
end

function runmc!(mcresults, nboot, P, R, k, α,
                oosstat!, bootmean!, truemean)
    y = Array(Float64, P+R)
    x = Array(Float64, P+R, k)
    yboot = similar(y)
    xboot = similar(x)
    f = Array(Float64,  P)
    ffull = Array(Float64, P+R)
    coef = Array(Float64, k)
    bootvalue = Array(Float64, nboot)

    for i in 1:length(mcresults)
        makedata!(y, x)
        oosvalue = oosstat!(f, coef, y, x) - truemean
        bootnull = bootmean!(ffull, coef, y, x)

        for j in 1:nboot
            blockboot!(yboot, xboot, y, x, 8)
            @inbounds bootvalue[j] = oosstat!(f, coef, yboot, xboot) - bootnull
        end
        crit = quantile(bootvalue, [α/2, 1 - α/2])
        @inbounds mcresults[i] = oosvalue < crit[1] || oosvalue > crit[2]
    end
end

mcresults1 = Array(Bool, 600)
@time runmc!(mcresults1, 499, 120, 240, 5, 0.05, oosstat!, bootmean!, 0)
mcresults2 = Array(Bool, 1500)
@time runmc!(mcresults2, 499, 120, 240, 5, 0.05, oosstat2!, bootmean2!, sqrt(2/π))

max_stars_repo_name: grayclhn-econ/oosbootstrap

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: grayclhn-econ/oosbootstrap

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: grayclhn-econ/oosbootstrap

max_forks_repo_licenses: MIT Unlicense


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

mutable struct IoK8sApiAuthorizationV1ResourceRule <: SwaggerModel
    apiGroups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiGroups
    resourceNames::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resourceNames
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiAuthorizationV1ResourceRule(;apiGroups=nothing, resourceNames=nothing, resources=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1ResourceRule, Symbol("apiGroups"), apiGroups)
        setfield!(o, Symbol("apiGroups"), apiGroups)
        validate_property(IoK8sApiAuthorizationV1ResourceRule, Symbol("resourceNames"), resourceNames)
        setfield!(o, Symbol("resourceNames"), resourceNames)
        validate_property(IoK8sApiAuthorizationV1ResourceRule, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiAuthorizationV1ResourceRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiAuthorizationV1ResourceRule

const _property_map_IoK8sApiAuthorizationV1ResourceRule = Dict{Symbol,Symbol}(Symbol("apiGroups")=>Symbol("apiGroups"), Symbol("resourceNames")=>Symbol("resourceNames"), Symbol("resources")=>Symbol("resources"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiAuthorizationV1ResourceRule = Dict{Symbol,String}(Symbol("apiGroups")=>"Vector{String}", Symbol("resourceNames")=>"Vector{String}", Symbol("resources")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1ResourceRule }) = collect(keys(_property_map_IoK8sApiAuthorizationV1ResourceRule))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1ResourceRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1ResourceRule[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1ResourceRule }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1ResourceRule[property_name]

function check_required(o::IoK8sApiAuthorizationV1ResourceRule)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

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


text:
export detderivates

using BandedMatrices
using StaticArrays

const SymRealHerm{T} = Union{Symmetric{T,<:BandedMatrix},Hermitian{T,<:BandedMatrix}}
struct Workspace{T,M<:SymRealHerm{T},TQ,TQp}
    A::M
    B::M
    Q::TQ
    Qp::TQp
end

"""
    detderivates(A, B, s::Real)

For a symmetric real or Hermtian matriices `A` and `B` calculate determinant of
`A - s * B` respectively related values:

`κ`: number of zeros of characteristic polynomial `p(x) = det(A - x *  B)` below `s`

`η`: `p'(s) / p(s)`

`ζ`: `-(p'(s) / p(s))' = (p'(s) / p(s))^2 - p''(s) / p(s) `

!!! note
    For a zero of `p` with multplicity `r`,
    the Laguerre iteration formula is `s -= n / ( η + sqrt((n * ζ - η^2) * (n-r)/r) )`.

    This iteration has far better area- and order of convergence than
    for example Newton's iteration `s -= 1 / η`.
"""
function detderivates(A::M, B::M, s::Real) where {T,M<:SymRealHerm{T}}
    ws = make_workspace(A, B)
    detderivates!(ws, s)
end
function detderivates!(ws::Workspace{T,M}, s, L::Union{Nothing,AbstractArray{T,N}}=nothing) where {T,M<:SymRealHerm{T},N}
    A, B, Q, Qp = ws.A, ws.B, ws.Q, ws.Qp
    R = real(T)
    Z = zero(R)
    n = size(A, 1)
    ϵ = T(eps(real(T)))
    if isinf(s)
        return ifelse(s < 0, 0, n), Z, Z, Z, Z, Z, Z, Z, Z, Z, Z
    end
    κ = 0
    dξp = dξpp = η = zero(R)
    ξ, ξp, ξpp = zero(R), zero(R), zero(R)
    ζ = zero(R)
    Q, Qp = initQ!(Q, Qp, A, B, s)
    ld = zero(R)
    for i = 1:n
        ξ, ξp, ξpp = R(real(Q[1,1])), real(Qp[1,1,1]), real(Qp[1,1,2])
        ld += log(abs(ξ))
        if iszero(ξ)
            ξ = ϵ^1.5
        end
        dξp = ξp / ξ
        dξpp = ξpp / ξ
        #println("i=$i ξ=$ξ dξp=$dξp dξpp = $dξpp")
        if L isa AbstractArray
            L[1,i,1] = ξ
            if N >= 3
                L[1,i,2] = dξp
            end
            for j = 2:size(L,1)
                L[j,i,1] = Q[j,1] / ξ
                if N >= 3
                    L[j,i,2] = (Qp[j,1,1] - Q[j,1] * dξp ) / ξ
                end
            end
        end
        κ += ξ < 0
        dξ2 = dξp^2
        if abs(η + dξp) < abs(η) * 1e-5
            # println("Δη small: i=$i $(abs(η + dξp) / abs(η))")
        end
        if abs(dξpp) > dξ2 * 1000
            #dξpp = zero(R)
        end
        dζ = dξ2 - dξpp
        ζ += dζ
        if ζ < -η^2 / i
            # println("discriminant negative: i=$i ξ=$ξ $ξp $ξpp  η=$η ζ=$ζ")
            ζ = R(NaN) # -η^2 / i
        end
        η += dξp
        if  i == n && ζ > 1.5e10 * η^2
            # println("discriminant too big: i=$i ξ=$ξ $ξp $ξpp  η=$η ζ=$ζ")
            # ζ = 1.5 * η^2
        end
        updateQ!(ξ, Q, Qp, i)
        stepQ!(Q, Qp, i, A, B, s)
    end
    λ = laguerre1(η, ζ, n, 1)
    κ, η, ζ, ld, λ, dξp, dξpp, ξ, ξp, ξpp
end

function initQ!(Q, Qp, A, B, s)
    k = min(size(Q, 1), size(A, 1))
    W = eltype(Q)
    for j = 1:k
        for i = j:k
            Q[i,j] = W(A[i,j]) - W(B[i,j]) * s
            Qp[i,j,1] = -B[i,j]
            Qp[i,j,2] = 0
        end
    end
    Q, Qp
end

function stepQ!(Q, Qp, k, A, B, s)
    q = size(Q, 1)
    n = size(A, 1)
    W = eltype(Q)
    for j = 1:q-1
        for i = j:q-1
            Q[i,j] = Q[i+1,j+1]
            Qp[i,j,1] = Qp[i+1,j+1,1]
            Qp[i,j,2] = Qp[i+1,j+1,2]
        end
        if k + q <= n
            Q[q,j] = W(A[k+q,k+j]) - W(B[k+q,k+j]) * s
            Qp[q,j,1] = -B[k+q,k+j]
            Qp[q,j,2] = 0
        end
    end
    if k + q <= n
        Q[q,q] = W(A[k+q,k+q]) - W(B[k+q,k+q]) * s
        Qp[q,q,1] = -B[k+q,k+q]
        Qp[q,q,2] = 0
    end
    nothing
end

function updateQ!(ξ, Q, Qp, k)
    T = eltype(Qp)
    R = real(T)
    q = size(Q, 1)
    aa = real(ξ)
    a = R(aa)
    ap = real(Qp[1,1,1])
    app = real(Qp[1,1,2])
    bb = inv(aa)
    b = R(bb)
    bp = -ap * b * b
    bpp = -(app * b + 2 * ap * bp) * b
    for j = 2:q
        aaj = Q[j,1]'
        aj = T(aaj)
        apj = Qp[j,1,1]'
        appj = Qp[j,1,2]'
        aajb = aaj * bb
        ajb = aj * b
        apjb = apj * b + aj * bp
        appjb = appj * b + 2 * apj * bp + aj * bpp
        for i = j:q
            aai = Q[i,1]
            ai = T(aai)
            api = Qp[i,1,1]
            appi = Qp[i,1,2]
            aajib = aajb * aai
            ajib = ajb * ai
            aij = Q[i,j]
            Q[i,j] -= aajib
            apjib = apjb * ai + ajb * api
            Qp[i,j,1] -= apjib
            appjib = appjb * ai + 2*apjb * api + ajb * appi
            appji = Qp[i,j,2]
            if abs(appji - appjib) < 1e-8*max(abs(appji), abs(appjib))
                # println("inaccurate at $k($i, $j): $(abs(appji) \ (appji - appjib))")
            end
            Qp[i,j,2] -= appjib
        end
    end
    nothing
end

function make_workspace(A::M, B::M) where {T,M<:SymRealHerm{T}}
    k = min(max(bandwidth(A), bandwidth(B)) + 1, size(A, 1))
    if isbitstype(T)
        Q = MMatrix{k,k}(zeros(T, k, k))
        Qp = MArray{Tuple{k,k,2}}(zeros(T, k, k, 2))
    else
        Q = zeros(T, k, k)
        Qp = zeros(T, k, k, 2)
    end
    Workspace(A, B, Q, Qp)
end


text:
module ibm2ieee

const IBM32_SIGN = (UInt32(0x80000000))
const IBM32_EXPT = (UInt32(0x7f000000))
const IBM32_FRAC = (UInt32(0x00ffffff))
const IBM32_TOP =  (UInt32(0x00f00000))
const TIES_TO_EVEN_MASK32 = (UInt32(0xfffffffd))

const IBM64_SIGN = (UInt64(0x8000000000000000))
const IBM64_EXPT = (UInt64(0x7f00000000000000))
const IBM64_FRAC = (UInt64(0x00ffffffffffffff))
const IBM64_TOP =  (UInt64(0x00f0000000000000))
const TIES_TO_EVEN_MASK64 = (UInt64(0xfffffffffffffffd))

# Masks used for 3-bit and 32-bit rounded right-shifts of a 64-bit quantity.
#   The masks comprise the parity bit and the trailing bits for the shift.
const TIES_TO_EVEN_RSHIFT3 =  (UInt64(0x000000000000000b))
const TIES_TO_EVEN_RSHIFT32 = (UInt64(0x000000017fffffff))

const IEEE32_MAXEXP = 254     # Maximum biased exponent for finite values.
const IEEE32_INFINITY = (UInt32(0x7f800000))

# Constant used to count number of leading bits in a nonzero hex digit
# via `(BITCOUNT_MAGIC >> (hex_digit*2)) & 3U`. */
const BITCOUNT_MAGIC = (UInt32(0x000055af))

"""
IBM single-precision bit pattern to IEEE single-precision bit pattern.
"""

function ibm32ieee32(ibm::UInt32)::UInt32
    # Overflow and underflow possible; rounding can only happen
    # in subnormal cases.
    # int ibm_expt, ieee_expt, leading_zeros;
    # npy_uint32 ibm_frac, top_digit;
    # npy_uint32 ieee_sign, ieee_frac;

    ieee_sign = ibm & IBM32_SIGN
    ibm_frac = ibm & IBM32_FRAC

    # Quick return for zeros.
    if ibm_frac == 0
        return ieee_sign
    end

    # Reduce shift by 2 to get a binary exponent from the hex exponent.
    ibm_expt = Int((ibm & IBM32_EXPT) >> 22)

    # Normalise significand, then count leading zeros in top hex digit.
    top_digit = ibm_frac & IBM32_TOP
    while top_digit == 0
        ibm_frac <<= 4
        ibm_expt -= 4
        top_digit = ibm_frac & IBM32_TOP
    end
    leading_zeros = Int32((BITCOUNT_MAGIC >> (top_digit >> 19)) & 0x3)
    ibm_frac <<= leading_zeros

    # Adjust exponents for the differing biases of the formats: the IBM bias
    # is 64 hex digits, or 256 bits. The IEEE bias is 127. The difference is
    # -129; we get an extra -1 from the different significand representations
    # (0.f for IBM versus 1.f for IEEE), and another -1 to compensate for an
    # evil trick that saves an operation: on the fast path: we don't remove the
    # hidden 1-bit from the IEEE significand, so in the final addition that
    # extra bit ends in incrementing the exponent by one. */
    ieee_expt = ibm_expt - 131 - leading_zeros

    if ieee_expt >= 0 && ieee_expt < IEEE32_MAXEXP
        # normal case; no shift needed
        ieee_frac = ibm_frac;
        return ieee_sign + (UInt32(ieee_expt) << 23) + ieee_frac;
    elseif ieee_expt >= IEEE32_MAXEXP
        # overflow
        return ieee_sign + IEEE32_INFINITY;
    elseif ieee_expt >= -32
        # possible subnormal result; shift significand right by -ieee_expt
        # bits, rounding the result with round-ties-to-even.
        #
        # The round-ties-to-even code deserves some explanation: out of the
        # bits we're shifting out, let's call the most significant bit the
        # "rounding bit", and the rest the "trailing bits". We'll call the
        # least significant bit that *isn't* shifted out the "parity bit".
        # So for an example 5-bit shift right, we'd label the bits as follows:
        #
        # Before the shift:
        #
        #         ...xxxprtttt
        #                    ^
        #    msb            lsb
        #
        # After the shift:
        #
        #              ...xxxp
        #                    ^
        #    msb            lsb
        #
        # with the result possibly incremented by one.
        #
        # For round-ties-to-even, we need to round up if both (a) the rounding
        # bit is 1, and (b) either the parity bit is 1, or at least one of the
        # trailing bits is 1. We construct a mask that has 1-bits in the
        # parity bit position and trailing bit positions, and use that to
        # check condition (b). So for example in the 5-bit shift right, the
        # mask looks like this:
        #
        #         ...000101111 : mask
        #         ...xxxprtttt : ibm_frac
        #                    ^
        #    msb            lsb
        #
        # We then shift right by (shift - 1), add 1 if (ibm & mask) is
        # nonzero, and then do a final shift by one to get the rounded
        # value. Note that this approach avoids the possibility of
        # trying to shift a width-32 value by 32, which would give
        # undefined behaviour (see C99 6.5.7p3).
        #
        mask = ~(TIES_TO_EVEN_MASK32 << (-1 - ieee_expt))
        round_up = (ibm_frac & mask) > 0x0
        ieee_frac = ((ibm_frac >> (-1 - ieee_expt)) + round_up) >> 1
        return ieee_sign + ieee_frac
    else
        # underflow to zero
        return ieee_sign;
    end
end

"""
IBM double-precision bit pattern to IEEE single-precision bit pattern.
"""

function ibm64ieee32( ibm::UInt32)::UInt32
    # Overflow and underflow possible; rounding can occur in both
    # normal and subnormal cases.
    # int ibm_expt, ieee_expt, leading_zeros;
    # npy_uint64 ibm_frac, top_digit;
    # npy_uint32 ieee_sign, ieee_frac;

    ieee_sign = (ibm & IBM64_SIGN) >> 32
    ibm_frac = ibm & IBM64_FRAC

    # Quick return for zeros.
    if ibm_frac==0
        return ieee_sign
    end

    # Reduce shift by 2 to get a binary exponent from the hex exponent.
    ibm_expt = Int((ibm & IBM64_EXPT) >> 54)

    # Normalise significand, then count leading zeros in top hex digit.
    top_digit = ibm_frac & IBM64_TOP
    while top_digit == 0
        ibm_frac <<= 4
        ibm_expt -= 4
        top_digit = ibm_frac & IBM64_TOP
    end
    leading_zeros = Int((BITCOUNT_MAGIC >> (top_digit >> 51)) & 0x3)

    ibm_frac <<= leading_zeros;
    ieee_expt = ibm_expt - 131 - leading_zeros;

    if ieee_expt >= 0 && ieee_expt < IEEE32_MAXEXP
        # normal case; shift right 32, with round-ties-to-even
        round_up = (ibm_frac & TIES_TO_EVEN_RSHIFT32) > 0x0
        ieee_frac = ((npy_uint32)(ibm_frac >> 31) + round_up) >> 1
        return ieee_sign + ((npy_uint32)ieee_expt << 23) + ieee_frac
    elseif ieee_expt >= IEEE32_MAXEXP
        # overflow
        return ieee_sign + IEEE32_INFINITY
    elseif ieee_expt >= -32
        # possible subnormal; shift right with round-ties-to-even
        mask = ~(TIES_TO_EVEN_MASK64 << (31 - ieee_expt))
        round_up = (ibm_frac & mask) > 0x0
        ieee_frac = UInt32((ibm_frac >> (31 - ieee_expt)) + round_up) >> 1
        return ieee_sign + ieee_frac
    else
        # underflow to zero
        return ieee_sign
    end
end

"""
IBM single-precision bit pattern to IEEE double-precision bit pattern.

This is the simplest of the four cases: there's no need to check for
overflow or underflow, no possibility of subnormal output, and never
any rounding.

"""

function ibm32ieee64(ibm::UInt32)::UInt64
    ieee_sign = UInt64((ibm & IBM32_SIGN)) << 32
    ibm_frac = ibm & IBM32_FRAC

    # Quick return for zeros.
    if ibm_frac == 0
        return ieee_sign
    end

    # Reduce shift by 2 to get a binary exponent from the hex exponent.
    ibm_expt = UInt64((ibm & IBM32_EXPT) >> 22)

    # Normalise significand, then count leading zeros in top hex digit.
    top_digit = ibm_frac & IBM32_TOP
    while top_digit == 0
        ibm_frac <<= 4
        ibm_expt -= 4
        top_digit = ibm_frac & IBM32_TOP
    end
    leading_zeros = Int32((BITCOUNT_MAGIC >> (top_digit >> 19)) & 0x3)
    ibm_frac <<= leading_zeros

    # Adjust exponents for the differing biases of the formats: the IBM bias
    # is 64 hex digits, or 256 bits. The IEEE bias is 1023. The difference is
    # 767; we get an extra -1 from the different significand representations
    # (0.f for IBM versus 1.f for IEEE), and another -1 to compensate for an
    # evil trick that saves an operation: we don't remove the hidden 1-bit
    # from the IEEE significand, so in the final addition that extra bit ends
    # in incrementing the exponent by one. */
    ieee_expt = ibm_expt + 765 - leading_zeros
    ieee_frac = UInt64(ibm_frac) << (29 + leading_zeros)
    ieee_frac = UInt64(ibm_frac) << (29)
    return ieee_sign + (UInt64(ieee_expt) << 52) + ieee_frac
end

"""
IBM double-precision bit pattern to IEEE double-precision bit pattern.
"""

function ibm64ieee64(ibm::UInt64)::UInt64
    # No overflow or underflow possible, but the precision of the
    # IBM double-precision format exceeds that of its IEEE counterpart,
    # so we'll frequently need to round. */
    # int ibm_expt, ieee_expt, leading_zeros;
    # npy_uint64 ibm_frac, top_digit;
    # npy_uint64 ieee_sign, ieee_frac, round_up;

    ieee_sign = ibm & IBM64_SIGN
    ibm_frac = ibm & IBM64_FRAC

    # Quick return for zeros.
    if ibm_frac == 0
        return ieee_sign;
    end

    # Reduce shift by 2 to get a binary exponent from the hex exponent.
    ibm_expt = Int((ibm & IBM64_EXPT) >> 54);

    # Normalise significand, then count leading zeros in top hex digit.
    top_digit = ibm_frac & IBM64_TOP
    while top_digit == 0
        ibm_frac <<= 4
        ibm_expt -= 4
        top_digit = ibm_frac & IBM64_TOP
    end
    leading_zeros = Int((BITCOUNT_MAGIC >> (top_digit >> 51)) & 0x3)

    ibm_frac <<= leading_zeros
    ieee_expt = ibm_expt + 765 - leading_zeros

    # Right-shift by 3 bits (the difference between the IBM and IEEE
    # significand lengths), rounding with round-ties-to-even.
    round_up = (ibm_frac & TIES_TO_EVEN_RSHIFT3) > 0x0;
    ieee_frac = ((ibm_frac >> 2) + round_up) >> 1
    return ieee_sign + (UInt64(ieee_expt) << 52) + ieee_frac;
end

end # module


text:
function has_quadratic(model::LinQuadOptimizer)
    return model.obj_type == QuadraticObjective ||
        length(cmap(model).q_less_than) > 0 ||
        length(cmap(model).q_greater_than) > 0 ||
        length(cmap(model).q_equal_to) > 0
end

#=
    Optimize the model
=#

function MOI.optimize!(model::LinQuadOptimizer)
    # reset storage
    fill!(model.variable_primal_solution, NaN)
    fill!(model.variable_dual_solution, NaN)
    fill!(model.constraint_primal_solution, NaN)
    fill!(model.constraint_dual_solution, NaN)
    model.primal_status = MOI.UnknownResultStatus
    model.dual_status   = MOI.UnknownResultStatus
    model.primal_result_count = 0
    model.dual_result_count = 0

    start_time = time()
    if has_integer(model)
        solve_mip_problem!(model)
    elseif has_quadratic(model)
        solve_quadratic_problem!(model)
    else
        solve_linear_problem!(model)
    end
    model.solvetime = time() - start_time

    # termination_status
    model.termination_status = get_termination_status(model)
    model.primal_status = get_primal_status(model)
    model.dual_status = get_dual_status(model)

    if model.primal_status in [MOI.FeasiblePoint, MOI.InfeasiblePoint]
        get_variable_primal_solution!(model, model.variable_primal_solution)
        get_linear_primal_solution!(model, model.constraint_primal_solution)
        if has_quadratic(model)
            get_quadratic_primal_solution!(model, model.qconstraint_primal_solution)
        end
        model.primal_result_count = 1
    elseif model.primal_status == MOI.InfeasibilityCertificate
        get_unbounded_ray!(model, model.variable_primal_solution)
        model.primal_result_count = 1
    end
    if model.dual_status in [MOI.FeasiblePoint, MOI.InfeasiblePoint]
        get_variable_dual_solution!(model, model.variable_dual_solution)
        get_linear_dual_solution!(model, model.constraint_dual_solution)
        if has_quadratic(model)
            get_quadratic_dual_solution!(model, model.qconstraint_dual_solution)
        end
        model.dual_result_count = 1
    elseif model.dual_status == MOI.InfeasibilityCertificate
        get_farkas_dual!(model, model.constraint_dual_solution)
        model.dual_result_count = 1
    end

    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense
        model.constraint_dual_solution *= -1
        model.variable_dual_solution *= -1
    end
end

#=
    Result Count
=#
MOI.canget(::LinQuadOptimizer, ::MOI.ResultCount) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ResultCount)
    return max(model.primal_result_count, model.dual_result_count)
end

#=
    Termination status
=#
MOI.canget(::LinQuadOptimizer, ::MOI.TerminationStatus) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.TerminationStatus)
    return model.termination_status
end

#=
    Primal status
=#
function MOI.canget(model::LinQuadOptimizer,  status::MOI.PrimalStatus)
    if status.N != 1
        error("Multiple primal statuses are not supported.")
    end
    return model.primal_result_count >= status.N
end
function MOI.get(model::LinQuadOptimizer, ::MOI.PrimalStatus)
    return model.primal_status
end

#=
    Dual status
=#
function MOI.canget(model::LinQuadOptimizer,  status::MOI.DualStatus)
    if status.N != 1
        error("Multiple dual statuses are not supported.")
    end
    return model.dual_result_count >= status.N
end
function MOI.get(model::LinQuadOptimizer, ::MOI.DualStatus)
    return model.dual_status
end

#=
    Objective Value
=#
function MOI.canget(::LinQuadOptimizer, attr::MOI.ObjectiveValue)
    return attr.resultindex == 1
end
function MOI.get(model::LinQuadOptimizer, attr::MOI.ObjectiveValue)
    if attr.resultindex == 1
        # Note: we add m.objective_constant here to account for any constant
        # term which was not passed to the solver itself (and which therefore
        # would not be accounted for in `get_objective_value(m)`. We do *not*
        # call `get_constant_objective(m)` because that would also pull any
        # constants which were passed to the solver, resulting those constants
        # being counted twice. This confusion will be alleviated when all LQOI
        # solvers implement `get_constant_objective()` and
        # `set_constant_objective!()` by actually passing the relevant constants
        # to the solvers, at which point we can just get rid of
        # m.objective_constant entirely.
        return get_objective_value(model) + model.objective_constant
    else
        error("Unable to access multiple objective values")
    end
end

#=
    Variable Primal solution
=#
MOI.canget(::LinQuadOptimizer, ::MOI.VariablePrimal, ::Type{VarInd}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.VariablePrimal, index::VarInd)
    column = get_column(model, index)
    return model.variable_primal_solution[column]
end

MOI.canget(::LinQuadOptimizer, ::MOI.VariablePrimal, ::Type{Vector{VarInd}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.VariablePrimal, indices::Vector{VarInd})
    MOI.get.(Ref(model), Ref(MOI.VariablePrimal()), indices)
end

#=
    Variable Dual solution
=#
"""
    is_binding(set, value::Float64)

Return true if `value` is an extreme point of the set `set`.
"""
is_binding(set::LE, value::Float64) = isapprox(set.upper, value)
is_binding(set::GE, value::Float64) = isapprox(set.lower, value)
is_binding(set::EQ, value::Float64) = isapprox(set.value, value)
is_binding(set::IV, value::Float64) = isapprox(set.lower, value) || isapprox(set.upper, value)

function MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintDual, ::Type{SVCI{S}}) where S <: LinSets
    return true
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintDual, index::SVCI{<: LinSets})
    column = get_column(model, model[index])
    # the variable reduced cost is only the constriant dual if the bound is active.
    set = MOI.get(model, MOI.ConstraintSet(), index)
    primal_value = model.variable_primal_solution[column]
    if is_binding(set, primal_value)
        return model.variable_dual_solution[column]
    else
        return 0.0
    end
end

function MOI.canget(s::LinQuadOptimizer, ::MOI.ConstraintDual,
                    ::Type{VVCI{S}}) where S <: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}
    return true
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintDual, index::VVCI{<: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}})
    return [model.constraint_dual_solution[row] for row in model[index]]
end

#=
    Variable Bound Primal solution
=#

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintPrimal, ::Type{SVCI{S}}) where S <: LinSets = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintPrimal, index::SVCI{<: LinSets})
    column = get_column(model, model[index])
    return model.variable_primal_solution[column]
end

function MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintPrimal,
                   ::Type{VVCI{S}}) where S <: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}
    return true
end

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintPrimal,
                 index::VVCI{<: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}})
    return [model.constraint_primal_solution[row] for row in model[index]]
end

#=
    Constraint Primal solution
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintPrimal, ::Type{LCI{S}}) where S <: LinSets = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintPrimal, index::LCI{<: LinSets})
    row = model[index]
    return model.constraint_primal_solution[row] + model.constraint_constant[row]
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintPrimal, ::Type{QCI{S}}) where S <: LinSets = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintPrimal, index::QCI{<: LinSets})
    row = model[index]
    return model.qconstraint_primal_solution[row]
end

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintPrimal, ::Type{VLCI{S}}) where S <: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives} = true

# vector valued constraint duals
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintPrimal, index::VLCI{<: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}})
    row = model[index]
    return model.constraint_primal_solution[row] + model.constraint_constant[row]
end

#=
    Constraint Dual solution
=#

__assert_dual_sense__(::LCI{LE}, dual) = @assert dual <= 0.0
__assert_dual_sense__(::LCI{GE}, dual) = @assert dual >= 0.0
__assert_dual_sense__(::LCI{IV}, dual) = nothing
__assert_dual_sense__(::LCI{EQ}, dual) = nothing

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintDual, ::Type{LCI{S}}) where S <: LinSets = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintDual, index::LCI{<: LinSets})
    row = model[index]
    dual = model.constraint_dual_solution[row]
    __assert_dual_sense__(index, dual)
    return dual
end

MOI.canget(model::LinQuadOptimizer, ::MOI.ConstraintDual, ::Type{QCI{S}}) where S <: LinSets = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintDual, index::QCI{<: LinSets})
    row = model[index]
    return model.qconstraint_dual_solution[row]
end

# vector valued constraint duals
MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintDual, ::Type{VLCI{S}}) where S <: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives} = true

function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintDual, index::VLCI{<: Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}})
    rows = model[index]
    return model.constraint_dual_solution[rows]
end

#=
    Solution Attributes
=#

MOI.supports(::LinQuadOptimizer, ::MOI.ObjectiveBound) = true
MOI.canget(::LinQuadOptimizer, ::MOI.ObjectiveBound) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ObjectiveBound)
    return get_objective_bound(model)
end

MOI.supports(::LinQuadOptimizer, ::MOI.RelativeGap) = true
MOI.canget(::LinQuadOptimizer, ::MOI.RelativeGap) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.RelativeGap)
    return get_relative_mip_gap(model)
end

MOI.supports(::LinQuadOptimizer, ::MOI.SolveTime) = true
MOI.canget(::LinQuadOptimizer, ::MOI.SolveTime) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.SolveTime)
    return model.solvetime
end

MOI.supports(::LinQuadOptimizer, ::MOI.SimplexIterations) = true
MOI.canget(::LinQuadOptimizer, ::MOI.SimplexIterations) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.SimplexIterations)
    return get_iteration_count(model)
end

MOI.supports(::LinQuadOptimizer, ::MOI.BarrierIterations) = true
MOI.canget(::LinQuadOptimizer, ::MOI.BarrierIterations) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.BarrierIterations)
    return get_barrier_iterations(model)
end

MOI.supports(::LinQuadOptimizer, ::MOI.NodeCount) = true
MOI.canget(::LinQuadOptimizer, ::MOI.NodeCount) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.NodeCount)
    return get_node_count(model)
end

MOI.supports(::LinQuadOptimizer, ::MOI.RawSolver) = true
MOI.canget(::LinQuadOptimizer, ::MOI.RawSolver) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.RawSolver)
    return model
end


text:
hessian_pos(V::SitePotential, at::AbstractAtoms) =
      _precon_or_hessian_pos(V, at, evaluate_dd!)

# implementation of a generic assembly of a global block-hessian from
# local site-hessians
function _precon_or_hessian_pos(V, at::AbstractAtoms{T}, hfun) where {T}
   nlist = neighbourlist(at, cutoff(V))
   maxN = maxneigs(nlist)
   hEs = zeros(JMat{T}, maxN, maxN)
   tmp = alloc_temp_dd(V, maxN)
   I, J, Z = Int[], Int[], JMat{T}[]
   # a crude size hint
   for C in (I, J, Z); sizehint!(C, npairs(nlist)); end
   for (i, neigs, R) in sites(nlist)
      nneigs = length(neigs)
      # [1] the "off-centre" component of the hessian:
      # h[a, b] = ∂_{Ra} ∂_{Rb} V     (this is a nneigs x nneigs block-matrix)
      fill!(hEs, zero(JMat{T}))
      hfun(hEs, tmp, V, R)
      for a = 1:nneigs, b = 1:nneigs
         # if norm(hEs[a,b], Inf) > 0
            push!(I, neigs[a])
            push!(J, neigs[b])
            push!(Z, hEs[a,b])
         # end
      end

      # [2] the ∂_{Ri} ∂_{Ra} terms
      # hib = ∂_{Ri} ∂_{Rb} V = - ∑_a ∂_{Ra} ∂_{Rb} V
      # also at the same time we pre-compute the centre-centre term:
      #    hii = ∂_{Ri} ∂_{Ri} V = - ∑_a ∂_{Ri} ∂_{Ra} V
      hii = zero(JMat{T})
      for b = 1:nneigs
         hib = -sum( hEs[a, b] for a = 1:nneigs )
         # if norm(hib, Inf) > 0
            hii -= hib
            append!(I, (i,         neigs[b] ))
            append!(J, (neigs[b],  i        ))
            append!(Z, (hib,       hib'     ))
         # end
      end

      # and finally add the  ∂_{Ri}^2 term, which is precomputed above
      # if norm(hii, Inf) > 0
         push!(I, i)
         push!(J, i)
         push!(Z, hii)
      # end
   end
   return sparse(I, J, Z, length(at), length(at))
end

# implementation of a generic assembly of a global block-hessian from
# local site-hessians
function _precon_or_hessian_pos(V::MSitePotential, at::AbstractAtoms{T}, hfun) where {T}
   nlist = neighbourlist(at, cutoff(V))
   maxN = maxneigs(nlist)
   hEs = zeros(JMat{T}, maxN, maxN)
   tmp = alloc_temp_dd(V, maxN)
   I, J, Z = Int[], Int[], JMat{T}[]
   # a crude size hint
   for C in (I, J, Z); sizehint!(C, npairs(nlist)); end
   for (i, neigs, R) in sites(nlist)
      z = at.Z[neigs]
      z0 = at.Z[i]
      nneigs = length(neigs)
      # [1] the "off-centre" component of the hessian:
      # h[a, b] = ∂_{Ra} ∂_{Rb} V     (this is a nneigs x nneigs block-matrix)
      fill!(hEs, zero(JMat{T}))
      hfun(hEs, tmp, V, R, z, z0)
      for a = 1:nneigs, b = 1:nneigs
         # if norm(hEs[a,b], Inf) > 0
            push!(I, neigs[a])
            push!(J, neigs[b])
            push!(Z, hEs[a,b])
         # end
      end

      # [2] the ∂_{Ri} ∂_{Ra} terms
      # hib = ∂_{Ri} ∂_{Rb} V = - ∑_a ∂_{Ra} ∂_{Rb} V
      # also at the same time we pre-compute the centre-centre term:
      #    hii = ∂_{Ri} ∂_{Ri} V = - ∑_a ∂_{Ri} ∂_{Ra} V
      hii = zero(JMat{T})
      for b = 1:nneigs
         hib = -sum( hEs[a, b] for a = 1:nneigs )
         # if norm(hib, Inf) > 0
            hii -= hib
            append!(I, (i,         neigs[b] ))
            append!(J, (neigs[b],  i        ))
            append!(Z, (hib,       hib'     ))
         # end
      end

      # and finally add the  ∂_{Ri}^2 term, which is precomputed above
      # if norm(hii, Inf) > 0
         push!(I, i)
         push!(J, i)
         push!(Z, hii)
      # end
   end
   return sparse(I, J, Z, length(at), length(at))
end

# ====== TODO: revisit the FD hessians =========

# """
# `fd_hessian(V, R, h) -> H`
#
# If `length(R) = N` and `length(R[i]) = d` then `H` is an N × N `Matrix{SMatrix}` with
# each element a d × d static array.
# """
# function fd_hessian(V::SitePotential, R::Vector{SVec{D,T}}, h) where {D,T}
#    N = length(R)
#    H = zeros(SMat{D, D, T}, N, N)
#    return fd_hessian!(H, V, R, h)
# end
#
# """
# `fd_hessian!(H, V, R, h) -> H`
#
# Fill `H` with the hessian entries; cf `fd_hessian`.
# """
# function fd_hessian!(H, V::SitePotential, R::Vector{SVec{D,T}}, h) where {D,T}
#    N = length(R)
#    # convert R into a long vector and H into a big matrix (same part of memory!)
#    Rvec = mat(R)[:]
#    Hmat = zeros(N*D, N*D)   # reinterpret(T, H, (N*D, N*D))
#    # now re-define ∇V as a function of a long vector (rather than a vector of SVecs)
#    dV(r) = (evaluate_d(V, r |> vecs) |> mat)[:]
#    # compute the hessian as a big matrix
#    for i = 1:N*D
#       Rvec[i] +=h
#       dVp = dV(Rvec)
#       Rvec[i] -= 2*h
#       dVm = dV(Rvec)
#       Hmat[:, i] = (dVp - dVm) / (2 * h)
#       Rvec[i] += h
#    end
#    Hmat = 0.5 * (Hmat + Hmat')
#    # convert to a block-matrix
#    for i = 1:N, j = 1:N
#       Ii = (i-1) * D + (1:D)
#       Ij = (j-1) * D + (1:D)
#       H[i, j] = SMat{D,D}(Hmat[Ii, Ij])
#    end
#    return H
# end
#
# function fd_hessian(calc::AbstractCalculator, at::AbstractAtoms{T}, h) where {T}
#    N = length(at)
#    H = zeros(JMat{T}, N, N)
#    return fd_hessian!(H, calc, at, h)
# end
#
#
# """
# `fd_hessian!{D,T}(H, calc, at, h) -> H`
#
# Fill `H` with the hessian entries; cf `fd_hessian`.
# """
# function fd_hessian!(H, calc::AbstractCalculator, at::AbstractAtoms, h)
#    D = 3
#    N = length(at)
#    X = positions(at) |> mat
#    x = X[:]
#    # convert R into a long vector and H into a big matrix (same part of memory!)
#    Hmat = zeros(N*D, N*D)
#    # now re-define ∇V as a function of a long vector (rather than a vector of SVecs)
#    dE(x_) = (site_energy_d(calc, set_positions!(at, reshape(x_, D, N)), 1) |> mat)[:]
#    # compute the hessian as a big matrix
#    for i = 1:N*D
#       x[i] += h
#       dEp = dE(x)
#       x[i] -= 2*h
#       dEm = dE(x)
#       Hmat[:, i] = (dEp - dEm) / (2 * h)
#       x[i] += h
#    end
#    Hmat = 0.5 * (Hmat + Hmat')
#    # convert to a block-matrix
#    for i = 1:N, j = 1:N
#       Ii = (i-1) * D + (1:D)
#       Ij = (j-1) * D + (1:D)
#       H[i, j] = SMat{D,D}(Hmat[Ii, Ij])
#    end
#    return H
# end


text:
module TravisTest

export MrFunc

#External func that divides by 2
include("MrsFunc.jl")

#Func defined in moudle script that multiplys by 2
function MrFunc(X)
Y=X*2
return(Y)
end

end #end of module.

max_stars_repo_name: Timmmdavis/JuliaTravisTest

max_issues_repo_name: Timmmdavis/JuliaTravisTest

max_forks_repo_name: Timmmdavis/JuliaTravisTest