text:
# Test all of the helper methods
@testset "Helper Methods" begin
    # Setup
    m = InfiniteModel();
    @infinite_parameter(m, t in [0,1])
    @infinite_parameter(m, x[1:2] in [0,5], num_supports = 2)
    @variable(m, T, Infinite(t, x))
    @variable(m, q, Infinite(t))
    oc = OrthogonalCollocation(3, GaussLobatto())
    bad_oc = OrthogonalCollocation{GaussLegendre}(3, GaussLegendre())
    tref = dispatch_variable_ref(t)
    # test generative_support_info
    @testset "generative_support_info" begin
        # test Fallback
        @test_throws ErrorException generative_support_info(TestGenMethod())
        # test NonGenerativeDerivativeMethod
        @test generative_support_info(TestMethod()) == NoGenerativeSupports()
        # test OrthogonalCollocation with GaussLobatto
        @test generative_support_info(oc) == UniformGenerativeInfo([0.5], InternalGaussLobatto)
        @test generative_support_info(OrthogonalCollocation(2)) == NoGenerativeSupports()
        # test bad OrthogonalCollocation
        @test_throws ErrorException generative_support_info(bad_oc)
    end
    # test support_label
    @testset "support_label" begin
        # test Fallback
        @test_throws ErrorException support_label(TestMethod())
        # test GenerativeDerivativeMethod
        @test support_label(oc) == InternalGaussLobatto
    end
    # test make_reduced_expr
    @testset "make_reduced_expr" begin
        # test MeasureIndex based
        meas = support_sum(T, x)
        expected = T(0, [5, 5]) + T(0, [0, 0])
        @test isequal_canonical(InfiniteOpt.make_reduced_expr(meas, t, 0.0, m), expected)
        # test InfiniteVariableIndex
        @test isequal(InfiniteOpt.make_reduced_expr(q, t, 0.0, m), q(0))
        @test isequal(InfiniteOpt.make_reduced_expr(T, t, 0.0, m), T(0, x))
        # test DerivativeIndex
        dT = deriv(T, x[1])
        dq = deriv(q, t)
        @test isequal(InfiniteOpt.make_reduced_expr(dq, t, 0.0, m), dq(0))
        v = dT(t, [x[1], 0])
        @test isequal(InfiniteOpt.make_reduced_expr(dT, x[2], 0.0, m), v)
        # test SemiInfiniteVariableIndex
        v2 = dT(t, [0, 0])
        @test isequal(InfiniteOpt.make_reduced_expr(v, x[1], 0.0, m), v2)
        @test isequal(InfiniteOpt.make_reduced_expr(v2, t, 0.0, m), dT(0, [0, 0]))
    end
end

# Test evaluate_derivative and its helpers
@testset "evaluate_derivative" begin
    # setup
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 10], num_supports = 3)
    @infinite_parameter(m, x[1:2] in [0, 1], num_supports = 2)
    @variable(m, y, Infinite(t))
    @variable(m, q, Infinite(t, x))
    meas = support_sum(q, t)
    d1 = deriv(y, t)
    d2 = deriv(q, t)
    d3 = deriv(q, x[2])
    d4 = deriv(meas, x[1])
    # test fallback
    @testset "Fallback" begin
        @test_throws ErrorException InfiniteOpt.evaluate_derivative(d1, TestMethod(), m)
    end
    # test _make_difference_expr for Forward
    @testset "_make_difference_expr (Forward)" begin
        supps = supports(t, label = All)
        expected = 5d1(0) - y(5) + y(0)
        @test isequal_canonical(InfiniteOpt._make_difference_expr(d1, y, t, 1, supps, m, Forward()), expected)
    end
    # test _make_difference_expr for Central
    @testset "_make_difference_expr (Central)" begin
        supps = supports(t, label = All)
        expected = 10d2(5, x) - q(10, x) + q(0, x)
        @test isequal_canonical(InfiniteOpt._make_difference_expr(d2, q, t, 2, supps, m, Central()), expected)
    end
    # test _make_difference_expr for FDBackward
    @testset "_make_difference_expr (Backward)" begin
        supps = sort(supports(x[2], label = All))
        expected = d3(t, [x[1], 1]) - q(t, [x[1], 1]) + q(t, [x[1], 0])
        @test isequal_canonical(InfiniteOpt._make_difference_expr(d3, q, x[2], 2, supps, m, Backward()), expected)
    end
    # test _make_difference_expr fallback
    @testset "_make_difference_expr (Fallback)" begin
        @test_throws ErrorException InfiniteOpt._make_difference_expr(d3, q, x[2], 2, [1], m, BadFiniteTech())
    end
    # test FiniteDifference with evaluate_derivative
    @testset "FiniteDifference" begin
        # test with independent parameter
        method = FiniteDifference(Central())
        exprs = [10d1(5) - y(10) + y(0)]
        @test isequal_canonical(InfiniteOpt.evaluate_derivative(d1, method, m), exprs )
        # test with dependent parameter
        method = FiniteDifference(Forward())
        exprs = [d4([0, x[2]]) - q(0, [1, x[2]]) - q(5, [1, x[2]]) -
                 q(10, [1, x[2]]) + q(0, [0, x[2]]) + q(5, [0, x[2]]) +
                 q(10, [0, x[2]])]
        @test isequal_canonical(InfiniteOpt.evaluate_derivative(d4, method, m), exprs)
        # test using Backward without boundary constraint
        method = FiniteDifference(Backward(), false)
        exprs = [5d1(5) - y(5) + y(0)]
        @test isequal_canonical(InfiniteOpt.evaluate_derivative(d1, method, m), exprs)
        # test Backward with boundary constraint
        method = FiniteDifference(Backward(), true)
        exprs = [5d1(5) - y(5) + y(0), 5d1(10) - y(10) + y(5)]
        @test isequal_canonical(InfiniteOpt.evaluate_derivative(d1, method, m), exprs)
        # test without supports
        delete_supports(x)
        @test_throws ErrorException InfiniteOpt.evaluate_derivative(d3, method, m)
        # test parameter function
        f = parameter_function(sin, t)
        df = deriv(f, t)
        exprs = [5df(5) - sin(5) + sin(0), 5df(10) - sin(10) + sin(5)]
        @test isequal_canonical(InfiniteOpt.evaluate_derivative(df, method, m), exprs)
    end
    # test OrthogonalCollocation with evaluate_derivative
    @testset "OrthogonalCollocation" begin
        # test first solve
        method = OrthogonalCollocation(3)
        set_derivative_method(t, method)
        Mt = [1. 5.; 1. 10]' \ [2.5 2.5^2; 5. 25.]'
        M = Mt'
        exprs = [@expression(m, M[1, 1] * d2(2.5, x) + M[1, 2] * d2(5, x) - q(2.5, x) + q(0, x)),
                 @expression(m, M[2, 1] * d2(2.5, x) + M[2, 2] * d2(2.5, x) - q(5, x) + q(0, x)),
                 @expression(m, M[1, 1] * d2(7.5, x) + M[1, 2] * d2(10, x) - q(7.5, x) + q(5, x)),
                 @expression(m, M[2, 1] * d2(7.5, x) + M[2, 2] * d2(7.5, x) - q(10, x) + q(5, x))]
        delete_supports(t, label = UserDefined)
        function rm_zeros(exs)
            for e in exs
                for (v, c) in e.terms
                    if abs(c) < 1e-15
                        delete!(e.terms, v)
                    end
                end
            end
            return exs
        end
        @test isequal_canonical(rm_zeros(evaluate_derivative(d2, method, m)), exprs)
        @test supports(t) == [0, 5, 10]
        @test supports(t, label = All) == [0, 2.5, 5, 7.5, 10]
        # test resolve
        @test isequal_canonical(rm_zeros(evaluate_derivative(d2, method, m)), exprs)
        @test supports(t) == [0, 5, 10]
        @test supports(t, label = All) == [0, 2.5, 5, 7.5, 10]
        @test has_generative_supports(t)
        @test has_internal_supports(t)
    end
end

# Test User Methods
@testset "User Methods" begin
    # setup
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 10], supports = [0, 3, 10])
    @infinite_parameter(m, x in [-1, 1], supports = [-1, 1],
                        derivative_method = OrthogonalCollocation(4))
    @variable(m, y, Infinite(t))
    @variable(m, q, Infinite(x))
    dy = @deriv(y, t)
    dq = @deriv(q, x)
    dy2 = @deriv(y, t^2)
    # test evaluate
    @testset "evaluate" begin
        # test evaluate first derivative
        @test derivative_method(t) isa FiniteDifference
        @test evaluate(dy) isa Nothing
        @test num_constraints(m) == 2
        @test !has_generative_supports(t)
        @test has_derivative_constraints(t)
        @test has_derivative_constraints(dy)
        @test length(derivative_constraints(dy)) == 2
        # test evaluate second derivative
        @test derivative_method(t) isa FiniteDifference
        @test evaluate(dy2) isa Nothing
        @test num_constraints(m) == 4
        @test !has_generative_supports(t)
        @test has_derivative_constraints(t)
        @test has_derivative_constraints(dy2)
        @test length(derivative_constraints(dy2)) == 2
    end
    # test evaluate_all_derivatives!
    @testset "evaluate_all_derivatives!" begin
        @test derivative_method(x) isa OrthogonalCollocation
        @test evaluate_all_derivatives!(m) isa Nothing
        @test num_constraints(m) == 7
        @test has_generative_supports(x)
        @test has_internal_supports(x)
        @test supports(x) == [-1, 1]
        @test num_supports(x, label = All) == 4
        @test supports(t) == [0, 3, 10]
        @test length(derivative_constraints(dy)) == 2
        @test length(derivative_constraints(dy2)) == 2
        @test length(derivative_constraints(dq)) == 3
        @test has_derivative_constraints(t)
        @test has_derivative_constraints(x)
        @test has_derivative_constraints(dy)
        @test has_derivative_constraints(dy2)
        @test has_derivative_constraints(dq)
    end
end


text:
module TraceUtils

export dif, find, pkpk

import Base: find

function dif(v::Vector{T}) where {T<:Real}
    return v[2]-v[1]
end

function find(v::Vector{T}, x::Number) where {T<:Real}
    return searchsortedlast(v, x)
end

function pkpk(v::Vector{T}) where {T<:Real}
    return maximum(v) - minimum(v)
end

end

ext: jl


text:
using Test
using Distances: Euclidean
using SparseArrays
using LinearAlgebra
using UMAP
using UMAP: fuzzy_simplicial_set, compute_membership_strengths, smooth_knn_dists, smooth_knn_dist, spectral_layout, optimize_embedding, pairwise_knn

