text:
# Adapted from
# https://github.com/JuliaData/AbstractTables.jl/blob/e5afc569504ecf08ec769fd52a78da4027eab35f/src/AbstractTable/show.jl

let
    local io = IOBuffer(Array(UInt8, 80), true, true)
    global ourstrwidth
    function ourstrwidth(x::Any)::Int
        truncate(io, 0)
        ourshowcompact(io, x)
        return position(io)
    end
    ourstrwidth(x::AbstractString) = strwidth(x) + 2 # -> Int
    # myconv = VERSION < v"0.4-" ? convert : Base.unsafe_convert
    ourstrwidth(s::Symbol) =
        @compat Int(
            ccall(:u8_strwidth, Csize_t, (Ptr{UInt8}, ),
            Base.unsafe_convert(Ptr{UInt8}, s))
        )
end

ourshowcompact(io::IO, x::Any)::Void = showcompact(io, x)
ourshowcompact(io::IO, x::AbstractString)::Void = showcompact(io, x)
ourshowcompact(io::IO, x::Symbol)::Void = print(io, x)

function getmaxwidths(r::Relation, rowlabel, limit, offset)::Vector{Int}
    attrs = attributes(r)
    ncols = length(attrs)
    widths = [ Vector{Int}() for j in 1:ncols ]
    maxwidths = Array{Int}(ncols + 1)
    undefstrwidth = ourstrwidth(Base.undef_ref_str)

    s = start(r)
    i = 1
    while (i <= offset) & (!done(r, s)) # skip through rows up to offset
        row, s = next(r, s)
        i += 1
    end
    while (i <= offset + limit) & (!done(r, s)) #
        row, s = next(r, s)
        for (j, v) in enumerate(row)
            try
                push!(widths[j], ourstrwidth(v))
            catch
                push!(widths[j], undefstrwidth)
            end
        end
        i += 1
    end
    for j in 1:ncols
        if isempty(widths[j]) # no rows
            maxwidths[j] = ourstrwidth(attrs[j])
        else # max width in column is max of max width of printed el and attribute width
            maxwidths[j] = max(maximum(widths[j]), ourstrwidth(attrs[j]))
        end
    end
    maxwidths[end] = max(ourstrwidth(rowlabel), ndigits(limit)+1)
    return maxwidths
end

function getprintedwidth(maxwidths::Vector{Int})::Int
    # Include length of line-initial |
    totalwidth = 1
    for i in 1:length(maxwidths)
        # Include length of field + 2 spaces + trailing |
        totalwidth += maxwidths[i] + 3
    end
    return totalwidth
end

pad(io, padding)::Void = print(io, " "^padding)

function print_bounding_line(io, maxwidths, j_left, j_right)::Void
    rowmaxwidth = maxwidths[end]
    print(io, '├', "─"^(rowmaxwidth + 2), '┼')
    for j in j_left:j_right
        print(io, "─"^(maxwidths[j] + 2))
        if j < j_right
            print(io, '┼')
        else
            print(io, '┤')
        end
    end
    print(io, '\n')
end

# NOTE: returns a Bool indicating whether or not there are more rows than what
# has printed
function print_rows(
  io, r, maxwidths, j_left, j_right, rowlabel, limit, offset)::Bool ##
    rowmaxwidth = maxwidths[end]
    attrs = attributes(r)
    p = project(r, attrs[j_left:j_right]...)
    s, i = start(p), 1
    while (i <= offset) & (!done(p, s)) # skip through offset
        row, s = next(p, s)
        i += 1
    end
    while (i <= offset + limit) & (!done(p, s))
        row, s = next(p, s)
        @printf(io, "│ %d", i)
        pad(io, rowmaxwidth - ndigits(i))
        # print(io, " │ ")
        print(io, "   ")
        # print entry
        for j in j_left:j_right
            v = getfield(row, j)
            strlen = ourstrwidth(v)
            ourshowcompact(io, v)
            pad(io, maxwidths[j] - strlen)
            if j == j_right
                if i == (limit+offset)
                    # print(io, " │")
                    print(io, "  ")
                else
                    # print(io, " │\n")
                    print(io, "  \n")
                end
            else
                print(io, "   ")
            end
        end
        i += 1
    end
    return (i > (limit+offset)) & (!done(p, s))
end

function print_footer(io, r, more_rows, j_right, limit, offset)::Void
    println(io, "\n⋮")
    attrs = attributes(r)
    ncols = length(attrs)
    if j_right < ncols
        @printf(io, "with more rows and %d more columns: ", ncols - j_right)
        for j in j_right+1:(ncol-1)
            attr = attrs[j]
            print(io, "$attr, ")
        end
        print(io, "$attr.")
    else
        print(io, "with more rows.")
    end
end

function print_footer{S<:Attributed}(
  io, r::Relation{S}, more_rows, j_right, limit, offset)::Void
    attrs = attributes(r)
    ncols, nrows = length(attrs), length(r)
    if more_rows
        println(io, "\n⋮")
        if j_right < ncols
            if offset > 0
                @printf(io, "with %d more rows (skipped the first %d rows) and %d more columns: ",
                        nrows-(limit+offset), offset, ncols-j_right)
            else
                @printf(io, "with %d more rows and %d more columns: ",
                        nrows-limit, ncols-j_right)
            end
            for j in j_right+1:(ncols-1)
                attr = attrs[j]; print(io, "$attr, ")
            end
            attr = attrs[ncols]; print(io, "$attr.")
        else
            if offset > 0
                @printf(io, "with %d more rows (skipped the first %d rows).",
                        nrows-(limit+offset), offset)
            else
                @printf(io, "with %d more rows.", nrows-limit)
            end
        end
    else
        if j_right < ncols
            println(io, "\n⋮")
            @printf(io, "with %d more columns: ", ncols-j_right)
            for j in j_right+1:(ncols-1)
                attr = attrs[j]
                print(io, "$attr, ")
            end
            attr = attrs[ncols]
            print(io, "$attr.")
        end
    end
    return
end

function print_header(io, r, maxwidths, j_left, j_right, rowlabel)::Void
    rowmaxwidth = maxwidths[end]
    attrs = attributes(r)
    @printf(io, "│ %s", rowlabel)
    pad(io, rowmaxwidth - ourstrwidth(rowlabel))
    # print(io, " │ ")
    print(io, "   ")
    for j in j_left:j_right
        attr = attrs[j]
        ourshowcompact(io, attr)
        pad(io, maxwidths[j] - ourstrwidth(attr))
        # j == j_right ? print(io, " │\n") : print(io, " │ ")
        j == j_right ? print(io, "  \n") : print(io, "   ")
    end
    print_bounding_line(io, maxwidths, j_left, j_right)
    return
end

function getchunkbounds(
  io, maxwidths::Vector{Int}, splitchunks::Bool,
  availablewidth::Int=displaysize(io)[2])::Vector{Int}
    ncols = length(maxwidths) - 1
    rowmaxwidth = maxwidths[end]
    if splitchunks
        chunkbounds = [0]
        # Include 2 spaces + 2 | characters for row/col label
        totalwidth = rowmaxwidth + 4
        for j in 1:ncols
            # Include 2 spaces + | character in per-column character count
            totalwidth += maxwidths[j] + 3
            if totalwidth > availablewidth
                push!(chunkbounds, j - 1)
                totalwidth = rowmaxwidth + 4 + maxwidths[j] + 3
            end
        end
        push!(chunkbounds, ncols)
    else
        chunkbounds = [0, ncols]
    end
    return chunkbounds
end

# 1 space for line-initial | + length of field + 2 spaces + trailing |
printedwidth(maxwidths)::Void = foldl((x,y)->x+y+3, 1, maxwidths)

function _show{S<:Attributed}(
  io, r::Relation{S}, rowlabel=Symbol("#"), displaysummary=true,
  splitchunks=true, allchunks=false, limit=10, offset=0)::Void ##
    ncols = length(attributes(r))
    if ncols > 0
        displaysummary && println(io, summary(r))
        # println("⋮")
        maxwidths = getmaxwidths(r, rowlabel, limit, offset)
        chunkbounds = getchunkbounds(io, maxwidths, splitchunks | !allchunks)
        if !allchunks # show only the first chunk
            j_left, j_right = 1, chunkbounds[2]
            print_header(io, r, maxwidths, j_left, j_right, rowlabel)
            more_rows = print_rows(
                io, r, maxwidths, j_left, j_right, rowlabel, limit, offset)
            print_footer(
                io, r, more_rows, j_right, limit, offset)
        else
            nchunks = length(chunkbounds) - 1
            more_rows = false
            for r in 1:nchunks
                j_left = chunkbounds[r] + 1
                j_right = chunkbounds[r + 1]
                print_tbl_header(io, r, maxwidths, j_left, j_right, rowlabel)
                more_rows |= print_rows(
                    io, r, maxwidths, j_left, j_right, rowlabel, limit, offset)
            end
            print_footer(
                io, r, more_rows, false, limit, offset)
        end
    else
        @printf(io, "An empty %s", typeof(tbl))
    end
    return
end

Base.summary(r::Relation) =
    @sprintf("Relation (source of type %s)", typeof(r.src))
Base.summary{S<:Attributed}(r::Relation{S}) = "Relation (in-memory Julia source)"
Base.show(io::IO, r::Relation) = _show(io, r)
Base.show(io::IO, ::MIME"text/plain", r::Relation) = _show(io, r)
Base.show(r::Relation; limit=10, offset=0) =
    _show(STDOUT, r, Symbol("#"), true, true, false, limit, offset)
# _show(r::Relation, limit, offset)::Void = _show(STDOUT, tbl, :Row, true, limit, offset)


text:
using WeaklyHard
"""
Example for comparing dominance and equivalence of WeaklyHard constraints
"""

# Constraints to test
lambda1 = AnyHitConstraint(48, 54)
lambda2 = AnyMissConstraint(12, 67)
lambda3 = AnyHitConstraint(12, 15)
lambda4 = RowMissConstraint(9)
lambda5 = RowHitConstraint(20, 47)
lambda6 = RowHitConstraint(5, 12)
lambda7 = RowMissConstraint(3)

lambdas = [lambda1, lambda2, lambda3, lambda4, lambda5, lambda6, lambda7]

# Compare individual constraints
for l1 in lambdas
    for l2 in lambdas
        if is_equivalent(l1, l2)
            println("$l1 \t===\t $l2")
        elseif is_dominant(l1, l2)
            println("$l1 \t<=\t $l2")
        else
            println("$l1 \t</=\t $l2")
        end # if
    end # for
    println(" ")
end # for

ext: jl


text:
using Plots
using DelimitedFiles

save = false

if save
    gr()
else
    plotly()
end

function LoadStylus(path, flattening)
    data = readdlm(path, ',', skipstart=34)[:, 1:2]

    correction = data[flattening, 2]
    for i in 1:length(data[:, 2])
        data[i, 2] = data[i, 2] - correction * i / flattening
    end

    data[:, 1], data[:, 2]
end

stylusfig = plot(title="Stylus measurements", xlabel="distance [μm]", ylabel="height [nm]", ticks=:native)
for sample in ["5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20"]
    xs, ys = LoadStylus("../NanoLab_process/20220405/S" * sample * ".csv", 16000)

    plot!(stylusfig, xs, ys, label=sample)
end
display(stylusfig)
if save savefig(stylusfig, "../../rapport/figures/profiles/stylus.pdf") end

stylusfig_contested = plot(title="Stylus measurements", xlabel="distance [μm]", ylabel="height [nm]", ticks=:native)
for sample in [
        "20220405/S5", "20220424/5A",
        "20220405/S8", "20220424/8A",
        "20220405/S9", "20220424/9A",
        "20220405/S10", "20220424/10A",
        "20220405/S12", "20220424/12A",
        "20220405/S13", "20220424/13A",
        "20220405/S14", "20220424/14A",
        "20220405/S18", "20220424/18A",
        "20220405/S19", "20220424/19A",
        "20220405/S20", "20220424/20A",
        ]
    xs, ys = LoadStylus("../NanoLab_process/" * sample * ".csv", 15500)

    plot!(stylusfig_contested, xs, ys, label=sample)
end
display(stylusfig_contested)
if save savefig(stylusfig_contested, "../../rapport/figures/profiles/stylus.pdf") end

