text:
function RegimeParameterDecomposition(parameter_values::Array{Float64,1},ergodic::Array{Float64,1})
    # outer constructor
    bar = dot(parameter_values,ergodic)
    hatarray = parameter_values - bar
    return RegimeParameterDecomposition(bar,hatarray)
end

function formatparametervalues(numregimes,numparameters,parameter_values,parameter_names)
    formatted = Array{Float64}(numregimes,numparameters)
    switching = BitArray(length(parameter_names))
    fill!(switching,false)
    for (p,value,name) in zip(1:numparameters,parameter_values,parameter_names)
        if !isa(value,Number)
            switching[p] = true
            if length(value)!=numregimes
                error("Regime-switching parameter $(name) has the wrong number of declared parameters values. It must have one for each regime.")
            else
                formatted[:,p]=collect(value)
            end
        else
            formatted[:,p]=value
        end
    end
    return formatted,switching
end

function createsymversions(vars::Array{String,1},suffix=""::String)
    out = Array{SymPy.Sym}(length(vars))
    for (i,var) in enumerate(vars)
        out[i] = Sym(string(var,suffix))
    end
    return out
end

function createsymparameters(names::Array{String,1},switches::BitArray{1})
    out = Array{SymPy.Sym}(2,length(names))
    for (i,name,switch) in zip(1:length(names),names,switches)
        if switch
            out[1,i] = Sym("$(name)_R")
            out[2,i] = Sym("$(name)_RP")
        else
            out[1,i] = Sym("$(name)")
            out[2,i] = out[1,i]
        end
    end
    return out
end

function createSSparameters(numparameters::Int,parameter_affectsSS::BitArray{1},barversions::Array{SymPy.Sym},
			     parameter_names::Array{String},parameters::Array{SymPy.Sym},switching::BitArray{1})
    parametersSS = Array{SymPy.Sym}(numparameters)
    for p = 1:numparameters
        if parameter_affectsSS[p]
            parametersSS[p] = barversions[p]
        elseif switching[p]
            # Define a new, non-state dependent symbol version
            parametersSS[p] = Sym(parameter_names[p])
        else
            # Just use the existing symbol
            parametersSS[p] = parameters[1,p]
        end
    end
    return parametersSS
end

function decomposeparameters(parameters,switching,affectsSS,parameter_values,ergodic,numparameters,numregimes)

    # Preallocate output arrays
    parameter_decomposition = Array{RegimeParameterDecomposition}(numparameters)
    parameter_bar_sym = Array{SymPy.Sym}(numparameters)
    parameter_regime_sym = Array{SymPy.Sym}(numregimes,numparameters)

    sympyzero = Sym(0.0)

    for (p,name) in enumerate(parameters)
        if affectsSS[p]
            parameter_decomposition[p] = RegimeParameterDecomposition(parameter_values[:,p],ergodic)
            parameter_bar_sym[p] = Sym("$(name)_BAR")
            for r = 1:numregimes
                parameter_regime_sym[r,p] = Sym("$(name)_HAT_R$(r)")
            end
        else
            parameter_decomposition[p] = RegimeParameterDecomposition(0.0,zeros(Float64,1))
            parameter_bar_sym[p] = sympyzero
            if switching[p]
                for r = 1:numregimes
                    parameter_regime_sym[r,p] = Sym("$(name)_R$(r)")
                end
            else
                parameter_regime_sym[:,p] = Sym("$(name)")
            end
        end
    end

    return parameter_decomposition, parameter_bar_sym, parameter_regime_sym
end

function creategenerichats(names,affectsSS)
    generics = Array{SymPy.Sym}(2,length(affectsSS))
    sympyzero = Sym(0.0)
    for (i,name) in enumerate(names)
        if affectsSS[i]
            generics[1,i] = Sym("$(name)_HAT_R")
            generics[2,i] = Sym("$(name)_HAT_RP")
        else
            generics[:,i] = sympyzero
        end
    end
    return generics
end

function Parameters!(meta,equations,vars,shocks,transmatrix,names,values,ss_guess)
    dictionary = Dict(zip(names,1:meta.numparameters))
    values,isswitching = formatparametervalues(meta.numregimes,meta.numparameters,values,names)
    generic_sym = createsymparameters(names,isswitching)

    affectsSS,ss_values = findparametersaffectSS(meta.numparameters,names,generic_sym,isswitching,
    values,equations.syms,vars.leads_sym,vars.states_sym,vars.contemps_sym,vars.states_indices,shocks.syms,ss_guess,transmatrix.ergodic)

    generic_hat = creategenerichats(names,affectsSS)

    ss_guess[:] = ss_values # Overwrite the ss_guess with the actual solved for steady state
    # Parameter decomposition and create the sym versions
    decomposition,bar_sym,specificregime_sym = decomposeparameters(names,isswitching,affectsSS,values,transmatrix.ergodic,meta.numparameters,meta.numregimes)
    # Create a values matrix that subs in bars
    valueswithbar = copy(values)
    for p in 1:meta.numparameters
        if affectsSS[p]
            valueswithbar[:,p] = decomposition[p].bar
        end
    end

    return Parameters(names,values,dictionary,isswitching,affectsSS,decomposition,valueswithbar,generic_sym,generic_hat,bar_sym,specificregime_sym)
end


text:
if VERSION >= v"0.7.0"
	import Random: randperm, shuffle!, shuffle, seed!
	import Printf.@printf
	import LinearAlgebra
	import LinearAlgebra: norm, Diagonal, dot
	import Statistics: var, mean, std
	using UnicodePlots
	import Base.minimum
end

using Requires

function __init__()
	@require Optim ="429524aa-4258-5aef-a3af-852621145aeb" @eval import Optim
end

ext: jl


text:
# Continuous quadratic problems

function qp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "QP1 - Quadratic objective" begin
        # simple quadratic objective
        # Min x^2 + xy + y^2 + yz + z^2
        # st  x + 2y + 3z >= 4 (c1)
        #     x +  y      >= 1 (c2)
        #     x,y \in R

        #@test MOI.supportsproblem(model, MOI.ScalarQuadraticFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        v = MOI.addvariables!(model, 3)
        @test MOI.get(model, MOI.NumberOfVariables()) == 3

        cf1 = MOI.ScalarAffineFunction(v, [1.0,2.0,3.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(cf1), MOI.GreaterThan{Float64})
        c1 = MOI.addconstraint!(model, cf1, MOI.GreaterThan(4.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

        @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
        c2 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([v[1],v[2]], [1.0,1.0], 0.0), MOI.GreaterThan(1.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 2

        obj = MOI.ScalarQuadraticFunction(MOI.VariableIndex[], Float64[], v[[1,1,2,2,3]], v[[1,2,2,3,3]], [2.0, 1.0, 2.0, 1.0, 2.0], 0.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

        if config.query
            @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
            @test obj ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())

            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c1))
            @test cf1 ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

            @test MOI.canget(model, MOI.ConstraintSet(), typeof(c1))
            @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 130/70 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.5714285714285715,0.4285714285714285,0.8571428571428572] atol=atol rtol=rtol
        end
    end
end

function qp2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "QP2" begin
        # same as QP0 but with duplicate terms
        # then change the objective and sense
        # simple quadratic objective
        # Min x^2 + xy + y^2 + yz + z^2
        # st  x + 2y + 3z >= 4 (c1)
        #     x +  y      >= 1 (c2)
        #     x,y \in R

        #@test MOI.supportsproblem(model, MOI.ScalarQuadraticFunction{Float64}, [(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        v = MOI.addvariables!(model, 3)
        @test MOI.get(model, MOI.NumberOfVariables()) == 3

        c1f = MOI.ScalarAffineFunction(v, [1.0,2.0,3.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(c1f), MOI.GreaterThan{Float64})
        c1 = MOI.addconstraint!(model, c1f, MOI.GreaterThan(4.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

        @test MOI.canaddconstraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
        c2 = MOI.addconstraint!(model, MOI.ScalarAffineFunction([v[1],v[2]], [1.0,1.0], 0.0), MOI.GreaterThan(1.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 2

        obj = MOI.ScalarQuadraticFunction(v, [0.0,0.0,0.0],[v[1], v[1], v[1], v[2], v[2], v[3], v[3]], [v[1], v[2], v[2], v[2], v[3], v[3], v[3]], [2.0, 0.5, 0.5, 2.0, 1.0, 1.0, 1.0], 0.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

        if config.query
            @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
            @test obj ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())

            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c1))
            @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

            @test MOI.canget(model, MOI.ConstraintSet(), typeof(c1))
            @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 130/70 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.5714285714285715,0.4285714285714285,0.8571428571428572] atol=atol rtol=rtol
        end

        # change objective to Max -2(x^2 + xy + y^2 + yz + z^2)
        obj2 = MOI.ScalarQuadraticFunction(v, [0.0,0.0,0.0],[v[1], v[1], v[1], v[2], v[2], v[3], v[3]], [v[1], v[2], v[2], v[2], v[3], v[3], v[3]], [-4.0, -1.0, -1.0, -4.0, -2.0, -2.0, -2.0], 0.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj2)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

        if config.query
            @test MOI.canget(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
            @test obj2 ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ -2*130/70 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [0.5714285714285715,0.4285714285714285,0.8571428571428572] atol=atol rtol=rtol
        end
    end
end

function qp3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "qp3test - Linear Quadratic objective" begin
        # simple quadratic objective
        #    minimize 2 x^2 + y^2 + xy + x + y + 1
        #       s.t.  x, y >= 0
        #             x + y = 1

        #@test MOI.supportsproblem(model, MOI.ScalarQuadraticFunction{Float64},
        #    [
        #        (MOI.SingleVariable,MOI.GreaterThan{Float64}),
        #        (MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64})
        #    ]
        #)

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x = MOI.addvariable!(model)
        MOI.canaddvariable(model)
        y = MOI.addvariable!(model)

        MOI.addconstraint!(model,
            MOI.ScalarAffineFunction([x,y], [1.0,1.0], 0.0),
            MOI.EqualTo(1.0)
        )

        @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
        MOI.addconstraint!(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
        MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

        obj = MOI.ScalarQuadraticFunction(
                [x,y], [1.0,1.0],
                [x,y,x], [x,y,y], [4.0, 2.0, 1.0],
                1.0
              )
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.875 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), [x,y]) ≈ [0.25, 0.75] atol=atol rtol=rtol
        end

        # change back to linear
        #        max 2x + y + 1
        #       s.t.  x, y >= 0
        #             x + y = 1
        # (x,y) = (1,0), obj = 3
        objf = MOI.ScalarAffineFunction([x,y], [2.0,1.0], 1.0)
        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol
            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), [x,y]) ≈ [1.0, 0.0] atol=atol rtol=rtol
        end

    end
end

function qptests(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    @testset "Quadratic Programs (quad. objective)" begin
        qp1test(model, atol=atol, rtol=rtol)
        qp2test(model, atol=atol, rtol=rtol)
        qp3test(model, atol=atol, rtol=rtol)
    end
end

#=
    Quadratically constrained programs
=#

function qcp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "qcp1" begin
        # quadratic constraint
        # Max x + y
        # st  - x + y >= 0 (c1[1])
        #       x + y >= 0 (c1[2])
        #     0.5x^2 + y <= 2 (c2)

        #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.VectorAffineFunction{Float64},MOI.Nonnegatives),(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x = MOI.addvariable!(model)
        MOI.canaddvariable(model)
        y = MOI.addvariable!(model)
        @test MOI.get(model, MOI.NumberOfVariables()) == 2

        @test MOI.canaddconstraint(model, MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)
        c1 = MOI.addconstraint!(model, MOI.VectorAffineFunction([1,1,2,2], [x,y,x,y], [-1.0,1.0,1.0,1.0], [0.0,0.0]), MOI.Nonnegatives(2))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}()) == 1

        c2f = MOI.ScalarQuadraticFunction([y],[1.0],[x],[x],[1.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(c2f), MOI.LessThan{Float64})
        c2 = MOI.addconstraint!(model, c2f, MOI.LessThan(2.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x,y], [1.0,1.0], 0.0))
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

        if config.query
            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c2))
            @test c2f ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.25 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), [x,y]) ≈ [0.5,1.75] atol=atol rtol=rtol
        end

        # try delete quadratic constraint and go back to linear

        # MOI.delete!(model, c2)
        #
        # MOI.optimize!(model)
        #
        # @test MOI.canget(model, MOI.TerminationStatus())
        # @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success
        #
        # @test MOI.canget(model, MOI.PrimalStatus())
        # @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint
        #
        # @test MOI.canget(model, MOI.ObjectiveValue())
        # @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
    end