include("umap_tests.jl")


text:
VERSION >= v"0.0.2" && __precompile__(true)

module FitnessForService

export

    Part5LTAFlawSizeLimitCriteria, Part5LTAFlawSizeLevel1Acceptance, sc, Part5LTALevel1, CreepRangeTemperature, DesignCodeCriteria, MaterialToughness, CyclicService, Part5ComponentType, Part5AsessmentApplicability

    include("applicability.jl")
    include("part5_LTA_functions.jl")
    include("data_grid_functions.jl")
    include("thickness_MAWP_stress_equations_for_FFS.jl")
    include("part5_local_metal_loss_assessment.jl")

end


text:
"""
    PeriodicStorage{T,scalar} <: AbstractStorage{T, scalar}

Storage for a (discontinuous) shift invariant, periodic function ``G(τ1, τ2) = G(τ1 - τ2)``, ``G(τ) = G(τ + β)``. \\
Represents data on the uniform grid ``{0, Δτ, 2Δτ, ..., (N-1)Δτ}`` where ``Δτ = β/(N-1)``. \\
Note don't generically have ``G(0) = G(β)`` so we need to store both endpoints. \\
Because of this, there is a discontinuity along the diagonal and the storage is not quite circulant. \\
Diagonal discontinuity is handled by `greater` argument to `getindex` method.
"""
struct PeriodicStorage{T,scalar} <: AbstractStorage{T, scalar}
  data::Array{T,3}
  norb::Int
  N::Int
  M::Int

  function PeriodicStorage{T,scalar}(data::AbstractArray{T,3}) where {T, scalar}
    norb = size(data,1)
    N = size(data,3)
    @assert size(data) == (norb, norb, N)
    @assert !(scalar && norb > 1)
    new(data,norb,N,N)
  end
end

function PeriodicStorage(data::AbstractArray{T,3}, scalar=false) where T
  PeriodicStorage{T,scalar}(data)
end

function PeriodicStorage(::Type{T}, N::Integer, norb=1, scalar=false) where T <: Number
  data = zeros(T, norb, norb, N)
  return PeriodicStorage(data, scalar)
end

PeriodicStorage(N::Integer, norb=1, scalar=false) = PeriodicStorage(ComplexF64, N, norb, scalar)

function Base.getindex(X::PeriodicStorage{T,scalar}, k, l, i::Int, j::Int, greater=true) where {T,scalar}
  greater = i == j ? greater : i > j
  return greater ? X.data[k,l,i-j+1] : X.data[k,l,i-j+X.N]
end
# need to define these methods to forward greater argument
Base.getindex(X::PeriodicStorage{T,true}, i, j, greater=true) where {T} = X[1,1,i,j, greater]
Base.getindex(X::PeriodicStorage{T,false}, i, j, greater=true) where {T} = X[:,:,i,j, greater]

function Base.setindex!(X::PeriodicStorage{T,scalar}, v, k, l, i::Int, j::Int) where {T,scalar}
  greater = i >= j
  if greater
    return X.data[k,l,i-j+1] = v
  else
    return X.data[k,l,i-j+X.N] = v
  end
end


text:
##### SGS TKE tendencies
export ShearProduction,
    BuoyancyProduction,
    Dissipation

#####
##### First order flux (Advective)
#####