stylusfig_extra = plot(title="Stylus measurements", xlabel="distance [μm]", ylabel="height [nm]", ticks=:native)
for sample in ["21", "22", "23", "24"]
    xs, ys = LoadStylus("../NanoLab_process/20220505/" * sample * "A.csv", 8000)

    plot!(stylusfig_extra, xs, ys, label=sample)
end
display(stylusfig_extra)
if save savefig(stylusfig_extra, "../../rapport/figures/profiles/stylus.pdf") end


text:
"""
Return default inscription value based on `PNTD`. Has meaning of unity, as in `one`.
"""
function default_inscription end
default_inscription(::PNTD) where {PNTD <: PnmlType} = one(Integer)
default_inscription(::PNTD) where {PNTD <: AbstractContinuousCore} = one(Float64)
default_inscription(pntd::PNTD) where {PNTD <: AbstractHLCore} = default_term(pntd)

#-------------------
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Labels an Arc.
"""
struct PTInscription{T<:Number}  <: Annotation
    value::T
    com::ObjectCommon
end

PTInscription() = PTInscription(one(Int))
PTInscription(value) = PTInscription(value, ObjectCommon())

"""
Evaluate a [`PTInscription`](@ref).
"""
(inscription::PTInscription)() = inscription.value

#
#-------------------
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Labels an Arc.
"""
struct HLInscription{TermType<:AbstractTerm} <: HLAnnotation
    text::Maybe{String}
    term::Maybe{TermType} # structure
    com::ObjectCommon
end

"""
$(TYPEDSIGNATURES)
"""

HLInscription() = HLInscription(nothing, Term(), ObjectCommon())
HLInscription(s::AbstractString) = HLInscription(s, Term())
HLInscription(t::Term) = HLInscription(nothing, t)
HLInscription(s::AbstractString, t::Term) = HLInscription(s, t, ObjectCommon())

"""
Evaluate a [`HLInscription`](@ref). Returns a value of the same sort as _TBD_.
"""
(hlm::HLInscription)() = _evaluate(hlm.term) #"HLInscription functor not implemented"


text:
# Raw memory management

import Base: copy!

export
    cualloc, cumemset, free

function cualloc{T, N<:Integer}(::Type{T}, len::N=1)
    if T.abstract || !T.isleaftype
        throw(ArgumentError("Cannot allocate pointer to abstract or non-leaf type"))
    end
    ptr_ref = Ref{Ptr{Void}}()
    nbytes = len * sizeof(T)
    if nbytes <= 0
        throw(ArgumentError("Cannot allocate $nbytes bytes of memory"))
    end
    @apicall(:cuMemAlloc, (Ptr{Ptr{Void}}, Csize_t), ptr_ref, nbytes)
    return unsafe_convert(DevicePtr{T}, reinterpret(Ptr{T}, ptr_ref[]))
end

cumemset(p::DevicePtr, value::Cuint, len::Integer) =
    @apicall(:cuMemsetD32, (Ptr{Void}, Cuint, Csize_t), p.inner, value, len)

free(p::DevicePtr) = @apicall(:cuMemFree, (Ptr{Void},), p.inner)

function copy!{T}(dst::DevicePtr{T}, src::T)
    if !Base.datatype_pointerfree(T)
        # TODO: recursive copy?
        throw(ArgumentError("Only pointer-free types can be copied"))
    end
    @apicall(:cuMemcpyHtoD, (Ptr{Void}, Ptr{Void}, Csize_t),
                            dst.inner, pointer_from_objref(src), sizeof(T))
    return nothing
end


text:
module Metaheuristics

using Distributions, Random, LinearAlgebra, Statistics
using Distributed

# Clearing
include(joinpath("clearing","struct.jl"))
include(joinpath("clearing","functions.jl"))
include(joinpath("clearing","optimize.jl"))
include(joinpath("clearing","utils.jl"))

end


text:
module ParametricOptInterface

using MathOptInterface

const MOI = MathOptInterface

const PARAMETER_INDEX_THRESHOLD = 1_000_000_000_000_000_000

@enum ConstraintsInterpretationCode ONLY_CONSTRAINTS ONLY_BOUNDS BOUNDS_AND_CONSTRAINTS

"""
    Parameter(val::Float64)

The `Parameter` structure stores the numerical value associated to a given
parameter.

## Example

```julia-repl
julia> ParametricOptInterface.Parameter(5)
ParametricOptInterface.Parameter(5)
```
"""
struct Parameter <: MOI.AbstractScalarSet
    val::Float64
end

# Utilities for using a CleverDict in Parameters
struct ParameterIndex
    index::Int64
end
function MOI.Utilities.CleverDicts.index_to_key(
    ::Type{ParameterIndex},
    index::Int64,
)
    return ParameterIndex(index)
end
function MOI.Utilities.CleverDicts.key_to_index(key::ParameterIndex)
    return key.index
end
function p_idx(vi::MOI.VariableIndex)
    return ParameterIndex(vi.value - PARAMETER_INDEX_THRESHOLD)
end

"""
    Optimizer{T, OT <: MOI.ModelLike} <: MOI.AbstractOptimizer

Declares a `Optimizer`, which allows the handling of parameters in a
optimization model.

## Keyword arguments

- `evaluate_duals::Bool`: If `true`, evaluates the dual of parameters. Users might want to set it to false
  to increase performance when the duals of parameters are not necessary. Defaults to `true`.

- `constraints_interpretation`: Decides how to interpret constraints with `ScalarAffineFunctions`. More details
  are in [`POI.ConstraintsInterpretation`](@ref) and [JuMP Example - Dealing with parametric expressions as variable bounds](@ref).
  Defaults to `ONLY_CONSTRAINTS`.

## Example

```julia-repl
julia> ParametricOptInterface.Optimizer(GLPK.Optimizer())
ParametricOptInterface.Optimizer{Float64,GLPK.Optimizer}
```
"""
mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT
    parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    }
    parameters_name::Dict{MOI.VariableIndex,String}
    # The updated_parameters dictionary has the same dimension of the
    # parameters dictionary and if the value stored is a NaN is means
    # that the parameter has not been updated.
    updated_parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    }
    variables::MOI.Utilities.CleverDicts.CleverDict{
        MOI.VariableIndex,
        MOI.VariableIndex,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    }
    last_variable_index_added::Int64
    last_parameter_index_added::Int64
    # Store the constraint function and set passed to POI by MOI
    original_constraint_function_and_set_cache::Dict{
        MOI.ConstraintIndex,
        Tuple{MOI.AbstractFunction,MOI.AbstractSet},
    }
    # Store the map for SAFs that might be transformed into VI
    affine_added_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        MOI.ConstraintIndex,
    }
    last_affine_added::Int64
    # Store reference to parameters of affine constraints with parameters: v + p
    affine_constraint_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarAffineTerm{Float64}},
    }
    # Store reference quadratic constraints with parameter * variable constraints: p * v
    quadratic_constraint_cache_pv::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarQuadraticTerm{Float64}},
    }
    # Store reference quadratic constraints with parameter * variable constraints: p * v
    quadratic_constraint_cache_pp::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarQuadraticTerm{Float64}},
    }
    # Store reference to constraints with quad_variable_term + affine_with_parameters: v * v + p
    quadratic_constraint_cache_pc::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarAffineTerm{Float64}},
    }
    # Store the reference to variables in the scalar affine part that are
    # multiplied by parameters in the scalar quadratic terms.
    # i.e.
    # If we have a constraint function with both scalar quadratic terms and
    # scalar affine terms such as p_1 * v_1 + 2.0 * v_1
    # When we need to update the constraint coefficient after updating the parameter
    # we must do (new_p_1 + 2.0) * v_1
    # This cache is storing the 2.0 * v_1 part.
    quadratic_constraint_variables_associated_to_parameters_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.ScalarAffineTerm{T}},
    }
    # Store the map for SQFs that might be transformed into SAF
    # for instance p*p + var -> ScalarAffine(var)
    quadratic_added_cache::Dict{MOI.ConstraintIndex,MOI.ConstraintIndex}
    last_quad_add_added::Int64
    vector_constraint_cache::MOI.Utilities.DoubleDicts.DoubleDict{
        Vector{MOI.VectorAffineTerm{Float64}},
    }
    # Ditto from the constraint caches but for the objective function
    original_objective_function::MOI.AbstractFunction
    affine_objective_cache::Vector{MOI.ScalarAffineTerm{T}}
    quadratic_objective_cache_pv::Vector{MOI.ScalarQuadraticTerm{T}}
    quadratic_objective_cache_pp::Vector{MOI.ScalarQuadraticTerm{T}}
    quadratic_objective_cache_pc::Vector{MOI.ScalarAffineTerm{T}}
    quadratic_objective_variables_associated_to_parameters_cache::Vector{
        MOI.ScalarAffineTerm{T},
    }
    multiplicative_parameters::Set{Int64}
    dual_value_of_parameters::Vector{Float64}
    evaluate_duals::Bool
    number_of_parameters_in_model::Int64
    constraints_interpretation::ConstraintsInterpretationCode
    function Optimizer(optimizer::OT; evaluate_duals::Bool = true) where {OT}
        return new{Float64,OT}(
            optimizer,
            MOI.Utilities.CleverDicts.CleverDict{ParameterIndex,Float64}(
                MOI.Utilities.CleverDicts.key_to_index,
                MOI.Utilities.CleverDicts.index_to_key,
            ),
            Dict{MOI.VariableIndex,String}(),
            MOI.Utilities.CleverDicts.CleverDict{ParameterIndex,Float64}(
                MOI.Utilities.CleverDicts.key_to_index,
                MOI.Utilities.CleverDicts.index_to_key,
            ),
            MOI.Utilities.CleverDicts.CleverDict{
                MOI.VariableIndex,
                MOI.VariableIndex,
            }(
                MOI.Utilities.CleverDicts.key_to_index,
                MOI.Utilities.CleverDicts.index_to_key,
            ),
            0,
            PARAMETER_INDEX_THRESHOLD,
            Dict{
                MOI.ConstraintIndex,
                Tuple{MOI.AbstractFunction,MOI.AbstractSet},
            }(),
            MOI.Utilities.DoubleDicts.DoubleDict{MOI.ConstraintIndex}(),
            0,
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarAffineTerm{Float64}},
            }(),
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarQuadraticTerm{Float64}},
            }(),
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarQuadraticTerm{Float64}},
            }(),
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarAffineTerm{Float64}},
            }(),
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.ScalarAffineTerm{Float64}},
            }(),
            Dict{MOI.ConstraintIndex,MOI.ConstraintIndex}(),
            0,
            MOI.Utilities.DoubleDicts.DoubleDict{
                Vector{MOI.VectorAffineTerm{Float64}},
            }(),
            MOI.VariableIndex(-1),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Vector{MOI.ScalarQuadraticTerm{Float64}}(),
            Vector{MOI.ScalarQuadraticTerm{Float64}}(),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Vector{MOI.ScalarAffineTerm{Float64}}(),
            Set{Int64}(),
            Vector{Float64}(),
            evaluate_duals,
            0,
            ONLY_CONSTRAINTS,
        )
    end
end

include("utils.jl")
include("duals.jl")
include("update_parameters.jl")

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           isempty(model.parameters) &&
           isempty(model.parameters_name) &&
           isempty(model.variables) &&
           isempty(model.updated_parameters) &&
           isempty(model.variables) &&
           model.last_variable_index_added == 0 &&
           model.last_parameter_index_added == PARAMETER_INDEX_THRESHOLD &&
           isempty(model.original_constraint_function_and_set_cache) &&
           isempty(model.affine_added_cache) &&
           model.last_affine_added == 0 &&
           isempty(model.affine_constraint_cache) &&
           isempty(model.quadratic_constraint_cache_pv) &&
           isempty(model.quadratic_constraint_cache_pp) &&
           isempty(model.quadratic_constraint_cache_pc) &&
           isempty(
               model.quadratic_constraint_variables_associated_to_parameters_cache,
           ) &&
           isempty(model.quadratic_added_cache) &&
           model.last_quad_add_added == 0 &&
           model.original_objective_function == MOI.VariableIndex(-1) &&
           isempty(model.affine_objective_cache) &&
           isempty(model.quadratic_objective_cache_pv) &&
           isempty(model.quadratic_objective_cache_pp) &&
           isempty(model.quadratic_objective_cache_pc) &&
           isempty(
               model.quadratic_objective_variables_associated_to_parameters_cache,
           ) &&
           isempty(model.dual_value_of_parameters) &&
           model.number_of_parameters_in_model == 0
