text:
export SplitExplicitLSRK2nSolver

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using StaticArrays
using ...SystemSolvers
using ...MPIStateArrays: array_device, realview
using ...GenericCallbacks

using ...ODESolvers:
    AbstractODESolver, LowStorageRungeKutta2N, update!, updatedt!, getdt
import ...ODESolvers: dostep!

using ...BalanceLaws:
    #   initialize_fast_state!,
    #   initialize_adjustment!,
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

LSRK2N = LowStorageRungeKutta2N

@doc """
    SplitExplicitLSRK2nSolver(slow_solver, fast_solver; dt, t0 = 0, coupled = true)

This is a time stepping object for explicitly time stepping the differential
equation given by the right-hand-side function `f` with the state `Q`, i.e.,

```math
  \\dot{Q_fast} = f_fast(Q_fast, Q_slow, t)
  \\dot{Q_slow} = f_slow(Q_slow, Q_fast, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

This method performs an operator splitting to timestep the Sea-Surface elevation
and vertically averaged horizontal velocity of the model at a faster rate than
the full model, using LowStorageRungeKutta2N time-stepping.

""" SplitExplicitLSRK2nSolver
mutable struct SplitExplicitLSRK2nSolver{SS, FS, RT, MSA} <: AbstractODESolver
    "slow solver"
    slow_solver::SS
    "fast solver"
    fast_solver::FS
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "storage for transfer tendency"
    dQ2fast::MSA

    function SplitExplicitLSRK2nSolver(
        slow_solver::LSRK2N,
        fast_solver,
        Q = nothing;
        dt = getdt(slow_solver),
        t0 = slow_solver.t,
    ) where {AT <: AbstractArray}
        SS = typeof(slow_solver)
        FS = typeof(fast_solver)
        RT = real(eltype(slow_solver.dQ))

        dQ2fast = similar(slow_solver.dQ)
        dQ2fast .= -0.0
        MSA = typeof(dQ2fast)
        return new{SS, FS, RT, MSA}(
            slow_solver,
            fast_solver,
            RT(dt),
            RT(t0),
            0,
            dQ2fast,
        )
    end
end

function dostep!(
    Qvec,
    split::SplitExplicitLSRK2nSolver{SS},
    param,
    time::Real,
) where {SS <: LSRK2N}
    slow = split.slow_solver
    fast = split.fast_solver

    Qslow = Qvec.slow
    Qfast = Qvec.fast

    dQslow = slow.dQ
    dQ2fast = split.dQ2fast

    slow_bl = slow.rhs!.balance_law
    fast_bl = fast.rhs!.balance_law

    groupsize = 256

    slow_dt = getdt(slow)
    fast_dt_in = getdt(fast)

    for slow_s in 1:length(slow.RKA)
        # Current slow state time
        slow_stage_time = time + slow.RKC[slow_s] * slow_dt

        # Fractional time for slow stage
        if slow_s == length(slow.RKA)
            fract_dt = (1 - slow.RKC[slow_s]) * slow_dt
        else
            fract_dt = (slow.RKC[slow_s + 1] - slow.RKC[slow_s]) * slow_dt
        end

        # Initialize fast model and set time-step and number of substeps we need
        # Note: to reproduce previous Fast output, 1) remove "firstStage = ..." line
        fast_steps = [0 0 0]
        FT = typeof(slow_dt)
        fast_time_rec = [fast_dt_in FT(0)]
        initialize_fast_state!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            fract_dt,
            fast_time_rec,
            fast_steps;
            firstStage = (slow_s == 1),
        )
        # Initialize tentency adjustment before evaluation of slow mode
        initialize_adjustment!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
        )

        # Evaluate the slow mode
        # --> save tendency for the fast
        slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false)

        # vertically integrate slow tendency to advance fast equation
        # and use vertical mean for slow model (negative source)
        # ---> work with dQ2fast as input
        tendency_from_slow_to_fast!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            dQ2fast,
        )

        # Compute (and RK update) slow tendency
        slow.rhs!(dQslow, Qslow, param, slow_stage_time, increment = true)

        # Update (RK-stage) slow state
        event = Event(array_device(Qslow))
        event = update!(array_device(Qslow), groupsize)(
            realview(dQslow),
            realview(Qslow),
            slow.RKA[slow_s % length(slow.RKA) + 1],
            slow.RKB[slow_s],
            slow_dt,
            nothing,
            nothing,
            nothing;
            ndrange = length(realview(Qslow)),
            dependencies = (event,),
        )
        wait(array_device(Qslow), event)

        # Determine number of substeps we need
        fast_dt = fast_time_rec[1]
        nsubsteps = fast_steps[3]

        updatedt!(fast, fast_dt)

        for substep in 1:nsubsteps
            fast_time = slow_stage_time + (substep - 1) * fast_dt
            dostep!(Qfast, fast, param, fast_time)

            # cumulate fast solution
            cummulate_fast_solution!(
                fast_bl,
                fast.rhs!,
                Qfast,
                fast_time,
                fast_dt,
                substep,
                fast_steps,
                fast_time_rec,
            )
        end

        # Reconcile slow equation using fast equation
        # Note: to reproduce previous Fast output, 2) remove "lastStage = ..." line
        reconcile_from_fast_to_slow!(
            slow_bl,
            fast_bl,
            slow.rhs!,
            fast.rhs!,
            Qslow,
            Qfast,
            fast_time_rec;
            lastStage = (slow_s == length(slow.RKA)),
        )

    end

    # reset fast time-step to original value
    updatedt!(fast, fast_dt_in)

    # now do implicit mixing step
    nImplSteps = slow_bl.numImplSteps
    if nImplSteps > 0
        # 1. get implicit mising model, model state variable array and solver handles
        ivdc_dg = slow.rhs!.modeldata.ivdc_dg
        ivdc_bl = ivdc_dg.balance_law
        ivdc_Q = slow.rhs!.modeldata.ivdc_Q
        ivdc_solver = slow.rhs!.modeldata.ivdc_bgm_solver
        # ivdc_solver_dt = getdt(ivdc_solver) # would work if solver time-step was set
        # FT = typeof(slow_dt)
        # ivdc_solver_dt = slow_dt / FT(nImplSteps) # just recompute time-step
        ivdc_solver_dt = ivdc_bl.parent_om.ivdc_dt
        # println("ivdc_solver_dt = ",ivdc_solver_dt )
        # 2. setup start RHS, initial guess and values for computing mixing coeff
        ivdc_Q.θ .= Qslow.θ
        ivdc_RHS = slow.rhs!.modeldata.ivdc_RHS
        ivdc_RHS.θ .= Qslow.θ
        ivdc_RHS.θ .= ivdc_RHS.θ ./ ivdc_solver_dt
        ivdc_dg.state_auxiliary.θ_init .= ivdc_Q.θ
        # 3. Invoke iterative solver

        println("BEFORE maximum(ivdc_Q.θ[:]): ", maximum(ivdc_Q.realdata[:]))
        println("BEFORE minimum(ivdc_Q.θ[:]): ", minimum(ivdc_Q.realdata[:]))

        lm!(y, x) = ivdc_dg(y, x, nothing, 0; increment = false)
        solve_tot = 0
        iter_tot = 0
        for i in 1:nImplSteps
            solve_time = @elapsed iters =
                linearsolve!(lm!, nothing, ivdc_solver, ivdc_Q, ivdc_RHS)
            solve_tot = solve_tot + solve_time
            iter_tot = iter_tot + iters
            # Set new RHS and initial values
            ivdc_RHS.θ .= ivdc_Q.θ ./ ivdc_solver_dt
            ivdc_dg.state_auxiliary.θ_init .= ivdc_Q.θ
        end
        println("solver iters, time: ", iter_tot, ", ", solve_tot)

        println("AFTER  maximum(ivdc_Q.θ[:]): ", maximum(ivdc_Q.realdata[:]))
        println("AFTER  minimum(ivdc_Q.θ[:]): ", minimum(ivdc_Q.realdata[:]))

        # exit()
        # Now update
        Qslow.θ .= ivdc_Q.θ

    end

    return nothing