function flux(::SGSTKE, ::Advect, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.sgstke.ρe_SGS
end

#####
##### Second order flux (Diffusive)
#####

#struct SGSDiffFlux <: TendencyDef{Flux{SecondOrder}} end
function flux(::SGSTKE, ::Diffusion, atmos, args)
    @unpack state, diffusive = args
    @unpack ν, D_t = args.precomputed.turbulence

    FT = eltype(state)
    K_m = (ν isa SDiagonal) ? diag(ν) : ν
    d_e_SGS = -FT(2) * K_m .* diffusive.sgstke.∇e_SGS
    return d_e_SGS * state.ρ
    # technically K_m instead of ν. Need to include this into TurbulenceClosures.
end

#####
##### Sources (Productions and dissipation)
#####

struct ShearProduction <: TendencyDef{Source} end
prognostic_vars(::ShearProduction) = (SGSTKE(),)
struct BuoyancyProduction <: TendencyDef{Source} end
prognostic_vars(::BuoyancyProduction) = (SGSTKE(),)

function source(::SGSTKE, ::ShearProduction, atmos, args)
    @unpack state, diffusive = args
    @unpack τ = args.precomputed.turbulence

    return - state.ρ * sum(τ .* diffusive.sgstke.∇u)
end

function source(::SGSTKE, ::BuoyancyProduction, atmos, args)
    @unpack state, aux, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    param_set = parameter_set(atmos)

    FT = eltype(state)
    g = FT(grav(param_set))

    K_h = D_t # could be FT or SVector{3, FT}
    flux_θ_liq_ice = state.ρ * (-K_h) .* diffusive.sgstke.∇θ_liq_ice
    flux_q_t = state.ρ * (-K_h) .* diffusive.moisture.∇q_tot

    # Deardorff 1980, for unsaturated air
    A = FT(1) + FT(0.61) * state.moisture.ρq_tot / state.ρ
    B = FT(0.61)
    flux_θ_vir = A * flux_θ_liq_ice[3] + B * flux_q_t[3]
    return (g / aux.ref_state.T ) * flux_θ_vir
end

# The effective geometric mean grid resolution at the point
struct Dissipation <: TendencyDef{Source} end
prognostic_vars(::Dissipation) = (SGSTKE(),)

function source(::SGSTKE, ::Dissipation, atmos, args)
    @unpack state, aux, diffusive = args
    @unpack D_t = args.precomputed.turbulence

    #turbulence_model(atmos) isa ConstantViscosity ? error("ConstantViscosity not supported") : nothing

    FT = eltype(state)

    param_set = parameter_set(atmos)
    g = FT(grav(param_set))

    Δs = aux.turbulence.Δ
    e = state.sgstke.ρe_SGS / state.ρ

    l_s = FT(0.76) * sqrt(e) / sqrt(abs((g / aux.ref_state.T ) * diffusive.sgstke.∇θ_liq_ice[3]))
    l = (l_s < Δs) ? l_s : Δs
    C = FT(0.19) + FT(0.51) * (l / Δs)
    return - state.ρ * C * e^FT(1.5) / l
end


text:
function add(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "add", RealMatrix, (RealMatrix,), arg0)
end

function add_to_entry(obj::RealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function copy(obj::RealMatrix)
    return jcall(obj, "copy", RealMatrix, ())
end

function copy_sub_matrix(obj::RealMatrix, arg0::Vector{jint}, arg1::Vector{jint}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "copySubMatrix", void, (Vector{jint}, Vector{jint}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function copy_sub_matrix(obj::RealMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::Vector{Vector{jdouble}})
    return jcall(obj, "copySubMatrix", void, (jint, jint, jint, jint, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3, arg4)
end

function create_matrix(obj::RealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", RealMatrix, (jint, jint), arg0, arg1)
end

function get_column(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getColumn", Vector{jdouble}, (jint,), arg0)
end

function get_column_matrix(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getColumnMatrix", RealMatrix, (jint,), arg0)
end

function get_column_vector(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getColumnVector", RealVector, (jint,), arg0)
end

function get_data(obj::RealMatrix)
    return jcall(obj, "getData", Vector{Vector{jdouble}}, ())
end

function get_entry(obj::RealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", jdouble, (jint, jint), arg0, arg1)
end

function get_frobenius_norm(obj::RealMatrix)
    return jcall(obj, "getFrobeniusNorm", jdouble, ())
end

function get_row(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getRow", Vector{jdouble}, (jint,), arg0)
end

function get_row_matrix(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getRowMatrix", RealMatrix, (jint,), arg0)
end

function get_row_vector(obj::RealMatrix, arg0::jint)
    return jcall(obj, "getRowVector", RealVector, (jint,), arg0)
end

function get_sub_matrix(obj::RealMatrix, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(obj, "getSubMatrix", RealMatrix, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function get_sub_matrix(obj::RealMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(obj, "getSubMatrix", RealMatrix, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function get_trace(obj::RealMatrix)
    return jcall(obj, "getTrace", jdouble, ())
end

function multiply(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiply", RealMatrix, (RealMatrix,), arg0)
end

function multiply_entry(obj::RealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "multiplyEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function operate(obj::RealMatrix, arg0::RealVector)
    return jcall(obj, "operate", RealVector, (RealVector,), arg0)
end

function operate(obj::RealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "operate", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function power(obj::RealMatrix, arg0::jint)
    return jcall(obj, "power", RealMatrix, (jint,), arg0)
end

function pre_multiply(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "preMultiply", RealMatrix, (RealMatrix,), arg0)
end

function pre_multiply(obj::RealMatrix, arg0::RealVector)
    return jcall(obj, "preMultiply", RealVector, (RealVector,), arg0)
end

function pre_multiply(obj::RealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "preMultiply", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function scalar_add(obj::RealMatrix, arg0::jdouble)
    return jcall(obj, "scalarAdd", RealMatrix, (jdouble,), arg0)
end

function scalar_multiply(obj::RealMatrix, arg0::jdouble)
    return jcall(obj, "scalarMultiply", RealMatrix, (jdouble,), arg0)
end

function set_column(obj::RealMatrix, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setColumn", void, (jint, Vector{jdouble}), arg0, arg1)
end

function set_column_matrix(obj::RealMatrix, arg0::jint, arg1::RealMatrix)
    return jcall(obj, "setColumnMatrix", void, (jint, RealMatrix), arg0, arg1)
end

function set_column_vector(obj::RealMatrix, arg0::jint, arg1::RealVector)
    return jcall(obj, "setColumnVector", void, (jint, RealVector), arg0, arg1)
end

function set_entry(obj::RealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "setEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function set_row(obj::RealMatrix, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setRow", void, (jint, Vector{jdouble}), arg0, arg1)
end

function set_row_matrix(obj::RealMatrix, arg0::jint, arg1::RealMatrix)
    return jcall(obj, "setRowMatrix", void, (jint, RealMatrix), arg0, arg1)
end

function set_row_vector(obj::RealMatrix, arg0::jint, arg1::RealVector)
    return jcall(obj, "setRowVector", void, (jint, RealVector), arg0, arg1)
end

function set_sub_matrix(obj::RealMatrix, arg0::Vector{Vector{jdouble}}, arg1::jint, arg2::jint)
    return jcall(obj, "setSubMatrix", void, (Vector{Vector{jdouble}}, jint, jint), arg0, arg1, arg2)
end

function subtract(obj::RealMatrix, arg0::RealMatrix)
    return jcall(obj, "subtract", RealMatrix, (RealMatrix,), arg0)
end

function transpose(obj::RealMatrix)
    return jcall(obj, "transpose", RealMatrix, ())
end

function walk_in_column_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_column_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_column_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_column_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_row_order(obj::RealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_row_order(obj::RealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end


text:
@doc raw"""
    struct SymmetryGroup{G<:Group, N::Integer}

Stores the charges and degeneracies?
TODO: docstring
"""
struct SymmetryGroup{G<:Group}
    charges::Array{<:AbstractIrrep{G},4}
    degeneracy::Dict{<:AbstractIrrep{G},Integer}
    basis::GradedSpace{<:AbstractIrrep{G}}
end

@doc raw"""
    SymmetryGroup{G}(g::TensorMap{ComplexSpace}) where {G<:Group}

Constructs a `SymmetryGroup` type given a group `G` and a group generator `g`.
"""
function SymmetryGroup(
    g::TensorMap{ComplexSpace}, basis::GradedSpace{<:AbstractIrrep{G}}
) where {G<:Group}
    return SymmetryGroup{G}(find_charges(G, g)..., basis)
end

struct Liouvillian{
    S<:IndexSpace,
    I<:Sector,
    A<:Union{TensorKit.SectorDict{I,<:DenseMatrix},<:DenseMatrix},
    F₁,
    F₂,
} <: AbstractPhysicalTensor{S,4,4,I,A,F₁,F₂}
    tensor::TensorMap{S,4,4,I,A,F₁,F₂}
    mapping
end

const SymmetricLiouvillian{S<:GradedSpace,I<:Sector,A<:TensorKit.SectorDict{I,<:DenseMatrix},F₁,F₂} = Liouvillian{S,I,A, F₁,F₂}

function SymmetricLiouvillian(t::TensorMap{S,4,4,I, A, F₁, F₂}, mapping ) where {S<:GradedSpace, I<:Sector, A<:TensorKit.SectorDict{I,<:DenseMatrix}, F₁,F₂}
    return SymmetricLiouvillian{S,I,A,F₁,F₂}(t, mapping)
end

const StandardLiouvillian{A<:DenseMatrix} = Liouvillian{ComplexSpace,Trivial,A,Nothing,Nothing}

function StandardLiouvillian(
    t::TensorMap{ComplexSpace,4,4,Trivial,A,Nothing,Nothing}, nothing
) where {A<:DenseMatrix}
    return StandardLiouvillian{A}(t, nothing)
end

@doc raw"""
    find_charges(G, g)

Given the vectorised form of a group generator, `g`, outputs a list of the charges and
their respective degeneracies. Assumes `g` is diagonalised (i.e. assumes that the basis used
is one such that `g` is diagonal), and that each diagonal element of `g` corresponds to a
valid irrep label of the group `G`.
"""
function find_charges(G::Type{<:TensorKit.Group}, g)
    chrgs = Array{Irrep[G]}(undef, 2, 2, 2, 2)
    degn = Dict{Irrep[G],Integer}([])
    for i in 1:2, j in 1:2, k in 1:2, l in 1:2
        chg = g[i, j, k, l, i, j, k, l]
        irrep = Irrep[G](Integer(chg))
        chrgs[i, j, k, l] = irrep # TODO: Convert `chg` to (half)-integer
        if haskey(degn, irrep)
            degn[irrep] = degn[irrep] + 1
        else
            degn[irrep] = 1
        end
    end
    return chrgs, degn
end

@doc raw"""
    find_blocks(chrgs, degn, L::TensorMap{ComplexSpace, 4, 4})

Outputs the blocks corresponding to each charge sector of a symmetric Liouvillian.
"""
function find_blocks(SG::SymmetryGroup{G}, L::StandardLiouvillian) where {G<:Group}
    L = L
    chrgs = SG.charges
    degn = SG.degeneracy
    blocks = Dict{Irrep[G],Matrix{ComplexF64}}([])
    cols = Dict{Irrep[G],Vector{Vector{ComplexF64}}}([])
    mapping = Dict{Irrep[G],Vector{CartesianIndex}}([])
    for key in keys(degn)
        cols[key] = []
        mapping[key] = []
    end
    # For each input state, find the linear combination of output states that this state
    # maps to. Output states should have the same charge as the input state.
    for i in 1:2, j in 1:2, k in 1:2, l in 1:2 # Input state
        chrg = chrgs[i, j, k, l]
        col = Vector{ComplexF64}()
        for a in 1:2, b in 1:2, c in 1:2, d in 1:2 # Output state
            if chrgs[a, b, c, d] == chrg
                push!(col, L[a, b, c, d, i, j, k, l])
            else
                # Check if symmetry-required non-zero elements are actually zero. If not
                # then L is not block diagonal so throw an error.
                if !(L[a, b, c, d, i, j, k, l] ≈ zero(L[a, b, c, d, i, j, k, l]))
                    throw(
                        ArgumentError(
                            "Liouvillian has non-zero elements at incompatible positions"
                        ),
                    )
                end
            end
        end
        push!(cols[chrg], col)
        push!(mapping[chrg], CartesianIndex(i, j, k, l))
    end
    for chrg in keys(degn)
        col = cols[chrg]
        if length(col) == 1
            ele = reshape(col[1], 1, 1)
            blocks[chrg] = ele
        else
            blocks[chrg] = reduce(hcat, col)
        end
    end
    return blocks, mapping
end

function symmetrise(SG::SymmetryGroup{<:Group}, L::StandardLiouvillian)
    blcks, mapping = find_blocks(SG, L)
    lin = TensorMap(zeros, ComplexF64, SG.basis^4, SG.basis^4)
    for i in keys(SG.degeneracy)
        blocks(lin)[i] = blcks[i]
    end
    return SymmetricLiouvillian(lin, mapping)
end

function unsymmetrise(L::SymmetricLiouvillian)
    tensor = L.tensor
    map = L.mapping
    blcks = blocks(tensor)
    cL = zeros(ComplexF64, 2, 2, 2, 2, 2, 2, 2, 2)
    # cL = TensorMap(zeros, ComplexF64, ℂ^4, ℂ^4)
    for (irrep, block) in blcks
        irrepmap = map[irrep]
        for i in CartesianIndices(block)
            row = i[1]
            col = i[2]
            cL[irrepmap[row], irrepmap[col]] = block[i]
        end
    end
    return StandardLiouvillian(
        TensorMap(cL, ℂ^2 * ℂ^2 * ℂ^2 * ℂ^2, ℂ^2 * ℂ^2 * ℂ^2 * ℂ^2), nothing
    )
end

#TODO: nothing -> identity mapping

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
struct IonImage

end


text:
import Base: @__doc__
using PowerDynamics: @DynamicNode
import PowerDynamics: construct_vertex, dimension, symbolsof, showdefinition, AbstractNode
using LinearAlgebra: Diagonal

begin
    @__doc__ struct CSImin_limit <: AbstractNode
            I_r
            imax
            Y_n
        end
    CSImin_limit(; I_r, imax, Y_n=0) = CSImin_limit(I_r, imax, Y_n)
    function construct_vertex(par::CSImin_limit)
        I_r = par.I_r
        imax = par.imax
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            i = total_current(e_s, e_d) + Y_n * u

            id_r = SmoothClamp(real(i), imax; scale=500.)
            iq_r = SmoothClamp(imag(i), sqrt(abs(imax^2 - id_r^2)); scale=500.)

            du = complex(id_r, iq_r) - I_r

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=2, mass_matrix=[0.0, 0.0], sym=Symbol[:u_r, :u_i])
    end
    symbolsof(::CSImin_limit) = begin
            [:u_r, :u_i]
        end
    dimension(::CSImin_limit) = begin
            2
        end
end

begin
    @__doc__ struct CSImin_FRT <: AbstractNode
            I_r
            Vac_ref
            V_dead
            imax
            k_FRT
            Y_n
        end
    CSImin_FRT(; I_r, Vac_ref, V_dead, imax, k_FRT, Y_n=0) = CSImin_FRT(I_r, Vac_ref, V_dead, imax, k_FRT, Y_n)
    function construct_vertex(par::CSImin_FRT)
        I_r = par.I_r
        Vac_ref = par.Vac_ref
        V_dead = par.V_dead
        imax = par.imax
        k_FRT = par.k_FRT
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            i = total_current(e_s, e_d) + Y_n * u

            Va_amp = abs(u)

            #err_state = SmoothStep(Va_amp, Vac_ref + V_dead, Vac_ref - V_dead; order=100.)
            err_state = Vac_ref - V_dead < Va_amp < Vac_ref + V_dead

            #fault err != 0
            id_r = SmoothClamp(real(i), imax; scale=200.)
            iq_r = SmoothClamp(imag(i), sqrt(abs(imax^2 - id_r^2)); scale=200.)

            iq_plus = 0.
            if Va_amp > Vac_ref + V_dead
                iq_plus = - k_FRT * (Va_amp - (Vac_ref + V_dead))
            elseif Va_amp < Vac_ref - V_dead
                iq_plus = - k_FRT * (Va_amp - (Vac_ref - V_dead))
            end
            iq_r_plus = SmoothClamp(iq_plus, imax; scale=200.)
            id_r_plus = SmoothClamp(iq_r_plus, sqrt(abs(imax^2 - iq_r^2)); scale=200.)

            #i_mix = err_state * complex(id_r, iq_r) + (1. - err_state) * complex(id_r_plus, iq_r_plus)
            i_mix = err_state ? complex(id_r, iq_r) : complex(id_r_plus, iq_r_plus)

            du = i_mix - I_r

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=2, mass_matrix=[0.0, 0.0], sym=Symbol[:u_r, :u_i])
    end
    symbolsof(::CSImin_FRT) = begin
            [:u_r, :u_i]
        end
    dimension(::CSImin_FRT) = begin
            2
        end
end

begin
    @__doc__ struct CSImin_PLL <: AbstractNode
            I_r
            K_pll
            Y_n
        end
    CSImin_PLL(; I_r, K_pll, Y_n=0) = CSImin_PLL(I_r, K_pll, Y_n)
    function construct_vertex(par::CSImin_PLL)
        I_r = par.I_r
        K_pll = par.K_pll
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            u = complex(x[1], x[2])
            i = total_current(e_s, e_d) + Y_n * u
            theta = x[3]

            dtheta = - K_pll * abs(u) * sin(theta - angle(u))
            du = i - I_r

            try
                dx[1] = real(du)
                dx[2] = imag(du)
                dx[3] = dtheta
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim=3, mass_matrix=Diagonal([0, 0, 1]), sym=Symbol[:u_r, :u_i, :θ])
    end
    symbolsof(::CSImin_PLL) = begin
            [:u_r, :u_i, :θ]
        end
    dimension(::CSImin_PLL) = begin
            3
        end
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
"""Check whether elements of a vector are increasing."""
function isincreasing(x::Vector)
    is_increasing = true
    for i in 2:length(x)
        is_increasing &= x[i] > x[i-1] ? true : false
    end
    return is_increasing
end


text:
facts("auth_token file exists") do
    context("auth_token is stubbed with message") do
        @fact length(readdlm(Pkg.dir("Quandl/src/token/auth_token.jl"))) => 28
    end
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, NESSie

const pages = [
    "Home" => "index.md",
    "Data formats" => [
        "formats/input.md",
        "formats/output.md"
    ],
    "Library" => [
        "lib/constants.md",
        "lib/electrostatics.md",
        "lib/models.md",
        "lib/quadrature.md",
        "lib/solvers.md",
        "lib/util.md"
    ],
    "Internals" => [
        "intern/base.md",
        "intern/bem.md",
        "intern/format.md",
        "intern/radon.md",
        "intern/rjasanow.md",
        "intern/testmodel.md"
    ]
]

makedocs(
    modules   = [NESSie, NESSie.BEM, NESSie.Format, NESSie.TestModel],
    clean     = true,
    doctest   = true,
    linkcheck = true,
    checkdocs = :all,
    pages     = pages,
    format    = Documenter.HTML(),
    sitename  = "NESSie.jl",
    repo      = "https://github.com/tkemmer/NESSie.jl/blob/master{path}"
)

deploydocs(
    repo = "github.com/tkemmer/NESSie.jl.git"
)


text:
module PiecewiseDeterministicMarkovProcesses
	using Random, LinearAlgebra, SparseArrays, Parameters
	using LSODA, Sundials, DiffEqJump, RecursiveArrayTools, DiffEqBase, SparseArrays
	using ForwardDiff
	using DiffEqBase, DiffEqJump
	import DiffEqBase: solve

	abstract type AbstractPDMPAlgorithm end
	abstract type AbstractCHV <: AbstractPDMPAlgorithm end
	abstract type AbstractCHVIterator <: AbstractCHV end
	abstract type AbstractRejection <: AbstractPDMPAlgorithm end
	abstract type AbstractRejectionExact <: AbstractRejection end
	abstract type AbstractRejectionIterator <: AbstractRejection end

	include("utilsforwarddiff.jl")
	include("rate.jl")
	include("problem.jl")
	include("utils.jl")
	include("chv.jl")
	include("chvdiffeq.jl")
	include("rejectiondiffeq.jl")
	include("rejection.jl")
	include("tau-leap.jl")
	include("diffeqwrap.jl")

	export ssa,
		chv!,chv,
		rejection!,
		rejection_exact,
		chv_diffeq!,
		rejection_diffeq!,
		pdmpArgs,
		pdmpResult,
		pdmp_data,
		ConstantRate, VariableRate, CompositeRate

	export PDMPProblem, CHV, Rejection, RejectionExact, solve
end # module


text:
# Quadratic decision function #
# ========================== #
struct QuadraticDecisionObjectiveBridge{T, LQ} <: MOIB.Objective.AbstractBridge
    decision_function::QuadraticDecisionFunction{T, LQ}
end

function MOIB.Objective.bridge_objective(::Type{QuadraticDecisionObjectiveBridge{T,LQ}}, model::MOI.ModelLike,
                                         f::QuadraticDecisionFunction{T, LQ}) where {T, LQ <: LinearPart{T}}
    # All decisions have been mapped to either the decision part constant or the variable part terms
    # at this point.
    F = MOI.ScalarAffineFunction{T}
    lq = f.linear_quadratic_terms
    # Set the bridged objective
    MOI.set(model, MOI.ObjectiveFunction{F}(),
            MOI.ScalarAffineFunction(
                lq.variable_part.terms,
                zero(T)))
    # Save decision function to allow modifications
    return QuadraticDecisionObjectiveBridge{T, LQ}(f)
end

function MOIB.Objective.bridge_objective(::Type{QuadraticDecisionObjectiveBridge{T,LQ}}, model::MOI.ModelLike,
                                         f::QuadraticDecisionFunction{T, LQ}) where {T, LQ <: QuadraticPart{T}}
    # All decisions have been mapped to either the decision part constant or the variable part terms
    # at this point.
    F = MOI.ScalarQuadraticFunction{T}
    lq = f.linear_quadratic_terms
    # Set the bridged objective
    MOI.set(model, MOI.ObjectiveFunction{F}(),
            MOI.ScalarQuadraticFunction(
                lq.variable_part.affine_terms,
                lq.variable_part.quadratic_terms,
                zero(T)))
    # Save decision function to allow modifications
    return QuadraticDecisionObjectiveBridge{T, LQ}(f)
end

function MOIB.Objective.supports_objective_function(
    ::Type{<:QuadraticDecisionObjectiveBridge}, ::Type{<:QuadraticDecisionFunction})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:QuadraticDecisionObjectiveBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:QuadraticDecisionObjectiveBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.Objective.concrete_bridge_type(::Type{<:QuadraticDecisionObjectiveBridge{T}},
                                             ::Type{QuadraticDecisionFunction{T, LinearPart{T}}}) where T
    return QuadraticDecisionObjectiveBridge{T, LinearPart{T}}
end
function MOIB.Objective.concrete_bridge_type(::Type{<:QuadraticDecisionObjectiveBridge{T}},
                                             ::Type{QuadraticDecisionFunction{T, QuadraticPart{T}}}) where T
    return QuadraticDecisionObjectiveBridge{T, QuadraticPart{T}}
end
function MOIB.set_objective_function_type(::Type{QuadraticDecisionObjectiveBridge{T,LinearPart{T}}}) where T
    return MOI.ScalarAffineFunction{T}
end
function MOIB.set_objective_function_type(::Type{QuadraticDecisionObjectiveBridge{T,QuadraticPart{T}}}) where T
    return MOI.ScalarQuadraticFunction{T}
end

function MOI.get(::QuadraticDecisionObjectiveBridge, ::MOI.NumberOfVariables)
    return 0
end

function MOI.get(::QuadraticDecisionObjectiveBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

function MOI.delete(::MOI.ModelLike, ::QuadraticDecisionObjectiveBridge)
    # Nothing to delete
    return nothing
end

function MOI.set(::MOI.ModelLike, ::MOI.ObjectiveSense,
                 ::QuadraticDecisionObjectiveBridge, ::MOI.OptimizationSense)
    # Nothing to handle if sense changes
    return nothing
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{QuadraticDecisionFunction{T, LinearPart{T}}},
                 bridge::QuadraticDecisionObjectiveBridge{T, LinearPart{T}}) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    G = MOI.ScalarAffineFunction{T}
    obj_val = MOI.get(model, MOIB.ObjectiveFunctionValue{G}(attr.result_index))
    # Calculate and add constant
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    return obj_val + constant
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{QuadraticDecisionFunction{T, QuadraticPart{T}}},
                 bridge::QuadraticDecisionObjectiveBridge{T, QuadraticPart{T}}) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    G = MOI.ScalarQuadraticFunction{T}
    obj_val = MOI.get(model, MOIB.ObjectiveFunctionValue{G}(attr.result_index))
    # Calculate and add constant
    constant = lq.variable_part.constant +
        lq.decision_part.constant +
        f.known_part.constant +
        f.known_decision_terms.constant
    return obj_val + constant
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOI.ObjectiveFunction{QuadraticDecisionFunction{T, LinearPart{T}}},
                 bridge::QuadraticDecisionObjectiveBridge{T, LinearPart{T}}) where T
    f = bridge.decision_function
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     f.decision_part.affine_terms)
        result |= any(t -> t.variable_index_1 == v || t.variable_index_2 == v,
                      f.decision_part.quadratic_terms)
        result |= any(t -> t.variable_index_1 == v,
                      f.cross_terms.affine_terms)
        result |= any(t -> t.variable_index == v,
                      f.known_decision_terms.affine_terms)
    end
    g = QuadraticDecisionFunction(
        LinearPart(
            MOIU.filter_variables(v -> !from_decision(v), f.variable_part),
            copy(f.decision_part)),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
    return g
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOI.ObjectiveFunction{QuadraticDecisionFunction{T, QuadraticPart{T}}},
                 bridge::QuadraticDecisionObjectiveBridge{T, QuadraticPart{T}}) where T
    f = bridge.decision_function
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     f.decision_part.affine_terms)
        result |= any(t -> t.variable_index_1 == v || t.variable_index_2 == v,
                      f.decision_part.quadratic_terms)
        result |= any(t -> t.variable_index_1 == v,
                      f.cross_terms.affine_terms)
        result |= any(t -> t.variable_index == v,
                      f.known_decision_terms.affine_terms)
    end
    g = QuadraticDecisionFunction(
        QuadraticPart(
            MOIU.filter_variables(v -> !from_decision(v), f.variable_part),
            copy(f.decision_part),
            copy(f.cross_terms)),
        copy(f.known_part),
        copy(f.known_variable_terms),
        copy(f.known_decision_terms))
    return g
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::MOI.ScalarConstantChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Modify constant of variable part
    lq.variable_part.constant = change.new_constant
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,LinearPart{T}}, change::MOI.ScalarCoefficientChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Update coefficient in variable part
    modify_coefficient!(lq.variable_part.terms, change.variable, change.new_coefficient)
    # Modify decision part if variable originated from a decision
    i = something(findfirst(t -> t.variable_index == change.variable,
                            lq.decision_part.terms), 0)
    if !iszero(i)
        modify_coefficient!(lq.decision_part.terms, change.variable, change.new_coefficient)
    end
    # Modify the variable part of the objective as usual
    F = MOI.ScalarAffineFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), change)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,QuadraticPart{T}}, change::MOI.ScalarCoefficientChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Update variable part
    modify_coefficient!(lq.variable_part.affine_terms, change.variable, change.new_coefficient)
    # Modify decision part if variable originated from a decision
    i = something(findfirst(t -> t.variable_index == change.variable,
                            lq.decision_part.affine_terms), 0)
    if !iszero(i)
        modify_coefficient!(lq.decision_part.affine_terms, change.variable, change.new_coefficient)
    end
    # Complete rebridge required
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::DecisionCoefficientChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    variable_terms = lq isa LinearPart ? lq.variable_part.terms : lq.variable_part.affine_terms
    decision_terms = lq isa LinearPart ? lq.decision_part.terms : lq.decision_part.affine_terms
    # Check if decision has been mapped
    i = something(findfirst(t -> t.variable_index == change.decision,
                            variable_terms), 0)
    if !iszero(i)
        # Update mapped variable through ScalarCoefficientChange
        F = QuadraticDecisionFunction{T}
        MOI.modify(model, MOI.ObjectiveFunction{F}(),
                   MOI.ScalarCoefficientChange(change.decision, change.new_coefficient))
    else
        i = something(findfirst(t -> t.variable_index == change.decision,
                                decision_terms), 0)
        # Query the fixed decision value
        unbridged = MOIB.unbridged_variable_function(model, change.decision)::MOI.SingleVariable
        decision_value = MOI.get(model, MOI.VariablePrimal(), unbridged.variable)
        # Update decision part of objective constant
        coefficient = iszero(i) ? zero(T) : decision_terms[i].coefficient
        lq.decision_part.constant +=
            (change.new_coefficient - coefficient) * decision_value
    end
    # Update coefficient in decision part
    modify_coefficient!(decision_terms, change.decision, change.new_coefficient)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::KnownCoefficientChange) where T
    f = bridge.decision_function
    i = something(findfirst(t -> t.variable_index == change.known,
                            f.known_part.affine_terms), 0)
    # Update known part of objective constant
    coefficient = iszero(i) ? zero(T) : f.known_part.affine_terms[i].coefficient
    known_value = MOI.get(model, MOI.VariablePrimal(), change.known)
    f.known_part.constant +=
        (change.new_coefficient - coefficient) * known_value
    modify_coefficient!(f.known_part.affine_terms, change.known, change.new_coefficient)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,LinearPart{T}}, ::DecisionStateChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    F = MOI.ScalarAffineFunction{T}
    # Switch on state transition
    if change.new_state == NotTaken
        # Check if decision is included in the decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.decision_part.terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = lq.decision_part.terms[i].coefficient
            # Update fixed decision value
            change.value_difference < 0 ||
                error("Decision value update should be negative when transitioning to NotTaken state.")
            lq.decision_part.constant +=
                coefficient * change.value_difference
            # Modify variable part
            push!(lq.variable_part.terms,
                  MOI.ScalarAffineTerm(coefficient, change.decision))
            # Modify the coefficient of the mapped variable
            MOI.modify(model, MOI.ObjectiveFunction{F},
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
        # Check if decision is included in the known decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                f.known_decision_terms.affine_terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = f.known_decision_terms.affine_terms[i].coefficient
            # Update fixed decision value
            change.value_difference < 0 ||
                error("Decision value update should be negative when transitioning to NotTaken state.")
            f.known_decision_terms.constant +=
                coefficient * change.value_difference
            # Modify variable part
            add_term!(lq.variable_part.terms,
                      MOI.ScalarAffineTerm(coefficient, change.decision))
            # Get current coefficient in variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                    lq.variable_part.terms), 0)
            coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
            # Modify the coefficient of the mapped variable
            MOI.modify(model, MOI.ObjectiveFunction{F},
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
    end
    if change.new_state == Taken
        # Check if decision is included in the decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.decision_part.terms), 0)
        if !iszero(i)
            # Update fixed decision value
            coefficient = lq.variable_part.terms[i].coefficient
            lq.decision_part.constant +=
                coefficient * change.value_difference
            # Modify variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                lq.variable_part.terms), 0)
            !iszero(i) && deleteat!(lq.variable_part.terms, i)
            # Modify the coefficient of the mapped variable
            MOI.modify(model, MOI.ObjectiveFunction{F},
                       MOI.ScalarCoefficientChange(change.decision, zero(T)))
        end
        # Check if decision is included in the known decision part
        i = something(findfirst(t -> t.variable_index == change.decision,
                                f.known_decision_terms.affine_terms), 0)
        if !iszero(i)
            # Get the decision coefficient
            coefficient = f.known_decision_terms.affine_terms[i].coefficient
            # Update fixed decision value
            f.known_decision_terms.constant +=
                coefficient * change.value_difference
            # Modify variable part
            remove_term!(lq.variable_part.terms,
                         MOI.ScalarAffineTerm(coefficient, change.decision))
            # Get current coefficient in variable part
            i = something(findfirst(t -> t.variable_index == change.decision,
                                    lq.variable_part.terms), 0)
            coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
            # Modify the coefficient of the mapped variable
            MOI.modify(model, MOI.ObjectiveFunction{F},
                       MOI.ScalarCoefficientChange(change.decision, coefficient))
        end
    end
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,QuadraticPart{T}}, ::DecisionStateChange) where T
    # Complete rebridge only option because quadratic coefficients cannot be modified
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,LinearPart{T}}, ::DecisionsStateChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Unbridge the objective function and cache
    # all occuring decisions.
    unbridged_func = zero(QuadraticDecisionFunction{T,LinearPart{T}})
    # Known/variable cross terms
    for (i, term) in enumerate(f.known_variable_terms.quadratic_terms)
        push!(unbridged_func.known_variable_terms.quadratic_terms, copy(term))
    end
    empty!(f.known_variable_terms.quadratic_terms)
    # Known/decision cross terms
    for (i, term) in enumerate(f.known_decision_terms.quadratic_terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_2)
        push!(unbridged_func.known_decision_terms.quadratic_terms,
              MOI.ScalarQuadraticTerm(term.coefficient,
                                      term.variable_index_1,
                                      unbridged.variable))
    end
    # Affine terms
    for (i,term) in enumerate(lq.decision_part.terms)
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index)
        push!(unbridged_func.linear_quadratic_terms.decision_part.terms,
              MOI.ScalarAffineTerm(term.coefficient, unbridged.variable))
    end
    empty!(lq.decision_part.affine_terms)
    # Rebridge objective
    F = QuadraticDecisionFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), unbridged_func)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,QuadraticPart{T}}, ::DecisionsStateChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     lq.decision_part.affine_terms)
        result |= any(t -> t.variable_index_1 == v || t.variable_index_2 == v,
                      lq.decision_part.quadratic_terms)
        result |= any(t -> t.variable_index_1 == v,
                      lq.cross_terms.quadratic_terms)
    end
    # Unbridge the objective function and cache
    # all occuring decisions.
    unbridged_func = QuadraticDecisionFunction(
        QuadraticPart(
            MOIU.filter_variables(v -> !from_decision(v), lq.variable_part),
            convert(MOI.ScalarQuadraticFunction{T}, zero(T)),
            convert(MOI.ScalarQuadraticFunction{T}, zero(T))),
        copy(f.known_part),
        copy(f.known_variable_terms),
        convert(MOI.ScalarQuadraticFunction{T}, zero(T)))
    # Affine terms
    unbridged_affine_terms = map(lq.decision_part.affine_terms) do term
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index)
        MOI.ScalarAffineTerm(term.coefficient, unbridged.variable)
    end
    append!(unbridged_func.linear_quadratic_terms.decision_part.affine_terms, unbridged_affine_terms)
    # Quadratic terms
    unbridged_quad_terms = map(lq.decision_part.quadratic_terms) do term
        unbridged_1 = MOIB.unbridged_variable_function(model, term.variable_index_1)
        unbridged_2 = MOIB.unbridged_variable_function(model, term.variable_index_2)
        MOI.ScalarQuadraticTerm(term.coefficient, unbridged_1.variable, unbridged_2.variable)
    end
    append!(unbridged_func.linear_quadratic_terms.decision_part.quadratic_terms, unbridged_quad_terms)
    # Cross terms
    for term in lq.cross_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(unbridged_func.linear_quadratic_terms.variable_part.affine_terms, term)
    end
    # Add unbridged decision/variable cross terms for rebridge
    unbridged_cross_terms = map(lq.cross_terms.quadratic_terms) do term
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_1)
        MOI.ScalarQuadraticTerm(term.coefficient, unbridged.variable, term.variable_index_2)
    end
    append!(unbridged_func.linear_quadratic_terms.cross_terms.quadratic_terms, unbridged_cross_terms)
    # Known/variable cross terms
    for term in f.known_variable_terms.affine_terms
        # Subtract any added variable/decision cross terms
        remove_term!(unbridged_func.variable_part.affine_terms, term)
    end
    # Known/decision cross terms
    unbridged_known_decision_affine_terms = map(f.known_decision_terms.affine_terms) do term
        # Add unbridged decision/variable cross terms for rebridge
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index)
        MOI.ScalarAffineTerm(term.coefficient, unbridged.variable)
    end
    append!(unbridged_func.known_decision_terms.affine_terms, unbridged_known_decision_affine_terms)
    # Add unbridged known/decision cross terms for rebridge
    unbridged_known_decision_quadratic_terms = map(f.known_decision_terms.quadratic_terms) do term
        unbridged = MOIB.unbridged_variable_function(model, term.variable_index_2)
        MOI.ScalarQuadraticTerm(term.coefficient, term.variable_index_1, unbridged.variable)
    end
    append!(unbridged_func.known_decision_terms.quadratic_terms, unbridged_known_decision_quadratic_terms)
    # Rebridge objective
    F = QuadraticDecisionFunction{T}
    MOI.set(model, MOI.ObjectiveFunction{F}(), unbridged_func)
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, ::KnownValueChange) where T
    # Complete rebridge only option because quadratic coefficients cannot be modified
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,LinearPart{T}}, change::KnownValuesChange) where T
    f = bridge.decision_function
    lq = f.linear_quadratic_terms
    F = MOI.ScalarAffineFunction{T}
    # Recalculate known constant
    known_val = zero(T)
    for term in f.known_part.affine_terms
        known_val += term.coefficient * MOI.get(model, MOI.VariablePrimal(), term.variable_index)
    end
    for term in f.known_part.quadratic_terms
        val_1 = MOI.get(model, MOI.VariablePrimal(), term.variable_index_1)
        val_2 = MOI.get(model, MOI.VariablePrimal(), term.variable_index_2)
        coeff = val_1 == val_2 ? term.coefficient / 2 : term.coefficient
        known_val += coeff * val_1 * val_2
    end
    # Update known part of objective
    f.known_part.constant = known_val
    # Handle cross terms
    for term in f.known_variable_terms.quadratic_terms
        coefficient = term.coefficient
        known_index = term.variable_index_1
        var_index = term.variable_index_2
        i = something(findfirst(t -> t.variable_index == var_index,
                                f.known_variable_terms.affine_terms), 0)
        old_coefficient = iszero(i) ? zero(T) : f.known_variable_terms.affine_terms[i].coefficient
        new_coefficient = coefficient * MOI.get(model, MOI.VariablePrimal(), known_index)
        # Update cross term
        add_term!(f.known_variable_terms.affine_terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, var_index))
        # Update variable part
        add_term!(lq.variable_part.terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, var_index))
        # Get current coefficient in variable part
        i = something(findfirst(t -> t.variable_index == var_index,
                                lq.variable_part.terms), 0)
        coefficient = iszero(i) ? zero(T) : lq.variable_part.terms[i].coefficient
        # Modify the coefficient of the mapped variable
        MOI.modify(model, MOI.ObjectiveFunction{F},
                   MOI.ScalarCoefficientChange(var_index, coefficient))
    end
    for term in f.known_decision_terms.quadratic_terms
        coefficient = term.coefficient
        known_index = term.variable_index_1
        dvar_index = term.variable_index_2
        i = something(findfirst(t -> t.variable_index == dvar_index,
                                f.known_decision_terms.affine_terms), 0)
        old_coefficient = iszero(i) ? zero(T) : f.known_decision_terms.affine_terms[i].coefficient
        new_coefficient = coefficient * MOI.get(model, MOI.VariablePrimal(), known_index)
        # Update cross term
        add_term!(f.known_decision_terms.affine_terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, dvar_index))
        # Update decision part
        add_term!(lq.decision_part.terms,
                  MOI.ScalarAffineTerm(new_coefficient - old_coefficient, dvar_index))
        # Get current coefficient in decision part
        i = something(findfirst(t -> t.variable_index == dvar_index,
                                lq.decision_part.terms), 0)
        coefficient = iszero(i) ? zero(T) : lq.decision_part.terms[i].coefficient
        # Proper decision modification
        unbridged = MOIB.unbridged_variable_function(model, dvar_index)::MOI.SingleVariable
        MOIB.modify_bridged_change(model, MOI.ObjectiveFunction{F},
                                   DecisionCoefficientChange(unbridged.variable, coefficient))
    end
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T,QuadraticPart{T}}, change::KnownValuesChange) where T
    # The objective must be rebridged to complete modification
    F = QuadraticDecisionFunction{T}
    MOI.modify(model, MOI.ObjectiveFunction{F}(), DecisionsStateChange())
    return nothing