end

function qcp2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "qcp2" begin
        # Max x
        # s.t. x^2 <= 2 (c)

        #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x = MOI.addvariable!(model)
        @test MOI.get(model, MOI.NumberOfVariables()) == 1

        cf = MOI.ScalarQuadraticFunction(MOI.VariableIndex[x],Float64[0.0],[x],[x],[1.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
        c = MOI.addconstraint!(model, cf, MOI.LessThan(2.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x], [1.0], 0.0))
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MaxSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MaxSense

        if config.query
            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c))
            @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            if config.duals
                @test MOI.canget(model, MOI.DualStatus())
                @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            end

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ sqrt(2) atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), x) ≈ sqrt(2) atol=atol rtol=rtol

            # TODO - duals
            # @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            # @test MOI.get(model, MOI.ConstraintDual(), c) ≈ 0.5/sqrt(2) atol=atol rtol=rtol
        end
    end
end

function qcp3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "qcp3" begin
        # Min -x
        # s.t. x^2 <= 2

        #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x = MOI.addvariable!(model)
        @test MOI.get(model, MOI.NumberOfVariables()) == 1

        cf = MOI.ScalarQuadraticFunction(MOI.VariableIndex[],Float64[],[x],[x],[1.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(cf), MOI.LessThan{Float64})
        c = MOI.addconstraint!(model, cf, MOI.LessThan(2.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([x], [-1.0], 0.0))
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

        if config.query
            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c))
            @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            if config.duals
                @test MOI.canget(model, MOI.DualStatus())
                @test MOI.get(model, MOI.DualStatus()) == MOI.FeasiblePoint
            end

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ -sqrt(2) atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), x) ≈ sqrt(2) atol=atol rtol=rtol

            # TODO - duals
            # @test MOI.canget(model, MOI.ConstraintDual(), typeof(c))
            # @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -0.5/sqrt(2) atol=atol rtol=rtol
        end
    end
end

function qcptests(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    @testset "Quadratic Constrainted Programs (quad. constraints only)" begin
        qcp1test(model, atol=atol, rtol=rtol)
        qcp2test(model, atol=atol, rtol=rtol)
        qcp3test(model, atol=atol, rtol=rtol)
    end
end

#=
    SOCP
=#

function socp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    @testset "socp1" begin
        # min t
        # s.t. x + y >= 1 (c1)
        #      x^2 + y^2 <= t^2 (c2)
        #      t >= 0 (bound)

        #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64}, [(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64}), (MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}), (MOI.SingleVariable,MOI.GreaterThan{Float64})])

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x = MOI.addvariable!(model)
        MOI.canaddvariable(model)
        y = MOI.addvariable!(model)
        MOI.canaddvariable(model)
        t = MOI.addvariable!(model)
        @test MOI.get(model, MOI.NumberOfVariables()) == 3

        c1f = MOI.ScalarAffineFunction([x,y],[1.0, 1.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(c1f), MOI.GreaterThan{Float64})
        c1 = MOI.addconstraint!(model, c1f, MOI.GreaterThan(1.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

        c2f = MOI.ScalarQuadraticFunction(MOI.VariableIndex[],Float64[],[x,y,t],[x,y,t],[1.0,1.0,-1.0], 0.0)
        @test MOI.canaddconstraint(model, typeof(c2f), MOI.LessThan{Float64})
        c2 = MOI.addconstraint!(model, c2f, MOI.LessThan(0.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

        @test MOI.canaddconstraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
        bound = MOI.addconstraint!(model, MOI.SingleVariable(t), MOI.GreaterThan(0.0))
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.GreaterThan{Float64}}()) == 1

        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([t], [1.0], 0.0))
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)
        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MinSense

        if config.query
            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c1))
            @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

            @test MOI.canget(model, MOI.ConstraintFunction(), typeof(c2))
            @test c2f ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
        end

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ sqrt(1/2) atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), [x,y,t]) ≈ [0.5,0.5,sqrt(1/2)] atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), [t,x,y,t]) ≈ [sqrt(1/2),0.5,0.5,sqrt(1/2)] atol=atol rtol=rtol
        end
    end
end

function socptests(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    @testset "Second Order Cone Programs" begin
        socp1test(model, atol=atol, rtol=rtol)
    end
end

const contquadratictests = Dict("quadratic1" => qp1test,
                                "quadratic2" => qp2test,
                                "quadratic3" => qp3test,
                                "quadratic4" => qcp1test,
                                "quadratic5" => qcp2test,
                                "quadratic6" => qcp3test,
                                "quadratic7" => socp1test)

@moitestset contquadratic


text:
function FcnNHCurlVec!(F,U,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
RhoPos=Param.RhoPos;
uPos=Param.uPos;
vPos=Param.vPos;
wPos=Param.wPos;
ThPos=Param.ThPos;
FCG = Param.FCG
FCG .= 0
KE = Param.KE
Pres = Param.Pres
RhoCG = Param.RhoCG
v1CG = Param.v1CG
v2CG = Param.v2CG
wCG = Param.wCG
wCCG = Param.wCCG
ThCG = Param.ThCG

RhoCG.= reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,RhoPos]
  ,OP,OP,NF,nz);
v1CG .= reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,uPos]
  ,OP,OP,NF,nz);
v2CG .= reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,vPos]
  ,OP,OP,NF,nz);
@views wCG[:,:,:,2:nz+1] .= reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,wPos]
  ,OP,OP,NF,nz);
@views wCCG = 0.5 .* (wCG[:,:,:,1:nz] .+ wCG[:,:,:,2:nz+1]);
@views wCCG[:,:,:,1] .= BoundaryW(v1CG,v2CG,CG,Param);
@views wCG[:,:,:,1] .= 2 .* wCCG[:,:,:,1] .- wCG[:,:,:,2];
ThCG .= reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,ThPos]
  ,OP,OP,NF,nz);

KE .= 0.5 .* (v1CG.*v1CG+v2CG.*v2CG+wCCG.*wCCG);

FDiv3Vec!(view(FCG,:,:,:,:,RhoPos),RhoCG,v1CG,v2CG,wCG,CG,Param);
if Param.RefProfile
# Param.Pres=Pressure(ThCG,RhoCG,KE,Param)-Param.pBGrd;
# FCG[:,:,:,:,uPos:wPos]=-FGrad3Vec(Param.Pres,CG,Param);
# FCG[:,:,:,:,uPos]=FCG[:,:,:,:,uPos]./RhoCG;
# FCG[:,:,:,:,vPos]=FCG[:,:,:,:,vPos]./RhoCG;
# if Param.Buoyancy
#   RhoF=0.5*(RhoCG[:,:,:,1:nz-1]+RhoCG[:,:,:,2:nz]);
#   FCG[:,:,:,1:nz-1,wPos]=(FCG[:,:,:,1:nz-1,wPos]-
#     Param.Grav*Param.JF[:,:,:,2:nz].*(RhoF-Param.RhoBGrdF)) ./ RhoF;
# end
else
  Pres .= Pressure(ThCG,RhoCG,KE,Param);
  FGrad3RhoVec!(FCG,Pres,RhoCG,CG,Param)
  if Param.Buoyancy
    @views FCG[:,:,:,1:nz-1,wPos] .-= Param.Grav*Param.JF[:,:,:,2:nz];
  end
end
FGrad3Vec!(FCG,KE,CG,Param)
FCurlNon3Vec!(FCG,v1CG,v2CG,wCG,wCCG,CG,Param);
if strcmp(Param.Thermo,"Energy")
else
  if Param.Upwind
    FCG[:,:,:,:,ThPos]=-FDiv3UpwindVec(ThCG,v1CG,v2CG,wCG,RhoCG,CG,Param);
  else
    FDiv3Vec!(view(FCG,:,:,:,:,ThPos),ThCG,v1CG,v2CG,wCG,CG,Param);
  end
end
if Param.HyperVisc
  if strcmp(Param.Thermo,"Energy")
  else
    HyperDiffusionVec!(FCG,v1CG,v2CG,wCG,ThCG,RhoCG,CG,Param)
  end
end

F .= 0
for iM=1:size(CG.FaceGlob,1)
  i = CG.FaceGlob[iM].Ind
  arr = reshape(CG.Glob[:,i,:],OP*OP*size(i,1))
  mat = reshape(FCG[:,:,i,:,:] ,OP*OP*size(i,1),nz,Param.NumV)
  @views F[arr,:,:] = F[arr,:,:] .+ mat;
end
@views F[:,:,RhoPos]=F[:,:,RhoPos]./CG.M;
@views F[:,:,uPos]=F[:,:,uPos]./CG.M;
@views F[:,:,vPos]=F[:,:,vPos]./CG.M;
@views F[:,1:nz-1,wPos]=F[:,1:nz-1,wPos]./CG.MW;
@views F[:,:,ThPos]=F[:,:,ThPos]./CG.M;

if Param.Damping
  @views F[:,1:nz-1,wPos]=F[:,1:nz-1,wPos]+Damping(U[:,:,wPos],Param);
end
if Param.Source
  Source!(F,U,CG,Param);
end
end

function FcnNHCurlVec(U,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
RhoPos=Param.RhoPos;
uPos=Param.uPos;
vPos=Param.vPos;
wPos=Param.wPos;
ThPos=Param.ThPos;
RhoCG=reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,RhoPos] ,OP,OP,NF,nz);
v1CG=reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,uPos] ,OP,OP,NF,nz);
v2CG=reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,vPos] ,OP,OP,NF,nz);
wCG=zeros(OP,OP,NF,nz+1);
wCG[:,:,:,2:nz+1]=reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,wPos] ,OP,OP,NF,nz);
wCCG=0.5*(wCG[:,:,:,1:nz]+wCG[:,:,:,2:nz+1]);
wCCG[:,:,:,1]=BoundaryW(v1CG,v2CG,CG,Param);
wCG[:,:,:,1]=2*wCCG[:,:,:,1]-wCG[:,:,:,2];
ThCG=reshape(U[reshape(CG.Glob,OP*OP*NF,1),:,ThPos]
  ,OP,OP,NF,nz);

FCG=zeros(OP,OP,NF,nz,Param.NumV);
KE=0.5*(v1CG.*v1CG+v2CG.*v2CG+wCCG.*wCCG);