end


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
theme(:wong)
gr(aspect_ratio=1,legend=:bottomright,axis=nothing,border=:none,ticks=nothing,legendfontsize=14,xlims=[1.5,4],ylims=[0.8,3],
   markerstrokewidth=0,markersize=1)
r1D, w1D = gauss_quad(0,0,2)
Nfp = length(r1D) # number of points per face
e = ones(Nfp) # vector of all ones
z = zeros(Nfp) # vector of all zeros
rf = [r1D; -r1D; -e];
sf = [-e; r1D; -r1D];

r_tri,s_tri,_ = Basis2DTri.quad_nodes_tri(4)
plot([-1;-1].+3,[-1;1].+2,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3,[1;-1].+2,seriescolor=:black,lw=3,label = nothing)
plot!([-1;1].+3,[-1;-1].+2,seriescolor=:black,lw=3,label = nothing)
#scatter!(r_tri.+3,s_tri.+2,markersize=8,markercolor=:deepskyblue,markerstrokecolor=:black,markerstrokewidth=4,label = nothing)
scatter!(rf.+3,sf.+2,markersize=8,markershape=:rect,markercolor=:darkorange1,markerstrokecolor=:black,markerstrokewidth=4,label = nothing)

savefig("triangle_facequad_nodes.png")

ext: jl

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
function Domain(
    boundary::Boundary{Interval},
    lattice::Lattice{1},
    dielectric::TDF = DielectricFunction(1),
    pump::TPF = PumpFunction(0),
    dispersion::Union{Tuple,AbstractDispersion} = ();
    type::Symbol = :generic,
    name::Symbol = :anonymous,
    fit::Bool = false
    ) where {TDF<:DielectricFunction, TPF<:PumpFunction}

    bnd = boundary
    lat = lattice
    if fit
        nx = ceil(Int,bnd.shape.a/lat.dx)
        dx = bnd.shape.a/nx
        lat = Lattice(dx;origin=dx/2)
    end

    i1 = latticeindex(lat,bnd.shape.start)
    i2 = latticeindex(lat,bnd.shape.stop)
    imin,imax = floor(Int,min(i1,i2)), floor(Int,max(i1,i2))
    x = lat[imin:imax]
    inds_keep = bnd.shape.(x)
    x = x[inds_keep]
    indices = CartesianIndex.(imin:imax)[inds_keep]
    ε = dielectric.(x)
    F = pump.(x)
    interior = trues(size(x))
    surface = falses(size(x))
    surface[1] = surface[end] = true
    bulk = interior .& .!surface
    corner = falses(size(x))
    nnm = (vcat(0,1:length(x)-1),)
    nnp = (vcat(2:length(x),0),)

    dispersions = (typeof(dispersion)<:AbstractDispersion ? (dispersion,) : dispersion)

    return Domain(type,name,bnd,lat,dielectric,pump,x,indices,ε,F,dispersions,interior,bulk,surface,corner,nnm,nnp)
end


text:
# You might not want to run this file completely, as the Pkg-commands can take a
# long time to complete.

# list all available packages:
Pkg.available()

# install one package (e.g. [Calculus](https://github.com/johnmyleswhite/Calculus.jl)) and all its dependencies:
Pkg.add("Calculus")

# to list all installed packages
Pkg.installed()

# to update all packages to their newest version
Pkg.update()

# to use a package:
using Calculus
# will import all functions of that package into the current namespace, so that
# it is possible to call
derivative(x -> sin(x), 1.0)
# without specifing the package it is included in.

import Calculus
# will enable you to specify which package the function is called from
Calculus.derivative(x -> cos(x), 1.0)

# Using `import` is especially useful if there are conflicts in function/type-names
# between packages.

max_stars_repo_name: inkydragon/JuliaByExample

max_issues_repo_name: inkydragon/JuliaByExample

max_forks_repo_name: inkydragon/JuliaByExample


text:
module Spk

using Statistics, FFTW, LinearAlgebra
using ..Histogram

const AbstractTS = AbstractVector{<:Real}

include("./PSTH.jl")
include("./Transform.jl")
include("./Xcorr.jl")
include("./Efficacy.jl")
#include("./Plot.jl")

end # end module


text:
@testset "validationplots.jl" begin
    @testset "ecdf(y)" begin
        n = 10
        data = -1 * collect(1:n)

        y, p = Extremes.ecdf(data)

        # y is length n and sorted
        @test length(y) == n
        @test issorted(y)

        # p is length n, sorted and between 0 and 1.
        @test length(p) == n
        @test issorted(p)
        @test p[1] >= 0 && p[end] <= 1

    end

    @testset "checkstationarity(m)" begin
        y = Variable("y", collect(1:10))

        sfm = MaximumLikelihoodEVA(BlockMaxima(y), [1.0, 1.0 ,0.0])
        nsfm = MaximumLikelihoodEVA(BlockMaxima(y, logscalecov = [y]), [1.0, 1.0, 1.0 ,0.0])

        # Model stationary and shouldbestationary no info
        @test_logs Extremes.checkstationarity(sfm.model)

        # Model non-stationary and shouldbestationary info
        @test_logs (:info, "The graph is optimized for stationary models and the model provided is not.") Extremes.checkstationarity(nsfm.model)

    end

    @testset "checknonstationarity(m)" begin
        y = Variable("y", collect(1:10))

        sfm = MaximumLikelihoodEVA(BlockMaxima(y), [1.0, 1.0 ,0.0])
        nsfm = MaximumLikelihoodEVA(BlockMaxima(y, logscalecov = [y]), [1.0, 1.0, 1.0 ,0.0])

        # Model stationary and !shouldbestationary info
        @test_logs (:info, "The graph is optimized for non-stationary models and the model provided is not.") Extremes.checknonstationarity(sfm.model)

        # Model non-stationary and !shouldbestationary no info
        @test_logs Extremes.checknonstationarity(nsfm.model)

    end

    include(joinpath("validationplots", "plots_std_tests.jl"))
    include(joinpath("validationplots", "plots_tests.jl"))

end


text:
module KRLS
using LinearAlgebra

include("kernel.jl")
include("filter.jl")
include("sim_env.jl")
end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# For curmod_*
include("testenv.jl")

using Test

@test_throws MethodError convert(Enum, 1.0)

macro macrocall(ex)
    @assert Meta.isexpr(ex, :macrocall)
    ex.head = :call
    for i in 2:length(ex.args)
        ex.args[i] = QuoteNode(ex.args[i])
    end
    insert!(ex.args, 3, __module__)
    return esc(ex)
end

@test_throws ArgumentError("no arguments given for Enum Foo") @macrocall(@enum Foo)

@enum Fruit apple orange kiwi
@test typeof(Fruit) == DataType
@test isbits(Fruit)
@test typeof(apple) <: Fruit <: Enum
@test Int(apple) == 0
@test Int(orange) == 1
@test Int(kiwi) == 2
@test Fruit(0) == apple
@test Fruit(1) == orange
@test Fruit(2) == kiwi
@test_throws ArgumentError Fruit(3)
@test_throws ArgumentError Fruit(-1)
@test Fruit(0x00) == apple
@test Fruit(big(0)) == apple
@test_throws MethodError Fruit(0.0)
@test typemin(Fruit) == apple
@test typemax(Fruit) == kiwi
@test Fruit(0) == apple
@test Fruit(1) == orange
@test Fruit(2) == kiwi
@test_throws ArgumentError Fruit(3)
@test_throws ArgumentError Fruit(-1)
@test UInt8(apple) === 0x00
@test UInt16(orange) === 0x0001
@test UInt128(kiwi) === 0x00000000000000000000000000000002
@test typeof(BigInt(apple)) <: BigInt
@test BigInt(apple) == 0
@test Bool(apple) == false
@test Bool(orange) == true
@test_throws InexactError Bool(kiwi)
@test instances(Fruit) == (apple, orange, kiwi)

f(x::Fruit) = "hey, I'm a Fruit"
@test f(apple) == "hey, I'm a Fruit"