end


text:
mutable struct Chord
    words::String
    modulation::Note
    notes::Vector{Note}
    # we need a separate property here for the list model
    notes_model::ListModel
end
export Chord

function Chord(words, modulation, notes)
    Chord(
        words,
        modulation,
        notes,
        property_model(notes, (:numerator, :denominator, :octave, :beats)),
    )
end

function Chord(;
    words = "",
    modulation = Note(),
    # need to allocate a new vector
    notes = Note[],
)
    Chord(words, modulation, notes)
end

const CHORD_REGEX = r"(?<modulation>.*): (?<notes>.*)"

function parse(::Type{Chord}, text::AbstractString; line_number, words = "")
    a_match = match(CHORD_REGEX, text)
    if a_match === nothing
        throw_parse_error(text, "chord", line_number)
    else
        notes_string = a_match["notes"]
        Chord(
            words = words,
            modulation = parse(Note, a_match["modulation"]; line_number = line_number),
            notes = if notes_string == ""
                Note[]
            else
                map(
                    let line_number = line_number
                        (sub_string -> parse(Note, sub_string; line_number = line_number))
                    end,
                    split(a_match["notes"], ", "),
                )
            end,
        )
    end
end

function print(io::IO, chord::Chord)
    words = chord.words
    if words != ""
        print(io, "# ")
        print(io, words)
        println(io)
    end
    print(io, chord.modulation)
    print(io, ": ")
    first_one = true
    for note in chord.notes
        if first_one
            first_one = false
        else
            print(io, ", ")
        end
        print(io, note)
    end
    println(io)