FDiv3Vec!(view(FCG,:,:,:,:,RhoPos),RhoCG,v1CG,v2CG,wCG,CG,Param);
if Param.RefProfile
# Param.Pres=Pressure(ThCG,RhoCG,KE,Param)-Param.pBGrd;
# FCG[:,:,:,:,uPos:wPos]=-FGrad3Vec(Param.Pres,CG,Param);
# FCG[:,:,:,:,uPos]=FCG[:,:,:,:,uPos]./RhoCG;
# FCG[:,:,:,:,vPos]=FCG[:,:,:,:,vPos]./RhoCG;
# if Param.Buoyancy
#   RhoF=0.5*(RhoCG[:,:,:,1:nz-1]+RhoCG[:,:,:,2:nz]);
#   FCG[:,:,:,1:nz-1,wPos]=(FCG[:,:,:,1:nz-1,wPos]-
#     Param.Grav*Param.JF[:,:,:,2:nz].*(RhoF-Param.RhoBGrdF)) ./ RhoF;
# end
else
  Param.Pres=Pressure(ThCG,RhoCG,KE,Param);
  FGrad3RhoVec!(FCG,Param.Pres,RhoCG,CG,Param)
  if Param.Buoyancy
    @views FCG[:,:,:,1:nz-1,wPos] .-= Param.Grav*Param.JF[:,:,:,2:nz];
  end
end
FGrad3Vec!(FCG,KE,CG,Param)
FCG[:,:,:,:,uPos:wPos]=FCG[:,:,:,:,uPos:wPos]+
  FCurlNon3Vec(v1CG,v2CG,wCG,wCCG,CG,Param);
if strcmp(Param.Thermo,"Energy")
else
  if Param.Upwind
    FCG[:,:,:,:,ThPos]=-FDiv3UpwindVec(ThCG,v1CG,v2CG,wCG,RhoCG,CG,Param);
  else
    FDiv3Vec!(view(FCG,:,:,:,:,ThPos),ThCG,v1CG,v2CG,wCG,CG,Param);
  end
end
if Param.HyperVisc
  if strcmp(Param.Thermo,"Energy")
  else
#   FCG[:,:,:,:,uPos:ThPos]=FCG[:,:,:,:,uPos:ThPos]+
#   HyperDiffusionVec(v1CG,v2CG,wCG,ThCG,RhoCG,CG,Param);
    HyperDiffusionVec!(FCG,v1CG,v2CG,wCG,ThCG,RhoCG,CG,Param)
  end
end

F=zeros(size(U));
for iM=1:size(CG.FaceGlob,1)
  i = CG.FaceGlob[iM].Ind
  arr = reshape(CG.Glob[:,i,:],OP*OP*size(i,1))
  mat = reshape(FCG[:,:,i,:,:] ,OP*OP*size(i,1),nz,Param.NumV)
  F[arr,:,:] = F[arr,:,:] .+ mat;
end
F[:,:,RhoPos]=F[:,:,RhoPos]./CG.M;
F[:,:,uPos]=F[:,:,uPos]./CG.M;
F[:,:,vPos]=F[:,:,vPos]./CG.M;
F[:,1:nz-1,wPos]=F[:,1:nz-1,wPos]./CG.MW;
F[:,:,ThPos]=F[:,:,ThPos]./CG.M;

if Param.Damping
  F[:,1:nz-1,wPos]=F[:,1:nz-1,wPos]+Damping(U[:,:,wPos],Param);
end
if Param.Source
  Source!(F,U,CG,Param);
end
return F
end

function Damping(W,Param)
F=zeros(size(W,1),size(W,2)-1);
for iz=Param.Grid.nz-1:-1:1
  zLoc=Param.Grid.z[iz];
  if zLoc>=Param.H-Param.StrideDamp
    Damp = Param.Relax*
      sin(0.5*pi*(1.0 - (Param.H - zLoc)/Param.StrideDamp))^2;
    F[:,iz]=-Damp*W[:,iz];
  else
    break
  end
end
return F
end


text:
"""
  punch_holes_nexus(xpoints, ypoints, zpoints, radius)

...
# Arguments

  - `xpoints::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `ypoints::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `zpoints::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `radius::Union{Float64,Vector{Float64}}`: the radius, or radii of the punch, if vector.
...

...
# Outputs

  - `absolute_indices::Vector{Int64}`: vector containing the indices of coordinates
  inside the punch

...
"""
function punch_holes_nexus(xpoints, ypoints, zpoints, radius)
    rad = (typeof(radius) <: Tuple) ? radius : (radius, radius, radius)
    radius_x, radius_y, radius_z = rad
    absolute_indices = Int64[]
    count = 1
    for i = 1:length(zpoints)
        ir = round(zpoints[i])
        for j = 1:length(ypoints)
            jr = round(ypoints[j])
            for h = 1:length(xpoints)
                hr = round(xpoints[h])
                if (((hr - xpoints[h])/radius_x)^2 + ((jr - ypoints[j])/radius_y)^2 + ((ir - zpoints[i])/radius_z)^2 <= 1.0)
                    append!(absolute_indices, count)
                end
                count=count+1
            end
        end
    end
    return absolute_indices
end

"""
  punch_holes_nexus_Cartesian(x, y, z, radius)

...
# Arguments

  - `x::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `y::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `z::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `radius::Union{Float64,Tuple{Float64}}`: the radius, or radii of the punch, if vector.
...

...
# Outputs

  - `inds::Vector{Int64}`: vector containing the indices of coordinates
  inside the punch

...
"""
function punch_holes_nexus_Cartesian(x, y, z, radius)
    radius_x, radius_y, radius_z = (typeof(radius) <: Tuple) ? radius : (radius, radius, radius)
    inds = filter(i -> (((x[i[1]] - round(x[i[1]])) / radius_x) ^2
                        + ((y[i[2]] - round(y[i[2]])) / radius_y) ^2
                        + ((z[i[3]] - round(z[i[3]])) / radius_z) ^2 <= 1.0),
                  CartesianIndices((1:length(x), 1:length(y), 1:length(z))))
    return inds
end

"""

  interp_1lu(x, y, z, imgg, punch_template, eps, m)

Interpolate around the bragg peaks in the image by tiling the punch_template

...
# Arguments
  - `x::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `y::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `z::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg::Array{Float64,3]`: the matrix containing the image
  - `punch_template::Array{Float64,3}`: one lattice unit template of the values to be filled
  - `eps::Float64 = 0.0`: Matern parameter eps
  - `m::Int64 = 1` : Matern parameter

# Outputs
  - array containing the interpolated image
...
"""
function interp_1lu(x, y, z, imgg, punch_template, eps, m)
  discard = findall(punch_template .> 0)
  res = interp(x, y, z, imgg, discard, eps, m)
  return res
end

"""

  interp_nexus(x, y, z, imgg, punch_template, eps, m)

Interpolate around the bragg peaks in the image by tiling the punch_template

...
# Arguments
  - `x::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `y::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `z::Vector{T} where T<:Real`: the vector containing the z coordinate
  - `imgg::Array{Float64,3]`: the matrix containing the image
  - `punch_template::Array{Float64,3}`: one lattice unit template of the values to be filled
  - `eps::Float64 = 0.0`: Matern parameter eps
  - `m::Int64 = 1` : Matern parameter

# Outputs
  - array containing the interpolated image
...

# Example

```<julia-repl>
x = y = z = collect(-0.5:0.1:10)
imgg = randn(length(x), length(y), length(z))
punch_template = zeros(10, 10, 10); punch_template[5, 5, 5] = 1
interp_nexus(x, y, z, imgg, punch_template)
```
"""
function interp_nexus(x, y, z, imgg, punch_template, eps = 0.0, m = 1)
  discard = findall(punch_template .> 0)
  x_unit, y_unit, z_unit = size(punch_template)
  Nx, Ny, Nz = size(imgg)
  Qx = ceil(x[end]) + 1
  Qy = ceil(y[end]) + 1
  Qz = ceil(z[end]) + 1
  corners_x = findall((x .- 0.5) .% 1 .== 0.0)
  corners_y = findall((y .- 0.5) .% 1 .== 0.0)
  corners_z = findall((z .- 0.5) .% 1 .== 0.0)
  corners = [CartesianIndex(xind, yind, zind) for xind in corners_x[1:end-1]
                                              for yind in corners_y[1:end-1]
                                              for zind in corners_y[1:end-1]]
  # punched = copy(imgg);
  # Threads.@threads for (i,c) in enumerate(corners)
  for (i,c) in enumerate(corners)
    tcx = minimum([Nx, Tuple(c)[1] + x_unit])
    tcy = minimum([Ny, Tuple(c)[2] + y_unit])
    tcz = minimum([Nz, Tuple(c)[3] + z_unit])
    tc = CartesianIndex(tcx, tcy, tcz)
    imgg[c:tc] = interp(x[c[1]:tc[1]], y[c[2]:tc[2]], z[c[3]:tc[3]],
                       imgg[c:tc], discard, Float64(eps), Int64(m))
    # punched[c:tc] .= 0.0
  end
  return imgg
end

max_stars_repo_name: rayosborn/nxpeaks

max_issues_repo_name: rayosborn/nxpeaks

max_forks_repo_name: axmas-anl/nxrefine