end

function MOI.supports_constraint(
    model::Optimizer,
    F::Union{
        Type{MOI.VariableIndex},
        Type{MOI.ScalarAffineFunction{T}},
        Type{MOI.VectorOfVariables},
        Type{MOI.VectorAffineFunction{T}},
    },
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(model.optimizer, F, S)
end

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{T}},
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(
        model.optimizer,
        MOI.ScalarAffineFunction{T},
        S,
    )
end

function MOI.supports_constraint(
    model::Optimizer,
    ::Type{MOI.VectorQuadraticFunction{T}},
    S::Type{<:MOI.AbstractSet},
) where {T}
    return MOI.supports_constraint(
        model.optimizer,
        MOI.VectorAffineFunction{T},
        S,
    )
end

function MOI.supports(
    model::Optimizer,
    attr::Union{
        MOI.ObjectiveSense,
        MOI.ObjectiveFunction{MOI.VariableIndex},
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
    },
) where {T}
    return MOI.supports(model.optimizer, attr)
end

function MOI.supports(model::Optimizer, ::MOI.NLPBlock)
    return MOI.supports(model.optimizer, MOI.NLPBlock())
end

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    return MOI.set(model.optimizer, MOI.NLPBlock(), nlp_data)
end

function MOI.supports(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}},
) where {T}
    return MOI.supports(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
    )
end

function MOI.supports_incremental_interface(model::Optimizer)
    return MOI.supports_incremental_interface(model.optimizer)
end
function MOI.supports(model::Optimizer, ::MOI.Name)
    return MOI.supports(model.optimizer, MOI.Name())
end
function MOI.get(model::Optimizer, ::MOI.ListOfModelAttributesSet)
    return MOI.get(model.optimizer, MOI.ListOfModelAttributesSet())
end
MOI.get(model::Optimizer, ::MOI.Name) = MOI.get(model.optimizer, MOI.Name())
function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    return MOI.set(model.optimizer, MOI.Name(), name)
end
function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return MOI.get(model.optimizer, MOI.ListOfVariableIndices())
end
function MOI.get(model::Optimizer, ::MOI.ListOfVariableAttributesSet)
    return MOI.get(model.optimizer, MOI.ListOfVariableAttributesSet())
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintAttributesSet{F,S},
) where {F,S}
    if F === MOI.ScalarQuadraticFunction
        error(
            "MOI.ListOfConstraintAttributesSet is not implemented for ScalarQuadraticFunction.",
        )
    end
    return MOI.get(model.optimizer, MOI.ListOfConstraintAttributesSet{F,S}())
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F,S},
    f::F,
) where {F,S}
    MOI.set(model.optimizer, MOI.ConstraintFunction(), c, f)
    return
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F,S},
    s::S,
) where {F,S}
    MOI.set(model.optimizer, MOI.ConstraintSet(), c, s)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
    chg::MOI.ScalarCoefficientChange{Float64},
) where {F,S}
    MOI.modify(model.optimizer, c, chg)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ObjectiveFunction{F},
    chg::MOI.ScalarCoefficientChange{Float64},
) where {F<:MathOptInterface.AbstractScalarFunction}
    MOI.modify(model.optimizer, c, chg)
    return
end
function MOI.modify(
    model::Optimizer,
    c::MOI.ObjectiveFunction{F},
    chg::MOI.ScalarConstantChange{Float64},
) where {F<:MathOptInterface.AbstractScalarFunction}
    MOI.modify(model.optimizer, c, chg)
    return
end

function MOI.empty!(model::Optimizer{T}) where {T}
    MOI.empty!(model.optimizer)
    empty!(model.parameters)
    empty!(model.parameters_name)
    empty!(model.updated_parameters)
    empty!(model.variables)
    model.last_variable_index_added = 0
    model.last_parameter_index_added = PARAMETER_INDEX_THRESHOLD
    empty!(model.original_constraint_function_and_set_cache)
    empty!(model.affine_added_cache)
    model.last_affine_added = 0
    empty!(model.affine_constraint_cache)
    empty!(model.quadratic_constraint_cache_pv)
    empty!(model.quadratic_constraint_cache_pp)
    empty!(model.quadratic_constraint_cache_pc)
    empty!(model.quadratic_constraint_variables_associated_to_parameters_cache)
    empty!(model.quadratic_added_cache)
    model.last_quad_add_added = 0
    model.original_objective_function = MOI.VariableIndex(-1)
    empty!(model.vector_constraint_cache)
    empty!(model.affine_objective_cache)
    empty!(model.quadratic_objective_cache_pv)
    empty!(model.quadratic_objective_cache_pp)
    empty!(model.quadratic_objective_cache_pc)
    empty!(model.quadratic_objective_variables_associated_to_parameters_cache)
    empty!(model.dual_value_of_parameters)
    model.number_of_parameters_in_model = 0
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.VariableName,
    v::MOI.VariableIndex,
    name::String,
)
    if is_parameter_in_model(model, v)
        model.parameters_name[v] = name
    else
        MOI.set(model.optimizer, attr, v, name)
    end
    return
end

function MOI.get(model::Optimizer, attr::MOI.VariableName, v::MOI.VariableIndex)
    if is_parameter_in_model(model, v)
        return get(model.parameters_name, v, "")
    else
        return MOI.get(model.optimizer, attr, v)
    end
end

function MOI.get(model::Optimizer, tp::Type{MOI.VariableIndex}, attr::String)
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveBound)
    return MOI.get(model.optimizer, attr)
end

function MOI.supports(
    model::Optimizer,
    attr::MOI.VariableName,
    tp::Type{MOI.VariableIndex},
)
    return MOI.supports(model.optimizer, attr, tp)
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ConstraintName,
    c::MOI.ConstraintIndex,
    name::String,
)
    MOI.set(model.optimizer, attr, c, name)
    return
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintName,
    c::MOI.ConstraintIndex,
)
    return MOI.get(model.optimizer, attr, c)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintName,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S},
) where {T,S}
    moi_ci = get(model.affine_added_cache, c, c)
    # This SAF constraint was transformed into variable bound
    if typeof(moi_ci) === MOI.ConstraintIndex{MOI.VariableIndex,S}
        v = MOI.get(model.optimizer, MOI.ConstraintFunction(), moi_ci)
        variable_name = MOI.get(model.optimizer, MOI.VariableName(), v)
        return "ParametricBound_$(S)_$(variable_name)"
    else
        return MOI.get(model.optimizer, attr, moi_ci)
    end
end

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return length(model.parameters) + length(model.variables)
end

function MOI.supports(
    model::Optimizer,
    attr::MOI.ConstraintName,
    tp::Type{<:MOI.ConstraintIndex},
)
    return MOI.supports(model.optimizer, attr, tp)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintFunction,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if haskey(model.original_constraint_function_and_set_cache, ci)
        return model.original_constraint_function_and_set_cache[ci][1]
    else
        return MOI.get(model.optimizer, attr, ci)
    end
end

function MOI.get(
    model::Optimizer,
    tp::Type{MOI.ConstraintIndex{F,S}},
    attr::String,
) where {F,S}
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.get(model::Optimizer, tp::Type{MOI.ConstraintIndex}, attr::String)
    return MOI.get(model.optimizer, tp, attr)
end

function MOI.is_valid(model::Optimizer, vi::MOI.VariableIndex)
    return MOI.is_valid(model.optimizer, vi)
end

function MOI.supports(model::Optimizer, ::MOI.NumberOfThreads)
    return MOI.supports(model.optimizer, MOI.NumberOfThreads())
end

function MOI.supports(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.supports(model.optimizer, MOI.TimeLimitSec())
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

function MOI.get(model::Optimizer, ::MOI.SolveTimeSec)
    return MOI.get(model.optimizer, MOI.SolveTimeSec())
end

function MOI.supports(model::Optimizer, ::MOI.Silent)
    return MOI.supports(model.optimizer, MOI.Silent())
end

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    MOI.set(model.optimizer, MOI.Silent(), value)
    return
end

MOI.get(model::Optimizer, ::MOI.Silent) = MOI.get(model.optimizer, MOI.Silent())

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return MOI.get(model.optimizer, MOI.RawStatusString())
end

function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F,S}
    return length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}()))
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if haskey(model.original_constraint_function_and_set_cache, ci)
        return model.original_constraint_function_and_set_cache[ci][2]
    else
        MOI.throw_if_not_valid(model, ci)
        return MOI.get(model.optimizer, attr, ci)
    end
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveSense)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveFunctionType)
    return typeof(model.original_objective_function)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ObjectiveFunction{F},
) where {
    F<:Union{
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{T},
        MOI.ScalarQuadraticFunction{T},
    },
} where {T}
    if !function_has_parameters(model, model.original_objective_function)
        return MOI.get(model.optimizer, attr)
    else
        if F === typeof(model.original_objective_function)
            return model.original_objective_function
        else
            throw(InexactError)
        end
    end
end

function MOI.get(model::Optimizer, attr::MOI.ResultCount)
    return MOI.get(model.optimizer, attr)
end

# In the AbstractBridgeOptimizer, we collect all the possible constraint types and them filter with NumberOfConstraints.
# If NumberOfConstraints is zero then we remove it from the list.
# Here, you can look over keys(quadratic_added_cache) and add the F-S types of all the keys in constraints.
# To implement NumberOfConstraints, you call NumberOfConstraints to the inner optimizer.
# Then you remove the number of constraints of that that in values(quadratic_added_cache)
function MOI.get(model::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    inner_ctrs = MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())
    if !has_quadratic_constraint_caches(model)
        return inner_ctrs
    end

    cache_keys = collect(keys(model.quadratic_added_cache))
    constraints = Set{Tuple{DataType,DataType}}()

    for (F, S) in inner_ctrs
        inner_index =
            MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
        cache_map_check =
            quadratic_constraint_cache_map_check.(model, inner_index)
        for type in typeof.(cache_keys[cache_map_check])
            push!(constraints, (type.parameters[1], type.parameters[2]))
        end
        # If not all the constraints are chached then also push the original type
        if !all(cache_map_check)
            push!(constraints, (F, S))
        end
    end

    return collect(constraints)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ListOfConstraintIndices{F,S},
) where {S,F<:Union{MOI.VectorOfVariables,MOI.VariableIndex}}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {
    S<:MOI.AbstractSet,
    F<:Union{MOI.ScalarAffineFunction{T},MOI.VectorAffineFunction{T}},
} where {T}
    inner_index = MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
    if !has_quadratic_constraint_caches(model)
        return inner_index
    end

    cache_map_check = quadratic_constraint_cache_map_check(mode, inner_index)
    return inner_index[cache_map_check]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {S<:MOI.AbstractSet,F<:MOI.ScalarQuadraticFunction{T}} where {T}
    inner_index = MOI.ConstraintIndex{F,S}[]
    if MOI.supports_constraint(model.optimizer, F, S)
        inner_index =
            MOI.get(model.optimizer, MOI.ListOfConstraintIndices{F,S}())
        if !has_quadratic_constraint_caches(model)
            return inner_index
        end
    end

    quadratic_caches = [
        :quadratic_constraint_cache_pc,
        :quadratic_constraint_cache_pp,
        :quadratic_constraint_cache_pv,
        # JD: Check if this applies here
        # :quadratic_constraint_variables_associated_to_parameters_cache
    ]

    for field in quadratic_caches
        cache = getfield(model, field)
        push!(inner_index, keys(cache)...)
    end
    return inner_index
end

function MOI.supports_add_constrained_variable(::Optimizer, ::Type{Parameter})
    return true
end

function MOI.supports_add_constrained_variables(
    model::Optimizer,
    ::Type{MOI.Reals},
)
    return MOI.supports_add_constrained_variables(model.optimizer, MOI.Reals)
end

function MOI.add_variable(model::Optimizer)
    next_variable_index!(model)
    return MOI.Utilities.CleverDicts.add_item(
        model.variables,
        MOI.add_variable(model.optimizer),
    )