end

# TODO: propertynames?

@inline function getproperty(chord::Chord, property_name::Symbol)
    if property_name === :numerator ||
       property_name === :denominator ||
       property_name === :denominator ||
       property_name === :octave ||
       property_name === :beats
        getproperty(chord.modulation, property_name)
    else
        getfield(chord, property_name)
    end
end

@inline function setproperty!(chord::Chord, property_name::Symbol, value)
    if property_name === :numerator ||
       property_name === :denominator ||
       property_name === :octave ||
       property_name === :beats
        setproperty!(chord.modulation, property_name, value)
    else
        setfield!(
            chord,
            property_name,
            convert(fieldtype(typeof(chord), property_name), value),
        )
    end
end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
using Random, Statistics

function execute_abrkga(instance::OBOPDataset;seed::Int64=nothing,local_search=true,debug=false,ls_type=1)

    if !isnothing(seed)
        Random.seed!(seed)
    end

    # Start timer
    time_total = time()

    # Initialize best solution
    best_solution = OBOPSolution(instance.total_itens,time_total)

    # Initialize parameters
    parameters = OBOPParameters(0.999)

    println("0.1.14")

    # Extended chromosome size
    chromosome_size_extended = instance.total_itens + 1

    # Population current
    pop_keys    = rand(chromosome_size_extended, parameters.population_size)
    pop_buckets = zeros(Int64, instance.total_itens, parameters.population_size)
    pop_fitness = zeros(Int64, parameters.population_size)

    # Decode and evaluate population fitness
    #Threads.@threads
    for item in 1:parameters.population_size
        decoder!(item, pop_buckets, pop_keys)
        objective!(item, pop_buckets, pop_fitness, instance)

        if pop_fitness[item] > best_solution.objective
            best_solution.id = item
            best_solution.total_time = time() - best_solution.start_time
            best_solution.generation = 0
            if debug
                println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
            end
            best_solution.objective = pop_fitness[item]
        end
    end

    # Best objective
    best_solution.bucket = deepcopy(pop_buckets[:,best_solution.id])

    # Sort pop by fitness
    index_order = sortperm(pop_fitness[1:parameters.population_size], rev=true)    # 999  akslaks

    # Get partial time
    time_partial = time() - time_total

    # Statistics
    statistics = OBOPStatistics([mean(pop_fitness[ index_order[1:parameters.elite_size] ])],
                                [ std(pop_fitness[ index_order[1:parameters.elite_size] ])],
                                [mean(pop_fitness)],
                                [ std(pop_fitness)],
                                [best_solution.objective],
                                [time_partial],
                                time_partial,
                                best_solution.total_time,
                                0,
                                0)

    # Temp population
    pop_keys_new = zeros(chromosome_size_extended, parameters.population_size)

    # Variables used inside main loop
    elite_id = non_elite_id = 0
    rhoe = 0.0
    offspring_begin = offspring_end = 0
    inicio_mutante = fim_mutante = 0

    # ========= #
    # Main loop #
    # ========= #
    generation::Int64 = 1

    while generation <= parameters.max_generations

        time_partial = time()

        # ========================== #
        # Start evolutionary process #
        # ========================== #

        # ============= #
        # Create elites #
        # ============= #
        for item in 1:parameters.elite_size
            pop_keys_new[:,item] = pop_keys[:,index_order[item]]
            #pop_fitness[item] = pop_fitness[index_order[item]]
            #pop_buckets[:,item] = pop_buckets[:,index_order[item]]

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end

        end

        # ================= #
        # Create offsprings #
        # ================= #
        offspring_begin = parameters.elite_size + 1
        offspring_end = parameters.population_size - parameters.mutation_size

        for item in offspring_begin:offspring_end
            # Select an elite parent
            elite_id = rand(collect(1:parameters.elite_size))
            elite_id = index_order[elite_id]

            # Select an non elite parent
            #non_elite_id = rand(collect(offspring_begin:offspring_end))                # don't include mutants
            non_elite_id = rand(collect(offspring_begin:parameters.population_size))    # including mutants
            non_elite_id = index_order[non_elite_id]

            # Get rhoe value from non elite last key
            rhoe = 0.65 + (pop_keys[chromosome_size_extended,non_elite_id]) * (0.80 - 0.65)

            # Mate: including rhoe
            for j in 1:chromosome_size_extended
                # copy allele of top chromosome to the new generation or not
                if rand() < rhoe
                    pop_keys_new[j,item] = pop_keys[j,elite_id]
                else
                    pop_keys_new[j,item] = pop_keys[j,non_elite_id]
                end
            end

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end
        end

        # ============== #
        # Create Mutants #
        # ============== #
        inicio_mutante = offspring_end + 1
        fim_mutante = parameters.population_size
        for item in inicio_mutante:fim_mutante
            pop_keys_new[:,item] = rand(chromosome_size_extended)

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end
        end

        # =================================== #
        # Set temp population to the original #
        # =================================== #
        pop_keys = pop_keys_new

        # =================== #
        # Sort pop by fitness #
        # =================== #
        index_order = sortperm(pop_fitness[1:parameters.population_size], rev=true)

        # ============ #
        # Local search #
        # ============ #
        if local_search
            clustering_search(parameters.elite_size,
                              index_order,
                              pop_keys,
                              pop_buckets,
                              pop_fitness,
                              instance,
                              best_solution,
                              generation,
                              statistics,ls_type)
        end

        # ================= #
        # Update parameters #
        # ================= #
        update_parameters!(generation, parameters)

        # ================================== #
        # Update partial time and statistics #
        # ================================== #
        time_partial = time() - time_partial

        push!(statistics.elite_mean, mean(pop_fitness[ index_order[1:parameters.elite_size] ]))
        push!(statistics.elite_std, std(pop_fitness[ index_order[1:parameters.elite_size] ]))
        push!(statistics.pop_mean, mean(pop_fitness))
        push!(statistics.pop_std, std(pop_fitness))
        push!(statistics.best_generations, best_solution.objective)
        push!(statistics.time_generations, time_partial)

        generation += 1
    end

    statistics.time_total = time() - time_total
    statistics.time_best = best_solution.total_time

    return best_solution.objective, create_bucket_order(best_solution.bucket), statistics