text:
"""
	gmtread(fname::String, data; kwargs...)

Read GMT object from file. The object is one of "grid" or "grd", "image" or "img",
"data" or "table", "cmap" or "cpt" and "ps" (for postscript).

Parameters
----------

Specify data type.  Choose among:
- **grd** | **grid** :: [Type => Any]

    Load a grid.
- **img** | **image** :: [Type => Any]

    Load an image.
- **cpt** | **cmap** :: [Type => Any]

    Load a GMT color palette.
- **dataset** | **table** :: [Type => Any]

    Load a dataset (a table of numbers).
- **ps** :: [Type => Any]

    Load a PostScript file
- **gdal** :: [Type => Any]

    Force reading the file via GDAL. Should only be used to read grids.
- **varname** :: [Type => Str]

    When netCDF files have more than one 2D (or higher) variables use *varname* to pick the wished variable.
    e.g. varname=:slp to read the variable named 'slp'. This option defaults data type to 'grid'
- **layer** | **band** :: [Type => Str, Number, Array]

    When files are multiband or nc files with 3D or 4D arrays, we access them via these keywords.
    layer=4 reads the fourth layer (or band) of the file. But the file can be a grid or an image. If it is a
    grid layer can be a scalar (to read 3D arrays) or an array of two elements (to read a 4D array).
    If file is an image 'layer' can be a 1 or a 1x3 array (to read a RGB image). Not that in this later case
    bands do not need to be contiguous. A band=[0,5,2] composes an RGB out of those bands. See more at
    $(GMTdoc)/GMT_Docs.html#modifiers-for-coards-compliant-netcdf-files)

- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_f)

Example: to read a nc called 'lixo.grd'

    G = gmtread("lixo.grd", grd=true);

to read a jpg image with the bands reversed (this example is currently broken in GMT5. Needs GMT6dev)

    I = gmtread("image.jpg", band=[2,1,0]);
"""
function gmtread(fname::String; kwargs...)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:R :V_params :f :i :h])
	cmd, opt_bi = parse_bi(cmd, d)

	# Process these first so they may take precedence over defaults set below
	opt_T = add_opt("", "Tg", d, [:grd :grid])
	if (opt_T != "")
		(haskey(d, :gdal)) && (fname = fname * "=gd")     # Force read via GDAL
	else
		opt_T = add_opt("", "Ti", d, [:img :image])
	end
	if (opt_T == "")  opt_T = add_opt("", "Td", d, [:dataset :table])  end
	if (opt_T == "")  opt_T = add_opt("", "Tc", d, [:cpt :cmap])  end
	if (opt_T == "")  opt_T = add_opt("", "Tp", d, [:ps])  end

	if (haskey(d, :varname))				# See if we have a nc varname / layer request
		if (isempty(opt_T))
			(haskey(d, :gdal)) && (fname = fname * "=gd")     # Force read via GDAL
			opt_T = " -Tg"
		end
		fname = fname * "?" * arg2str(d[:varname])
		if ((val = find_in_dict(d, [:layer :band])[1]) !== nothing)
			if (isa(val, Number))     fname *= @sprintf("[%d]", val)
			elseif (isa(val, Array))  fname *= @sprintf("[%d,%d]", val[1], val[2])	# A 4D array
			end
		end
	else									# See if we have a bands request
		if ((val = find_in_dict(d, [:layer :band])[1]) !== nothing)
			fname = fname * "+b"
			if (isa(val, String) || isa(val, Symbol) || isa(val, Number))
				fname = string(fname, val)
			elseif (isa(val, Array) || isa(val, Tuple))
				if (length(val) == 3)
					fname = fname * @sprintf("%d,%d,%d", val[1], val[2], val[3])
				else
					error("Number of bands in the 'band' option can only be 1 or 3")
				end
			end
			if (isempty(opt_T))	opt_T = " -Ti"	end
		end
	end

	if (opt_T == "")
		error("Must select one input data type (grid, image, dataset, cmap or ps)")
	else
		opt_T = opt_T[1:4]      				# Remove whatever was given as argument to type kwarg
	end

	if (opt_T == " -Ti" || opt_T == " -Tg")		# See if we have a mem layout request
		if (haskey(d, :layout))
			t = arg2str(d[:layout])
			cmd = (opt_T == " -Ti") ? cmd * " -%" * t : cmd * " -&" * t
		end
	end

	if (opt_T == " -Td" && !isempty(opt_bi))  cmd *= opt_bi  end		# Read from binary file
	cmd *= opt_T
	(haskey(d, :Vd)) && println(@sprintf("\tread %s %s", fname, cmd))

	O = gmt("read " * fname * cmd)
end

"""
	gmtwrite(fname::String, data; kwargs...)

Write a GMT object to file. The object is one of "grid" or "grd", "image" or "img",
"dataset" or "table", "cmap" or "cpt" and "ps" (for postscript).

When saving grids we have a large panoply of formats at our disposal.

Parameters
----------

- **id** ::  [Type => Str]

    Use an ``id`` code when not not saving a grid into a standard COARDS-compliant netCDF grid. This ``id``
    is made up of two characters like ``ef`` to save in ESRI Arc/Info ASCII Grid Interchange format (ASCII float).
    See the full list of ids at $(GMTdoc)grdconvert.html#format-identifier.

    ($(GMTdoc)grdconvert.html#g)
- **scale** | **offset** :: [Type => Number]

    You may optionally ask to scale the data and then offset them with the specified amounts.
    These modifiers are particularly practical when storing the data as integers, by first removing an offset
    and then scaling down the values.
- **nan** | **novalue** | **invalid** | **missing** :: [Type => Number]

    Lets you supply a value that represents an invalid grid entry, i.e., ‘Not-a-Number’.
- **gdal** :: [Type => Bool]

    Force the use of the GDAL library to write the grid (to be used only with grids).
    ($(GMTdoc)GMT_Docs.html#grid-file-format-specifications)
- **driver** :: [Type => Str]

    When saving in other than the netCDF format we must tell the GDAL library what is wished format.
    That is done by specifying the driver name used by GDAL itself (e.g., netCDF, GTiFF, etc...).
- **datatype** :: [Type => Str] 		``Arg = u8|u16|i16|u32|i32|float32``

    When saving with GDAL we can specify the data type from u8|u16|i16|u32|i32|float32 where ‘i’ and ‘u’ denote
    signed and unsigned integers respectively.
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_bo)
- $(GMT.opt_f)

Example: write the GMTgrid 'G' object into a nc file called 'lixo.grd'

	gmtwrite("lixo.grd", G);
"""
function gmtwrite(fname::String, data; kwargs...)

	if (fname == "")
		error("First argument cannot be empty. It must contain the file name to write.")
	end

	d = KW(kwargs)
	cmd, opt_R = parse_R("", d)
	cmd = parse_V_params(cmd, d)

	if (isa(data, GMTgrid))
		opt_T = " -Tg"
		fname = fname * parse_grd_format(d)		# If we have format requests
		cmd, = parse_f(cmd, d)
	elseif (isa(data, GMTimage))
		opt_T = " -Ti"
		fname *= parse_grd_format(d)		# If we have format requests
	elseif (isa(data, GMTdataset))
		opt_T = " -Td"
		cmd, = parse_bo(cmd, d)					# Write to binary file
	elseif (isa(data, GMTcpt))
		opt_T = " -Tc"
	elseif (isa(data, GMTps))
		opt_T = " -Tp"
	elseif (isa(data, Array{UInt8}))
		fmt = parse_grd_format(d)				# See if we have format requests
		if (fmt == "")							# If no format, write a dataset
			opt_T = " -Td"
			cmd, = parse_bo(cmd, d)				# Write to binary file
		else
			data = mat2img(data)
			fname *= fmt
			opt_T = " -Ti"
		end
	elseif (isa(data, AbstractArray))
		fmt = parse_grd_format(d)				# See if we have format requests
		if (fmt == "")							# If no format, write a dataset
			opt_T = " -Td"
			cmd, = parse_bo(cmd, d)				# Write to binary file
		else
			data = mat2grid(data)
			fname *= fmt
			opt_T = " -Tg"
		end
	end
	cmd = cmd * opt_T

	if (opt_T == " -Ti" || opt_T == " -Tg")		# See if we have a mem layout request
		if (haskey(d, :layout))
			t = arg2str(d[:layout])
			cmd = (opt_T == " -Ti") ? cmd * " -%" * t : cmd * " -&" * t
		end
	end

	(haskey(d, :Vd)) && println(@sprintf("\twrite %s %s", fname, cmd))

	gmt("write " * fname * cmd, data)
	return nothing
end

# -----------------------------------------------------------------------------------------------
function parse_grd_format(d::Dict)
# Scan options to fill any of  [=<id>][+s<scale>][+o<offset>][+n<nan>][:<driver>[/<dataType>]
# that control the grid/image output format
	out = ""
	for sym in [:id :gdal :driver]
		if (haskey(d, sym))
			if (sym == :gdal || sym == :driver)
				out = "=gd"
			else
				t = arg2str(d[sym])
				if (length(t) != 2)
					error(@sprintf("Format code MUST have 2 characters and not %s", t))
				end
				out = "=" * t
			end
			break
		end
	end
	if (haskey(d, :scale))  out *= "+s" * arg2str(d[:scale])  end
	if (haskey(d, :offset)) out *= "+o" * arg2str(d[:offset]) end
	if ((val = find_in_dict(d, [:nan :novalue :invalid :missing])[1]) !== nothing)
		out *= "+n" * arg2str(val)
	end
	if (haskey(d, :driver))
		out *= ":" * arg2str(d[:driver])
		if (haskey(d, :datatype))  out *= "/" * arg2str(d[:datatype])  end
	end
	return out
end


text:
abstract FullStrategy

type SeparateTurnover <: FullStrategy
    OptStrat::InitialStrategy
    TurnOvStrat::TOFilter

    function SeparateTurnover(strat::InitialStrategy,
                              toStrat::TOFilter)
        return new(strat, toStrat)
    end
end

##################
## optimizeWgts ##
##################

## without history
##----------------

## without investment history (initial investment)
function optimizeWgts(strat::AssetMgmt.SeparateTurnover,
                      univ::AssetMgmt.MuSigmaUniverse)

    ## first get ideal weights without transaction costs
    idealWgts = optimizeWgts(strat.OptStrat, univ)

    return idealWgts
end

## regular case
##-------------

## case: ThresHoldDeviance Filter

## Note: any given TOFilter should work for all initial optimization
## strategies -> parametric type
function optimizeWgts(strat::AssetMgmt.SeparateTurnover,
                      univ::AssetMgmt.MuSigmaUniverse,
                      invHistory::Investments)

    ## first get ideal weights without transaction costs
    idealWgts = optimizeWgts(strat, univ)

    ## get last morning weights
    ## TODO: invHistory[end, :] returns DataFrame!!
    lastDate = idx(invHistory)[end]
    invTm = convert(Timematr, invHistory)
    lastWeights = asArr(invTm[lastDate, :], Float64, NaN)

    ## get corresponding return data
    rets = asArr(univ.Data[lastDate, :], Float64, NaN)

    ## calculate current weights after daily evolution
    eveningWgts = evolWgtsCore(lastWeights, rets)

    ## check deviance between real weights and desired weights
    wgtsDeviance = sum(abs(idealWgts[:] .- eveningWgts[:]))
    display(wgtsDeviance)
    if wgtsDeviance > strat.TurnOvStrat.thres
        ## println("rebalancing wgts:")
        ## display(idealWgts)
        return idealWgts
    else
        ## println("evolving wgts:")
        ## display(eveningWgts)
        return eveningWgts
    end
end

###################
## applyStrategy ##
###################

function applyStrategy(strat::AssetMgmt.FullStrategy,
                       estimator::Type{AssetMgmt.ExpWeighted},
                       data::Timematr)
    ## indicator for initial investment
    initialInv = true

    ## iterate over dates
    allDats = idx(data)
    nDats = length(allDats)

    ## preallocation
    nObs, nAss = size(data)
    allWgts = Array(Float64, nObs, nAss)
    pfMoments = Array(Float64, nObs, 2)
    modelEstimated = falses(nObs)

    for ii=1:nDats
        thisDate = allDats[ii]
        univ = AssetMgmt.estimate(estimator, data, thisDate)
        if AssetMgmt.isDef(univ)
            if initialInv
                wgts = AssetMgmt.optimizeWgts(strat, univ)
                initialInv = false
            else
                ## get history
                wgtsHist = allWgts[modelEstimated, :]

                dateHist = allDats[modelEstimated]
                wgtsHistDf = composeDataFrame(wgtsHist, names(data))
                invHistory = AssetMgmt.Investments(wgtsHistDf,
                                         dateHist)

                ## optimize with given history
                wgts = AssetMgmt.optimizeWgts(strat, univ, invHistory)

            end

            ## set estimation indicator to true
            modelEstimated[ii] = true

            allWgts[ii, :] = wgts

            ## calculate portfolio moments
            pfMu = AssetMgmt.getPMean(wgts[:], univ.Universe.mu)
            pfSigma = sqrt(AssetMgmt.getPVar(wgts[:], univ.Universe.sigma))
            pfMoments[ii, :] = [pfMu pfSigma]
        end
    end

    ## remove non-estimatable dates
    allWgts = allWgts[modelEstimated, :]
    pfMoments = pfMoments[modelEstimated, :]
    allDats = allDats[modelEstimated]

    ## DataFrame
    wgtsDf = composeDataFrame(allWgts, names(data))
    pfMomentsDf = Timematr(pfMoments, [:mu, :sigma], allDats)

    invs = Investments(wgtsDf, allDats)

    return invs, pfMomentsDf