d = Dict(apple=>"apple",orange=>"orange",kiwi=>"kiwi")
@test d[apple] == "apple"
@test d[orange] == "orange"
@test d[kiwi] == "kiwi"
vals = [apple,orange,kiwi]
for (i,enum) in enumerate(instances(Fruit))
    @test enum == vals[i]
end

@enum(QualityofFrenchFood, ReallyGood)
@test length(instances(QualityofFrenchFood)) == 1
@test typeof(ReallyGood) <: QualityofFrenchFood <: Enum
@test Int(ReallyGood) == 0

@enum Binary _zero=0 _one=1 _two=10 _three=11
@test Int(_zero) === 0
@test Int(_one) === 1
@test Int(_two) === 10
@test Int(_three) === 11
@enum Negative _neg1=-1 _neg2=-2
@test Int(_neg1) === -1
@test Int(_neg2) === -2
@test_throws InexactError UInt8(_neg1)
@enum Negative2 _neg5=-5 _neg4 _neg3
@test Int(_neg5) === -5
@test Int(_neg4) === -4
@test Int(_neg3) === -3

@test_throws ArgumentError("invalid value for Enum Test1, _zerofp = 0.0=0.0; values must be integers") @macrocall(@enum Test1 _zerofp=0.0)
@test_throws ArgumentError("invalid value for Enum Test11, _zerofp2 = 0.5=0.5; values must be integers") @macrocall(@enum Test11 _zerofp2=0.5)
@enum Test111 _zerobi=BigInt(1)
@test Integer(_zerobi) == 1

# can't use non-identifiers as enum members
@test_throws ArgumentError("""invalid argument for Enum Test2: if x
                                  1
                              else
                                  2
                              end""") @macrocall(@enum Test2  x ? 1 : 2)
@test_throws ArgumentError("invalid argument for Enum Test22: 1 = 2") @macrocall(@enum Test22 1=2)

# other Integer types of enum members
@enum Test3::UInt8 _one_Test3=0x01 _two_Test3=0x02 _three_Test3=0x03
@test Test3.size == 1
@test UInt8(_one_Test3) === 0x01
@test length(instances(Test3)) == 3

@enum Test4::UInt16 _one_Test4=0x01 _two_Test4=0x0002 _three_Test4=0x03
@test Test4.size == 2

@enum Test5::UInt32 _one_Test5=0x01 _two_Test5=0x00000002 _three_Test5=0x00000003
@test Test5.size == 4

@enum Test6::UInt128 _one_Test6=0x00000000000000000000000000000001 _two_Test6=0x00000000000000000000000000000002
@test Test6.size == 16
@test typeof(Integer(_one_Test6)) == UInt128

# enum values must be integers
@test_throws ArgumentError("invalid value for Enum Test7, _zero = \"zero\"=zero; values must be integers") @macrocall(@enum Test7 _zero="zero")
@test_throws ArgumentError("invalid value for Enum Test8, _zero = '0'=0; values must be integers") @macrocall(@enum Test8 _zero='0')
@test_throws ArgumentError("invalid value for Enum Test9, _zero = 0.5=0.5; values must be integers") @macrocall(@enum Test9 _zero=0.5)

# test macro handles keyword arguments
@enum(Test11, _zero_Test11=2,
              _one_Test11,
              _two_Test11=5,
              _three_Test11)

@test Int(_zero_Test11) == 2
@test Int(_one_Test11) == 3
@test Int(_two_Test11) == 5
@test Int(_three_Test11) == 6

# don't allow enum value to overflow
@test_throws ArgumentError("overflow in value \"y\" of Enum EnumOvf") @macrocall(@enum EnumOvf x=typemax(Int32) y)

# test for unique Enum values
@test_throws ArgumentError("values for Enum Test14 are not unique") @macrocall(@enum(Test14, _zero_Test14, _one_Test14, _two_Test14=0))

@test repr(apple) == "apple::$(string(Fruit)) = 0"
@test string(apple) == "apple"

@test reprmime("text/plain", Fruit) == "Enum $(string(Fruit)):\napple = 0\norange = 1\nkiwi = 2"
@test reprmime("text/plain", orange) == "orange::$(curmod_prefix)Fruit = 1"

@enum LogLevel DEBUG INFO WARN ERROR CRITICAL
@test DEBUG < CRITICAL

# serialization
let b = IOBuffer()
    serialize(b, apple)
    seekstart(b)
    @test deserialize(b) === apple
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
using  DataFrames
using LinearAlgebra

mutable struct HiddenLayer

    number_of_neurons::Integer
    weight_matrix::Matrix{Float64}
    weight_initializer::Function
    bias_vector::Vector{Float64}
    activation_function::Function

end

function return_name(x)
    return String(Symbol(x))
end