end

function execute_abrkga_vnd(instance::OBOPDataset;seed::Int64=nothing,local_search=true,debug=false,ls_type=1)

    if !isnothing(seed)
        Random.seed!(seed)
    end

    # Start timer
    time_total = time()

    # Initialize best solution
    best_solution = OBOPSolution(instance.total_itens,time_total)

    # Initialize parameters
    parameters = OBOPParameters(0.999)

    println("0.1.14")

    # Extended chromosome size
    chromosome_size_extended = instance.total_itens + 1

    # Population current
    pop_keys    = rand(chromosome_size_extended, parameters.population_size)
    pop_buckets = zeros(Int64, instance.total_itens, parameters.population_size)
    pop_fitness = zeros(Int64, parameters.population_size)

    # Decode and evaluate population fitness
    #Threads.@threads
    for item in 1:parameters.population_size
        decoder!(item, pop_buckets, pop_keys)
        objective!(item, pop_buckets, pop_fitness, instance)

        if pop_fitness[item] > best_solution.objective
            best_solution.id = item
            best_solution.total_time = time() - best_solution.start_time
            best_solution.generation = 0
            if debug
                println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
            end
            best_solution.objective = pop_fitness[item]
        end
    end

    # Best objective
    best_solution.bucket = deepcopy(pop_buckets[:,best_solution.id])

    # Sort pop by fitness
    index_order = sortperm(pop_fitness[1:parameters.population_size], rev=true)    # 999  akslaks

    # Get partial time
    time_partial = time() - time_total

    # Statistics
    statistics = OBOPStatistics([mean(pop_fitness[ index_order[1:parameters.elite_size] ])],
                                [ std(pop_fitness[ index_order[1:parameters.elite_size] ])],
                                [mean(pop_fitness)],
                                [ std(pop_fitness)],
                                [best_solution.objective],
                                [time_partial],
                                time_partial,
                                best_solution.total_time,
                                0,
                                0)

    # Temp population
    pop_keys_new = zeros(chromosome_size_extended, parameters.population_size)

    # Variables used inside main loop
    elite_id = non_elite_id = 0
    rhoe = 0.0
    offspring_begin = offspring_end = 0
    inicio_mutante = fim_mutante = 0

    return best_solution.objective, create_bucket_order(best_solution.bucket), statistics

    # ========= #
    # Main loop #
    # ========= #
    generation::Int64 = 1

    while generation <= parameters.max_generations

        time_partial = time()

        # ========================== #
        # Start evolutionary process #
        # ========================== #

        # ============= #
        # Create elites #
        # ============= #
        for item in 1:parameters.elite_size
            pop_keys_new[:,item] = pop_keys[:,index_order[item]]
            #pop_fitness[item] = pop_fitness[index_order[item]]
            #pop_buckets[:,item] = pop_buckets[:,index_order[item]]

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end

        end

        # ================= #
        # Create offsprings #
        # ================= #
        offspring_begin = parameters.elite_size + 1
        offspring_end = parameters.population_size - parameters.mutation_size

        for item in offspring_begin:offspring_end
            # Select an elite parent
            elite_id = rand(collect(1:parameters.elite_size))
            elite_id = index_order[elite_id]

            # Select an non elite parent
            #non_elite_id = rand(collect(offspring_begin:offspring_end))                # don't include mutants
            non_elite_id = rand(collect(offspring_begin:parameters.population_size))    # including mutants
            non_elite_id = index_order[non_elite_id]

            # Get rhoe value from non elite last key
            rhoe = 0.65 + (pop_keys[chromosome_size_extended,non_elite_id]) * (0.80 - 0.65)

            # Mate: including rhoe
            for j in 1:chromosome_size_extended
                # copy allele of top chromosome to the new generation or not
                if rand() < rhoe
                    pop_keys_new[j,item] = pop_keys[j,elite_id]
                else
                    pop_keys_new[j,item] = pop_keys[j,non_elite_id]
                end
            end

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end
        end

        # ============== #
        # Create Mutants #
        # ============== #
        inicio_mutante = offspring_end + 1
        fim_mutante = parameters.population_size
        for item in inicio_mutante:fim_mutante
            pop_keys_new[:,item] = rand(chromosome_size_extended)

            decoder!(item, pop_buckets, pop_keys_new)
            objective!(item, pop_buckets, pop_fitness, instance)

            if pop_fitness[item] > best_solution.objective
                best_solution.id = item
                best_solution.total_time = time() - best_solution.start_time
                best_solution.generation = generation
                best_solution.local_search = false
                if debug
                    println("Generation 0 $(pop_fitness[item]) > $(best_solution.objective) $(best_solution.total_time)")
                end
                best_solution.objective = pop_fitness[item]
            end
        end

        # =================================== #
        # Set temp population to the original #
        # =================================== #
        pop_keys = pop_keys_new

        # =================== #
        # Sort pop by fitness #
        # =================== #
        index_order = sortperm(pop_fitness[1:parameters.population_size], rev=true)

        # ============ #
        # Local search #
        # ============ #
        if local_search
            clustering_search_new(parameters.elite_size,
                              index_order,
                              pop_keys,
                              pop_buckets,
                              pop_fitness,
                              instance,
                              best_solution,
                              generation,
                              statistics)
        end

        # ================= #
        # Update parameters #
        # ================= #
        update_parameters!(generation, parameters)

        # ================================== #
        # Update partial time and statistics #
        # ================================== #
        time_partial = time() - time_partial

        push!(statistics.elite_mean, mean(pop_fitness[ index_order[1:parameters.elite_size] ]))
        push!(statistics.elite_std, std(pop_fitness[ index_order[1:parameters.elite_size] ]))
        push!(statistics.pop_mean, mean(pop_fitness))
        push!(statistics.pop_std, std(pop_fitness))
        push!(statistics.best_generations, best_solution.objective)
        push!(statistics.time_generations, time_partial)

        # Stop if elite std is low
        if statistics.elite_std[end] < 0.0001
            break
        end

        generation += 1

    end

    statistics.time_total = time() - time_total
    statistics.time_best = best_solution.total_time

    return best_solution.objective, create_bucket_order(best_solution.bucket), statistics