end

## example:
## investment strategy could be a certain MaxMu plain Markowitz
## strategy together with slow moving turnover constraint

## Turnover filter:
## - each second period only
## - convex combination between current and desired portfolio
## - threshold: only if too far off
## 	- in terms of weight differences
##    - in terms of mu-sigma differences

##############
## old code ##
##############

## function randHolding(discRet::Timematr)
##     ## create random weights and hold them without re-balancing

## end

## function randInvestments(discRet::Timematr)
##     ## create random investments for given returns
##     nams = names(discRet)

##     (nObs, nAss) = size(discRet)

##     wgts = randWgts(nObs, nAss)
##     df = convert(DataFrame, wgts)
##     rename!(df, names(df), nams)

##     return AssetMgmt.Investments(df, idx(discRet))
## end

## function equWgtInvestments(discRet::Timematr)
##     ## create equal weights investments for given returns
##     nams = names(discRet)

##     (nObs, nAss) = size(discRet)

##     wgts = 1/nAss * ones(nObs, nAss)
##     df = convert(DataFrame, wgts)
##     rename!(df, names(df), nams)

##     return AssetMgmt.Investments(df, idx(discRet))
## end

## function gmvSampleMoments(minDays, increasingSample)

##     ## estimation framework:
##     ## - minimum number of days?
##     ## - increasing sample size?
##     ##

## end


text:
"""
    snr = calculateSNR(hwfft, f, nsig=1)

Estimate the signal-to-noise ratio, given the in-band bins of a (Hann-windowed)
fft and the location of the input signal (`f` > 0). For `nsig` = 1, the input tone is
contained in `hwfft`[`f`:`f`+2]; this range is appropriate for a Hann-windowed fft.
Each increment in `nsig` adds a bin to either side. The `snr` is expressed in dB.
"""
function calculateSNR(hwfft, f, nsig=1)
    signalBins = (f-nsig+1):(f+nsig+1)
    signalBins = signalBins[signalBins .> 0]
    signalBins = signalBins[signalBins .<= length(hwfft)]
    s = norm(hwfft[signalBins])

    noiseBins = collect(1:length(hwfft))
    deleteat!(noiseBins, signalBins)
    n = norm(hwfft[noiseBins])

    if n == 0
        snr = Inf
    else
        snr = dbv(s/n)
    end

    return snr
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
using ChainRules, Test, FDM, LinearAlgebra, LinearAlgebra.BLAS, Random, Statistics
using Base.Broadcast: broadcastable
import LinearAlgebra: dot

# For testing purposes we use a lot of
using ChainRulesCore: add, cast, extern, accumulate, accumulate!, store!, @scalar_rule,
    Wirtinger, wirtinger_primal, wirtinger_conjugate, add_wirtinger, mul_wirtinger,
    Zero, add_zero, mul_zero, One, add_one, mul_one, Casted, cast, add_casted, mul_casted,
    DNE, Thunk, Casted, DNERule

include("test_util.jl")

@testset "ChainRules" begin
    include("helper_functions.jl")
    @testset "rulesets" begin
        @testset "Base" begin
            include(joinpath("rulesets", "Base", "base.jl"))
            include(joinpath("rulesets", "Base", "array.jl"))
            include(joinpath("rulesets", "Base", "mapreduce.jl"))
            include(joinpath("rulesets", "Base", "broadcast.jl"))
        end

        @testset "LinearAlgebra" begin
            include(joinpath("rulesets", "LinearAlgebra", "dense.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "structured.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "factorization.jl"))
            include(joinpath("rulesets", "LinearAlgebra", "blas.jl"))
        end

        @testset "packages" begin
            include(joinpath("rulesets", "packages", "NaNMath.jl"))
            include(joinpath("rulesets", "packages", "SpecialFunctions.jl"))
        end
    end
end


text:
using Espresso: @simple_rule

struct Rule
    left
    right
    condition
end

# These symbols must always correspond to numbers
const NUM_PHS = [:a, :n]
integralTable = []

addMatch!(x::Symbol) = Expr(:ref, :match, QuoteNode(x))
function addMatch!(ex::Expr)
    for (i, x) ∈ enumerate(ex.args)
        if i == 1 && ex.head == :call
            continue
        elseif x isa Symbolic
            ex.args[i] = addMatch!(x)
        end
    end
    ex
end

addMatch(ex::Expr) = Expr(:->, :match, addMatch!(ex))

macro integration_rule(pat, rpat, condition = nothing)
    conditions = Any[m -> m[p] isa Number for p ∈ NUM_PHS if contains(p, pat)]
    if condition ≠ nothing
        condition = condition |> addMatch |> eval
        push!(conditions, condition)
    end

    fullCondition = _ -> true
    if !isempty(conditions)
        fullCondition = match -> all(f -> f(match), conditions)
    end

    rule = Rule(pat, rpat, fullCondition)
    push!(integralTable, rule)
    nothing
end

@simple_rule    (+)(x)     x
@simple_rule    (*)(x)     x
@simple_rule    x - -y     x + y
@simple_rule    x + -y     x - y

@integration_rule   x       x^2 / 2
@integration_rule   1/x     log(x)
@integration_rule   x^n     x^(n+1) / (n+1)       n ≠ -1
@integration_rule   a^x     a^x / log(a)
@integration_rule   exp(x)  e^x
@integration_rule   e^x     e^x
@integration_rule   exp(-x) -e^(-x)
@integration_rule   e^(-x)  -e^(-x)

@integration_rule   sin(x)  -cos(x)
@integration_rule   cos(x)  sin(x)
@integration_rule   1 / sin(x)^2      -cot(x)
@integration_rule   1 / cos(x)^2      tan(x)

@integration_rule   log(x)            x*log(x) - x
@integration_rule   1/(x*log(x))      log(log(x))
@integration_rule   x * log(x)        x^2 * (log(x)/2 - 1/4)
@integration_rule   x^n * log(x)      x^(n+1) * (log(x)/(n+1) - 1/(n+1)^2)


text:
using LightGraphs
using DataFrames
using Random
using Statistics
using JLD
using JLD2
using StatsBase
using CSV
using GraphIO

for script in readdir("src")
    if script != "simulation.jl" && script != "03runeval.jl"
        include(script)
    end
end

mutable struct Simulation

    name::String
    runnr::Int64
    repnr::Int64
    rng::MersenneTwister
    config::Config
    init_state::Any
    final_state::Any
    agent_log::DataFrame
    post_log::Any
    graph_list::Array{AbstractGraph}

    function Simulation(config=Config())
        new(
            "",
            0,
            0,
            MersenneTwister(),
            config,
            (nothing, nothing),
            (nothing, nothing),
            DataFrame(),
            DataFrame(),
            Array{AbstractGraph, 1}(undef, 0)
        )
    end

    function Simulation(sim_transferred::Tuple{String,Int64,Int64,MersenneTwister,Config,Tuple{SimpleDiGraph{Int64},Array{Agent,1}},Tuple{SimpleDiGraph{Int64},Array{Agent,1}},DataFrame,DataFrame,Array{AbstractGraph,1}})
        new(
            sim_transferred[1],
            sim_transferred[2],
            sim_transferred[3],
            sim_transferred[4],
            sim_transferred[5],
            sim_transferred[6],
            sim_transferred[7],
            sim_transferred[8],
            sim_transferred[9],
            sim_transferred[10]
        )
    end

end

function tick!(
    state::Tuple{AbstractGraph, AbstractArray},
    tick_nr::Int64, config::Config, rng::MersenneTwister
)
    agent_list = state[2]
    if config.simulation.agent_logging
        prev_state = deepcopy(state)
    end

    for agent_idx in shuffle(rng, 1:length(agent_list))
        this_agent = agent_list[agent_idx]
        update_feed!(state, agent_idx, config)
        update_perceiv_publ_opinion!(state, agent_idx)
        update_opinion!(state, agent_idx, config)

        drop_friends!(state, agent_idx, config)
        if indegree(state[1], agent_idx) < config.network.m0
            if config.simulation.addfriends == "neighborsofneighbors"
                add_friends_neighbors_of_neighbors!(
                    state, agent_idx,
                    config, config.network.new_follows, rng
                )
            elseif config.simulation.addfriends == "random"
                add_friends_random!(
                    state, agent_idx,
                    config, config.network.new_follows
                )
            else
                add_friends_neighbors_of_neighbors!(
                    state, agent_idx,
                    config, floor(Int,config.network.new_follows/2), rng
                )
                add_friends_random!(
                    state, agent_idx,
                    config, floor(Int,config.network.new_follows/2)
                )
            end
        end

        inclin_interact = deepcopy(this_agent.inclin_interact)
        while inclin_interact > 0
            if rand() < inclin_interact
                publish_post!(state, agent_idx, rng, tick_nr)
            end
            inclin_interact -= 1.0
        end
    end

    if config.simulation.agent_logging
        return log_network(state, prev_state, tick_nr)
    else
        return state
    end
end

function run!(
    simulation::Simulation;
    name::String = "result"
)
    if name * ".jld2" in readdir("results")
        raw = load(joinpath("results", name * ".jld2"))
        rep = raw[first(keys(raw))]
    else
        rep = Simulation[]
    end

    config = simulation.config
    simulation.name = name

    for current_rep in 1:config.simulation.repcount

        simulation.repnr = length(rep) + 1

        simulation.rng = Random.seed!(sum(codeunits(name)) + simulation.repnr)

        graph = SimpleDiGraph(barabasi_albert(
                        config.network.agent_count,
                        config.network.m0))
        simulation.init_state = (graph, create_agents(graph, simulation.rng))
        state = deepcopy(simulation.init_state)

        if config.simulation.agent_logging
            simulation.agent_log = DataFrame(
            TickNr = 0,
            AgentID = [agent.id for agent in state[2]],
            Opinion = [agent.opinion for agent in state[2]]
            )
        end

        for i in 1:config.simulation.ticks

            if config.simulation.agent_logging
                append!(simulation.agent_log, tick!(state, i, config, simulation.rng))
            else
                tick!(state, i, config, simulation.rng)
            end

            if i % ceil(config.simulation.ticks / 10) == 0
                if name != "result"
                    print(".")
                end
            end
        end

        simulation.final_state = state

        push!(rep, deepcopy(simulation))
    end

    if !in("results", readdir())
        mkdir("results")
    end
    save(joinpath("results", name * ".jld2"), name, rep)
    # rm(joinpath("tmp", name * ".jld2"))

    print("\n---\nFinished simulation run with the following specifications:\n $config\n---\n")

    return rep
end

function run_batch(
    configlist::Array{Config, 1};
    resume_at::Int64=1,
    stop_at::Int64=length(configlist),
    batch_name::String = ""
)

    for i in resume_at:stop_at
        run_nr = lpad(string(i),length(string(length(configlist))),"0")
        current_sim = Simulation(configlist[i])
        current_sim.runnr = i
        run!(
            current_sim,
            name = (batch_name * "_run$run_nr")
        )
    end
end