function activations(layer::HiddenLayer,x::Matrix{Float64})
    # Calculates the activations of the hidden layer neurons
    #
    # Parameters
    # ----------
    # `layer` is the HiddenLayer with neurons
    # `x` is the input matrix
    #
    # Returns
    # -------
    # Activation matrix after passing through hidden layer

    n_observations = size(x)[1]
    activation_matrix = zeros(layer.n_of_neurons, n_observations)

        for i = 1:n_observations
            activation_matrix[:, i] = layer.activation_func(layer.weight_matrix * x[i, :]' + layer.bias_vector)
        end

    activation_matrix
 end

 mutable struct ExtremeLearningMachine
    # Extreme Learning Machine
    #
    # Properties
    # ----------
    # `hidden_layer` is the hidden layer inside this ELM
    # `output_weights` is the weight matrix calculated during training

    hidden_layer::HiddenLayer
    output_weights::Matrix{Float64}

    function ExtremeLearningMachine(n_hidden_neurons::Integer)
	this = new()

	weight_matrix =Array{Float64}(undef,0,0)
	bias_vector = rand(n_hidden_neurons)
	this.hidden_layer = HiddenLayer(n_hidden_neurons, weight_matrix,glorot_uniform, bias_vector, sigmoid)
	println("######################### ELM STRUCTURE ##############################")
    println("# Number of Hidden Neurons: ",this.hidden_layer.number_of_neurons)
    println("# Activation Function: ",return_name(this.hidden_layer.activation_function))
    println("######################################################################")
    this

    end

end

function compile(elm::ExtremeLearningMachine,input_dim::Integer)
    elm.hidden_layer.weight_matrix=elm.hidden_layer.weight_initializer(input_dim,elm.hidden_layer.number_of_neurons)
    println("compiled")
end

function find_activations(layer::HiddenLayer,x::Matrix{Float64})
# Calculates the activations of the hidden layer neurons
#
# Parameters
# ----------
# `layer` is the HiddenLayer with neurons
# `x` is the input matrix
#
# Returns
# -------
# Activation matrix after passing through hidden layer

n_observations = size(x)[1]
act_matrix = zeros(layer.number_of_neurons, n_observations)
transposedx=x'
for i = 1:n_observations
    act_matrix[:, i] =layer.weight_matrix * transposedx[ :,i] + layer.bias_vector
end

map((x)->layer.activation_function(x),act_matrix)

end

function fit!(elm::ExtremeLearningMachine,x,y)
    # Trains the elm using the given training data
    #
    # Parameters
    # ----------
    # `elm` the ELM to train
    # `x` input data
    # `y` output data

    if typeof(x) == DataFrame
    x_mat = float(array(x))
    else
    x_mat = x
    end

    if typeof(y) != Vector{Float64}
    y_vec = float(array(y))
    else
    y_vec = y
    end

    n_observations, n_inputs = size(x_mat)
    weight_matrix = rand(elm.hidden_layer.number_of_neurons, n_inputs) * 2 .- 1
    elm.hidden_layer.weight_matrix = weight_matrix

    act_matrix = find_activations(elm.hidden_layer, x_mat)

    output_weights = y_vec' * pinv(act_matrix)
    elm.output_weights = output_weights
    end

    function fitSVD!(elm::ExtremeLearningMachine,x,y)
        # Trains the elm using the given training data
        #
        # Parameters
        # ----------
        # `elm` the ELM to train
        # `x` input data
        # `y` output data

        if typeof(x) == DataFrame
        x_mat = float(array(x))
        else
        x_mat = x
        end

        if typeof(y) != Vector{Float64}
        y_vec = float(array(y))
        else
        y_vec = y
        end

        n_observations, n_inputs = size(x_mat)
        weight_matrix = rand(elm.hidden_layer.number_of_neurons, n_inputs) * 2 .- 1
        elm.hidden_layer.weight_matrix = weight_matrix

        act_matrix = find_activations(elm.hidden_layer, x_mat)

        mySvd =  SVDCM(act_matrix)
        println(size(mySvd.V))
        println(size(((mySvd.S_inv))))
        c=mySvd.V*mySvd.S_inv
        d=c*mySvd.U'
        output_weights = y_vec' *d
        elm.output_weights = output_weights
        end

function predict(elm::ExtremeLearningMachine,x)
# Predicts the output
#
# Parameters
# ----------
# `elm` the trained ELM
# `x` input data to predict (Matrix)

if typeof(x) == DataFrame
x_mat = float(array(x))
else
x_mat = x
end

act_matrix = find_activations(elm.hidden_layer, x_mat)
vec(elm.output_weights * act_matrix)
end


text:
# src/secondorder.jl contains function for solving second order ODEs
# can use either numerical integration or closed for formulas
# secondorder(a,b,c) generates the function you pass to ODE.jl RK methods
# solution(a,b,c) is the closed form solver
# ivp solves the initial value problem.
using Optim
using NLsolve

""" secondorder(α,β,γ) yields a function f
    which calculates the ode rule

        αx'' + βx' + γx = 0

    using the substitution

        y = [x; x']

        y[1]' = y[2]

        y[2]' = -α^-1 (βy[2] + γy[1])

"""
function secondorder(α, β, γ)
        return (t, y) -> [y[2]; -(β*y[2] + γ*y[1])/α]
end

function secondorder(x::Vector)
    return secondorder(x...)
end

"""descriminant(a,b,c): quadtratic equation descriminant of ax^2 + bx + c = 0
"""
descriminant(a,b,c) = √(complex(b^2-4a*c))

"""roots(a,b,c): solve the quadratic equation ax^2 + bx + c = 0 """
function roots(a,b,c)
    d = descriminant(a,b,c);
    return (-b+d)/2a, (-b-d)/2a
end

doc"""solution(a,b,c,y₀,t₀) general solution of second order linear homogenous

Arguments:
   - a,b,c::Real the coeffients of the ODE

Returns:
   - func(t, c₁, c₂): the solution at t when the coefficients are c₁, c₂
   - deriv(t, c₁, c₂): the derivative at t when the coefficients are c₁, c₂

Notes: func is derived from the textbook definition and deriv is based on
 applying calculus rules to that (a symbolic algebra tool did this).
"""
function solution(a,b,c)
    @show d = descriminant(a,b,c)
    @show r₁, r₂ = roots(a,b,c)
    if real(d) == d
        info("characteristic polynomial has real roots")
        return (t, c₁, c₂) -> c₁ * exp(r₁*t) + c₂*exp(r₂*t),
               (t, c₁, c₂) -> c₁*r₁*exp(r₁*t) + c₂*r₂*exp(r₂*t)

    end
    if imag(d) != 0.0
        info("characteristic polynomial has complex roots")
        @show α, β = -b/2a, imag(d/2a)
        func(t, c₁, c) = begin
            return (c₁*exp(α*t)*cos(β*t) + c*exp(α*t)*sin(β*t))
        end
        deriv(t,c₁,c₂) = begin
            expα = exp(α*t)
            cosβ = cos(β*t)
            sinβ = sin(β*t)
            return (-β*c₁*expα*sinβ +
                    α*c₂*expα*sinβ +
                    α*c₁*expα*cosβ +
                    β*c₂*expα*cosβ
                    )
            end
        return func, deriv
    end
end

"""ivp(a,b,c,y₀,t₀) solves second order linear homogenous initial value problem

Arguments:
   - a,b,c::Real the coeffients of the ODE
   - y₀ the initial value as a vector
   - t₀ the initial time value

Returns:
   - result : the Optim.jl result object of solving the equations
   - general: the general form solution (see `solution`)
   - deriv:   the general form derivative (see `solution`)

Note: result.zero is the coefficients to plug into the general solution
to get the correct particular solution.
"""
function ivp(a,b,c, y₀, t₀)
    x₀ = y₀[1]
    v₀ = y₀[2]
    general, deriv = solution(a,b,c)

    function f!(params, fvec)
        fvec[1] = x₀ - general(t₀, params...)
        fvec[2] = v₀ - deriv(t₀, params...)
    end
    result = nlsolve(f!,
                     [0.0,0.0],
                     autodiff=true,
                     method=:trust_region,
                     store_trace=true)
    return result, general, deriv
end

max_stars_repo_name: jpfairbanks/ODEOpt

max_issues_repo_name: jpfairbanks/ODEOpt

max_forks_repo_name: jpfairbanks/ODEOpt


text:
using KeywordCalls
using Test

f(nt::NamedTuple{(:c,:b,:a)}) = nt.a^3 + nt.b^2 + nt.c
@kwcall f(c,b,a)

struct Foo{N,T}
    nt::NamedTuple{N,T}
end

Foo(nt::NamedTuple{(:a,:b),T}) where {T} = Foo{(:a,:b), T}(nt)

@kwcall Foo(a,b)

@testset "KeywordCalls.jl" begin
    @testset "Functions" begin
        # @test @inferred f(1, 2, 3) == 32
        @test @inferred f(a=1, b=2, c=3) == 8
        @test @inferred f((a=1, b=2, c=3)) == 8
    end

    @testset "Constructors" begin
        @test @inferred Foo((b=1,a=2)).nt == (a=2,b=1)
    end
end


text:
using ReservoirComputing

#model parameters
const approx_res_size = 30
const radius = 1.2
const degree = 6
const first_activation = tanh
const second_activation = identity
const first_lambda = 0.45
const second_lambda = 0.3
const sigma = 0.3
const beta = 0.0
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 3
const out_size = 3
const extended_states = false
const h_steps = 2

W_in = init_input_layer(approx_res_size, in_size, sigma)
W = init_reservoir_givendeg(approx_res_size, radius, degree)

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#constructor 1
esn = dafESN(approx_res_size,
            train,
            degree,
            radius,
            first_lambda,
            second_lambda,
            first_activation = first_activation,
            second_activation = second_activation,
            sigma = sigma,
            alpha = alpha,
            nla_type = nla_type,
            extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(first_activation, esn.first_activation)
@test isequal(second_activation, esn.second_activation)
@test isequal(first_lambda, esn.first_lambda)
@test isequal(second_lambda, esn.second_lambda)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = dafESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

#constructor 2
esn = dafESN(W,
            train,
            first_lambda,
            second_lambda,
            first_activation = first_activation,
            second_activation = second_activation,
            sigma = sigma,
            alpha = alpha,
            nla_type = nla_type,
            extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(first_activation, esn.first_activation)
@test isequal(second_activation, esn.second_activation)
@test isequal(first_lambda, esn.first_lambda)
@test isequal(second_lambda, esn.second_lambda)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = dafESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

#constructor 3
esn = dafESN(approx_res_size,
            train,
            degree,
            radius,
            first_lambda,
            second_lambda,
            W_in,
            first_activation = first_activation,
            second_activation = second_activation,
            alpha = alpha,
            nla_type = nla_type,
            extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(first_activation, esn.first_activation)
@test isequal(second_activation, esn.second_activation)
@test isequal(first_lambda, esn.first_lambda)
@test isequal(second_lambda, esn.second_lambda)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = dafESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

#constructor 4
esn = dafESN(W,
            train,
            first_lambda,
            second_lambda,
            W_in,
            first_activation = first_activation,
            second_activation = second_activation,
            alpha = alpha,
            nla_type = nla_type,
            extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), esn.res_size)
@test isequal(train, esn.train_data)
@test isequal(alpha, esn.alpha)
@test isequal(first_activation, esn.first_activation)
@test isequal(second_activation, esn.second_activation)
@test isequal(first_lambda, esn.first_lambda)
@test isequal(second_lambda, esn.second_lambda)
@test isequal(nla_type, esn.nla_type)
@test size(esn.W) == (esn.res_size, esn.res_size)
@test size(esn.W_in) == (esn.res_size, esn.in_size)
@test size(esn.states) == (esn.res_size, train_len)

#test train
W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size)
#test predict
output = dafESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)

output = dafESNpredict_h_steps(esn, predict_len, h_steps, test, W_out)
@test size(output) == (out_size, predict_len)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -88.40110862821636 <= q <= 80.25263273980681)

                     add_NL_constraint(m, :(gelu(-0.3402051397109518 + -0.28456924845213605*gelu(-0.42610478134649465 + 0.43223652991075534*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.6888193490874222*gelu(-0.17528834888704425 + -0.1425247358809414*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.026837640444652866*gelu(-0.005280183583000753 + 0.6172883711874637*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.08739699596173622*gelu(0.26121381489212014 + -0.018200625348021138*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.2360509887003217*gelu(-0.09690410604182231 + 0.6473593808554208*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + gelu(0.2351041904461102 + 0.11979579421065711*gelu(-0.42610478134649465 + 0.43223652991075534*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.012124046222294904*gelu(-0.17528834888704425 + -0.1425247358809414*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.704332995743346*gelu(-0.005280183583000753 + 0.6172883711874637*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.7162932434639226*gelu(0.26121381489212014 + -0.018200625348021138*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.3351546298978829*gelu(-0.09690410604182231 + 0.6473593808554208*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + gelu(-0.4731175976624935 + -0.5144698215629462*gelu(-0.42610478134649465 + 0.43223652991075534*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.005430014157027596*gelu(-0.17528834888704425 + -0.1425247358809414*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.05201762086802919*gelu(-0.005280183583000753 + 0.6172883711874637*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.06728299071980626*gelu(0.26121381489212014 + -0.018200625348021138*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.8893971601390014*gelu(-0.09690410604182231 + 0.6473593808554208*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + gelu(-0.8380258290664089 + 0.49770714260784343*gelu(-0.42610478134649465 + 0.43223652991075534*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.732919146629087*gelu(-0.17528834888704425 + -0.1425247358809414*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.2983691588163282*gelu(-0.005280183583000753 + 0.6172883711874637*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.6200389337758314*gelu(0.26121381489212014 + -0.018200625348021138*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.03199113104777407*gelu(-0.09690410604182231 + 0.6473593808554208*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) + gelu(0.7168755118620993 + -0.38856981611164576*gelu(-0.42610478134649465 + 0.43223652991075534*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.5575457569666615*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.7868073625875285*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.5078849457963623*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.2883931154790984*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.9818525806719265*gelu(-0.17528834888704425 + -0.1425247358809414*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + 0.07932201643104397*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.9443705506849667*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.7094815218028745*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + -0.7352206438426769*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.9400532013818674*gelu(-0.005280183583000753 + 0.6172883711874637*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.07641401576689111*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.42701726221355685*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + -0.5680942985776078*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.9418726475599537*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + -0.928831604242299*gelu(0.26121381489212014 + -0.018200625348021138*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.03849358234625422*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + -0.856604860683484*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.2143578715735437*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.706261759527687*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5]))) + 0.6913287958041021*gelu(-0.09690410604182231 + 0.6473593808554208*gelu(-0.9360387382859474 + 0.27356363334073697*$(x[1]) + 0.5623428270051978*$(x[2]) + -0.9519487490584724*$(x[3]) + 0.5639775172803199*$(x[4]) + 0.8360124736609604*$(x[5])) + -0.6752519993458765*gelu(0.3895701371443887 + 0.5298392176190725*$(x[1]) + -0.9900818182438851*$(x[2]) + 0.9005096697487907*$(x[3]) + 0.008235253700211498*$(x[4]) + 0.03120241694369552*$(x[5])) + 0.771479525571356*gelu(0.7802421455771191 + -0.3899284089358539*$(x[1]) + -0.510799262552962*$(x[2]) + -0.16925880893067768*$(x[3]) + -0.7702701325797845*$(x[4]) + -0.21353502867805174*$(x[5])) + 0.8906247556095672*gelu(-0.4643333519515682 + -0.743154946296293*$(x[1]) + 0.08142803382004216*$(x[2]) + 0.8067509325852682*$(x[3]) + -0.6642228946665827*$(x[4]) + -0.03791136441363552*$(x[5])) + 0.560146720685823*gelu(0.6488423858989929 + -0.3655334329180211*$(x[1]) + 0.8449912783652063*$(x[2]) + -0.7635680273216123*$(x[3]) + 0.17014039681297044*$(x[4]) + -0.14496676502055283*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "BareExpansion" begin

MAX_ORDER=2
diag=Tethys.Diag(0,1.0,MAX_ORDER,[0.1,0.2,1.0 , 0.3,0.4,1.0], [0,0,0,0,0,0])

H=Tethys.FrohlichHamiltonian(α=1.0, μ=-1.2)

GF=Tethys.BareExpansion(diag, H)
@test GF ≈ -0.7444 atol=0.1

println("Splish splash: ",diag," ==♓> GF = ", GF)

GF=Tethys.Monte!(diag, H)
println("Monte: ",diag," ==♓>  GF= ", GF)

@time Tethys.Monte!(diag,H)

fout=open("Monte.dat","w")
for i in 1:100
    GF=Tethys.Monte!(diag, H)
    @printf(fout,"%d %f \n",i,GF)
end
close(fout)

end


text:
# Test the mut_evolve_subgoal() function from src/Evolve.jl

#include("../src/CGP.jl")
include("../src/Pop_evolution.jl")
using Main.CGP
MyInt = Main.CGP.MyInt

#=
# This is a test function but is included here for convenience
function test_pop_evolve(p::Parameters,popsize,funcs,max_pop_steps::Int64,max_indiv_steps::Int64,ngoals::Int64;
    robust_sel::Bool=false)
  print_parameters(p)
  println("ngoals: ",ngoals)
  gl = randgoallist(ngoals,p.numinputs,p.numoutputs)
  (max_fit_all_gens, pop) = pop_evolve( p, popsize, max_pop_steps, max_indiv_steps, gl, funcs, robust_sel=robust_sel)
end
=#

numinputs = 2
numoutputs = 2
popsize = 3
use_robust=true
funcs = default_funcs(numinputs)
println("funcs: ",funcs)
p = Parameters( numinputs=numinputs, numoutputs=numoutputs, numinteriors=4, numlevelsback=5 )
context = construct_context(numinputs)
max_indiv_steps = 20
max_pop_steps = 12
ngoals = 4   # number goals in goallist
gl = randgoallist(ngoals,p.numinputs,p.numoutputs)
robust_sel=true
all_max_sel = true
#(steps,worse,same,c,output,goallist,matched_goals,matched_goals_list)= test_pop_evolve(p,funcs,maxsteps,ngoals)
#result = test_pop_evolve(p,popsize,funcs,max_pop_steps,max_indiv_steps,ngoals,robust_sel=use_robust)
#(max_fit_all_gens, pop) = pop_evolve( p, popsize, max_pop_steps, max_indiv_steps, gl, funcs, robust_sel=robust_sel,all_max_sel=all_max_sel)
result = pop_evolve( p, popsize, max_pop_steps, max_indiv_steps, gl, funcs, robust_sel=robust_sel,all_max_sel=all_max_sel);
nothing
#println("result: ",result)


text:
struct CIFAR10 <: ColorImage{MLImage} end
getModule(::CIFAR10) = MLDatasets.CIFAR10
size(::CIFAR10) = (50000, 0, 10000)
problem(::CIFAR10) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
using Test

function find_smallest(arr)
    smallest = arr[1]
    smallest_index = 1
    for i = 1:length(arr)
        if arr[i] < smallest
            smallest = arr[i]
            smallest_index = i
        end
    end
    return smallest_index
end

function selection_sort(arr)
    new_arr = Array{Int64}(undef,0)
    for i = 1:length(arr)
        smallest = find_smallest(arr)
        append!(new_arr,arr[smallest])
        deleteat!(arr,smallest)
    end
    return new_arr
end

arr = [5,3,6,2,10]
@test selection_sort(arr) == [2,3,5,6,10]


text:
module test_stdlib_statistics

using Statistics, Test

X = rand(5)

@test cor(X, 2X) == 1.0
@test cov(X, 2X) - cov(X, X) == var(X)

end # module test_stdlib_statistics


text:
"""
This module contains functions and types necessary for performing
nonadiabatic molecular dynamics.

Dynamics is performed using [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/).
As such, this module is centered around the implementation of the functions
necessary to integrate the dynamics.

For deterministic Hamiltonian methods, the central function is [`DynamicsMethods.motion!`](@ref),
which is the inplace form of the function
to be integrated by [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/).

Further, methods that have discontinuities, such as surface hopping, use the
[callback interface](https://diffeq.sciml.ai/stable/features/callback_functions/#callbacks)
provided by `DifferentialEquations.jl`.
"""
module DynamicsMethods

using DiffEqBase: DiffEqBase
using TypedTables: TypedTables
using UnitfulAtomic: austrip
using Reexport: @reexport
using OrdinaryDiffEq: OrdinaryDiffEq
using ComponentArrays: ComponentVector

using NQCDynamics: AbstractSimulation, Simulation, RingPolymerSimulation,
    DynamicsOutputs
using NQCBase: austrip_kwargs

export DynamicsVariables
export run_trajectory

"""
Each type of dynamics subtypes `Method` which is passed to
the `AbstractSimulation` as a parameter to determine the type of
dynamics desired.
"""
abstract type Method end

"""
    motion!(du, u, sim, t)

As per `DifferentialEquations.jl`, this function is implemented for
each method and defines the time derivatives of the `DynamicalVariables`.

We require that each implementation ensures `du` and `u` are subtypes
of `DynamicalVariables` and `sim` subtypes `AbstractSimulation`.
"""
function motion! end

"Provides the DEProblem for each type of simulation."
create_problem(u0, tspan, sim) =
    OrdinaryDiffEq.ODEProblem(motion!, u0, tspan, sim; callback=get_callbacks(sim))

"Choose a default algorithm for solving the differential equation."
select_algorithm(::AbstractSimulation) = OrdinaryDiffEq.VCABM5()

"Select the default callbacks for this simulation type."
get_callbacks(::AbstractSimulation) = nothing

"""
    DynamicsVariables(::AbstractSimulation, args...)

For each dynamics method this function is implemented to provide the variables for the
dynamics in the appropriate format.

By default, `DynamicsVariables` is set up for the classical case and takes
`sim`, `v`, `r` as arguments and returns a `ComponentVector(v=v, r=r)`
which is used as a container for the velocities and positions during
classical dynamics.
"""
DynamicsVariables(::AbstractSimulation, v, r) = ComponentVector(v=v, r=r)

"""
    run_trajectory(u0, tspan::Tuple, sim::AbstractSimulation;
                   output=(:u,), saveat=[], algorithm=select_algorithm(sim), kwargs...)

Solve a single trajectory starting from `u0` with a timespan `tspan`
for the simulation `sim`.

# Keyword arguments

`output` specifies the quantities that should be saved during the dynamics simulation.
The options for this keyword are any of the functions found in
`src/DynamicsMethods/output.jl`.

The rest of the keywords are the usual arguments for the `solve` function from
`DifferentialEquations.jl`.
It is possible to use `Unitful` quantities for any of the arguments since these are
automatically converted to atomic units internally.

# Output

The function returns a `Table` from `TypedTables.jl` with columns for time `t` and
every quantity specified in the `output` tuple.
"""
function run_trajectory(u0, tspan::Tuple, sim::AbstractSimulation;
        output=(:u,), saveat=[], callback=nothing, algorithm=select_algorithm(sim),
        kwargs...)

    if !(output isa Tuple)
        output = (output,)
    end

    problem = create_problem(u0, austrip.(tspan), sim)

    saving_callback, vals = DynamicsOutputs.create_saving_callback(output; saveat=austrip.(saveat))
    callback_set = DiffEqBase.CallbackSet(callback, saving_callback, get_callbacks(sim))
    problem = DiffEqBase.remake(problem, callback=callback_set)

    stripped_kwargs = austrip_kwargs(;kwargs...)
    @info "Performing a single trajectory."
    stats = @timed DiffEqBase.solve(problem, algorithm; stripped_kwargs...)
    @info "Finished after $(stats.time) seconds."
    out = [(;zip(output, val)...) for val in vals.saveval]
    TypedTables.Table(t=vals.t, out)
end

include("ClassicalMethods/ClassicalMethods.jl")
@reexport using .ClassicalMethods

include("SurfaceHoppingMethods/SurfaceHoppingMethods.jl")
@reexport using .SurfaceHoppingMethods

include("MappingVariableMethods/MappingVariableMethods.jl")
@reexport using .MappingVariableMethods

include("EhrenfestMethods/EhrenfestMethods.jl")
@reexport using .EhrenfestMethods

include("IntegrationAlgorithms/IntegrationAlgorithms.jl")

end # module


text:
import Base: cat
export partition

mutable struct AllocateArray{T,N} <: ArrayOp{T,N}
    eltype::Type{T}
    f::Function
    domain::ArrayDomain{N}
    domainchunks
end
size(a::AllocateArray) = size(a.domain)

function _cumlength(len, step)
    nice_pieces = div(len, step)
    extra = rem(len, step)
    ps = [step for i=1:nice_pieces]
    cumsum(extra > 0 ? vcat(ps, extra) : ps)
end

function partition(p::Blocks, dom::ArrayDomain)
    DomainBlocks(map(first, indexes(dom)),
        map(_cumlength, map(length, indexes(dom)), p.blocksize))
end

function stage(ctx, a::AllocateArray)
    alloc(idx, sz) = a.f(idx, a.eltype, sz)
    thunks = [delayed(alloc)(i, size(x)) for (i, x) in enumerate(a.domainchunks)]
    DArray(a.eltype,a.domain, a.domainchunks, thunks)
end

function Base.rand(p::Blocks, eltype::Type, dims)
    s = rand(UInt)
    f = function (idx, x...)
        rand(MersenneTwister(s+idx), x...)
    end
    d = ArrayDomain(map(x->1:x, dims))
    AllocateArray(eltype, f, d, partition(p, d))
end

Base.rand(p::Blocks, t::Type, dims::Integer...) = rand(p, t, dims)
Base.rand(p::Blocks, dims::Integer...) = rand(p, Float64, dims)
Base.rand(p::Blocks, dims::Tuple) = rand(p, Float64, dims)

function Base.randn(p::Blocks, dims)
    s = rand(UInt)
    f = function (idx, x...)
        randn(MersenneTwister(s+idx), x...)
    end
    d = ArrayDomain(map(x->1:x, dims))
    AllocateArray(Float64, f, d, partition(p, d))
end
Base.randn(p::Blocks, dims::Integer...) = randn(p, dims)

function Base.ones(p::Blocks, eltype::Type, dims)
    d = ArrayDomain(map(x->1:x, dims))
    AllocateArray(eltype, (_, x...) -> ones(x...), d, partition(p, d))
end
Base.ones(p::Blocks, t::Type, dims::Integer...) = ones(p, t, dims)
Base.ones(p::Blocks, dims::Integer...) = ones(p, Float64, dims)
Base.ones(p::Blocks, dims::Tuple) = ones(p, Float64, dims)

function Base.zeros(p::Blocks, eltype::Type, dims)
    AllocateArray(eltype, (_, x...) -> zeros(x...), ArrayDomain(map(x->1:x, dims)), p)
end
Base.zeros(p::Blocks, t::Type, dims::Integer...) = zeros(p, t, dims)
Base.zeros(p::Blocks, dims::Integer...) = zeros(p, Float64, dims)
Base.zeros(p::Blocks, dims::Tuple) = zeros(p, Float64, dims)

function Base.sprand(p::Blocks, m::Integer, n::Integer, sparsity::Real)
    s = rand(UInt)
    f = function (idx, t,sz)
        sprand(MersenneTwister(s+idx), sz...,sparsity)
    end
    d = ArrayDomain((1:m, 1:n))
    AllocateArray(Float64, f, d, partition(p, d))
end

function Base.sprand(p::Blocks, n::Integer, sparsity::Real)
    s = rand(UInt)
    f = function (idx,t,sz)
        sprand(MersenneTwister(s+idx), sz...,sparsity)
    end
    AllocateArray(Float64, f, d, partition(p, ArrayDomain((1:n,))))
end


text:
abstract type IndexFunction end

(f::IndexFunction)(A::BaseOperator) = BaseOperator(A.t,A.name,f.(A.inds))
(f::IndexFunction)(d::δ) = δ(f(d.iA),f(d.iB))
(f::IndexFunction)(par::Param) = Param(par.name,par.state,f.(par.inds))
(f::IndexFunction)(ops::BaseOpProduct) = BaseOpProduct(f.(ops.v))
(f::IndexFunction)(ev::T) where T<:Union{ExpVal,Corr} = T(f(ev.ops))
(f::IndexFunction)(A::QuTerm,nsuminds=A.nsuminds) = QuTerm(nsuminds,f.(A.δs),f.(A.params),f.(A.expvals),f.(A.corrs),f(A.bares))
# some index functions have to be reset for each new term in a sum
(f::IndexFunction)(A::QuExpr) = _map_quexpr_ops(x->(adapt_for_new_term!(f); f(x)),A)
# by default, nothing to do
adapt_for_new_term!(f::IndexFunction) = f

function _canon_ind(n)
    i,c = divrem(n-1,15)
    i == 0 ? QuIndex('i'+c) : QuIndex('i'+c,i+1)
end
_canon_ind_to_num(ind::QuIndex) = ind.sym-'i'+1 + (ind.num==typemin(ind.num) ? 0 : (ind.num-1)*15)

mutable struct canon_inds <: IndexFunction
    icurr::IndexInt
    canon_inds() = new(0)
end
(f::canon_inds)(ind::QuIndex) = isnoindex(ind) ? ind : _canon_ind(f.icurr+=one(f.icurr))
adapt_for_new_term!(f::canon_inds) = (f.icurr = 0; f)

mutable struct canon_inds_remember <: IndexFunction
    icurr::IndexInt
    replacements::Dict{QuIndex,QuIndex}
    canon_inds_remember() = new(0,Dict{QuIndex,QuIndex}())
end
function (f::canon_inds_remember)(ind::QuIndex)
    isnoindex(ind) && return ind
    i = _canon_ind(f.icurr+=one(f.icurr))
    f.replacements[i] = ind
    i
end
adapt_for_new_term!(f::canon_inds_remember) = error("canon_inds_remember should not be applied on QuExpr.")

struct shift_sumind <: IndexFunction
    n::IndexInt
end
(f::shift_sumind)(ind::QuIndex) = issumindex(ind) ? sumindex(ind.num+f.n) : ind

struct replace_inds <: IndexFunction
    replacements::Dict{QuIndex,QuIndex}
    replace_inds(args...) = new(Dict{QuIndex,QuIndex}(args...))
end
(f::replace_inds)(ind::QuIndex) = get(f.replacements,ind,ind)

mutable struct reorder_suminds <: IndexFunction
    icurr::IndexInt
    replacements::Dict{QuIndex,QuIndex}
    reorder_suminds() = new(0,Dict{QuIndex,QuIndex}())
end
(f::reorder_suminds)(ind::QuIndex) = issumindex(ind) ? get!(()->sumindex(f.icurr += one(f.icurr)), f.replacements, ind) : ind
adapt_for_new_term!(f::reorder_suminds) = (f.icurr=0; empty!(f.replacements); f)

mutable struct popped_out_suminds <: IndexFunction
    conserved_ind::QuIndex
    icurr::IndexInt
    suminds::Vector{QuIndex}
    popped_out_suminds(n) = new(n,0,QuIndex[])
end
function (f::popped_out_suminds)(ind::QuIndex)
    if isnoindex(ind)
        ind
    elseif ind==f.conserved_ind
        sumindex(1)
    else
        push!(f.suminds,ind)
        _canon_ind(f.icurr+=one(f.icurr))
    end
end
adapt_for_new_term!(f::popped_out_suminds) = error("popped_out_suminds should not be applied on QuExpr.")

hasind(ind::QuIndex,d::δ) = ind in (d.iA,d.iB)
hasind(ind::QuIndex,A::Union{BaseOperator,Param}) = ind in A.inds
hasind(ind::QuIndex,ev::Union{ExpVal,Corr}) = hasind(ind,ev.ops)
hasind(ind::QuIndex,A::BaseOpProduct) = any(hasind.((ind,),A.v))

split_by_ind(ind::QuIndex,v) = (iwiths = hasind.((ind,),v); (v[.!iwiths], v[iwiths]))

indices(d::δ)::Vector{QuIndex} = [d.iA,d.iB]
indices(A::Union{BaseOperator,Param})::Vector{QuIndex} = collect(assignedinds(A.inds)) # collect to guarantee copy
indices(v::AbstractVector)::Vector{QuIndex} = mapreduce(indices, vcat, v; init=QuIndex[])
indices(A::BaseOpProduct)::Vector{QuIndex} = indices(A.v)
indices(A::Union{ExpVal,Corr})::Vector{QuIndex} = indices(A.ops)
indices(A::QuTerm, withδs::Bool=true)::Vector{QuIndex} = begin
    args = withδs ? (A.δs,A.params,A.expvals,A.corrs,A.bares) : (A.params,A.expvals,A.corrs,A.bares)
    mapreduce(indices, vcat, args; init=QuIndex[])
end
indices(A::QuExpr)::Vector{QuIndex} = mapreduce(indices, vcat, sort!(collect(keys(A.terms))); init=QuIndex[])

"`extindices(A)` return externally visible indices of an expression"
extindices(A) = filter(!issumindex,indices(A))

"`symmetric_index_nums(A)` return sequence of numbers of exchange-symmetric indices"
symmetric_index_nums(A) = symmetric_index_nums(QuExpr(A))
function symmetric_index_nums(A::QuExpr)
    inds = extindices(A)
    Nsyms = [1]
    for ii=2:length(inds)
        i1,i2 = inds[ii-1], inds[ii]
        Aexc = normal_form(replace_inds(i1=>i2,i2=>i1)(A))
        if A == Aexc
            Nsyms[end] += 1
        else
            push!(Nsyms,1)
        end
    end
    Nsyms
end


text:
using DelayEmbeddings
using PerronFrobenius
using Simplices
using Test

@testset "Grid estimators" begin
    n = 200
    x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
    push!(x, 5)
    τs = (0, -2, -3)
    js = (1, 1, 1)
    D = genembed(x, τs, js)

    @testset "SingleGrid" begin
        @test SingleGrid(RectangularBinning(5)) isa SingleGrid
        @test SingleGrid(RectangularBinning(5), "circular") isa SingleGrid

        method = SingleGrid(RectangularBinning(5))
        tog = transopergenerator(D, method)
        TO = tog()
        E = TransferOperatorGenerator{SingleGrid{RectangularBinning}}
        @test TO isa TransferOperatorApproximation{<:E}
        M = TO.M
        @test all([M[i, :] |> sum .≈ 1.0 for i = 1:size(M, 2)])

        @test transferoperator(D, method) isa TransferOperatorApproximation{<:E}
    end
end
@testset "Triangulation estimators" begin
    n = 20
    x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
    push!(x, 5)
    τs = (0, -2, -3)
    js = (1, 1, 1)
    D = genembed(x, τs, js)

    @testset "SimplexExact" begin
        tog = transopergenerator(D, SimplexExact())
        TO = tog()
        E = TransferOperatorGenerator{SimplexExact}
        @test TO isa TransferOperatorApproximation{<:E}
        M = TO.M
        @test all(sum(M, dims = 2) .≈ 1.0)
    end

    @testset "SimplexApprox" begin
        tog = transopergenerator(D, SimplexPoint())
        TO = tog(n = 50)
        E = TransferOperatorGenerator{SimplexPoint}
        @test TO isa TransferOperatorApproximation{<:E}
        M = TO.M
        @test all(sum(M, dims = 2) .≈ 1.0)
    end
end

@testset "Invariant measures" begin

n = 200
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
push!(x, 5)
τs = (0, -2, -3)
js = (1, 1, 1)
D = genembed(x, τs, js)

@testset "Grid estimators" begin
    method = SingleGrid(RectangularBinning(5))
    to = transferoperator(D, method)
    iv = invariantmeasure(to)
    @test iv isa InvariantDistribution
end

@testset "Triangulation estimators" begin
n = 20
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
push!(x, 5)
τs = (0, -2, -3)
js = (1, 1, 1)
D = genembed(x, τs, js)

    @testset "SimplexPoint" begin
        method = SimplexPoint()
        to = transferoperator(D, method)
        iv = invariantmeasure(to)
        @test invariantmeasure(to) isa InvariantDistribution
        @test invariantmeasure(D, method) isa InvariantDistribution
    end

    @testset "SimplexExact" begin
        method = SimplexExact()
        to = transferoperator(D, method)
        iv = invariantmeasure(to)
        @test invariantmeasure(to) isa InvariantDistribution
        @test invariantmeasure(D, method) isa InvariantDistribution
    end

end

end


text:
using Complementarity, Test

@testset "mcp test 1 with PATHSolver" begin

    m = MCPModel()

    @variable(m, x3 >= 0)
    @variable(m, x4 >= 0)
    @variable(m, x1 >= 0)
    @variable(m, x2 >= 0)

    @mapping(m, F2, x3-2x4 +2)
    @mapping(m, F3, x1-x2+2x3-2x4 -2)
    @mapping(m, F4, x1+2x2-2x3+4x4 -6)
    @mapping(m, F1, -x3-x4 +2)

    @complementarity(m, F4, x4)
    @complementarity(m, F1, x1)
    @complementarity(m, F3, x3)
    @complementarity(m, F2, x2)

    # This will show deprecation warning.
    status = solveLCP(m)
    z = [result_value(x1), result_value(x2), result_value(x3), result_value(x4)]
    @show z
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

    # This will show deprecation warning.
    status = solveMCP(m; linear=true)
    z = [result_value(x1), result_value(x2), result_value(x3), result_value(x4)]
    @show z
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

    status = solveMCP(m; convergence_tolerance=1e-8, output="yes", time_limit=3600)
    z = [result_value(x1), result_value(x2), result_value(x3), result_value(x4)]
    @show z
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

    status = solveMCP(m)
    z = [result_value(x1), result_value(x2), result_value(x3), result_value(x4)]
    @show z
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

end

println("------------------------------------------------------------------")

@testset "mcp test 2 with PATHSolver" begin

    m = MCPModel()

    @variable(m, x3 >= 0)
    @variable(m, x4 >= 0)
    @variable(m, x1 >= 0)
    @variable(m, x2 >= 0)

    @mapping(m, F2, x3-2x4 +2)
    @mapping(m, F3, x1-x2+2x3-2x4 -2)
    @mapping(m, F4, x1+2x2-2x3+4x4 -6)
    @mapping(m, F1, -x3-x4 +2)

    @complementarity(m, F2, x2)
    @complementarity(m, F3, x3)
    @complementarity(m, F1, x1)
    @complementarity(m, F4, x4)

    status = solveMCP(m, convergence_tolerance=1e-8, output="yes", time_limit=3600)
    z = [result_value(x1), result_value(x2), result_value(x3), result_value(x4)]
    @show z
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

println("------------------------------------------------------------------")

@testset "mcp test 3 with PATHSolver" begin

    m = MCPModel()

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    lb = zeros(4)
    ub = Inf*ones(4)

    @variable(m, lb[i] <= myvariablename[i in 1:4] <= ub[i])
    @mapping(m, myconst[i=1:4], sum(M[i,j]*myvariablename[j] for j in 1:4) + q[i])
    @complementarity(m, myconst, myvariablename)

    status = solveMCP(m; convergence_tolerance=1e-8, output="yes", time_limit=3600)
    z = result_value.(myvariablename)
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

println("------------------------------------------------------------------")

@testset "mcp test 4 with PATHSolver" begin

    m = MCPModel()

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    lb = zeros(4)
    ub = Inf*ones(4)

    items = 1:4

    @variable(m, lb[i] <= x[i in items] <= ub[i])
    # @variable(m, x[i in items] >= 0)
    @mapping(m, F[i in items], sum(M[i,j]*x[j] for j in items) + q[i])
    @complementarity(m, F, x)

    status = solveMCP(m; convergence_tolerance=1e-8, output="no", time_limit=3600)
    z = result_value.(x)
    @test isapprox(z.data, [2.8, 0.0, 0.8, 1.2])

end

println("------------------------------------------------------------------")

@testset "mcp test 5 with PATHSolver" begin

    m = nothing
    m = MCPModel()

    lb = zeros(4)
    ub = Inf*ones(4)
    items = 1:4
    @variable(m, lb[i] <= x[i in items] <= ub[i])

    @mapping(m, F1, 3*x[1]^2 + 2*x[1]*x[2] + 2*x[2]^2 + x[3] + 3*x[4] -6)
    @mapping(m, F2, 2*x[1]^2 + x[1] + x[2]^2 + 3*x[3] + 2*x[4] -2)
    @mapping(m, F3, 3*x[1]^2 + x[1]*x[2] + 2*x[2]^2 + 2*x[3] + 3*x[4] -1)
    @mapping(m, F4, x[1]^2 + 3*x[2]^2 + 2*x[3] + 3*x[4] - 3)

    @complementarity(m, F1, x[1])
    @complementarity(m, F2, x[2])
    @complementarity(m, F3, x[3])
    @complementarity(m, F4, x[4])

    set_start_value(x[1], 1.25)
    set_start_value(x[2], 0.)
    set_start_value(x[3], 0.)
    set_start_value(x[4], 0.5)

    status = solveMCP(m)
    @show status

    z = result_value.(x)
    @show z

    @test isapprox(z.data, [1.22474, 0.0, 0.0, 0.5]; atol=1e-4)

    set_start_value(x[1], 1.19)
    set_start_value(x[2], 0.2)
    set_start_value(x[3], 0.1)
    set_start_value(x[4], 0.49)

    status = solveMCP(m)
    @show status

    z = result_value.(x)

    @show z

    @test isapprox(z.data, [1.22474, 0.0, 0.0, 0.5]; atol=1e-4)

end

println("------------------------------------------------------------------")


text:
"""
    DoubleWell(mass=1, ω=1, γ=1, Δ=1)

Two state double well, also called the one-dimensional spin-boson model.
See: [J. Chem. Phys. 150, 244102 (2019)](https://doi.org/10.1063/1.5096276)
"""
Parameters.@with_kw struct DoubleWell{M,W,Y,D} <: DiabaticModel
    mass::M = 1
    ω::W = 1
    γ::Y = 1
    Δ::D = 1
end

NQCModels.ndofs(::DoubleWell) = 1
NQCModels.nstates(::DoubleWell) = 2

function NQCModels.potential(model::DoubleWell, R::Real)

    V0(R) = 0.5 * model.mass * model.ω^2 * R^2

    V₀ = V0(R)
    v = sqrt(2)*model.γ*R
    V11 = V₀ + v
    V22 = V₀ - v
    V12 = model.Δ/2

    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::DoubleWell, R::Real)

    D0(R) = model.mass * model.ω^2 * R

    D₀ = D0(R[1])
    v = sqrt(2)*model.γ
    D11 = D₀ + v
    D22 = D₀ - v
    return Hermitian(SMatrix{2,2}(D11, 0, 0, D22))
end