end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
using ReverseDiffExpressions: Model, Variable, Func, getvar!, returns!, uses!, onevar, targetvar, addfunc!, isref
using LoopVectorization: Instruction

read_model(q::Expr, mod::Module) = read_model!(Model(mod), q)

function read_model!(m::Model, q::Expr)
    for arg ∈ q.args
        arg isa Expr || continue
        ex = arg::Expr
        if ex.head === :for
            add_loopset!(m, q)
        elseif ex.head === :block
            read_model!(m, ex)
        else
            read_line!(m, ex)
        # else
            # throw("Don't know how to handle block $ex")
        end
    end
    for (varid, v) ∈ enumerate(m.vars)
        vid = v.varid
        # @show v
        @assert varid == vid + 1
        # if ReverseDiffExpressions.isinitialized(v) && !isref(v)
        #     # If it must already be initialized, yet it isn't a ref to a constant
        #     push!(m.inputvars, vid)
        # end
    end
    m
end

function add_loopset!(m::Model, ex)
    throw("Not yet implemented")
    ls = LoopSet(stage1distsub!(ex), Symbol(m.mod))
end

function read_line!(m::Model, ex::Expr)
    if ex.head === :(=)
        LHS = first(ex.args)::Union{Symbol,Expr}
        # If LHS isa Expr, we assume it is a tuple being unpacked
        RHS = ex.args[2]
        read_call!(m, RHS, LHS)
    elseif ex.head === :(.=)
        read_broadcast!(m, ex)
    elseif ex.head === :call
        read_call!(m, ex)
    end
end

# Reads function arguments, returning a Variable
read_argument!(m::Model, s::Symbol)::Variable = getvar!(m, s)
function read_argument!(m::Model, ex::Expr)::Variable
    if ex.head === :call
        # To unnest the expression, we assign to LHS
        read_call!(m, ex, gensym(:LHS))
    elseif ex.head === :(.)
        read_broadcast_tempalloc!(m, ex, gensym(:LHS))
    elseif ex.head === :ref
        read_ref!(m, ex, gensym(:LHS))
    elseif ex.head === Symbol("'")
        read_call!(m, Instruction(:adjoint), ex.args, gensym(:LHS))
    else
        println(ex)
        @show ex.head === Symbol("'")
        @show ex.head, Symbol("'")
        throw("Expression not found.")
    end
end
function read_argument!(m::Model, x)::Variable
    xv = getvar!(m, gensym())
    xv.ref = x
    # xv.initialized = true
    xv
end
ReverseDiffExpressions.uses!(func::Func, m::Model, x) = uses!(func, read_argument!(m, x))

function read_ref!(m::Model, ex::Expr, LHS)
    retv = getvar!(m, LHS)
    func = Func(Instruction(:Base,:view), false)
    returns!(func, retv)
    foreach(arg -> uses!(func, m, arg), ex.args)
    addfunc!(m, func)
end

function read_broadcast!(m::Model, ex::Expr)
    read_broadcast!(m, first(Meta.lower(m.mod, ex).args)::Core.CodeInfo)
end
function read_broadcast!(m::Model, ci::Core.CodeInfo)

end
instr_from_expr(ex::Expr) = Instruction(first(ex.args)::Union{Symbol,Expr})::Instruction
function LoopVectorization.Instruction(ex::Expr)
    ex.head === :(.) || throw("Could not parse instruction $ex.")
    Instruction((ex.args[1])::Symbol, (((ex.args[2])::QuoteNode).value)::Symbol)
end

function read_sampling_statement!(m::Model, f::Instruction, ex::Expr)
    arg1 = (ex.args[2])::Symbol
    call = (ex.args[3])::Expr
    @assert call.head === :call
    v = getvar!(m, arg1)
    f = instr_from_expr(call)
    scale = if f == Instruction(:^) # Then it is weighted
        call = (call.args[2])::Expr
        f = instr_from_expr(call)
        read_argument!(m, call.args[3])
    else
        onevar(m)
    end
    func = Func(f, true)
    target = targetvar(m)
    returns!(func, target)
    uses!(func, v)
    for i ∈ 2:length(call.args)
        uses!(func, m, call.args[i])
    end
    uses!(func, scale)
    addfunc!(m, func)
end

# Reads a call without a return
function read_call!(m::Model, ex::Expr, ::Nothing = nothing)
    f = instr_from_expr(ex)
    if f.instr === :~
        return read_sampling_statement!(m, f, ex)
    else
        @show ex
        throw("Currently only `~` and `=` type calls are supported.")
    end
end
# Reads a call with a return
function read_call!(m::Model, ex::Expr, LHS::Symbol)
    broadcasts = (:(.+), :(.*), :(.-), :(./), :(.*ˡ), :(.÷))
    f = instr_from_expr(ex)
    # @show ex
    if f.instr ∈ broadcasts
        return read_broadcast!(m, ex, LHS)
    # elseif f.instr === :~
        # return read_sampling_statement!(m, ex)
    elseif f.instr === :(:)
        return read_range_expr!(m, ex, LHS)
    end
    read_call!(m, f, @view(ex.args[2:end]), LHS)