function rerun_single(
    simulation_imported::Simulation;
    logging::Bool = true
    )
    simulation = deepcopy(simulation_imported)

    simulation.config = Config(
        network = simulation.config.network,
        simulation = cfg_sim(
            ticks = simulation.config.simulation.ticks,
            addfriends = simulation.config.simulation.addfriends,
            repcount = simulation.config.simulation.repcount,
            agent_logging = logging
        ),
        opinion_threshs = simulation.config.opinion_threshs,
        agent_props = simulation.config.agent_props
    )

    config = simulation.config

    simulation.rng = Random.seed!(simulation.rng.seed)

    graph = SimpleDiGraph(barabasi_albert(
                    config.network.agent_count,
                    config.network.m0))
    simulation.init_state = (graph, create_agents(graph, simulation.rng))
    state = deepcopy(simulation.init_state)

    if config.simulation.agent_logging
        simulation.agent_log = DataFrame(
        TickNr = 0,
        AgentID = [agent.id for agent in state[2]],
        Opinion = [agent.opinion for agent in state[2]],
        OpinionChange = 0.0,
        PerceivPublOpinion = [agent.opinion for agent in state[2]]
        )
    end

    for i in 1:config.simulation.ticks

        if config.simulation.agent_logging
            append!(simulation.agent_log, tick!(state, i, config, simulation.rng))
        else
            tick!(state, i, config, simulation.rng)
        end

        if i % ceil(config.simulation.ticks / 10) == 0
            print(".")
        end
    end

    communities = DataFrame(
        AgentID = 1:length(state[2]),
        Community = label_propagation(state[1])[1]
    )

    simulation.final_state = state[1], state[2], communities

    if !in("results", readdir())
        mkdir("results")
    end
    save(joinpath("results", simulation.name * "_singlerun.jld2"), simulation.name, simulation)
    # rm(joinpath("tmp", name * ".jld2"))

    print("\n---\nFinished simulation run with the following specifications:\n $config\n---\n")

    return simulation
end

# suppress output of include()
;


text:
#########################Initial Condition Updating#########################################
# TODO: Consider when more than one UC model is used for the stages that the counts need
# to be scaled.
function calculate_ic_quantity(
    ::ICKey{TimeDurationON, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.Component}
    cache = get_cache(simulation_cache, ic.cache_type, T)
    name = device_name(ic)
    time_cache = cache_value(cache, name)

    current_counter = time_cache[:count]
    last_status = time_cache[:status]
    var_status = isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
    @debug last_status, var_status, abs(last_status - var_status)
    @assert abs(last_status - var_status) < ABSOLUTE_TOLERANCE

    return last_status >= 1.0 ? current_counter : 0.0
end

function calculate_ic_quantity(
    ::ICKey{TimeDurationOFF, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.Component}
    cache = get_cache(simulation_cache, ic.cache_type, T)
    name = device_name(ic)
    time_cache = cache_value(cache, name)

    current_counter = time_cache[:count]
    last_status = time_cache[:status]
    var_status = isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
    @debug last_status, var_status, abs(last_status - var_status)
    @assert abs(last_status - var_status) < ABSOLUTE_TOLERANCE

    return last_status >= 1.0 ? 0.0 : current_counter
end

function calculate_ic_quantity(
    ::ICKey{DeviceStatus, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.Component}
    current_status = isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
    return current_status
end

function calculate_ic_quantity(
    ::ICKey{DevicePower, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.ThermalGen}
    cache = get_cache(simulation_cache, TimeStatusChange, T)
    # This code determines if there is a status change in the generators. Takes into account TimeStatusChange for the presence of UC stages.
    dev = get_device(ic)
    min_power = PSY.get_active_power_limits(dev).min
    if isnothing(cache)
        # Transitions can't be calculated without cache
        status_change_to_on =
            get_condition(ic) <= min_power && var_value >= ABSOLUTE_TOLERANCE
        status_change_to_off =
            get_condition(ic) >= min_power && var_value <= ABSOLUTE_TOLERANCE
        status_remains_off =
            get_condition(ic) <= min_power && var_value <= ABSOLUTE_TOLERANCE
        status_remains_on =
            get_condition(ic) >= min_power && var_value >= ABSOLUTE_TOLERANCE
    else
        # If the min is 0.0 this calculation doesn't matter
        name = device_name(ic)
        time_cache = cache_value(cache, name)
        series = time_cache[:series]
        if min_power > 0.0
            current = min(time_cache[:current], length(series))
            current_status = isapprox(series[current], 1.0; atol = ABSOLUTE_TOLERANCE)
            previous_status = isapprox(series[current - 1], 1.0; atol = ABSOLUTE_TOLERANCE)
            status_change_to_on = current_status && !previous_status
            status_change_to_off = !current_status && previous_status
            status_remains_on = current_status && previous_status
            status_remains_off = !current_status && !previous_status
        else
            status_remains_on = true
            status_remains_off = false
            status_change_to_off = false
            status_change_to_on = false
        end
        time_cache[:current] += 1
    end

    if status_remains_off
        return 0.0
    elseif status_change_to_off
        return 0.0
    elseif status_change_to_on
        return min_power
    elseif status_remains_on
        return var_value
    else
        @assert false
    end
end

function calculate_ic_quantity(
    ::ICKey{DevicePower, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.Device}
    return var_value
end

function calculate_ic_quantity(
    ::ICKey{EnergyLevel, T},
    ic::InitialCondition,
    var_value::Float64,
    simulation_cache::Dict{<:CacheKey, AbstractCache},
) where {T <: PSY.Device}
    cache = get_cache(simulation_cache, ic.cache_type, T)
    name = device_name(ic)
    energy_cache = cache_value(cache, name)
    if energy_cache != var_value
        return var_value
    end
    return energy_cache
end

############################# Initial Conditions Initialization ############################
function _make_initial_conditions!(
    psi_container::PSIContainer,
    devices::Union{IS.FlattenIteratorWrapper{T}, Vector{T}},
    key::ICKey,
    make_ic_func::Function, # Function to make the initial condition object
    get_val_func::Function, # Function to get the value from the device to intialize
    cache = nothing,
) where {T <: PSY.Component}
    length_devices = length(devices)
    parameters = model_has_parameters(psi_container)
    ic_container = get_initial_conditions(psi_container)
    if !has_initial_conditions(ic_container, key)
        @debug "Setting $(key.ic_type) initial conditions for all devices $(T) based on system data"
        ini_conds = Vector{InitialCondition}(undef, length_devices)
        set_initial_conditions!(ic_container, key, ini_conds)
        for (ix, dev) in enumerate(devices)
            val_ = get_val_func(dev, key)
            val = parameters ? PJ.add_parameter(psi_container.JuMPmodel, val_) : val_
            ic = make_ic_func(ic_container, dev, val, cache)
            ini_conds[ix] = ic
            @debug "set initial condition" key ic val_
        end
    else
        ini_conds = get_initial_conditions(ic_container, key)
        ic_devices = Set((IS.get_uuid(ic.device) for ic in ini_conds))
        for dev in devices
            IS.get_uuid(dev) in ic_devices && continue
            @debug "Setting $(key.ic_type) initial conditions device $(PSY.get_name(dev)) based on system data"
            val_ = get_val_func(dev, key)
            val = parameters ? PJ.add_parameter(psi_container.JuMPmodel, val_) : val_
            ic = make_ic_func(ic_container, dev, val, cache)
            push!(ini_conds, ic)
            @debug "set initial condition" key ic val_
        end
    end

    @assert length(ini_conds) == length_devices
    return
end

######################### Initialize Functions for ThermalGen ##############################
"""
Status Init is always calculated based on the Power Output of the device
This is to make it easier to calculate when the previous model doesn't
contain binaries. For instance, looking back on an ED model to find the
IC of the UC model
"""
function status_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    _make_initial_conditions!(
        psi_container,
        devices,
        ICKey(DeviceStatus, T),
        _make_initial_condition_active_power,
        _get_status_value,
    )

    return
end

function output_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    _make_initial_conditions!(
        psi_container,
        devices,
        ICKey(DevicePower, T),
        _make_initial_condition_active_power,
        _get_active_power_output_value,
    )
    return
end

function output_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.ThermalMultiStart},
)
    _make_initial_conditions!(
        psi_container,
        devices,
        ICKey(DevicePower, PSY.ThermalMultiStart),
        _make_initial_condition_active_power,
        _get_active_power_output_above_min_value,
    )
end

function duration_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    for key in (ICKey(TimeDurationON, T), ICKey(TimeDurationOFF, T))
        _make_initial_conditions!(
            psi_container,
            devices,
            key,
            _make_initial_condition_active_power,
            _get_duration_value,
            TimeStatusChange,
        )
    end

    return
end

######################### Initialize Functions for Storage #################################
# TODO: This IC needs a cache for Simulation over long periods of tim
function storage_energy_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.Storage}
    key = ICKey(EnergyLevel, T)
    _make_initial_conditions!(
        psi_container,
        devices,
        key,
        _make_initial_condition_energy,
        _get_initial_energy_value,
    )

    return
end

######################### Initialize Functions for Hydro #################################
function status_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    _make_initial_conditions!(
        psi_container,
        devices,
        ICKey(DeviceStatus, T),
        _make_initial_condition_active_power,
        _get_status_value,
        # Doesn't require Cache
    )
end

function output_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    _make_initial_conditions!(
        psi_container,
        devices,
        ICKey(DevicePower, T),
        _make_initial_condition_active_power,
        _get_active_power_output_value,
        # Doesn't require Cache
    )

    return
end

function duration_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    for key in (ICKey(TimeDurationON, T), ICKey(TimeDurationOFF, T))
        _make_initial_conditions!(
            psi_container,
            devices,
            key,
            _make_initial_condition_active_power,
            _get_duration_value,
            TimeStatusChange,
        )
    end

    return
end

function storage_energy_init(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    key = ICKey(EnergyLevel, T)
    _make_initial_conditions!(
        psi_container,
        devices,
        key,
        _make_initial_condition_reservoir_energy,
        _get_reservoir_energy_value,
        StoredEnergy,
    )

    return
end

function area_control_init(psi_container::PSIContainer, services::Vector{PSY.AGC})
    key = ICKey(AreaControlError, PSY.AGC)
    _make_initial_conditions!(
        psi_container,
        services,
        key,
        _make_initial_condition_area_control,
        _get_ace_error,
        # Doesn't require Cache
    )

    return
end

function _make_initial_condition_active_power(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_active_power(T, container), value, cache)
end

function _make_initial_condition_energy(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_energy(T, container), value, cache)
end

function _make_initial_condition_reservoir_energy(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_reservoir_energy(T, container), value, cache)
end

function _make_initial_condition_area_control(
    container,
    device::PSY.AGC,
    value,
    cache = nothing,
)
    return InitialCondition(device, _get_ref_ace_error(PSY.AGC, container), value, cache)
end

function _get_status_value(device, key)
    return PSY.get_status(device) ? 1.0 : 0.0
end

function _get_active_power_output_value(device, key)
    if !PSY.get_status(device)
        return 0.0
    end
    return PSY.get_active_power(device)
end

function _get_active_power_output_above_min_value(device, key)
    if !PSY.get_status(device)
        return 0.0
    end
    power_above_min = PSY.get_active_power(device) - PSY.get_active_power_limits(device).min
    @assert power_above_min >= -ABSOLUTE_TOLERANCE
    return power_above_min
end

function _get_initial_energy_value(device, key)
    return PSY.get_initial_energy(device)
end

function _get_reservoir_energy_value(device, key)
    return PSY.get_initial_storage(device)
end

function _get_ace_error(device, key)
    return PSY.get_initial_ace(device)
end

function _get_duration_value(dev, key)
    if key.ic_type == TimeDurationON
        value = PSY.get_status(dev) ? PSY.get_time_at_status(dev) : 0.0
    else
        @assert key.ic_type == TimeDurationOFF
        value = !PSY.get_status(dev) ? PSY.get_time_at_status(dev) : 0.0
    end

    return value
end

function _get_ref_active_power(
    ::Type{T},
    container::InitialConditions,
) where {T <: PSY.Component}
    if get_use_parameters(container)
        return UpdateRef{JuMP.VariableRef}(T, ACTIVE_POWER)
    else
        return UpdateRef{T}(ACTIVE_POWER, "get_active_power")
    end
end

function _get_ref_energy(::Type{T}, container::InitialConditions) where {T <: PSY.Component}
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, ENERGY) :
           UpdateRef{T}(ENERGY, "get_initial_energy")
end

function _get_ref_reservoir_energy(
    ::Type{T},
    container::InitialConditions,
) where {T <: PSY.Component}
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, ENERGY) :
           UpdateRef{T}(ENERGY, "get_hydro_budget")
end
function _get_ref_ace_error(::Type{PSY.AGC}, container::InitialConditions)
    T = PSY.AGC
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, "ACE") :
           UpdateRef{T}("ACE", "get_initial_ace")