end

function MOI.add_constrained_variable(model::Optimizer, set::Parameter)
    next_parameter_index!(model)
    p = MOI.VariableIndex(model.last_parameter_index_added)
    MOI.Utilities.CleverDicts.add_item(model.parameters, set.val)
    cp = MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(
        model.last_parameter_index_added,
    )
    MOI.Utilities.CleverDicts.add_item(model.updated_parameters, NaN)
    update_number_of_parameters!(model)
    return p, cp
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    set::MOI.AbstractScalarSet,
)
    if is_parameter_in_model(model, f)
        error("Cannot constrain a parameter")
    elseif !is_variable_in_model(model, f)
        error("Variable not in the model")
    end
    return MOI.add_constraint(model.optimizer, f, set)
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    set::S,
) where {T,S}
    vars, params, param_constant =
        separate_possible_terms_and_calculate_parameter_constant(model, f.terms)
    model.last_affine_added += 1
    if model.constraints_interpretation == ONLY_BOUNDS
        if (length(vars) == 1) && isone(MOI.coefficient(vars[1]))
            poi_ci =
                add_vi_constraint(model, vars, params, param_constant, f, set)
        else
            error(
                "It was not possible to interpret this constraint as a variable bound.",
            )
        end
    elseif model.constraints_interpretation == ONLY_CONSTRAINTS
        poi_ci = add_saf_constraint(model, vars, params, param_constant, f, set)
    elseif model.constraints_interpretation == BOUNDS_AND_CONSTRAINTS
        if (length(vars) == 1) && isone(MOI.coefficient(vars[1]))
            poi_ci =
                add_vi_constraint(model, vars, params, param_constant, f, set)
        else
            poi_ci =
                add_saf_constraint(model, vars, params, param_constant, f, set)
        end
    end
    model.original_constraint_function_and_set_cache[poi_ci] = (f, set)
    return poi_ci
end

function add_saf_constraint(
    model::Optimizer,
    vars::Vector{MOI.ScalarAffineTerm{T}},
    params::Vector{MOI.ScalarAffineTerm{T}},
    param_constant::T,
    f::MOI.ScalarAffineFunction{T},
    set::S,
) where {T,S}
    moi_ci = MOI.Utilities.normalize_and_add_constraint(
        model.optimizer,
        MOI.ScalarAffineFunction(vars, f.constant + param_constant),
        set,
    )
    poi_ci = create_new_poi_ci_and_save_affine_caches(model, params, moi_ci)
    return poi_ci
end

function add_vi_constraint(
    model::Optimizer,
    vars::Vector{MOI.ScalarAffineTerm{T}},
    params::Vector{MOI.ScalarAffineTerm{T}},
    param_constant::T,
    f::MOI.ScalarAffineFunction{T},
    set::S,
) where {T,S}
    moi_ci = MOI.Utilities.normalize_and_add_constraint(
        model.optimizer,
        MOI.VariableIndex(vars[1].variable.value),
        update_constant!(set, f.constant + param_constant),
    )
    poi_ci = create_new_poi_ci_and_save_affine_caches(model, params, moi_ci)
    return poi_ci
end

function create_new_poi_ci_and_save_affine_caches(
    model::Optimizer,
    params::Vector{MOI.ScalarAffineTerm{T}},
    moi_ci::MOI.ConstraintIndex{F,S},
) where {T,F,S}
    poi_ci = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S}(
        model.last_affine_added,
    )
    model.affine_constraint_cache[poi_ci] = params
    model.affine_added_cache[poi_ci] = moi_ci
    return poi_ci
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
    set::MOI.AbstractScalarSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    v::MOI.VariableIndex,
)
    if is_parameter_in_model(model, v)
        return model.parameters[p_idx(v)]
    elseif is_variable_in_model(model, v)
        return MOI.get(model.optimizer, attr, model.variables[v])
    else
        error("Variable not in the model")
    end
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariableBasisStatus,
    v::MOI.VariableIndex,
)
    if is_variable_in_model(model, v)
        return MOI.get(model.optimizer, attr, model.variables[v])
    else
        error("VariableBasisStatus is not supported for parameters.")
    end
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    cp::MOI.ConstraintIndex{MOI.VariableIndex,Parameter},
    set::Parameter,
)
    p = MOI.VariableIndex(cp.value)
    if !is_parameter_in_model(model, p)
        error("Parameter not in the model")
    end
    return model.updated_parameters[p_idx(p)] = set.val
end

"""
    ParameterValue <: MOI.AbstractVariableAttribute

Attribute defined to set and get parameter values

# Example

```julia
MOI.set(model, POI.ParameterValue(), p, 2.0)
MOI.get(model, POI.ParameterValue(), p)
```
"""
struct ParameterValue <: MOI.AbstractVariableAttribute end

# We need a CachingOptimizer fallback to
# get ParameterValue working correctly on JuMP
# TODO: Think of a better solution for this

function MOI.set(
    opt::MOI.Utilities.CachingOptimizer,
    ::ParameterValue,
    var::MOI.VariableIndex,
    val::Float64,
)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(var.value)
    set = MOI.set(opt, MOI.ConstraintSet(), ci, Parameter(val))
    return nothing
end

function MOI.set(
    model::Optimizer,
    ::ParameterValue,
    var::MOI.VariableIndex,
    val::Float64,
)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(var.value)
    set = MOI.set(model, MOI.ConstraintSet(), ci, Parameter(val))
    return nothing
end

function MOI.set(
    opt::MOI.Utilities.CachingOptimizer,
    ::ParameterValue,
    vi::MOI.VariableIndex,
    val::Real,
)
    return MOI.set(opt, ParameterValue(), vi, convert(Float64, val))
end

function MOI.set(
    model::Optimizer,
    ::ParameterValue,
    vi::MOI.VariableIndex,
    val::Real,
)
    return MOI.set(model, ParameterValue(), vi, convert(Float64, val))
end

function MOI.get(
    opt::MOI.Utilities.CachingOptimizer,
    ::ParameterValue,
    var::MOI.VariableIndex,
)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(var.value)
    set = MOI.get(opt, MOI.ConstraintSet(), ci)
    return set.val
end

function MOI.get(model::Optimizer, ::ParameterValue, var::MOI.VariableIndex)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,Parameter}(var.value)
    set = MOI.get(model, MOI.ConstraintSet(), ci)
    return set.val
end

"""
    ConstraintsInterpretation <: MOI.AbstractOptimizerAttribute

Attribute to define how [`POI.Optimizer`](@ref) should interpret constraints.

- `POI.ONLY_CONSTRAINTS`: Only interpret `ScalarAffineFunction` constraints as linear constraints
  If an expression such as `x >= p1 + p2` appears it will be trated like a new constraint.
  **This is the default behaviour of [`POI.Optimizer`](@ref)**

- `POI.ONLY_BOUNDS`: Only interpret `ScalarAffineFunction` constraints as a variable bound.
  This is valid for constraints such as `x >= p` or `x >= p1 + p2`. If a constraint `x1 + x2 >= p` appears,
  which is not a valid variable bound it will throw an error.

- `POI.BOUNDS_AND_CONSTRAINTS`: Interpret `ScalarAffineFunction` constraints as a variable bound if they
  are a valid variable bound, i.e., `x >= p` or `x >= p1 + p2` and interpret them as linear constraints
  otherwise.

# Example

```julia
MOI.set(model, POI.InterpretConstraintsAsBounds(), POI.ONLY_BOUNDS)
MOI.set(model, POI.InterpretConstraintsAsBounds(), POI.ONLY_CONSTRAINTS)
MOI.set(model, POI.InterpretConstraintsAsBounds(), POI.BOUNDS_AND_CONSTRAINTS)
```
"""
struct ConstraintsInterpretation <: MOI.AbstractOptimizerAttribute end

function MOI.set(
    model::Optimizer,
    ::ConstraintsInterpretation,
    value::ConstraintsInterpretationCode,
)
    return model.constraints_interpretation = value
end

function empty_objective_function_caches!(model::Optimizer)
    empty!(model.affine_objective_cache)
    empty!(model.quadratic_objective_cache_pv)
    empty!(model.quadratic_objective_cache_pp)
    empty!(model.quadratic_objective_cache_pc)
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    f::MOI.ScalarAffineFunction{T},
) where {T}
    # clear previously defined objetive function cache
    empty_objective_function_caches!(model)
    if !function_has_parameters(model, f)
        MOI.set(model.optimizer, attr, f)
    else
        vars, params, param_constant =
            separate_possible_terms_and_calculate_parameter_constant(
                model,
                f.terms,
            )
        MOI.set(
            model.optimizer,
            attr,
            MOI.ScalarAffineFunction(vars, f.constant + param_constant),
        )
        model.affine_objective_cache = params
    end
    model.original_objective_function = f
    return
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    v::MOI.VariableIndex,
)
    if is_parameter_in_model(model, v)
        error("Cannot use a parameter as objective function alone")
    elseif !is_variable_in_model(model, v)
        error("Variable not in the model")
    end
    model.original_objective_function = v
    return MOI.set(model.optimizer, attr, model.variables[v])
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    MOI.set(model.optimizer, attr, sense)
    return
end

function MOI.get(
    model::Optimizer,
    attr::T,
) where {
    T<:Union{
        MOI.TerminationStatus,
        MOI.ObjectiveValue,
        MOI.DualObjectiveValue,
        MOI.PrimalStatus,
        MOI.DualStatus,
    },
}
    return MOI.get(model.optimizer, attr)
end

function MOI.get(
    model::Optimizer,
    attr::T,
    c::MOI.ConstraintIndex,
) where {
    T<:Union{MOI.ConstraintPrimal,MOI.ConstraintDual,MOI.ConstraintBasisStatus},
}
    return MOI.get(model.optimizer, attr, c)
end

function MOI.get(
    model::Optimizer,
    attr::AT,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S},
) where {
    AT<:Union{
        MOI.ConstraintPrimal,
        MOI.ConstraintDual,
        MOI.ConstraintBasisStatus,
    },
    T,
    S<:MOI.AbstractScalarSet,
}
    moi_ci = get(model.affine_added_cache, c, c)
    return MOI.get(model.optimizer, attr, moi_ci)
end

function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, val::Any)
    MOI.set(model.optimizer, attr, val)
    return
end

function MOI.get(model::Optimizer, ::MOI.SolverName)
    name = MOI.get(model.optimizer, MOI.SolverName())
    return "Parametric Optimizer with $(name) attached"
end

function MOI.get(model::Optimizer, ::MOI.SolverVersion)
    return MOI.get(model.optimizer, MOI.SolverVersion())
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VectorOfVariables,
    set::MOI.AbstractVectorSet,
) where {T}
    if function_has_parameters(model, f)
        error("VectorOfVariables does not allow parameters")
    end
    return MOI.add_constraint(model.optimizer, f, set)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    set::MOI.AbstractVectorSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
    set::MOI.AbstractVectorSet,
) where {T}
    vars, params, param_constants =
        separate_possible_terms_and_calculate_parameter_constant(model, f, set)
    ci = MOI.add_constraint(
        model.optimizer,
        MOI.VectorAffineFunction(vars, f.constants + param_constants),
        set,
    )
    model.vector_constraint_cache[ci] = params
    model.original_constraint_function_and_set_cache[ci] = (f, set)
    return ci
end