end
function read_call!(m::Model, f::Instruction, args, LHS)
    func = Func(f, false)
    retv = getvar!(m, LHS)
    returns!(func, retv)
    foreach(arg -> uses!(func, m, arg), args)
    addfunc!(m, func)
end

function read_range_expr!(m::Model, ex::Expr, LHS)
    if length(ex.args) == 3
        read_range_args!(m, ex.args[2], ex.args[3], LHS)
    else
        @assert length(ex.args) == 4
        read_range_args!(m, ex.args[2], ex.args[3], ex.args[4], LHS)
    end
end
function read_range_args!(m::Model, l::Number, u::Number, LHS::Symbol)
    retv = getvar!(m, LHS)
    retv.ref = Expr(:call, Expr(:curly, :StaticUnitRange, l, u))
    # retv.initialized = true
    retv
end
function read_range_args!(m::Model, l, u, LHS::Symbol)
    lin = l isa Number; uin = u isa Number
    f = if lin
        :StaticUpperUnitRange
    elseif uin
        :StaticLowerUnitRange
    else
        :(:)
    end
    func = Func(Instruction(f), false)
    retv = getvar!(m, LHS)
    returns!(func, retv)
    uses!(func, m, l)
    uses!(func, m, u)
    addfunc!(m, func)
end
# function read_range_args!(m::Model, l, s, u, LHS::Symbol)
# end


text:
# # Noncommutative variables

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Noncommutative and Hermitian/noncommutative_variables.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Noncommutative and Hermitian/noncommutative_variables.ipynb)
# **Adapted from**: Examples 2.11 and 2.2 of [BKP16]
#
# [BKP16] Sabine Burgdorf, Igor Klep, and Janez Povh.
# *Optimization of polynomials in non-commuting variables*.
# Berlin: Springer, 2016.

# ## Example 2.11
#
# We consider the Example 2.11 of [BKP16] in which the polynomial with noncommutative variables
# $(x * y + x^2)^2 = x^4 + x^3y + xyx^2 + xyxy$ is tested to be sum-of-squares.

using Test #src
using DynamicPolynomials
@ncpolyvar x y
p = (x * y + x^2)^2

# We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.

using SumOfSquares
import CSDP
optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = Model(optimizer_constructor)
con_ref = @constraint(model, p in SOSCone())

optimize!(model)

# We see that both the monomials `xy` and `yx` are considered separately, this is a difference with the commutative version.

certificate_basis(con_ref)

# We see that the solution correctly uses the monomial `xy` instead of `yx`. We also identify that only the monomials `x^2` and `xy` would be needed. This would be dectected by the Newton chip method of [Section 2.3, BKP16].

gram_matrix(con_ref).Q

# When asking for the SOS decomposition, the numerically small entries makes the solution less readable.

@test length(sos_decomposition(con_ref).ps) == 2 #src
sos_decomposition(con_ref)

# They are however easily discarded by using a nonzero tolerance:

dec = sos_decomposition(con_ref, 1e-6)                    #src
@test length(dec.ps) == 1                                 #src
@test sign(first(coefficients(dec.ps[1]))) * dec.ps[1] ≈ x * y + x^2 rtol=1e-5 atol=1e-5 #src
sos_decomposition(con_ref, 1e-6)

# ## Example 2.2
#
# We consider now the Example 2.2 of [BKP16] in which the polynomial with noncommutative variables
# $(x + x^{10}y^{20}x^{10})^2$ is tested to be sum-of-squares.

using DynamicPolynomials
@ncpolyvar x y
n = 10
p = (x + x^n * y^(2n) * x^n)^2

using SumOfSquares
model = Model(optimizer_constructor)
con_ref = @constraint(model, p in SOSCone())

optimize!(model)

# Only two monomials were considered for the basis of the gram matrix thanks to the Augmented Newton chip method detailed in [Section 2.4, BKP16].

certificate_basis(con_ref)

gram_matrix(con_ref).Q

sos_decomposition(con_ref, 1e-6)

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
"""
```
dsgevecm_likelihood(m::AbstractDSGEVECMModel{S}, data::Matrix{S};
    apply_altpolicy::Bool = false) where {S<:Real}
```
evaluates the likelihood of a DSGE-VECM. It is assumed that
`get_λ(m)` retrieves the prior weight `λ` on the DSGE.

The matrix `data` is assumed an nobs x T+lags matrix,
where the lags indicate we cut off the data for presampling.

In the future, this function may be combined with `dsgevar_likelihood`
because the two functions are almost exactly the same. We currently have
them implemented separately so that the naming convention does not
create any confusion. One can think of DSGE-VECMs as DSGE-VARs
with an additional regressor that corrects for errors to account for
cointegration, hence it would not be "wrong" per se to add
a couple `if-else` conditions inside `dsgevar_likelihood` and let this
function cover both DSGE-VARs and DSGE-VECMS. But the current decision
to have separate functions makes it clear that the `dsgevar_likelihood`
specifically operates on DSGE-VARs without an error correction term, and
`dsgevecm_likelihood` operates on DSGE-VECMs specifically.
"""
function dsgevecm_likelihood(m::AbstractDSGEVECMModel{S}, data::Matrix{S};
                            apply_altpolicy::Bool = false,
                            use_intercept::Bool = true) where {S<:Real}

    # Set up. Default to non-regime switching if these settings don't exist.
    lags = n_lags(m)
    n_regimes        = haskey(get_settings(m), :n_regimes) ?
        get_setting(m, :n_regimes) : 1
    regime_switching = haskey(get_settings(m), :regime_switching) && n_regimes > 1 ?
        get_setting(m, :regime_switching) : false # if n_regimes == 1, then no switching needed

    # Get population moments
    YYYY, XXYY, XXXX = compute_var_population_moments(data, lags; use_intercept = use_intercept)

    if regime_switching
        error("Regime switching has not been implemented yet.")
        yyyyd = OrderedDict{Int, Matrix{S}}()
        xxyyd = OrderedDict{Int, Matrix{S}}()
        xxxxd = OrderedDict{Int, Matrix{S}}()

        for reg = 1:n_regimes
            yyyyd[reg], xxyyd[reg], xxxxd[reg] =
                compute_system(m; apply_altpolicy = apply_altpolicy,
                               regime_switching = regime_switching,
                               regime = reg, n_regimes = n_regimes,
                               get_population_moments = true,
                               use_intercept = use_intercept)
        end
    else
        yyyyd, xxyyd, xxxxd = compute_system(m; apply_altpolicy = apply_altpolicy,
                                             get_population_moments = true,
                                             use_intercept = use_intercept)
    end

    return dsgevecm_likelihood(YYYY, XXYY, XXXX, yyyyd, xxyyd, xxxxd,
                      size(data, 2) - lags, get_λ(m); n_cointadd = n_cointegrating_add(m))
end

"""
```
dsgevecm_likelihood(YYYY::Matrix{S}, XXYY::Matrix{S},
    XXXX::Matrix{S}, YYYYD::Matrix{S}, XXYYD::Matrix{S},
    XXXXD::Matrix{S}, T::Int, λ::S,
    lags::Int, n_obs::Int,
    coint_inds::Union{Vector{Int}, UnitRange{Int}}) where {S<:Real}
```
evaluates the likelihood of a DSGE-VECM given the
population moments of the raw data (`YYYY`, `XXYY`, `XXXX`)
and the population moments implied by the DSGE prior (`YYYYD`, `XXYYD`, `XXXXD`).

### Other Inputs
* `T`: number of time periods in data
* `λ`: weight placed on the DSGE prior
* `lags`: number of lags in the VECM
* `n_obs`: number of observables (e.g. number of time series).
* `coint_inds`: indices of the cointegrated observables
"""
function dsgevecm_likelihood(YYYY::Matrix{S}, XXYY::Matrix{S},
                             XXXX::Matrix{S}, YYYYD::Matrix{S}, XXYYD::Matrix{S},
                             XXXXD::Matrix{S}, T::Int, λ::S;
                             n_cointadd::Int = 0) where {S <: Real}

    if isinf(λ)
        if n_cointadd > 0
            inds = 1 + n_cointadd:size(XXXX, 1)
            XXXXD = XXXXD[inds, inds]
            XXYYD = XXYYD[inds, :]
            XXXX  = XXXXD[inds, inds]
            XXYY  = XXYY[inds, :]
        end
        β = XXXXD \ XXYYD
        Σ = YYYYD - XXYYD' * β
        n_obs  = size(YYYY, 1)
        loglik = -(n_obs * T / 2) * log(2 * π) - (T / 2) * log(det(Σ)) -
            .5 * sum(diag(Σ \ (YYYY - β' * XXYY - XXYY' * β + β' * XXXX * β)))
        return real(loglik)
    else
        # Compute DSGEVECM components for weighted population moments
        if n_cointadd > 0
            Λ⁻¹ = Diagonal{S}(I, size(XXXX, 1))    # these lines are the only ones which distinguish
            Λ⁻¹[1:n_cointadd, 1:n_cointadd] ./= λT # this likelihood from dsgevar_likelihood
            XXYY = Λ⁻¹ * XXYY
            XXXX = Λ⁻¹ * XXXX * Λ⁻¹
        end

        return dsgevar_likelihood(YYYY, XXYY, XXXX, YYYYD, XXYYD, XXXXD, T, λ)
    end
end


text:
"""
    GB2Luxor(p)

Transforms a `Point` defined in `GeometryBasics` into one defined
in `Luxor`.
"""
GB2Luxor(p) = Point(p...)

"""
    draw_node(args; center, radius, action)

Function used to draw a circle. Very similar to `Javis.JCircle` with the only
difference that it does not use `sethue` so animating a color change is easier.
"""
function draw_node(args...; center, radius, action)
    (args...; center = center, radius = radius, action = action) -> begin
        circle(center, radius, action)
        return center
    end
end

"""
    _adapt_to_lims(points, width, height)

Rescales a vector of points to fit in height and width to keep it simple
it uses the smallest among x and y needed variation to rescale both.
"""
function _adapt_to_lims(points, width, height)
    xs = [p[1] for p in points]
    ys = [p[2] for p in points]

    xmin, xmax = abs(minimum(xs)), abs(maximum(xs))
    ymin, ymax = abs(minimum(ys)), abs(maximum(ys))

    xscaling = min(width / 2xmin, width / 2xmax)
    yscaling = min(height / 2ymin, height / 2ymax)

    scaling = min(xscaling, yscaling)

    return [Point(p[1] * scaling, p[2] * scaling) for p in points]
end