end


text:
using Test
using LinearAlgebra

using BEAST, CompScienceMeshes, SauterSchwabQuadrature, StaticArrays

t1 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.0, -0.980785, -0.19509],
    @SVector[0.0, -0.92388, -0.382683])
t2 = simplex(
    @SVector[0.373086, -0.881524, -0.289348],
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.294908, -0.944921, -0.141962])

# s1 = simplex(
#     @SVector[0.180878, -0.941848, -0.283207],
#     @SVector[0.0, -0.980785, -0.19509],
#     @SVector[0.0, -0.92388, -0.382683])
# s2 = simplex(
#     @SVector[0.180878, -0.941848, -0.283207],
#     @SVector[0.373086, -0.881524, -0.289348],
#     @SVector[0.294908, -0.944921, -0.141962])

# Common face case:
rt = BEAST.RTRefSpace{Float64}()
tqd = BEAST.quadpoints(rt, [t1,t2], (12,))
bqd = BEAST.quadpoints(rt, [t1,t2], (13,))

op1 = BEAST.MWSingleLayer3D(1.0, 250.0, 1.0)
op2 = BEAST.MWSingleLayer3D(1.0)
op3 = BEAST.MWDoubleLayer3D(0.0)

SE_strategy = BEAST.WiltonSERule(
  tqd[1,1],
  BEAST.DoubleQuadRule(
	tqd[1,1],
	bqd[1,1],
  ),
)
SS_strategy = SauterSchwabQuadrature.CommonFace(BEAST._legendre(8,0.0,1.0))

z_se = zeros(3,3)
z_ss = zeros(3,3)

BEAST.momintegrals!(op1, rt, rt, t1, t1, z_se, SE_strategy)
BEAST.momintegrals!(op1, rt, rt, t1, t1, z_ss, SS_strategy)

@test z_se ≈ z_ss atol=1e-4

z_se2 = zeros(3,3)
z_ss2 = zeros(3,3)

BEAST.momintegrals!(op2, rt, rt, t1, t1, z_se2, SE_strategy)
BEAST.momintegrals!(op2, rt, rt, t1, t1, z_ss2, SS_strategy)

@test z_se2 ≈ z_ss2 atol=1e-4

## Common Vertex Case
SE_strategy = BEAST.WiltonSERule(
  tqd[1,1],
  BEAST.DoubleQuadRule(
	tqd[1,1],
	bqd[1,2]))
SS_strategy = SauterSchwabQuadrature.CommonVertex(BEAST._legendre(12,0.0,1.0))

z_cv_se_1 = zeros(3,3)
z_cv_ss_1 = zeros(3,3)
z_cv_se_2 = zeros(3,3)
z_cv_ss_2 = zeros(3,3)
z_cv_se_3 = zeros(3,3)
z_cv_ss_3 = zeros(3,3)

BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_se_1, SE_strategy)
BEAST.momintegrals!(op1, rt, rt, t1, t2, z_cv_ss_1, SS_strategy)

@test z_cv_se_1 ≈ z_cv_ss_1 atol=1e-7

BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_se_2, SE_strategy)
BEAST.momintegrals!(op2, rt, rt, t1, t2, z_cv_ss_2, SS_strategy)

@test z_cv_se_2 ≈ z_cv_ss_2 atol=1e-7

BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_se_3, SE_strategy)
BEAST.momintegrals!(op3, rt, rt, t1, t2, z_cv_ss_3, SS_strategy)

@show z_cv_se_3
@show z_cv_ss_3
@test z_cv_se_3 ≈ z_cv_ss_3 atol=1e-7

## Common Edge Case:
t1 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.0, -0.980785, -0.19509],
    @SVector[0.0, -0.92388, -0.382683],
    )
t2 = simplex(
    @SVector[0.180878, -0.941848, -0.283207],
    @SVector[0.158174, -0.881178, -0.44554],
    @SVector[0.0, -0.92388, -0.382683],
    )

tqd = BEAST.quadpoints(rt, [t1,t2], (12,))
bqd = BEAST.quadpoints(rt, [t1,t2], (13,))

SE_strategy = BEAST.WiltonSERule(
    tqd[1,1],
    BEAST.DoubleQuadRule(
        tqd[1,1],
        bqd[1,2]))
SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(12,0.0,1.0))

z_ce_se_1 = zeros(3,3)
z_ce_ss_1 = zeros(3,3)
z_ce_se_2 = zeros(3,3)
z_ce_ss_2 = zeros(3,3)
z_ce_se_3 = zeros(3,3)
z_ce_ss_3 = zeros(3,3)

BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_se_1, SE_strategy)
BEAST.momintegrals!(op1, rt, rt, t1, t2, z_ce_ss_1, SS_strategy)

@show norm(z_ce_se_1 - z_ce_ss_1)
@test z_ce_se_1 ≈ z_ce_ss_1 atol=1e-5

BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_se_2, SE_strategy)
BEAST.momintegrals!(op2, rt, rt, t1, t2, z_ce_ss_2, SS_strategy)

@show norm(z_ce_se_2 - z_ce_ss_2)
@test z_ce_se_2 ≈ z_ce_ss_2 atol=1e-5

BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_se_3, SE_strategy)
BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3, SS_strategy)
@show norm(z_ce_se_3 - z_ce_ss_3)
@test z_ce_se_3 ≈ z_ce_ss_3 atol=1e-5

SS_strategy = SauterSchwabQuadrature.CommonEdge(BEAST._legendre(18,0.0,1.0))
z_ce_ss_3_18 = zeros(3,3)
BEAST.momintegrals!(op3, rt, rt, t1, t2, z_ce_ss_3_18, SS_strategy)
@show norm(z_ce_ss_3 - z_ce_ss_3_18)
@test z_ce_ss_3 ≈ z_ce_ss_3_18 atol=1e-14


text:
using RecipesBase
using Colors
using .Plots

struct NgSpiceGraphs end

const graph = NgSpiceGraphs()

# `fort` is frequency or time vector. It is the `X` axis and should always be real.
@recipe function ng(::NgSpiceGraphs, vectype, fort, veclist, title, xlims=nothing)
    ft = getrealvec(fort)
    label := permutedims(veclist)
    grid := true
    title := "$title"
    background_color --> Colors.RGB(0.0, 0.0, 0.0)
    xguide := fort
    guide := length(veclist) > 1 ? "signals" : "signal"
    color_palette --> :default
    seriestype := :path
    overwrite_figure := false
    vec = vectype.(veclist)
    #isa(xlims, Nothing) || xlims := xlims
    ft, vec
end

# source function for sp files needs to specially process plot commands
function source_sp(netpath)
    netlines = readlines(netpath)
    circpath = joinpath(pwd(), "net.cir")
    for l in range(1, stop=length(netlines))
        if netlines[l] == ".control"
            open(circpath, "a") do f
                for line in netlines[1:l-1]
                    write(f, line*"\n")
                end
            end
            try source(circpath) catch; rm(circpath) end
            temp = l
            while netlines[temp+1] != ".endc"
                temp += 1
                try
                    if occursin("plot", netlines[temp])
                        _plotswitch(netlines[temp])
                    else
                        cmd(netlines[temp])
                    end
                catch
                    isfile(circpath) && rm(circpath)
                end
            end
        end
    end
    isfile(circpath) && rm(circpath)
end

function _plotswitch(pltstr)
    params = split(pltstr[6:end], " ")
    xlims = nothing
    xlim_index = findall(occursin.("xlim", params))
    length(xlim_index) == 0 || (xlim_index = xlim_index[1]; xlims = parse.(Float64, (params[xlim_index+1], params[xlim_index+2]));
        deleteat!(params, xlim_index:1:xlim_index+2))
    fort = "frequency" ∈  listcurvecs() ? "frequency" : "time"
    pushfirst!(params, fort)
    if !occursin("-", params[1]) plt = plot(graph, getrealvec, params[1], params[2:end], "Real plot", xlims)
    elseif params[1] ∈ ("--real", "-r") plt = plot(graph, getrealvec, params[2], params[3:end], "Real plot", xlims)
    elseif params[1] ∈ ("--imaginary", "-i") plt = plot(graph, getimagvec, params[2], params[3:end], "Imaginary plot", xlims)
    elseif params[1] ∈ ("--magnitude", "-m") plt = plot(graph, getmagnitudevec, params[2], params[3:end], "Magnitude plot", xlims)
    elseif params[1] ∈ ("--phase", "-p") plt = plot(graph, getphasevec, params[2], params[3:end], "Phase plot", xlims)
    else throw("$(params[1]) is not a valid plot type")
    end
    Base.display(plt)
end

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
# Check literals.
_lit(v::Variable) = false
_lit(v::Type) = false
_lit(v) = true

# Lift literals to types.
_lift(t::Type) = t
_lift(t) = typeof(t)

# Lower a function signature directly to IR.
function lower_to_ir(call, argtypes...)
    sig = length(argtypes) == 1 && argtypes[1] == Tuple{} ? begin
        Tuple{typeof(call)}
    end : Tuple{typeof(call), argtypes...}
    m = meta(sig)
    ir = IR(m)
    return ir
end

# Returns the signature of the matching method, specializing types as necessary.
function signature(fn_type::Type, arg_types::Type...)
    meta = IRTools.meta(Tuple{fn_type, arg_types...})
    if isnothing(meta) return nothing end
    sig, typevars = unparameterize(meta.method.sig)
    sig_types = collect(sig.parameters[2:end])
    for (i, a) in enumerate(arg_types)
        if sig_types[i] isa TypeVar
            sig_types[i] = a
        else
            sig_types[i] = reparameterize(sig_types[i], typevars)
        end
    end
    return sig_types
end

# Unwraps UnionAll types, returning the innermost body with a list of TypeVars.
function unparameterize(@nospecialize(T::Type))
    vars = TypeVar[]
    while T isa UnionAll
        pushfirst!(vars, T.var)
        T = T.body
    end
    return T, vars
end