function add_constraint_with_parameters_on_function(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
    set::S,
) where {T,S<:MOI.AbstractScalarSet}
    (
        quad_vars,
        quad_aff_vars,
        quad_params,
        aff_terms,
        variables_associated_to_parameters,
        quad_param_constant,
    ) = separate_possible_terms_and_calculate_parameter_constant(
        model,
        f.quadratic_terms,
    )

    (
        aff_vars,
        aff_params,
        terms_with_variables_associated_to_parameters,
        aff_param_constant,
    ) = separate_possible_terms_and_calculate_parameter_constant(
        model,
        f.affine_terms,
        variables_associated_to_parameters,
    )

    aff_terms = vcat(aff_terms, aff_vars)
    const_term = f.constant + aff_param_constant + quad_param_constant
    quad_terms = quad_vars
    f_quad = if !isempty(quad_vars)
        MOI.ScalarQuadraticFunction(quad_terms, aff_terms, const_term)
    else
        MOI.ScalarAffineFunction(aff_terms, const_term)
    end
    model.last_quad_add_added += 1
    ci =
        MOI.Utilities.normalize_and_add_constraint(model.optimizer, f_quad, set)
    # This part is used to remember that ci came from a quadratic function
    # It is particularly useful because sometimes the constraint mutates
    new_ci = MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{T},S}(
        model.last_quad_add_added,
    )
    model.quadratic_added_cache[new_ci] = ci
    model.original_constraint_function_and_set_cache[new_ci] = (f, set)
    fill_quadratic_constraint_caches!(
        model,
        new_ci,
        quad_aff_vars,
        quad_params,
        aff_params,
        terms_with_variables_associated_to_parameters,
    )
    return new_ci
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
    set::MOI.AbstractScalarSet,
) where {T}
    if !function_has_parameters(model, f)
        return MOI.add_constraint(model.optimizer, f, set)
    else
        return add_constraint_with_parameters_on_function(model, f, set)
    end
end

function MOI.delete(model::Optimizer, v::MOI.VariableIndex)
    delete!(model.variables, v)
    MOI.delete(model.optimizer, v)
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {F<:MOI.ScalarAffineFunction,S<:MOI.AbstractSet}
    if haskey(model.affine_constraint_cache, c)
        delete!(model.affine_constraint_cache, c)
    end
    MOI.delete(model.optimizer, c)
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{MOI.VariableIndex,MOI.VectorOfVariables,MOI.VectorAffineFunction},
    S<:MOI.AbstractSet,
}
    MOI.delete(model.optimizer, c)
    return
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{MOI.VariableIndex,MOI.VectorOfVariables,MOI.VectorAffineFunction},
    S<:MOI.AbstractSet,
}
    return MOI.is_valid(model.optimizer, c)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{F,S},
) where {F<:MOI.ScalarAffineFunction,S<:MOI.AbstractSet}
    return MOI.is_valid(model.optimizer, c)
end