# Rewraps parametric types with the provided TypeVars.
function reparameterize(@nospecialize(T::Type), vars::Vector{TypeVar})
    if T isa UnionAll
        T, origvars = unparameterize(T)
        vars = unique!([origvars; vars])
    end
    parameters = Set(T.parameters)
    for v in vars
        if (v in parameters) T = UnionAll(v, T) end
    end
    return T
end

# Returns a Vector of Tuple{Vararg{Int}} representing all possible combinations of flow of control.
function get_control_flow_paths(blk::Int, v::Vector{Vector{Int}})
    isempty(v[blk]) && return [(blk, )]
    paths = []
    for tar in filter(b -> b > blk, v[blk]) # prevent looping
        append!(paths, map(get_control_flow_paths(tar, v)) do p
                    (blk, p...)
                end)
    end
    paths
end
function get_control_flow_paths(cfg::CFG)
    graph = cfg.graph
    paths = get_control_flow_paths(1, graph)
    paths
end
@inline get_control_flow_paths(ir::IR) = get_control_flow_paths(CFG(ir))

# Check for control flow in IR - if multiple basic blocks, transfer of control is present.
@inline control_flow_check(ir) = !(length(ir.blocks) > 1)


text:
using Primes

# These primes are good example primes, as they have p = ±4 (mod 9), so taking cube roots is easy in F_p^2
plist = [5, 13, 23, 31, 41, 59, 67]
bigplist = [big(2)^107 - 1, big(2)^521 - 1]

# The prime used in the OG SIDH paper. Does not have above property^
OGp = big(2)^63*big(3)^41*11-1

# If p = 2 (mod 3) we have an easy supersingular curve
# E = Affine_EC_TwiHes(F, inv(4*one(F)), zero(F))
# Or twisted to
# E = Affine_EC_TwiHes(F, sqrt(F(b))*inv(4*one(F)), zero(F))
# This is isomorphic to the Weierstrass curve y^2 = x^3 + 1

# Similarly (slightly worse though) If p = 3 (mod 4) we have an easy supersingular curve
# y^2 = x^3 - x
# Translating this to twisted Hessian form will be future work

# If p = 1 (mod 12), just run away, dont look back.

#Scratch all this, much easier than that...
#If p = 2 (mod 3), then 0 is a supersingular j-invariant
#If p = 3 (mod 4), then 1728 is a supersingular j-invariant

#Example:
OGp = big(2)^63*big(3)^41*11-1
F = Fpsquare(OGp)
E = TwiHesEC(F, inv(4*one(F)), zero(F))
#w = F([1850222081870264162797200520152908562431, 2673953603825912900316077551548884232036])
#MAYBE: all points of order two are of the form (X : 1 : 1)
#Want to find roots of x^3-a = 0, where a = 3700444163740528325594401040305817124855
#One lies in F_p, the other two are the roots of x^2 + 3700444163740528325594401040305817124861*x + 4
#Roots are r1 = 3700444163740528325594401040305817124862 + 1647463043911297475037754062791951339209*i
#and r2 = 3700444163740528325594401040305817124862 + 2052981119829230850556646977513865785654*i
#r1 = F([3700444163740528325594401040305817124862, 1647463043911297475037754062791951339209])
#P = E(-r1, one(F), one(F))
#E0 = Isogeny(E,P).codomain
#c = F([3398654792751657841130893587213062545659])*w
#Q = E0(zero(F), -w, one(F))
#E1 = Isogeny(E0, Q).codomain

#Simple example with isogeny over F_p
F = Zmod(31)
E = TwiHesEC(F, F(1), F(7))
listP = allPointsNaive(E)
phi = Isogeny(E, E(F(13), F(1), F(1)))

#Simple example with isogeny over F_p^2
F = Fpsquare(31)
E = TwiHesEC(F, F([4, 5]), F([1,2]))
listP = allPointsNaive(E)
phi = Isogeny(E, E(F([9,19]), F([1]), F([1])))


text:
"""
Provides the [`crossref`](@ref) function used to automatically calculate link URLs.
"""
module CrossReferences

import ..Documenter:
    Anchors,
    Builder,
    Documents,
    Expanders,
    Documenter,
    Utilities

using DocStringExtensions
import Markdown

"""
$(SIGNATURES)

Traverses a [`Documents.Document`](@ref) and replaces links containg `@ref` URLs with
their real URLs.
"""
function crossref(doc::Documents.Document)
    for (src, page) in doc.blueprint.pages
        empty!(page.globals.meta)
        for element in page.elements
            crossref(page.mapping[element], page, doc)
        end
    end
end

function crossref(elem, page, doc)
    Documents.walk(page.globals.meta, elem) do link
        xref(link, page.globals.meta, page, doc)
    end
end

# Dispatch to `namedxref` / `docsxref`.
# -------------------------------------

const NAMED_XREF = r"^@ref (.+)$"

function xref(link::Markdown.Link, meta, page, doc)
    link.url == "@ref"             ? basicxref(link, meta, page, doc) :
    occursin(NAMED_XREF, link.url) ? namedxref(link, meta, page, doc) : nothing
    return false # Stop `walk`ing down this `link` element.
end
xref(other, meta, page, doc) = true # Continue to `walk` through element `other`.

function basicxref(link::Markdown.Link, meta, page, doc)
    if length(link.text) === 1 && isa(link.text[1], Markdown.Code)
        docsxref(link, link.text[1].code, meta, page, doc)
    elseif isa(link.text, Vector)
        # No `name` was provided, since given a `@ref`, so slugify the `.text` instead.
        text = strip(sprint(Markdown.plain, Markdown.Paragraph(link.text)))
        if occursin(r"#[0-9]+", text)
            issue_xref(link, lstrip(text, '#'), meta, page, doc)
        else
            name = Utilities.slugify(text)
            namedxref(link, name, meta, page, doc)
        end
    end
end

# Cross referencing headers.
# --------------------------

function namedxref(link::Markdown.Link, meta, page, doc)
    # Extract the `name` from the `(@ref ...)`.
    slug = match(NAMED_XREF, link.url)[1]
    if isempty(slug)
        text = sprint(Markdown.plaininline, link)
        push!(doc.internal.errors, :cross_references)
        @warn "'$text' missing a name after '#' in $(Utilities.locrepr(page.source))."
    else
        if Anchors.exists(doc.internal.headers, slug)
            namedxref(link, slug, meta, page, doc)
        elseif length(link.text) === 1 && isa(link.text[1], Markdown.Code)
            docsxref(link, slug, meta, page, doc)
        else
            namedxref(link, slug, meta, page, doc)
        end
    end
end

function namedxref(link::Markdown.Link, slug, meta, page, doc)
    headers = doc.internal.headers
    # Add the link to list of local uncheck links.
    doc.internal.locallinks[link] = link.url
    # Error checking: `slug` should exist and be unique.
    # TODO: handle non-unique slugs.
    if Anchors.exists(headers, slug)
        if Anchors.isunique(headers, slug)
            # Replace the `@ref` url with a path to the referenced header.
            anchor   = Anchors.anchor(headers, slug)
            path     = relpath(anchor.file, dirname(page.build))
            link.url = string(path, '#', slug, '-', anchor.nth)
        else
            push!(doc.internal.errors, :cross_references)
            @warn "'$slug' is not unique in $(Utilities.locrepr(page.source))."
        end
    else
        push!(doc.internal.errors, :cross_references)
        @warn "reference for '$slug' could not be found in $(Utilities.locrepr(page.source))."
    end
end

# Cross referencing docstrings.
# -----------------------------

function docsxref(link::Markdown.Link, code, meta, page, doc)
    # Add the link to list of local uncheck links.
    doc.internal.locallinks[link] = link.url
    # Parse the link text and find current module.
    keyword = Symbol(strip(code))
    local ex
    if haskey(Docs.keywords, keyword)
        ex = QuoteNode(keyword)
    else
        try
            ex = Meta.parse(code)
        catch err
            !isa(err, Meta.ParseError) && rethrow(err)
            push!(doc.internal.errors, :cross_references)
            @warn "unable to parse the reference '[`$code`](@ref)' in $(Utilities.locrepr(page.source))."
            return
        end
    end
    mod = get(meta, :CurrentModule, Main)

    # Find binding and type signature associated with the link.
    local binding
    try
        binding = Documenter.DocSystem.binding(mod, ex)
    catch err
        push!(doc.internal.errors, :cross_references)
        @warn "unable to get the binding for '[`$code`](@ref)' in $(Utilities.locrepr(page.source)) from expression '$(repr(ex))' in module $(mod)" exception = err
        return
    end

    local typesig
    try
        typesig = Core.eval(mod, Documenter.DocSystem.signature(ex, rstrip(code)))
    catch err
        push!(doc.internal.errors, :cross_references)
        @warn "unable to evaluate the type signature for '[`$code`](@ref)' in $(Utilities.locrepr(page.source)) from expression '$(repr(ex))' in module $(mod)" exception = err
        return
    end

    # Try to find a valid object that we can cross-reference.
    object = find_object(doc, binding, typesig)
    if object !== nothing
        # Replace the `@ref` url with a path to the referenced docs.
        docsnode = doc.internal.objects[object]
        path     = relpath(docsnode.page.build, dirname(page.build))
        slug     = Utilities.slugify(object)
        link.url = string(path, '#', slug)
    else
        push!(doc.internal.errors, :cross_references)
        @warn "no doc found for reference '[`$code`](@ref)' in $(Utilities.locrepr(page.source))."
    end
end

"""
$(SIGNATURES)

Find the included `Object` in the `doc` matching `binding` and `typesig`. The matching
heuristic isn't too picky about what matches and will only fail when no `Binding`s matching
`binding` have been included.
"""
function find_object(doc::Documents.Document, binding, typesig)
    object = Utilities.Object(binding, typesig)
    if haskey(doc.internal.objects, object)
        # Exact object matching the requested one.
        return object
    else
        objects = get(doc.internal.bindings, binding, Utilities.Object[])
        if isempty(objects)
            # No bindings match the requested object == FAILED.
            return nothing
        elseif length(objects) == 1
            # Only one possible choice. Use it even if the signature doesn't match.
            return objects[1]
        else
            candidate = find_object(binding, typesig)
            if candidate in objects
                # We've found an actual match out of the possible choices! Use it.
                return candidate
            else
                # No match in the possible choices. Use the one that was first included.
                return objects[1]
            end
        end
    end
end
function find_object(binding, typesig)
    if Documenter.DocSystem.defined(binding)
        local λ = Documenter.DocSystem.resolve(binding)
        return find_object(λ, binding, typesig)
    else
        return Utilities.Object(binding, typesig)
    end
end
function find_object(λ::Union{Function, DataType}, binding, typesig)
    if hasmethod(λ, typesig)
        signature = getsig(λ, typesig)
        return Utilities.Object(binding, signature)
    else
        return Utilities.Object(binding, typesig)
    end
end
find_object(::Union{Function, DataType}, binding, ::Union{Union,Type{Union{}}}) = Utilities.Object(binding, Union{})
find_object(other, binding, typesig) = Utilities.Object(binding, typesig)

getsig(λ::Union{Function, DataType}, typesig) = Base.tuple_type_tail(which(λ, typesig).sig)

# Issues/PRs cross referencing.
# -----------------------------

function issue_xref(link::Markdown.Link, num, meta, page, doc)
    link.url = isempty(doc.internal.remote) ? link.url :
        "https://github.com/$(doc.internal.remote)/issues/$num"
end

end


text: using FranklinCLI; FranklinCLI.comonicon_install()