function MOI.set(
    model::Optimizer,
    attr::MOI.ObjectiveFunction,
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    # clear previously defined objetive function cache
    empty_objective_function_caches!(model)
    if !function_has_parameters(model, f)
        MOI.set(model.optimizer, attr, f)
        return
    end
    (
        quad_vars,
        quad_aff_vars,
        quad_params,
        aff_terms,
        variables_associated_to_parameters,
        quad_param_constant,
    ) = separate_possible_terms_and_calculate_parameter_constant(
        model,
        f.quadratic_terms,
    )

    (
        aff_vars,
        aff_params,
        terms_with_variables_associated_to_parameters,
        aff_param_constant,
    ) = separate_possible_terms_and_calculate_parameter_constant(
        model,
        f.affine_terms,
        variables_associated_to_parameters,
    )

    aff_terms = vcat(aff_terms, aff_vars)
    const_term = f.constant + aff_param_constant + quad_param_constant

    if !isempty(quad_vars)
        MOI.set(
            model.optimizer,
            attr,
            MOI.ScalarQuadraticFunction(quad_vars, aff_terms, const_term),
        )
    else
        MOI.set(
            model.optimizer,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
            MOI.ScalarAffineFunction(aff_terms, const_term),
        )
    end

    model.original_objective_function = f
    model.quadratic_objective_cache_pv = quad_aff_vars
    model.quadratic_objective_cache_pp = quad_params
    model.quadratic_objective_cache_pc = aff_params
    model.quadratic_objective_variables_associated_to_parameters_cache =
        terms_with_variables_associated_to_parameters

    return
end

function _poi_default_copy_to(dest::T, src::MOI.ModelLike) where {T}
    if !MOI.supports_incremental_interface(dest)
        error("Model $(typeof(dest)) does not support copy_to.")
    end
    MOI.empty!(dest)
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    index_map = MOI.IndexMap()
    # The `NLPBlock` assumes that the order of variables does not change (#849)
    # Therefore, all VariableIndex and VectorOfVariable constraints are added
    # seprately, and no variables constrained-on-creation are added.

    # This is not valid for NLPs with Parameters, they should enter
    has_nlp = MOI.NLPBlock() in MOI.get(src, MOI.ListOfModelAttributesSet())
    constraints_not_added = if has_nlp
        vcat(
            Any[
                MOI.get(src, MOI.ListOfConstraintIndices{F,S}()) for
                (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent()) if
                MOI.Utilities._is_variable_function(F) && S != Parameter
            ],
            Any[MOI.Utilities._try_constrain_variables_on_creation(
                dest,
                src,
                index_map,
                Parameter,
            )],
        )
    else
        Any[
            MOI.Utilities._try_constrain_variables_on_creation(
                dest,
                src,
                index_map,
                S,
            ) for S in MOI.Utilities.sorted_variable_sets_by_cost(dest, src)
        ]
    end
    MOI.Utilities._copy_free_variables(dest, index_map, vis_src)
    # Copy variable attributes
    MOI.Utilities.pass_attributes(dest, src, index_map, vis_src)
    # Copy model attributes
    MOI.Utilities.pass_attributes(dest, src, index_map)
    # Copy constraints
    MOI.Utilities._pass_constraints(dest, src, index_map, constraints_not_added)
    MOI.Utilities.final_touch(dest, index_map)
    return index_map
end

function MOI.Utilities.default_copy_to(
    dest::MOI.Bridges.LazyBridgeOptimizer{Optimizer{Float64,T}},
    src::MOI.ModelLike,
) where {T}
    return _poi_default_copy_to(dest, src)
end

function MOI.Utilities.default_copy_to(
    dest::Optimizer{Float64,T},
    src::MOI.ModelLike,
) where {T}
    return _poi_default_copy_to(dest, src)
end

function MOI.optimize!(model::Optimizer)
    if !isempty(model.updated_parameters)
        update_parameters!(model)
    end
    MOI.optimize!(model.optimizer)
    if model.evaluate_duals
        calculate_dual_of_parameters(model)
    end
    return
end

end # module


text:
#=
Illustrates preimages of functions

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 07/09/2014
=#
using PyPlot

f(x) = 0.6 .* cos(4.0 .* x) .+ 1.3

xmin, xmax = -1.0, 1.0
x = linspace(xmin, xmax, 160)
y = f(x)
ya, yb = minimum(y), maximum(y)

fig, axes = subplots(2, 1, figsize=(8, 8))

for ax in axes
    for spine in ["left", "bottom"]
        ax[:spines][spine][:set_position]("zero")
    end

    for spine in ["right", "top"]
        ax[:spines][spine][:set_color]("none")
    end

    ax[:set_xlim](xmin, xmax)
    ax[:set_ylim](-0.6, 3.2)
    ax[:set_xticks]([])
    ax[:set_yticks]([])

    ax[:plot](x, y, "k-", lw=2, label=L"$f$")
    ax[:fill_between](x, ya, yb, facecolor="blue", alpha=0.05)
    ax[:vlines]([0], ya, yb, lw=3, color="blue", label=L"range of $f$")
    ax[:text](0.04, -0.3, L"$0$", fontsize=16)
end

ax = axes[1]
ax[:legend](loc="upper right", frameon=false)
ybar = 1.5
ax[:plot](x, x .* 0 .+ ybar, "k--", alpha=0.5)
ax[:text](0.05, 0.8 * ybar, L"$y$", fontsize=16)

for (i, z) in enumerate([-0.35, 0.35])
    ax[:vlines](z, 0, f(z), linestyle="--", alpha=0.5)
    ax[:text](z, -0.2, LaTeXString("\$x_$i\$"), fontsize=16)
end

ax = axes[2]

ybar = 2.6
ax[:plot](x, x * 0 + ybar, "k--", alpha=0.5)
ax[:text](0.04, 0.91 * ybar, L"$y$", fontsize=16)


text:
function initial_state_ctqw(::Type{T}, size::Int) where T<:Number
   fill(T(1. /sqrt(size)), size)
end

function initial_state(qws::QWSearch{<:AbstractCTQW})
   hamiltonian = qws.parameters[:hamiltonian]
   initial_state_ctqw(eltype(hamiltonian), size(hamiltonian,1))
end

function evolve(qwd::QWDynamics{<:AbstractCTQW},
                state::Vector{<:Number},
                runtime::Real)
   hamiltonian_evolution(qwd.parameters[:hamiltonian], state, runtime)
end

function measure(::QWDynamics{<:AbstractCTQW}, state::Vector{<:Number})
   abs.(state).^2
end

function measure(::QWDynamics{<:AbstractCTQW}, state::Vector{<:Number}, vertices::Vector{Int})
   abs.(state[vertices]).^2
end


text:
"""
    surrogate_model(plan::AbstractArray{T,2}, samples::AbstractArray{T,2}; options=Options()) where T

Returns a surrogate model function based on an optimized Radial Basis Function
interpolant. Depending on the options, the kernel, kernel width and scaling of
input data is optimized.

...
# Arguments
- `plan::AbstractArray{T,2}`:
    sample locations where each column corresponds to the location of one point.
    `size(plan) = (num_dimensions,num_samples)`.
- `samples::AbstractArray{T,2}`:
    function value at each sample location. each column contains one value from
    the corresponding plan location. `size(samples) = (1,num_samples)`.
- `options=Options()`:
    all options available to customize the surrogate optimization.
...
"""
function surrogate_model(plan::AbstractArray{T,2}, samples::AbstractArray{T,2}; options::Options=Options()) where T

    @unpack num_interpolants, trace, parallel_surrogate = options

    trace && println("Creating optimized surrogate model ...")

    (length(samples) != size(plan,2)) && error("plan and samples do not have the correct sizes")

    #Optimize RBF hypers for the ensamble of interpolants
    if parallel_surrogate
        optres = pmap(  (x)->rbf_hypers_opt(samples, plan, options),
                        1:num_interpolants)
    else
        optres = map(  (x)->rbf_hypers_opt(samples, plan, options),
                        1:num_interpolants)
    end

    # Sample order for ScatteredInterpolation
    samples_vec = vec(samples)

    # Interpolation object based on the optimization results
    itp = Array{ScatteredInterpolation.RBFInterpolant,1}(undef,num_interpolants)
    for i = 1:num_interpolants
        # Preprocess the points based on the settings used.
        preprocessed_point = preprocess_point(plan,optres[i],base_scale=plan)

        itp[i] = interpolate(optres[i].kernelFunc, preprocessed_point,
                            samples_vec, smooth=optres[i].smooth)
    end

    # Pre allocate and calculate to reduce cost of evaluating surrogate model.
    preprocessed_est_point = Array{Float64,2}(undef,size(plan,1),1)
    old_min = minimum(plan,dims=2)
    old_max = maximum(plan,dims=2)

    sm_func = (estimation_point) ->  surrogate_evaluate.(Ref(preprocessed_est_point),
                                        Ref(estimation_point),Tuple(itp),Tuple(optres),
                                        Ref(old_min),Ref(old_max))

    return sm_func, optres
end


text:
module WRDSMerger

##############################################################################
##
## Dependencies
##
##############################################################################
using DataFrames
using Dates
using BusinessDays
using Statistics
using LibPQ
using DBInterface
using AbstractTrees
using ShiftedArrays: lead
using InteractiveUtils
using StatsBase
using Statistics
using LinearAlgebra

##############################################################################
##
## Exported methods and types
##
##############################################################################

# identifiers and linking items
export link_identifiers, Permno, Cusip, NCusip,
    GVKey, CIK, Ticker, IbesTicker, LinkTable,
    link_table

# downloads and WRDS exploration functions
export comp_data, crsp_data, crsp_market, crsp_stocknames,
    crsp_adjust, crsp_delist, list_libraries, list_tables,
    describe_table, get_table, raw_sql, ff_data

# extra utilities
export range_join, BDay, Conditions

##############################################################################
##
## Load files
##
##############################################################################

include("utils.jl")
include(joinpath("utils", "dateFunctions.jl"))
include(joinpath("utils", "identifierTypes.jl"))
include(joinpath("utils", "linkTree.jl"))
include(joinpath("utils", "utils.jl"))

include("crspFunctions.jl")
include("compFunctions.jl")
include("mergeFunctions.jl")
include("exploreDB.jl")
include("ffData.jl")

global default_tables = Dict{String, String}(
    "comp_funda" => "compa.funda",
    "comp_fundq" => "compa.fundq",
    "crsp_stocknames" => "crsp.stocknames",
    "crsp_index" => "crsp.dsi",
    "crsp_stock_data" => "crsp.dsf",
    "crsp_delist" => "crsp.dsedelist",
    "crsp_a_ccm_ccmxpf_lnkhist" => "crsp_a_ccm.ccmxpf_lnkhist",
    "wrdsapps_ibcrsphist" => "wrdsapps.ibcrsphist",
    "comp_company" => "comp.company",
    "ff_factors" => "ff.factors_daily",
)

end


text:
function ij_pair(N)
  npairs = div(N*(N-1),2)
  ij = Matrix{Int}(undef,npairs,2)
  for ipair in 1:npairs
    ij[ipair,1], ij[ipair,2] = linear_index(ipair)
  end
  return ij
end
export ij_pair

# energy calculation using linked lists method
function utotal_parallel2(box, DATA, first_atom, next_atom, nc)

  nthreads = Threads.nthreads()        # number of threads available
  @unpack N, side, cutoff = DATA       # extraction some variables from DATA
  npairs = div(N * (N - 1), 2)         # Number of pairs
  Ut = zeros(nthreads)                 # vector with nthreads positions
  n_per_thread = div(N, nthreads)  # number of pairs per thread

  if mod(npairs,nthreads) != 0
    error(" mod(npairs,nthreads) != 0 ")   # return error if npairs % nthreads != 0
  end

  Threads.@threads for id in 1:nthreads
    ifirst = (id-1)*n_per_thread + 1
    ilast = ifirst - 1 + n_per_thread

    for iat in ifirst:ilast

  #    ip = ij[ipair,1]
  #    jp = ij[ipair,2]

      icell = trunc(Int64,box[iat][1]/cutoff) + 1
      jcell = trunc(Int64,box[iat][2]/cutoff) + 1

 #     println("atribuicao de particulas em celulas")

      for i in icell-1:icell+1
        for j in jcell-1:jcell+1

  #        println("loop pelas celulas vizinhas")

          iw, jw = wrapcell(nc, i, j)
          jat = first_atom[iw,jw]

   #       println("acessando listas ligadas")

    #      println(iat,"   ",jat)
          while jat > 0
            if jat > iat

              rij = pbcseparation(box[jat],box[iat],side)
     #         println("calculando distancia entre particulas")
              if rij <= 2.

      #          println("calculando")
                Ut[id] += upair(rij,DATA)
              end
            end
            jat = next_atom[jat]
          end
        end
      end
    end
  end
  return sum(Ut)
end


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

mutable struct NameAvailabilityInformation <: SwaggerModel
    nameAvailable::Any # spec type: Union{ Nothing, Bool } # spec name: nameAvailable
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    message::Any # spec type: Union{ Nothing, String } # spec name: message

    function NameAvailabilityInformation(;nameAvailable=nothing, reason=nothing, message=nothing)
        o = new()
        validate_property(NameAvailabilityInformation, Symbol("nameAvailable"), nameAvailable)
        setfield!(o, Symbol("nameAvailable"), nameAvailable)
        validate_property(NameAvailabilityInformation, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(NameAvailabilityInformation, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        o
    end
end # type NameAvailabilityInformation

const _property_map_NameAvailabilityInformation = Dict{Symbol,Symbol}(Symbol("nameAvailable")=>Symbol("nameAvailable"), Symbol("reason")=>Symbol("reason"), Symbol("message")=>Symbol("message"))
const _property_types_NameAvailabilityInformation = Dict{Symbol,String}(Symbol("nameAvailable")=>"Bool", Symbol("reason")=>"String", Symbol("message")=>"String")
Base.propertynames(::Type{ NameAvailabilityInformation }) = collect(keys(_property_map_NameAvailabilityInformation))
Swagger.property_type(::Type{ NameAvailabilityInformation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_NameAvailabilityInformation[name]))}
Swagger.field_name(::Type{ NameAvailabilityInformation }, property_name::Symbol) =  _property_map_NameAvailabilityInformation[property_name]

function check_required(o::NameAvailabilityInformation)
    true
end

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


text:
## Exercise 4-8
## Enter the code in this chapter in a notebook.

## 1. Draw a stack diagram that shows the state of the program while executing circle(🐢, radius).
## You can do the arithmetic by hand or add print statements to the code.
println("Ans 1: ")

println(" turtle --> Turtle")
println("     turtle --> Turtle")
println("     radius --> 100")
println("     circumference --> 628.318...")
println("     n --> 212")
println("     len --> 2")
println("         turtle --> Turtle")
println("         nsides --> 212")
println("         len --> 2")
println("         i --> 1:212")
println(" return")

## 2. The version of arc in Refactoring is not very accurate because the linear approximation of the circle
## is always outside the true circle. As a result, the turtle ends up a few pixels away from the correct destination.
## My solution shows a way to reduce the effect of this error. Read the code and see if it makes sense to you.
## If you draw a diagram, you might see how it works.
using ThinkJulia

function polyline(t, n, len, angle)
    for i in 1:n
        forward(t, len)
        turn(t, -angle)
    end
end

"""
arc(t, r, angle)

Draws an arc with the given radius and angle:

    t: turtle
    r: radius
    angle: angle subtended by the arc, in degrees
"""
function arc(t, r, angle)
    arc_len = 2 * π * r * abs(angle) / 360
    n = trunc(arc_len / 4) + 3
    step_len = arc_len / n
    step_angle = angle / n

    # making a slight left turn before starting reduces
    # the error caused by the linear approximation of the arc
    turn(t, -step_angle/2)
    polyline(t, n, step_len, step_angle)
    turn(t, step_angle/2)
end

println("Ans 2: ")

@svg begin
    turtle = Turtle()
    # forward(turtle, 30)  # for checking the position of the turtle
    arc(turtle, 100, 125)
end

##

println("End.")


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats, Statistics, HypothesisTests

d = 4
η = abs(3randn())
η₀ = 1

G = LKJ(d, η)
F = LKJ(d, η₀)

@testset "LKJ construction errors" begin
    @test_throws ArgumentError LKJ(-1, η)
    @test_throws ArgumentError LKJ(d, -1)
    @test LKJ(-1, η, check_args = false) isa LKJ{typeof(η), typeof(-1)}
end

@testset "LKJ params" begin
    η̃ = params(G)
    η̃₀ = params(F)
    @test η̃ == η
    @test η̃₀ == η₀
end

@testset "LKJ dim" begin
    @test dim(G) == d
    @test dim(F) == d
end

@testset "LKJ size" begin
    @test size(G) == (d, d)
    @test size(F) == (d, d)
end

@testset "LKJ rank" begin
    @test rank(G) == d
    @test rank(F) == d
    @test rank(G) == rank(rand(G))
    @test rank(F) == rank(rand(F))
end

@testset "LKJ insupport" begin
    @test insupport(G, rand(G))
    @test insupport(F, rand(F))

    @test !insupport(G, rand(G) + rand(G) * im)
    @test !insupport(G, randn(d, d + 1))
    @test !insupport(G, randn(d, d))
end

@testset "LKJ mode" begin
    @test mode(LKJ(5, 1.5)) == mean(LKJ(5, 1.5))
    @test_throws ArgumentError mode( LKJ(5, 0.5) )
end

@testset "LKJ sample moments" begin
    @test isapprox(mean(rand(G, 100000)), mean(G) , atol = 0.1)
    @test isapprox(var(rand(G, 100000)), var(G) , atol = 0.1)
end

@testset "LKJ marginals" begin
    M = 10000
    α = 0.05
    L = sum(1:(d - 1))
    ρ = Distributions._marginal(G)
    mymats = zeros(d, d, M)
    for m in 1:M
        mymats[:, :, m] = rand(G)
    end
    for i in 1:d
        for j in 1:i-1
            kstest = ExactOneSampleKSTest(mymats[i, j, :], ρ)
            @test pvalue(kstest) >= α / L  #  multiple comparisons
        end
    end
end

@testset "LKJ conversion" for elty in (Float32, Float64, BigFloat)
    Gel1 = convert(LKJ{elty}, G)
    Gel2 = convert(LKJ{elty}, G.d, G.η, G.logc0)

    @test Gel1 isa LKJ{elty, typeof(d)}
    @test Gel2 isa LKJ{elty, typeof(d)}
    @test partype(Gel1) == elty
    @test partype(Gel2) == elty
end

@testset "check d = 1 edge case" begin
    lkj = LKJ(1, 2abs(randn()))
    @test var(lkj) == zeros(1, 1)
    @test rand(lkj) == ones(1, 1)
end

@testset "LKJ integrating constant" begin
    #  =============
    #  odd non-uniform
    #  =============
    d = 5
    η = 2.3
    lkj = LKJ(d, η)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
    @test lkj.logc0 == -Distributions.lkj_onion_loginvconst(d, η)
    #  =============
    #  odd uniform
    #  =============
    d = 5
    η = 1.0
    lkj = LKJ(d, η)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst_uniform_odd(d)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_vine_loginvconst_uniform(d)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.corr_logvolume(d)
    @test lkj.logc0 == -Distributions.lkj_onion_loginvconst_uniform_odd(d)
    #  =============
    #  even non-uniform
    #  =============
    d = 6
    η = 2.3
    lkj = LKJ(d, η)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
    @test lkj.logc0 == -Distributions.lkj_onion_loginvconst(d, η)
    #  =============
    #  even uniform
    #  =============
    d = 6
    η = 1.0
    lkj = LKJ(d, η)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_onion_loginvconst_uniform_even(d)
    @test Distributions.lkj_vine_loginvconst(d, η) ≈ Distributions.lkj_vine_loginvconst_uniform(d)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.lkj_loginvconst_alt(d, η)
    @test Distributions.lkj_onion_loginvconst(d, η) ≈ Distributions.corr_logvolume(d)
    @test lkj.logc0 == -Distributions.lkj_onion_loginvconst_uniform_even(d)
end

@testset "check integrating constant as a volume" begin
    #  d = 2: Lebesgue measure of the set of correlation matrices is 2.
    volume2D = 2
    @test volume2D ≈ exp( Distributions.lkj_onion_loginvconst(2, 1) )
    @test 1 / volume2D ≈ exp( LKJ(2, 1).logc0 )
    #  d = 3: Lebesgue measure of the set of correlation matrices is π²/2.
    #  See here: https://www.jstor.org/stable/2684832
    volume3D = 0.5π^2
    @test volume3D ≈ exp( Distributions.lkj_onion_loginvconst(3, 1) )
    @test 1 / volume3D ≈ exp( LKJ(3, 1).logc0 )
    #  d = 4: Lebesgue measure of the set of correlation matrices is (32/27)π².
    #  See here: https://doi.org/10.4169/amer.math.monthly.123.9.909
    volume4D = (32 / 27)*π^2
    @test volume4D ≈ exp( Distributions.lkj_onion_loginvconst(4, 1) )
    @test 1 / volume4D ≈ exp( LKJ(4, 1).logc0 )
end

@testset "check logpdf against archived Stan output" begin
    #  Compare to archived output from Stan's lkj_corr_lpdf function.
    #  https://mc-stan.org/docs/2_22/functions-reference/lkj-correlation.html
    #  https://mc-stan.org/math/db/d4f/lkj__corr__lpdf_8hpp_source.html
    R = [1 0.962395133838894 -0.436307195544856;
         0.962395133838894 1 -0.301102833786894;
        -0.436307195544856 -0.301102833786894 1]
    n = size(R, 1)
    @test isapprox(logpdf(LKJ(n, 0.5), R), -0.9874823, atol = 1e-6)
    @test isapprox(logpdf(LKJ(n, 1),   R), -1.596313, atol = 1e-6)
    @test isapprox(logpdf(LKJ(n, 3.4), R), -7.253798, atol = 1e-6)
end

@testset "importance sampling check" begin
    d = 3
    M = 20000
    f = LKJ(d, 2)
    g = LKJ(d, 1)
    h = mean(logdet.(rand(f, M)))
    ĥ = mean(logdet(R) * pdf(f, R) / pdf(g, R) for R in (rand(g) for i in 1:M))
    @test isapprox(h, ĥ, atol = 0.1)
end


text:
struct Name{name} end

"""
    name_str(name)

Create a typed [`Name`](@ref).

```jldoctest
julia> using LightQuery

julia> name"a"
name"a"
```
"""
macro name_str(name)
    esc(Name{Symbol(name)}())
end
export @name_str

"""
    Name(name)

Create a typed `Name`. Inverse of [`unname`](@ref). See also [`@name_str`](@ref).

```jldoctest name
julia> using LightQuery

julia> Name(:a)
name"a"
```

`Name`s can be used as selector functions.

```jldoctest name
julia> using Test: @inferred

julia> data = (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
(a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)

julia> @inferred name"c"(data)
1
```

Multiple names can be used as selector functions

```jldoctest name
julia> @inferred (name"c", name"f")(data)
(c = 1, f = 1.0)
```

A final use for names can be as a way to construct NamedTuples from pairs. Using `Name`s
instead of `Symbol`s gives type stability.

```jldoctest name
julia> @inferred NamedTuple((name"a", 1), (name"b", 2))
(a = 1, b = 2)
```
"""
@pure function Name(name)
    Name{name}()
end
export Name

"""
    unname(::Name{name}) where name

Inverse of [`Name`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred unname(Name(:a))
:a
```
"""
function unname(::Name{name}) where {name}
    name
end
export unname

function show(output::IO, ::Name{name}) where {name}
    print(output, "name\"", name, '"')
end

const MyNamedTuple = SomeOf{Tuple{Name,Any}}

function recursive_get(initial, name)
    throw(BoundsError(initial, (name,)))
end
function recursive_get(initial, name, (check_name, value), rest...)
    if name === check_name
        value
    else
        recursive_get(initial, name, rest...)
    end
end

function (::Name{name})(object) where {name}
    getproperty(object, name)
end
function (name::Name)(data::MyNamedTuple)
    # TODO: avoid recursion for large tuples
    recursive_get(data, name, data...)
end

function get_pair(data, name)
    name, name(data)
end

function (some_names::SomeOf{Name})(data::MyNamedTuple)
    partial_map(get_pair, data, some_names)
end
# TODO: a version which works on all types?
function (some_names::SomeOf{Name})(data::NamedTuple)
    NamedTuple(some_names(named_tuple(data))...)
end

function isless(::Name{name1}, ::Name{name2}) where {name1,name2}
    isless(name1, name2)
end

# to override recusion limit on constant propagation
@pure function to_Names(some_names::SomeOf{Symbol})
    my_map(Name, some_names)
end
function to_Names(them)
    my_map(Name, Tuple(them))
end

function NamedTuple(data::Tuple{Name,Any}...)
    NamedTuple{my_map((function ((key, value),)
        unname(key)
    end), data)}(my_map(value, data))
end

function named_tuple(data)
    partial_map(get_pair, data, to_Names(propertynames(data)))
end
# arghhh im a pirate
@pure function propertynames(data::NamedTuple{names, <:Any}) where {names}
    names
end
export MyNamedTuple

"""
    remove(data, old_names...)

Remove `old_names` from `data`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0);

julia> @inferred remove(data, name"c", name"f")
(a = 1, b = 1.0, d = 1.0, e = 1)
```
"""
function remove(data::MyNamedTuple, old_names...)
    my_setdiff(my_map(key, data), old_names)(data)
end
function remove(data, old_names...)
    NamedTuple(remove(named_tuple(data), old_names...)...)
end

export remove

"""
    transform(data, assignments...)

Merge `assignments` into `data`, overwriting old values.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0);

julia> @inferred transform(data, c = 2.0, f = 2, g = 1, h = 1.0)
(a = 1, b = 1.0, d = 1.0, e = 1, c = 2.0, f = 2, g = 1, h = 1.0)
```
"""
function transform(data::MyNamedTuple, assignments...)
    remove(data, my_map(key, assignments)...)..., assignments...
end
function transform(data; assignments...)
    NamedTuple(transform(named_tuple(data), named_tuple(assignments.data)...)...)
end
export transform

maybe_named_tuple(something) = something
maybe_named_tuple(something::MyNamedTuple) = NamedTuple(something...)

function rename_one(data, (new_name, old_name))
    new_name, maybe_named_tuple(old_name(data))
end
"""
    rename(data, new_name_old_names...)

Rename `data`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0);

julia> @inferred rename(data, c2 = name"c", f2 = name"f")
(a = 1, b = 1.0, d = 1.0, e = 1, c2 = 1, f2 = 1.0)
```
"""
function rename(data::MyNamedTuple, new_name_old_names...)
    remove(data, my_map(value, new_name_old_names)...)...,
    partial_map(rename_one, data, new_name_old_names)...
end
function rename(data; new_name_old_names...)
    NamedTuple(rename(named_tuple(data), named_tuple(new_name_old_names.data)...)...)
end
export rename

"""
    gather(data, new_name_old_names...)

For each `new_name, old_names` pair in `new_name_old_names`, gather the `old_names` into a single `new_name`. Inverse of [`spread`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0);

julia> @inferred gather(data, g = (name"b", name"e"), h = (name"c", name"f"))
(a = 1, d = 1.0, g = (b = 1.0, e = 1), h = (c = 1, f = 1.0))
```
"""
function gather(data::MyNamedTuple, new_name_old_names...)
    remove(data, my_flatten(my_map(value, new_name_old_names))...)...,
    partial_map(rename_one, data, new_name_old_names)...
end
function gather(data; new_name_old_names...)
    NamedTuple(gather(named_tuple(data), named_tuple(new_name_old_names.data)...)...)
end
export gather

"""
    spread(data, some_names...)

Unnest nested named tuples. Inverse of [`gather`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> gathered = (a = 1, d = 1.0, g = (b = 1.0, e = 1), h = (c = 1, f = 1.0));

julia> @inferred spread(gathered, name"g", name"h")
(a = 1, d = 1.0, b = 1.0, e = 1, c = 1, f = 1.0)
```
"""
function spread(data::MyNamedTuple, some_names...)
    remove(data, some_names...)...,
    my_flatten(my_map((function ((name, value),)
        named_tuple(value)
    end), some_names(data)))...
end
function spread(data, some_names...)
    NamedTuple(spread(named_tuple(data), some_names...)...)
end

export spread


text:
#using HDF5

import MPI

# Initialize MPI environment before loading QCMaterial
MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

include("uccgsd_hubbard.jl")

nsite = 4
ham = generate_ham(nsite)
n_electron = 4
cost_history, circuit, exact_gs_ene, opt = solve(ham, n_electron, comm=comm)

#if rank == 0
#    println("cost_history", cost_history)
#    println(opt["x"])
#
#    h5open("opt.h5", "w") do file
#        write(file, "theta_list", opt["x"])
#        write(file, "cost_history", cost_history)
#        write(file, "exact_gs_ene", exact_gs_ene)
#        write(file, "n_electron", n_electron)
#        write(file, "nsite", nsite)
#    end
#end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
using StaticArrays
using Revise
using LinearAlgebra
using CSV
using DataFrames

global G = 6.67e-11 # m^3/(kg*s^2)

global THETA = 0.7

global timesInSec = 0.00

global frameRate = 1/24

global M = 2e30 # 1 = 1 solar mass

global timeScalar = 60*60*60*24*365*5e6 #1 frame

global distanceRate = 3.3e16 # 1 = 1 parsec = 3.3 light year = 3.3e16 m

global strengthOfInteraction = 1

global rotatingPlane = [1.0,-1.0,-1.0]

global minimumDis = 5e33 # m

global epsilon = 3.0e45

mutable struct Body
    id::Int64
    mass::Float64
    position::SVector{3, Float64}
    acceleration::SVector{3, Float64}
    velocity::SVector{3, Float64}
end
Body() = Body(0, 0.0, SVector{3, Float64}(zeros(Float64,3)),
SVector{3, Float64}(zeros(Float64,3)), SVector{3, Float64}(zeros(Float64,3)))

mutable struct Node
    hasChildren::Bool
    children::Array{Node,1}
    numOfChild::Int64
    body::Body
    mass::Float64
    centerOfMass::SVector{3, Float64}
    minBounds::SVector{3, Float64}
    maxBounds::SVector{3, Float64}
end

Node() = Node(
false, Array{Node,1}(undef, 8),0,
Body(), 0.0, SVector{3, Float64}(zeros(3)),
SVector{3, Float64}(zeros(3)),SVector{3, Float64}(zeros(3))
)

#this function returns which octants will the position belong to in the area defined by minBounds and maxBounds
# this function will output an integer: 1--> top, left, back ; 1--> top, left, back; 2--> top, right, back
#   3--> top, left, front ; 4--> top, right, front; 5--> down, left, back; 6--> down, right, back ;
# 7--> down, left, front; 1--> down, right, front
function determineWhichoctants(position::SArray{Tuple{3},Float64,1,3},
                                minBounds::SArray{Tuple{3},Float64,1,3},
                                maxBounds::SArray{Tuple{3},Float64,1,3})

    x = position[1]
    y = position[2]
    z = position[3]
    midX = (maxBounds[1] - minBounds[1])/2 + minBounds[1]
    midY = (maxBounds[2] - minBounds[2])/2 + minBounds[2]
    midZ = (maxBounds[3] - minBounds[3])/2 + minBounds[3]
    maxX = maxBounds[1]
    maxY = maxBounds[2]
    maxZ = maxBounds[3]
    minX = minBounds[1]
    minY = minBounds[2]
    minZ = minBounds[3]

    output = 0
    if minX <= x < midX && minY <= y < midY && minZ <= z < midZ
        output =  5
    elseif minX <= x < midX && minY <= y < midY && maxZ >= z >= midZ
        output =  7
    elseif minX <= x < midX && maxY >= y >= midY && minZ <= z < midZ
        output =  1
    elseif minX <= x < midX &&  maxY >= y >= midY && maxZ>= z >= midZ
        output =  3
    elseif maxX >= x >= midX && minY <= y < midY && minZ <= z < midZ
        output =  6
    elseif maxX >= x >= midX && minY <= y < midY && maxZ >= z >= midZ
        output =  8
    elseif maxX >= x >= midX && maxY >= y >= midY && minZ <= z < midZ
        output =  2
    elseif maxX >= x >= midX && maxY >= y >= midY && maxZ >= z >= midZ
        output =  4
    end

    return output
end

# this function returns two lists: a list of minBounds for 8 octantss, and a list of maxBounds for 8 octantss
function createSubBounds(minBounds::SVector{3, Float64},maxBounds::SVector{3, Float64})
    # ::Array{SArray{Tuple{3},Float64,1,3},1},::Array{SArray{Tuple{3},Float64,1,3},1}
    minX = minBounds[1]
    maxX = maxBounds[1]
    midX = minX + (maxX-minX)/2
    minY = minBounds[2]
    maxY = maxBounds[2]
    midY = minY + (maxY-minY)/2
    minZ = minBounds[3]
    maxZ = maxBounds[3]
    midZ = minZ + (maxZ-minZ)/2
    quadMinBounds = [[minX,midY,minZ],
    [midX,midY,minZ],
    [minX,midY,midZ],
    [midX,midY,midZ],
    [minX,minY,minZ],
    [midX,minY,minZ],
    [minX,minY,midZ],
    [midX,minY,midZ]]
    quadMaxBounds = [[midX,maxY,midZ],
    [maxX,maxY,midZ],
    [midX,maxY,maxZ],
    [maxX,maxY,maxZ],
    [midX,midY,midZ],
    [maxX,midY,midZ],
    [midX,midY,maxZ],
    [maxX,midY,maxZ]]
    return quadMinBounds, quadMaxBounds
end

# this function divides a space into 8 octantss and makes a new node for each octants
function divideToOctants(node::Node)
    quadMinBounds, quadMaxBound = createSubBounds(node.minBounds, node.maxBounds)
    for i in 1:8
        newNode = Node(
        false, Array{Node,1}(undef, 8),0,
        Body(), 0.0, SVector{3, Float64}(zeros(Float64,3)),
        quadMinBounds[i],quadMaxBound[i]
        )
        node.children[i] = newNode
    end
end

# this function merges two bodies that is considered to be close to each other
function mergeBodies(b1::Body, b2::Body)
    totalMass = b1.mass + b2.mass
    b1.position = (b1.mass .* b1.position + b2.mass .* b2.position)/totalMass
    b1.velocity = (b1.mass .* b1.velocity + b2.mass .* b2.velocity)/totalMass
    b1.acceleration = (b1.mass .* b1.acceleration + b2.mass .* b2.acceleration)/totalMass
    b1.mass = b1.mass + b2.mass
    return b1

end
# this functions simulates adding a body into a node, subjecting to the quadtree structure
function insertBody(node::Node, body::Body)

    # which octants does the original body locates
    octantsOriginal = determineWhichoctants(node.centerOfMass, node.minBounds, node.maxBounds)
    # which octants does the new body locates
    octantsNew = determineWhichoctants(body.position, node.minBounds, node.maxBounds)
    # get the bounds for different Octants
    quadMinBounds, quadMaxBound = createSubBounds(node.minBounds, node.maxBounds)

    # update mass, centerOfMass
    totalMass = node.mass + body.mass

    node.centerOfMass = node.mass/totalMass.*node.centerOfMass + body.mass/totalMass.*body.position

    node.mass = totalMass
    # where to insert the body?
    if node.hasChildren && node.numOfChild == 1 # we need to make new subnodes
        if dis(node.centerOfMass, body.position)*distanceRate < minimumDis
            # if the body is too closed to the node, consider them one big node
            node.body = mergeBodies(node.body, body)
            return
        else
            node.numOfChild += 1
            divideToOctants(node)
            insertBody(node.children[octantsOriginal], node.body)
            insertBody(node.children[octantsNew],body)
        end
    elseif node.numOfChild > 1 # we already have subnodes
        insertBody(node.children[octantsNew], body)
        node.numOfChild += 1
    else # this is a previously empty node
        node.numOfChild += 1
        node.body = body
        node.centerOfMass = body.position
        node.mass = body.mass
        node.hasChildren = true
    end
end

# this function returns the norm given two coordinates
function dis(x::SVector{3, Float64}, y::SVector{3, Float64})
    return distanceRate*sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2 + (x[3]-y[3])^2)
end

# this function returns the direction subject to the plane of rotation (x-y)
function getDirectionOfVelocity(x::SVector{3, Float64}, y::SVector{3, Float64})
    displacement = distanceRate.*[y[1]-x[1], y[2]-x[2],y[3]-x[3]]
    velocityDir = cross(rotatingPlane, displacement)
    return velocityDir/dis(x,y)
end

# this function calculate the acceleration of a body acted by a node
function calculateAcceleration(node::Node, body::Body)
    currentAcceleration =  calculateAccelerationHelper(node, body)
    body.acceleration = strengthOfInteraction*currentAcceleration
end

function calculateAccelerationHelper(node::Node, body::Body)
    currentAcceleration = [0.0,0.0,0.0]
    if node.hasChildren && node.numOfChild > 1
        # calculate the s/d ratio
        d = dis(node.centerOfMass,body.position)
        delx = (node.maxBounds[1] - node.minBounds[1])*distanceRate
        dely = (node.maxBounds[2] - node.minBounds[2])*distanceRate
        delz = (node.maxBounds[3] - node.minBounds[3])*distanceRate
        s = max(delx,dely,delz)
        if s/d > THETA
            # this means that the node is considered close to the body
            for i in 1:8
                currentAcceleration = currentAcceleration .+ calculateAccelerationHelper(node.children[i], body)
            end
        else
            # this means that the node is far away from the body
            currentAcceleration = currentAcceleration .+ G*node.mass./((dis(node.centerOfMass,body.position)*distanceRate)^2+epsilon*distanceRate) .*(node.centerOfMass .- body.position)/dis(node.centerOfMass, body.position)
        end
    elseif node.hasChildren && node.numOfChild == 1 && node.body.id!=body.id
        if node.body.position == body.position
            println("crap")
        end
        currentAcceleration = currentAcceleration .+ G*node.mass./((dis(node.centerOfMass,body.position)*distanceRate)^2+epsilon*distanceRate) .*(node.body.position .- body.position)/dis(node.centerOfMass,body.position)# try switching direction of acceleration
    end
    return currentAcceleration
end

# this function first calls apply acceleration to get updates on the new net acceleration acted on the body,
# and then call creepHelper to update the new position and velocity for each body
function creep(node::Node)
    applyAcceleration(node)
    creepHelper(node)
end

function creepHelper(node::Node)
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            creepHelper(node.children[i])
        end
    elseif node.hasChildren && node.numOfChild == 1
        acceleration = node.body.acceleration
        time = frameRate*timeScalar
        newPosition = node.body.position*distanceRate + time .* node.body.velocity + 0.5*time^2 .*acceleration
        node.body.position = newPosition/distanceRate
        velocity = node.body.velocity .+ time .* acceleration
        speed = sqrt(velocity[1]^2+velocity[2]^2+velocity[3]^2)
        velocityDir = velocity / speed
        if speed >= 3e8
            velocity = 3e8 .* velocityDir
            println("speed of light!")
        end
        node.body.velocity = velocity
    end
end

# this function traverse the node and calculate and apply the new acceleration of each body
function applyAcceleration(node::Node)
    applyAccelerationHelper(node, node)

end

# this function is a recursive helper for applyAcceleration
function applyAccelerationHelper(originalNode::Node, node::Node)
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            applyAccelerationHelper(originalNode, node.children[i])
        end
    elseif node.hasChildren && node.numOfChild == 1
        calculateAcceleration(originalNode, node.body)
    end
end

# this function takes a node write the positions of the children nodes into a file
function writeStatsToFrame(node::Node)
    df = DataFrame( mass = Float64[], positionX = Float64[], positionY = Float64[],positionZ = Float64[],
    velocityX = Float64[], velocityY = Float64[], velocityZ = Float64[],
    accelerationX = Float64[], accelerationY = Float64[], accelerationZ = Float64[] )
    # first entry is the center of mass
    push!(df, [node.mass, node.centerOfMass[1],node.centerOfMass[2],node.centerOfMass[3],
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    writeStatsToFrameHelper(node, df)
    # sort!(df)
    return df
end

# this function is a recursive helper function for writeStatsToFile
function writeStatsToFrameHelper(node::Node, df::DataFrame)
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            writeStatsToFrameHelper(node.children[i],df)
        end
    elseif node.hasChildren && node.numOfChild == 1
        push!(df, [node.body.mass, node.body.position[1],node.body.position[2],node.body.position[3],
         node.body.velocity[1],node.body.velocity[2],node.body.velocity[3],
          node.body.acceleration[1],node.body.acceleration[2],node.body.acceleration[3]])
    end
end

# this funciton makes a body given the data from the matrix from the DataFrame from the CVS file
function makeBodyFromStats(id,stats::Array{Float64,1})
    body = Body(id, stats[1],
    [stats[2],stats[3], stats[4]],
    [stats[8],stats[9], stats[10]],
    [stats[5],stats[6], stats[7]])
    return body
end

function detectBounds(stats::Matrix{Float64})
    minimumX = Inf
    minimumY = Inf
    minimumZ = Inf
    maximumX = 0
    maximumY = 0
    maximumZ = 0

    row = size(stats)[1]
    for i in 1:row
        if stats[i,2] < minimumX
            minimumX = stats[i,2]
        end
        if stats[i,2] > maximumX
            maximumX = stats[i,2]
        end
        if stats[i,3] < minimumY
            minimumY = stats[i,3]
        end
        if stats[i,3] > maximumY
            maximumY = stats[i,3]
        end
        if stats[i,4] < minimumZ
            minimumZ = stats[i,4]
        end
        if stats[i,4] > maximumZ
            maximumZ = stats[i,4]
        end
    end
    return minimumX-100, minimumY-100, minimumZ-100, maximumX+100, maximumY+100,maximumZ+100
end

# this function make a tree from the csv file
function makeTreeFromFrame(df)
    stats = convert(Matrix{Float64}, df)
    numOfBodies = size(stats)[1]
    # make the first body (essential to initialize the tree)
    firstBody = makeBodyFromStats(0,stats[1,1:10])
    minimumX, minimumY, minimumZ, maximumX, maximumY,maximumZ = detectBounds(stats)

    # initialize the tree
    C= Node(
       false, # hasChildren
       Array{Node,1}(undef, 8), # children
       0, # numOfChild
       Body(), # body
       0.0, # mass
       SVector{3, Float64}(zeros(3)), # centerOfMass
       SVector{3, Float64}(minimumX,minimumY,minimumZ), # minBounds
       SVector{3, Float64}(maximumX,maximumY,maximumZ)) # maxBounds

    # add bodies, skip the center of mass column
    for i in 2:numOfBodies
        body = makeBodyFromStats(i,stats[i,1:10])
        insertBody(C, body)
    end

    return C
end

# this function update the tree of bodies with stable orbit velocity
function caluculateInitialSpeed(node::Node)
    caluculateInitialSpeedHelper(node, node.centerOfMass, node.mass)
end

function caluculateInitialSpeedHelper(node::Node, centerOfMass::SVector{3, Float64}, totalMass::Float64)
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            caluculateInitialSpeedHelper(node.children[i], centerOfMass, totalMass)
        end
    elseif node.hasChildren && node.numOfChild == 1 &&node.body.mass < 1e34
        node.body.velocity = sqrt(G*totalMass/dis(centerOfMass,node.body.position)/distanceRate) .*getDirectionOfVelocity(centerOfMass,node.body.position)
    end
end

# this function update the tree of bodies with random velocity
function caluculateInitialSpeedRandom(node::Node)
    caluculateInitialSpeedRandomHelper(node)
end

function caluculateInitialSpeedRandomHelper(node::Node)
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            caluculateInitialSpeedRandomHelper(node.children[i])
        end
    elseif node.hasChildren && node.numOfChild == 1
        if node.body.mass < 1e34

            node.body.velocity = [rand(1:3e7-1)+rand(),rand(1:3e7-1)+rand(),rand(1:3e7-1)+rand()]
        else

        end
    end
end

# this function update the tree of bodies with stable orbit velocity
function caluculateInitialSpeedPlus(node::Node, centerOfMassVelocity::SVector{3, Float64})
    caluculateInitialSpeedPlusHelper(node, node.centerOfMass, node.mass, centerOfMassVelocity)
end

function caluculateInitialSpeedPlusHelper(node::Node, centerOfMass::SVector{3, Float64}, totalMass::Float64, centerOfMassVelocity::SVector{3, Float64})
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            caluculateInitialSpeedPlusHelper(node.children[i], centerOfMass, totalMass,centerOfMassVelocity )
        end
    elseif node.hasChildren && node.numOfChild == 1&&node.body.mass < 1e33
        node.body.velocity = centerOfMassVelocity + sqrt(G*totalMass/dis(centerOfMass,node.body.position)/distanceRate) .*getDirectionOfVelocity(centerOfMass,node.body.position)
    elseif node.hasChildren && node.numOfChild == 1&&node.body.mass >= 1e33
        node.body.velocity = centerOfMassVelocity
    end
end

# this function update the tree of bodies with stable orbit velocity
function caluculateInitialSpeedForEnvironment(node::Node, vector::Array{Float64,1})
    caluculateInitialSpeedForEnvironmentHelper(node, vector)
end

function caluculateInitialSpeedForEnvironmentHelper(node::Node, vector::Array{Float64,1})
    if node.hasChildren && node.numOfChild > 1
        for i in 1:8
            caluculateInitialSpeedForEnvironmentHelper(node.children[i], vector)
        end
    elseif node.hasChildren && node.numOfChild == 1 && node.body.velocity == [0.0,0.0,0.0]
        x = rand()
        if x>0.5

            node.body.velocity =  (rand(0:1000) + rand(Float64)) .* vector
        end
    end
end

# this function merges two nodes into one nodes
function mergeTwoNodes(C1::Node, C2::Node)

    df = writeStatsToFrame(C1)

    writeStatsToFrameHelper(C2, df)

    C = makeTreeFromFrame(df)

    return C
end

max_stars_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_issues_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator

max_forks_repo_name: KeKeJin/Barnes-Hut-Galaxy-Simulator