text:
using Revise, WaltonForaging
## Leaving times threshold rate
gr(size = (600,600))
patch = ["poor", "medium","rich"]
travel = ["long", "short"]
stim = [false, true]
Thr_SimulatedStim  =  DataFrame(patch = String[], travel = String[], stim = Bool[], time = Float64[],
    color = Symbol[], line = Symbol[], group = String[])
for p in patch
    pv = p == "rich" ? 3 : p == "medium" ? 2 : 1
    for t in travel
        tv = t == "short" ? -1 : 1
        l = t == "short" ? :solid : :dash
        tgroup = t =="short" ? "Short Travel" : "Long Travel"
        for s in stim
            if s
                pre = t == "short" ? 2.5 : 1
                sv = 1.3 - pre*(0.4)
                c = :blue
                sgroup = " Stim"
            else
                sv = 0
                c = :black
                sgroup = " No stim"
            end
            value = 3 * pv + tv + sv
            label = tgroup * sgroup
            push!(Thr_SimulatedStim, (patch = p, travel = t, stim = s, time = value, color = c, line = l, group = label))
        end
    end
end
@df Thr_SimulatedStim plot(:patch,:time, group = :group,
    legend = :topleft, color = :color, linestyle = :line,
    ylabel = "Patch leaving time (simulated)",
    xlabel = "Patch Richness", ylims = (0,11))
##
folder = "/Users/dariosarra/Dropbox (Personal)/2022_BBSRC_Foraging/Figures/Illustrative_Serotonin_effects"
savefig(joinpath(folder,"Bhv_trheshold1.pdf"))
##
thr_df1 = combine(groupby(Thr_SimulatedStim,[:patch, :stim]), :time => mean)
change_thr_patch=unstack(thr_df1, :stim, :time_mean)
change_thr_patch[!,:change] = change_thr_patch[:,3] - change_thr_patch[:,2]
change_thr_patch[!,:color] .= :cyan
change_thr_patch[!,:Hypothesis] .= "Global"
@df change_thr_patch bar(:patch, :change,
    ylabel = "Stimulated change in patch leving time", xlabel = "Patch Richness",
    color = :blue, legend = false, ylims = (0,1))
    # @df change_Stim_patch plot!(:patch, :change,linestyle = :dot, color = :grey)
##
savefig(joinpath(folder,"Bhv_trheshold2.pdf"))
##
thr_df2 = combine(groupby(Thr_SimulatedStim,[:travel, :stim]), :time => mean)
change_thr_trav=unstack(thr_df2, :stim, :time_mean)
change_thr_trav[!,:change] = change_thr_trav[:,3] - change_thr_trav[:,2]
change_thr_trav[!,:color] .= :cyan
change_thr_trav[!,:Hypothesis] .= "Global"
@df change_thr_trav bar(:travel, :change,
    ylabel = "Stimulated change in patch leving time",
    xlabel = "Travel distance", color = :blue, legend = false, ylims = (0,1))
##
savefig(joinpath(folder,"Bhv_trheshold3.pdf"))
##

open_html_table(change_thr_trav)
open_html_table(change_lr_trav)

##
patch = ["poor","rich"]
time = [0,1,2,3]
stim = [false, true]
SimulatedPhoto  =  DataFrame(time = Float64[], patch = String[], stim = Bool[], RPE = Float64[],
    color = Symbol[], line = Symbol[], group = String[])
for t in time
    for p in patch
        pv = p == "poor" ? 1.2 : 0.8
        c = p == "poor" ? :brown : :green
        pgroup = p == "poor" ? "Poor patch" : "Rich patch"
        for s in stim
            sv = s ? 0.85 : 1
            l = s ? :dash : :solid
            sgroup = s ? " Stim" : " No stim"
            value = (pv*sv)*t + 1
            label = pgroup * sgroup
            push!(SimulatedPhoto,(time = t, patch = p, stim = s, RPE = value,
                color = c, line = l, group = label))
        end
    end
end
@df SimulatedPhoto plot(:time, :RPE, group = :group,
    legend = :topleft, color = :color, linestyle = :line,
    ylabel = "RPE response at reward (simulated)",
    xlabel = "Elapsed time")
#"u expect reward to have not reduce as much"
##
savefig(joinpath(folder,"DA_trheshold1.pdf"))
##
df1 = combine(groupby(SimulatedPhoto,[:patch, :stim]), :RPE => mean)
change_Stim_patch = unstack(df1, :stim, :RPE_mean)
change_Stim_patch[!,:change] = change_Stim_patch[:,3] - change_Stim_patch[:,2]
@df change_Stim_patch bar(:patch, abs.(:change),
    ylabel = "Stimulated RPE absolute magnitude change", xlabel = "Patch Richness",
    color = :blue, legend = false)
    # @df change_Stim_patch plot!(:patch, :change,linestyle = :dot, color = :grey)
##
savefig(joinpath(folder,"DA_trheshold2.pdf"))
##
Trav_change = vcat(change_thr_trav,change_lr_trav)
Trav_change.color = [x == "long" ? :grey : :lightgrey for x in Trav_change.travel]
@df Trav_change groupedbar(:Hypothesis, :change, group = :travel,
    color = :color, ylims=(0,1), bar_width=0.4)
##
Patch_change = vcat(change_thr_patch,change_lr_patch)
posdic = Dict("poor" => 1, "medium" => 2, "rich" => 3)
colordic = Dict("poor" => :red, "medium" => :orange, "rich" => :yellow)
Patch_change.color = [get(colordic,x,0) for x in Patch_change.patch]
Patch_change[!,:pos] = [get(posdic,x,0) for x in Patch_change.patch]
@df Patch_change groupedbar(:Hypothesis, :change, group = :pos,
    color = :color, ylims=(0,1), bar_width=0.9)
    # xticks = ([1,2,3],["poor","medium", "rich"]))


text:
module CLE

using JuLIP: AbstractAtoms, AbstractCalculator, calculator,
             stress, cell, set_cell!, apply_defm!

using StaticArrays, Einsum

using JuLIPMaterials: Vec3, Mat3, Ten33, Ten43,
         MVec3, MMat3, MTen33, MTen43, ForceConstantMatrix1

# TODO: get rid of this?
const Tensor{T} = Array{T, 4}

# extend `angle` to avoid going via ℂ
Base.angle(x, y) = atan(y, x)

"convert normalised vector to spherical coordinates"
spherical(x) = angle(x[1], x[2]), angle(norm(x[1:2]), x[3])

"convert spherical to euclidean coordinates on the unit sphere"
euclidean(φ, ψ) = Vec3(cos(ψ) * cos(φ), cos(ψ) * sin(φ), sin(ψ))

"given a vector x ∈ ℝ³, return `z0, z1` where `(x/norm(x),z0,z1)` form a right--handed ONB."
function onb3D(x::Vec3{T}) where {T}
   x /= norm(x)
   φ, ψ = spherical(x)
   return Vec3{T}( sin(ψ)*cos(φ), sin(ψ)*sin(φ), -cos(ψ) ),
          Vec3{T}(-sin(φ), cos(φ), zero(T) )
end

"""
* `elastic_moduli(at::AbstractAtoms)`
* `elastic_moduli(calc::AbstractCalculator, at::AbstractAtoms)`
* `elastic_moduli(C::Matrix)` : convert Voigt moduli to 4th order tensor

computes the 3 x 3 x 3 x 3 elastic moduli tensor

*Notes:* this is a naive implementation that does not exploit
any symmetries at all; this means it performs 9 centered finite-differences
on the stress. The error should be in the range 1e-10
"""
elastic_moduli(at::AbstractAtoms) = elastic_moduli(calculator(at), at)

function elastic_moduli(calc::AbstractCalculator, at::AbstractAtoms)
   F0 = cell(at)' |> Matrix
   Ih = Matrix(1.0*I, 3,3)
   h = eps()^(1/3)
   C = zeros(3,3,3,3)
   for i = 1:3, a = 1:3
      Ih[i,a] += h
      apply_defm!(at, Ih)
      Sp = stress(calc, at)
      Ih[i,a] -= 2*h
      apply_defm!(at, inv(Ih))
      Sm = stress(calc, at)
      C[i, a, :, :] = (Sp - Sm) / (2*h)
      Ih[i,a] += h
   end
   # symmetrise it - major symmetries C_{iajb} = C_{jbia}
   for i = 1:3, a = 1:3, j=1:3, b=1:3
      C[i,a,j,b] = C[j,b,i,a] = 0.5 * (C[i,a,j,b] + C[j,b,i,a])
   end
   # minor symmetries - C_{iajb} = C_{iabj}
   for i = 1:3, a = 1:3, j=1:3, b=1:3
      C[i,a,j,b] = C[i,a,b,j] = 0.5 * (C[i,a,j,b] + C[i,a,b,j])
   end
   return C
end

function elastic_moduli(FCM::ForceConstantMatrix1)
   ℂ = zeros(3,3,3,3)
   for i=1:length(FCM.R)
       ρ = FCM.R[i]
       Hmat = FCM.H[i]
       @einsum ℂ[a,b,c,d] = Hmat[a,c] * ρ[b] * ρ[d]
   end
   return ℂ
end

"""
`voigt_moduli`: compute elastic moduli in the format of Voigt moduli.

Methods:
* `voigt_moduli(at)`
* `voigt_moduli(calc, at)`
* `voigt_moduli(C)`
"""
voigt_moduli(at::AbstractAtoms) = voigt_moduli(calculator(at), at)

voigt_moduli(calc::AbstractCalculator, at::AbstractAtoms) =
   voigt_moduli(elastic_moduli(calc, at))

const voigtinds = [1, 5, 9, 4, 7, 8]

voigt_moduli(C::Array{T,4}) where {T} = reshape(C, 9, 9)[voigtinds, voigtinds]

function elastic_moduli(Cv::AbstractMatrix{T}) where {T}
   @assert size(Cv) == (6,6)
   C = zeros(T, 9,9)
   C[voigtinds, voigtinds] = Cv
   C = reshape(C, 3,3,3,3)
   # now apply all the symmetries to recover C
   for i = 1:3, a = 1:3, j = 1:3, b = 1:3
      if C[i,a,j,b] != 0
         C[a,i,j,b] = C[i,a,b,j] = C[a,i,b,j] = C[j,b,i,a] =
            C[b,j,i,a] = C[j,b,a,i] = C[b,j,a,i] = C[i,a,j,b]
      end
   end
   return C
end

"""
`isotropic_moduli(λ, μ)`: compute 4th order tensor of elastic moduli
corresponding to the Lame parameters λ, μ.
"""
function isotropic_moduli(λ, μ)
   K = λ + μ * 2 / 3
   C = @SArray [ K * I[i,j] * I[k,l] + μ * (I[i,k]*I[j,l] + I[i,l]*I[j,k] - 2/3*I[i,j]*I[k,l])
         for i = 1:3, j = 1:3, k = 1:3, l = 1:3 ]
   return C
end

"""
`isotropic_moduli2D(λ, μ)`: compute 4th order tensor of elastic moduli
(in 2D planar elasticity) corresponding to the Lame parameters λ, μ.
"""
function isotropic_moduli2D(λ, μ)
   K = λ + μ * 2 / 3
   C = @SArray [ K * I[i,j] * I[k,l] + μ * (I[i,k]*I[j,l] + I[i,l]*I[j,k] - 2/3*I[i,j]*I[k,l])
         for i = 1:2, j = 1:2, k = 1:2, l = 1:2 ]
   return C
end

function zener_anisotropy_index(C::Tensor)
    Cv = voigt_moduli(C)
    A = 2*Cv[4,4]/(Cv[1,1] - Cv[1,2])
    return A
end

zener_anisotropy_index(at::AbstractAtoms) =
                  zener_anisotropy_index(elastic_moduli(at))

"""
compute the Lame parameters for an elasticity tensor C or throw
an error if the material is not isotropic.
"""
function lame_parameters(C::Tensor; aniso_threshold=1e-3)
    A = zener_anisotropy_index(C)
    @assert abs(A - 1.0) < aniso_threshold
    Cv = voigt_moduli(C)
    μ = Cv[1,2]
    λ = Cv[4,4]
    return λ, μ
end
lame_parameters(at::AbstractAtoms) = lame_parameters(elastic_moduli(at))

poisson_ratio(λ, μ) = 0.5 * λ / (λ + μ)
poisson_ratio(C::Tensor) = poisson_ratio(lame_parameters(C)...)
poisson_ratio(at::AbstractAtoms) = poisson_ratio(elastic_moduli(at))

function youngs_modulus(λ, μ)
    ν = poisson_ratio(λ, μ)
    return λ*(1 + ν)*(1 - 2ν) / ν
end

youngs_modulus(C::Tensor) = youngs_modulus(lame_parameters(C)...)
youngs_modulus(at::AbstractAtoms) = youngs_modulus(elastic_moduli(at))

"""
check whether the elasticity tensor is isotropic; return true/false
"""
function is_isotropic(C::Tensor)
   try
      lame_parameters(C)
      return true
   catch
      return false
   end
end

# """
# `module GreensFunctions`
#
# Implements some CLE Green's functions, both in analytic form or
# semi-analytic using the formulas from BBS79.
# """

include("cle_greenfunctions.jl")

include("cle_gfcorrectors3D.jl")

include("cle_dislocations.jl")

end

ext: jl


text:
# Make a DAG (Directed Acyclic Graph) by storing references to each symbol

"""
Structure representing a DAG.
Maintains a `symbol_map` giving the currently-known symbols and the corresponding
vertex number in the graph.
"""
struct DirectedAcyclicGraph <: LabelledDiGraph
    g::DiGraph
    labels::Vector{Any}
    symbol_map::Dict{Symbol, Int}
end

DirectedAcyclicGraph() = DirectedAcyclicGraph(DiGraph(), Symbol[], Dict())

"""
Adds a symbol to the DAG if it doesn't already exist.
Returns the vertex number
# Make numbers unique:
"""
function add_symbol!(dag::DirectedAcyclicGraph, s)  # number
    vertex = add_numbered_vertex!(dag.g)
    push!(dag.labels, s)
    return vertex
end

function lookup!(dag::DirectedAcyclicGraph, s)
    add_symbol!(dag, s)
end

"""
Look up a symbol to see if it has already been seen.
"""
function lookup!(dag::DirectedAcyclicGraph, s::Symbol)
    if haskey(dag.symbol_map, s)
        return dag.symbol_map[s]

    else  # make new one:
        vertex = add_numbered_vertex!(dag.g)
        push!(dag.labels, s)
        dag.symbol_map[s] = vertex
        return vertex
    end
end

make_dag!(dag::DirectedAcyclicGraph, s) = lookup!(dag, s)

"""
Update a Directed Acyclic Graph with the result of traversing the given `Expr`ession.
"""
function make_dag!(dag::DirectedAcyclicGraph, ex::Expr)

    local top

    if ex.head == :block
        for arg in ex.args
            make_dag!(dag, arg)
        end
        return -1

    elseif ex.head == :(=)  # treat assignment as just giving pointers to the tree
        local_var = ex.args[1]

        top = make_dag!(dag, ex.args[2])

        dag.symbol_map[local_var] = top  # add an alias to the corresponding tree node

        return top

    end

    where_start = 1  # which argument to start with

    if ex.head == :call
        f = ex.args[1]   # the function name
        top = add_symbol!(dag, f)

        where_start = 2   # drop "call" from tree

    else
        @show ex.head
        top = add_symbol!(dag, ex.head)
    end

    # @show top

    for arg in ex.args[where_start:end]

        # @show arg, typeof(arg)

        if isa(arg, Expr)

            child = make_dag!(dag, arg)
            # @show "Expr", top, child
            add_edge!(dag.g, top, child)

        else
            child = lookup!(dag, arg)
            # @show top, child
            add_edge!(dag.g, top, child)

        end
    end

    return top

end

"""
Make a Directed Acyclic Graph (DAG) from a Julia expression.
"""
function make_dag(ex::Expr)

    dag = DirectedAcyclicGraph()

    make_dag!(dag, MacroTools.striplines(ex))

    return dag

end

"""
Make a Directed Acyclic Graph (DAG) from a Julia expression.
"""
macro dag(ex::Expr)
    make_dag(ex)
end

"""
Perform common subexpression elimination on a Julia `Expr`ession,
and make a Directed Acyclic Graph (DAG) of the result.
"""
macro dag_cse(ex::Expr)
    make_dag(cse(ex))  # common subexpression elimination
end

import Base.show
function show(io::IO, mime::MIME"image/svg+xml", dag::DirectedAcyclicGraph)
    p = tikz_representation(dag)  # TikzPicture object
    show(io, mime, p)
end


text:
function lookup_gene_ec_mapping_record(kegg_gene_id::String)::Union{Set{String},Nothing}

    # initialize -
    record_set = Set{String}()

    # call the KEGG API to get the ec number -
    list_of_ec_numbers = get_ec_number_for_gene(kegg_gene_id)
    if list_of_ec_numbers !== nothing

        # create a record, add to the mapping buffer -
        for ec_number in list_of_ec_numbers
            push!(record_set,ec_number)
        end

        # return -
        return record_set
    end

    return nothing
end

function load_cobra_model_file(path_to_cobra_mat_file::String, model_name::String)::Dict{String,Any}

    # open the cobra file -
    file = matopen(path_to_cobra_mat_file)
    cobra_dictionary = read(file,model_name)
    close(file)

    # return -
    return cobra_dictionary
end

function export_reaction_mapping_file(cobra_dictionary::Dict{String,Any}, path_to_mapping_file::String)

    # Check the dictionary -
    is_cobra_dictionary_ok(cobra_dictionary)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// ------------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// Reaction->reaction name map")
    +(mapping_buffer,"// record: reacton_tag=reaction_name")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// reaction_name: from the rxnNames field of the COBRA mat file")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// The order of the reactions in this file is used to order the")
    +(mapping_buffer,"// columns of the stoichiometric matrix.")
    +(mapping_buffer,"// ------------------------------------------------------------ //")

    list_of_reactions = cobra_dictionary["rxns"]
    list_of_reaction_names = cobra_dictionary["rxnNames"]

    # process each gene -
    for (index,reaction_tag) in enumerate(list_of_reactions)

        # init -
        line = ""

        # get content -
        reaction_name = list_of_reaction_names[index]

        # build line -
        line *= "$(reaction_tag)=$(reaction_name)"

        # cache -
        +(mapping_buffer, line)
    end

    # Write out the mapping file
    write_file_to_path(path_to_mapping_file,mapping_buffer)
end

function export_gene_order_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # initalize -
    mapping_buffer = String[]
    list_of_genes = cobra_dictionary["genes"]

    # process each gene -
    for gene_id in list_of_genes

        # make the KEGG gene id -
        kegg_gene_id = "$(kegg_organism_code):$(gene_id)"

        # cache -
        +(mapping_buffer, kegg_gene_id)
    end

    # Write out the vff file -
    open("$(path_to_mapping_file)", "w") do f

        for line_item in mapping_buffer
            write(f,"$(line_item)\n")
        end
    end
end

function export_reaction_tag_to_gene_mapping_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// --------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"")
    +(mapping_buffer,"// Reaction->gene mapping - ")
    +(mapping_buffer,"// record: reacton_tag={gene_id}")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// gene_id: (kegg organism code):(gene location code)")
    +(mapping_buffer,"// --------------------------------------------------------- //")

    # get the list of reaction tags, and genes from the cobra dictionary -
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_genes = cobra_dictionary["genes"]

    # get the rxn gene mapping matrix -
    RGM = Matrix(cobra_dictionary["rxnGeneMat"])

    # What is the size of the system?
    (number_of_reactions, number_of_genes) = size(RGM)

    # main loop -
    for reaction_index = 1:number_of_reactions

        # what is the tag for this reaction?
        reaction_tag = list_of_reaction_tags[reaction_index]

        # initialize -
        record = ""
        record *= "$(reaction_tag)="

        # does this reaction have associated genes?
        idx_genes = findall(x->x==1,RGM[reaction_index,:])
        local_gene_set = Set{String}()
        if (!isempty(idx_genes))

            # ok, we have genes, grab them -
            gene_id_array = list_of_genes[idx_genes]

            # process each gene -
            for gene_id in gene_id_array

                # make the KEGG gene id -
                kegg_gene_id = gene_id
                if (occursin(".",gene_id) == true)
                    # need to cutoff the trailing *.1
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id[1:end-2])"
                else
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id)"
                end

                # cache -
                push!(local_gene_set,kegg_gene_id)
            end
        end

        if (isempty(local_gene_set) == false)

            # make an ec record -
            gene_record = ""
            for gene_number in local_gene_set
                gene_record*="$(gene_number),"
            end

            # remove trailing ,
            gene_record = gene_record[1:end-1]

            # add the record -
            record*="$(gene_record)"

            # push into buffer and go around again -
            +(mapping_buffer,record)
        end
    end

    # Write out the file -
    write_file_to_path(path_to_mapping_file, mapping_buffer)
end

function export_reaction_tag_to_ec_mapping_file(cobra_dictionary::Dict{String,Any}, kegg_organism_code::String, path_to_mapping_file::String)

    # initalize -
    mapping_buffer = String[]

    # add header text to mapping file -
    filename = splitdir(path_to_mapping_file)[2]
    +(mapping_buffer,"// --------------------------------------------------------- //")
    +(mapping_buffer,"// $(filename)")
    +(mapping_buffer,"// GENERATED BY: CBModelTools")
    +(mapping_buffer,"// GENERATED ON: $(Dates.now())")
    +(mapping_buffer,"// SOURCE: https://github.com/varnerlab/CBModelTools")
    +(mapping_buffer,"//")
    +(mapping_buffer,"// Reaction->ecnumber mapping - ")
    +(mapping_buffer,"// record: reacton_tag={ecnumber}")
    +(mapping_buffer,"// reaction_tag: from the rxns field of the COBRA mat file")
    +(mapping_buffer,"// ec_number: possible ecnumber estimated from KEGG")
    +(mapping_buffer,"// --------------------------------------------------------- //")

    # get the list of reaction tags, and genes from the cobra dictionary -
    list_of_reaction_tags = cobra_dictionary["rxns"]
    list_of_genes = cobra_dictionary["genes"]

    # get the rxn gene mapping matrix -
    RGM = Matrix(cobra_dictionary["rxnGeneMat"])

    # What is the size of the system?
    (number_of_reactions, number_of_genes) = size(RGM)

    # main loop -
    for reaction_index = 1:number_of_reactions

        # what is the tag for this reaction?
        reaction_tag = list_of_reaction_tags[reaction_index]

        # user message -
        msg = "Starting $(reaction_tag) ($(reaction_index) of $(number_of_reactions))"

        # initialize -
        record = ""
        record *= "$(reaction_tag)="

        # does this reaction have associated genes?
        idx_genes = findall(x->x==1,RGM[reaction_index,:])

        # init -
        ec_record_set = Set{String}()
        if (!isempty(idx_genes))

            # ok, we have genes, grab them -
            gene_id_array = list_of_genes[idx_genes]

            # process each gene -
            for gene_id in gene_id_array

                # make the KEGG gene id -
                kegg_gene_id = gene_id
                if (occursin(".",gene_id) == true)
                    # need to cutoff the trailing *.1
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id[1:end-2])"
                else
                    kegg_gene_id = "$(kegg_organism_code):$(gene_id)"
                end

                # lookup -
                local_ec_record_set = lookup_gene_ec_mapping_record(kegg_gene_id)
                if local_ec_record_set != nothing

                    # push into ec_record_set -
                    for ec_number in local_ec_record_set
                        push!(ec_record_set,ec_number)
                    end
                end
            end
        end

        if (isempty(ec_record_set) == false)

            # make an ec record -
            ec_record = ""
            for ec_number in ec_record_set
                ec_record*="$(ec_number),"
            end

            # remove trailing ,
            ec_record = ec_record[1:end-1]

            # add the record -
            record*="$(ec_record)"

            # push into buffer and go around again -
            +(mapping_buffer,record)
        end
    end

    # Write out the file -
    write_file_to_path(path_to_mapping_file,mapping_buffer)
end

# --- PUBLIC METHODS ------------------------------------------------------------------------------------------ #
function extract_files_from_cobra_model(path_to_cobra_mat_file::String, path_to_exported_files::String, modelName::String;
    kegg_organism_code_symbol::Symbol = :eco)

    # checks - do the paths exist?
    is_file_path_ok(path_to_cobra_mat_file)
    is_dir_path_ok(path_to_exported_files)

    # load the cobra dictionary -
    cobra_dictionary = load_cobra_model_file(path_to_cobra_mat_file, modelName)

    # check the cobra dictionary - is it legit?
    is_cobra_dictionary_ok(cobra_dictionary)

    # ok, looks ok if we get here - write out a bunch of files that we'll need -
    kegg_organism_code = String(kegg_organism_code_symbol)

    # reaction mapping -
    path_to_reaction_mapping_file = joinpath(path_to_exported_files,"ReactionNameMap.dat")
    export_reaction_mapping_file(cobra_dictionary,path_to_reaction_mapping_file)

    # gene order -
    path_to_gene_order_file = joinpath(path_to_exported_files,"GeneOrder.dat")
    export_gene_order_file(cobra_dictionary, kegg_organism_code, path_to_gene_order_file)

    # reaction tags to gene mapping -
    path_reaction_gene_mapping_file = joinpath(path_to_exported_files,"ReactionGeneMap.dat")
    export_reaction_tag_to_gene_mapping_file(cobra_dictionary, kegg_organism_code, path_reaction_gene_mapping_file)

    # reaction tag to ec mapping -
    path_reaction_to_ec_number_mapping = joinpath(path_to_exported_files, "ReactionECNumberMap.dat")
    export_reaction_tag_to_ec_mapping_file(cobra_dictionary, kegg_organism_code, path_reaction_to_ec_number_mapping)
end
# ------------------------------------------------------------------------------------------------------------- #

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
module FluxModels

using DocStringExtensions
using Flux
using Flux: @functor
using Parameters
using ModelUtils

abstract type ModuleSpec end

include("./OutputSizes.jl")
using .OutputSizes

include("./activations.jl")
include("./layers.jl")
include("./blocks.jl")
include("./heads.jl")

include("./efficientnet.jl")
include("./mobilenetv3.jl")
include("./xresnet.jl")

export
    outputsize,
    SqueezeExcitation,
    PixelShuffle,

    efficientnetb0,
    mobilenetv3_small,
    mobilenetv3_large,
    xresnet18,
    xresnet50

end # module


text:
module TestSplice

include("preamble.jl")

@testset begin
    @test splice!!((1, 2, 3), 2) === ((1, 3), 2)
    @test splice!!([1, 2, 3], 2) == ([1, 3], 2)
    @test splice!!([1, 2, 3], 2, 20) == ([1, 20, 3], 2)
    @test splice!!((1, 2, 3), 2, 20) === ((1, 20, 3), 2)
    @test splice!!([1, 2, 3], 2:3) == ([1], [2, 3])
    @test splice!!([1, 2, 3], 2:3, [20, 30]) == ([1, 20, 30], [2, 3])
end

end  # module


text:
function nested_containers()::Bool
    CODES = Dict("KEY" => [1, 3])
    return 1 ∈ CODES["KEY"]
end

if abspath(PROGRAM_FILE) == @__FILE__
    if nested_containers()
        println("OK")
    end
end


text:
ScriptVersion = "0.0.1alpha"
# Very unstable alpha Version!

# Locally include and use packages
if isdefined(@__MODULE__, :LookUpTable)
    LookUpTable isa Module || error("LookUpTable is present and it is not a Module")
else
    include("LookUpTable.jl")
end
using .LookUpTable

if isdefined(@__MODULE__, :FileCodeOutput)
    FileCodeOutput isa Module || error("FileCodeOutput is present and it is not a Module")
else
    include("FileCodeOutput.jl")
end
using .FileCodeOutput

if isdefined(@__MODULE__, :JSONPhraser)
    JSONPhraser isa Module || error("JSONPhraser is present and it is not a Module")
else
    include("JSONPhraser.jl")
end
using .JSONPhraser

#=======================User configured Section==========================#
# Your Version Number
YourVersion= "0.1"

# C(suitable for C and C++),Julia
codeStyle = "Julia"

#=======================User configured Section==========================#

M1=ReliAffixECDConverter();
M2=ReliLevelECDConverter();

open("output.txt", "w") do file
    PrintScriptVersion(file, ScriptVersion, YourVersion;codeStyle=codeStyle);
    Print3DMatrix(file, M1 ; varName="reliSubstat",codeStyle=codeStyle);
    Print3DMatrix(file, M2 ; varName="reliMainstat",codeStyle=codeStyle);
end


text:
function triangle(adj_mtx_file)
    if ~isfile( adj_mtx_file )
        println("unable to open input file");
        return (-1);
    end

    # load input data
    t_read_adj=@elapsed ii = readdlm( adj_mtx_file, '\t', Int64);
    println("adjacency matrix read time : ", t_read_adj);

    t_create_adj=@elapsed A = sparse( ii[:,1], ii[:,2], ii[:,3] );
    println("sparse adj. matrix creation time : ", t_create_adj);

    tic;
    nt = sum( A^2 .*  A )/6;
    t = toc;
    println("triangle count time : ", t);
    println("number of triangles : ", nt);

    return(nt)
end

#######################################################
# Graph Challenge Benchmark
# Developer : Dr. Siddharth Samsi (ssamsi@mit.edu)
#
# MIT
########################################################
# (c) <2017> Massachusetts Institute of Technology
########################################################

max_stars_repo_name: bkmgit/GraphChallenge

max_issues_repo_name: bkmgit/GraphChallenge

max_forks_repo_name: bkmgit/GraphChallenge


text:
using NumericalMethodsforEngineers, IterativeSolvers

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

upper = chol(a)
lower = upper'
y = lower \ b
c = upper \ y

tol = 1.0e-10
iters = 20
x,ch = cg(a, b;tol=tol,maxiter=iters);

println("\nSolution Vector: $x")
println("Converged: $(ch.isconverged)")
println("Number of iters: $(ch.mvps)")

@test round(x, 9) == round(c, 9)


text:
"""
    function get_thermal_slack(sys::System)

Gets the slack generator from a system and returns the generator, location in the vector
of generators, and name.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)

# Throws
- `NotFoundError`:                          No Slack bus found in the system
"""
function get_thermal_slack(sys::System)
    # Define elements
    gen_loc = 0
    ngen_slack = 0 #number of gen slacks
    has_slack = false
    gen_thermal_slack = 0
    gen_thermal_slack_loc = 0
    gen_thermal_slack_name = "Empty"

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)

    # Identify Slack Generator and save its information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus)== BusTypes.REF
            has_slack = true
            gen_thermal_slack = gen_thermal
            gen_thermal_slack_loc = gen_loc
            gen_thermal_slack_name = gen_thermal.name
        end
    end
    if !has_slack
        error("NotFoundError: No Slack found in the system")
    end
    # Get the Slack Generator
    gen_thermal_slack = get_component(ThermalStandard, sys, gen_thermal_slack_name)

    return gen_thermal_slack,gen_thermal_slack_loc,gen_thermal_slack_name
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
function _n_node!(num::Int, p::AbstractNode, q::AbstractNode)
	for link in q.links
		link.to == p && continue
		num = _n_node!(num, q, link.to)
	end

	return num + 1
end

"""
	n_node(tree)
	n_node(p [, q])

Count the nodes in a tree or subtree denoted by the nodes `p` and `q`.
"""
n_node(p::AbstractNode, q::AbstractNode) = _n_node!(0, p, q)
n_node(p::AbstractNode) = n_node(p, p)
n_node(tree::AbstractTree) = n_node(tree.origin)

function _n_tip!(num::Int, p::AbstractNode, q::AbstractNode)
	for link in q.links
		link.to == p && continue
		num = _n_tip!(num, q, link.to)
	end

	return ifelse(istip(q), num + 1, num)
end

"""
	n_tip(tree)
	n_tip(p [, q])

Count the tip nodes in a tree or subtree denoted by the nodes `p` and `q`.
"""
n_tip(p::AbstractNode, q::AbstractNode) = _n_tip!(0, p, q)
n_tip(p::AbstractNode) = n_tip(p, p)
n_tip(tree::AbstractTree) = n_tip(tree.origin)

"""
	n_branch(tree)
	n_branch(p [, q])

Count the branches in a `tree` or subtree denoted by the nodes `p` and `q`.
"""
n_branch(p::AbstractNode, q::AbstractNode) = n_node(p, q) - 1
n_branch(p::AbstractNode) = n_branch(p, p)
n_branch(tree::AbstractTree) = n_branch(tree.origin)

function _n_species!(num::Int, p::AbstractNode, q::AbstractNode)
	for link in q.links
		link.to == p && continue
		num = _n_species!(num, q, link.to)
	end

	return ifelse(q.species > 0, num + 1, num)
end

"""
	n_species(tree)
	n_species(p [, q])

Count the species in a `tree` or subtree denoted by the nodes `p` and `q`.
"""
n_species(p::AbstractNode, q::AbstractNode) = _n_species!(0, p, q)
n_species(p::AbstractNode) = n_species(p, p)
n_species(tree::AbstractTree) = n_species(tree.origin)

function _nodelabels!(labs::Vector{String}, p::AbstractNode, q::AbstractNode)
	push!(labs, label(q))
	for link in q.links
		link.to == p && continue
		_nodelabels!(labs, q, link.to)
	end
	return labs
end

"""
	nodelabels(tree)
	nodelabels(p [, q])

Get the labels all the nodes (in preorder) in a `tree` or subtree denoted by the nodes `p` and `q`.
"""
nodelabels(p::AbstractNode, q::AbstractNode) = _nodelabels!(String[], p, q)
nodelabels(p::AbstractNode) = nodelabels(p, p)
nodelabels(tree::AbstractTree) = nodelabels(tree.origin)

function _tiplabels!(labs::Vector{String}, p::AbstractNode, q::AbstractNode)
	istip(q) && push!(labs, label(q))
	for link in q.links
		link.to == p && continue
		_tiplabels!(labs, q, link.to)
	end

	return labs
end

"""
	tiplabels(tree)
	tiplabels(p [, q])

Get the labels all the tip nodes (in preorder) in a `tree` or subtree denoted by the nodes `p` and `q`.
"""
tiplabels(p::AbstractNode, q::AbstractNode) = _tiplabels!(String[], p, q)
tiplabels(p::AbstractNode) = tiplabels(p, p)
tiplabels(tree::AbstractTree) = tiplabels(tree.origin)

function _bipartitions(list::Vector{Bipartition}, p::AbstractNode, q::AbstractNode)
	for link in q.links
		link.to == p && continue
		push!(list, link.branch.bipart)
		_bipartitions(list, q, link.to)
	end

	return list
end

function _tips!(list::Vector{T}, p::T, q::T) where T <: AbstractNode
	istip(q) && push!(list, q)
	for link in q.links
		link.to == p && continue
		_tips!(list, q, link.to)
	end

	return list
end

"""
	tips(tree)
	tips(p [, q])

Get a list of all the tip nodes (in preorder) in a `tree` or subtree denoted by the nodes `p` and `q`.
"""
tips(p::T, q::T) where T <: AbstractNode = _tips!(T[], p, q)
tips(p::AbstractNode) = tips(p, p)
tips(tree::AbstractTree) = tips(tree.origin)

"""
	bipartitions(tree)
	bipartitions(p [, q])

Get all the bipartitions (in preorder) of a `tree` or subtree denoted by the nodes `p` and `q`.
"""
bipartitions(p::AbstractNode, q::AbstractNode) =
	_bipartitions(Bipartition[], p, q)
bipartitions(p::AbstractNode) = bipartitions(p, p)
bipartitions(tree::AbstractTree) = bipartitions(tree.origin)

function _getspecies!(list::Vector{Int}, p::AbstractNode, q::AbstractNode)
	push!(list, q.species)
	for link in q.links
		link.to == p && continue
		_getspecies!(list, q, link.to)
	end

	return list
end

"""
	getspecies(tree)
	getspecies(p , q)

Get the species numbers of all the nodes (in preorder) of a `tree` or subtree denoted by the nodes `p` and `q`.
"""
getspecies(p::AbstractNode, q::AbstractNode) = _getspecies!(Int[], p, q)
getspecies(tree::AbstractTree) = getspecies(tree.origin, tree.origin)

function _find_nonsplitting!(list::NodeVector, p::AbstractNode, q::AbstractNode)
	n_neighbour(q) == 2 && push!(list, q)
	for link in q.links
		link.to == p && continue
		list = _find_nonsplitting!(list, q, link.to)
	end

	return list
end

"""
	find_nonsplitting(tree)
	find_nonsplitting(p [, q])

Get a list of non-splitting internal nodes in a `tree` or subtree denoted by the nodes `p` and `q`. Non-splitting internal nodes are nodes with only two neighbours.
"""
find_nonsplitting(p::AbstractNode, q::AbstractNode) =
	_find_nonsplitting!(Vector{typeof(p)}(), p, q)
find_nonsplitting(p::AbstractNode) = find_nonsplitting(p, p)
find_nonsplitting(tree::AbstractTree) = find_nonsplitting(tree.origin)


text:
module BHAtp

using DataFrames
using Interpolations
using Distributed
using CSV

# Types and derived types

include("Types.jl")
include("Parameters.jl")
include("Dataframes.jl")

# Exported methods

include("Bha.jl")
include("TheoreticalPerformance.jl")
include("ShowFunctions.jl")

# Base level methods

include("InertiaForces.jl")
include("CurvatureForces.jl")
include("InitialRelease.jl")
include("ExceedanceAdjustment.jl")
include("FinalRelease.jl")
include("Fem.jl")
include("TpRunSetup.jl")
include("Finalinit.jl")
include("CreateMaterialTable.jl")
include("CreateMediaTable.jl")
include("CheckInput.jl")
include("CreateMesh.jl")
include("Interpolate.jl")
include("Xyinit.jl")
include("DalphaInit.jl")
include("Stiffness.jl")
include("GeometricMatrix.jl")
include("Weightforces.jl")
include("Endforces.jl")

export
  BHAJ,
  tp!,
  show_solution,
  show_tp,
  create_final_tp_df,
  create_element_df,
  create_node_df

end # module


text:
"""
This function is a wrapper for running expected utility and first year profit calculations
for making investment decisions.
"""
function calculate_expected_utility(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               risk_preference::R,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes) where {P <: Project{<: BuildPhase}, R <: RiskPreference}

    finance_data = get_finance_data(project)

    update_project_utility!(project,
                              scenario_data,
                              market_prices,
                              risk_preference,
                              rep_hour_weight,
                              iteration_year,
                              yearly_horizon,
                              queue_cost,
                              capacity_forward_years,
                              solver)

        annual_revenue = 0.0

        for scenario in scenario_data
            annual_revenue += get_probability(scenario) *
                             sum(get_scenario_profit(finance_data)[get_name(scenario)][iteration_year][:, get_construction_year(project)])
        end

        annual_profit = annual_revenue - get_fixed_OM_cost(finance_data)
        project_utility = get_expected_utility(finance_data)[iteration_year]

        return annual_profit, project_utility

end

"""
This function adds options of each technology type (also incorporating capital cost multiplier)
to the list of profitable options which can be sent to the queue.
Returns the updated list of profitable options.
"""
function add_profitable_option(projects::Vector{Project},
                                  max_new_options::Dict{String, Int64},
                                  profitable_options::Vector{Project},
                                  scenario_data::Vector{Scenario},
                                  market_prices::MarketPrices,
                                  risk_preference::R,
                                  capital_cost_multiplier::Float64,
                                  investor_name::String,
                                  iteration_year::Int64,
                                  yearly_horizon::Int64,
                                  simulation_years::Int64,
                                  rep_hour_weight::Vector{Float64},
                                  capacity_forward_years::Int64,
                                  solver::JuMP.MOI.OptimizerWithAttributes) where {R <: RiskPreference}

    if length(projects) >= 1
        if  get_construction_year(projects[1]) <= simulation_years

            total_counter = 1
            counter_by_zone = AxisArrays.AxisArray(ones(length(projects)), get_zone.(get_tech.(projects)))
            profitable_type_options = Project[]
            base_inv_cost = get_effective_investment_cost(get_finance_data(projects[1]))
            inv_cost_array_length = length(get_investment_cost(get_finance_data(projects[1])))
            queue_cost = get_queue_cost(get_finance_data(projects[1]))

            for project in projects

                annual_profit, project_utility = calculate_expected_utility(project,
                                                                    scenario_data,
                                                                    market_prices,
                                                                    risk_preference,
                                                                    rep_hour_weight,
                                                                    iteration_year,
                                                                    yearly_horizon,
                                                                    queue_cost,
                                                                    capacity_forward_years,
                                                                    solver)

                if project_utility >= 0 && annual_profit >= -1 && counter_by_zone[get_zone(get_tech(project))] <= max_new_options[get_name(project)]
                    new_option = deepcopy(project)
                    push!(profitable_type_options, new_option)
                end

            end

            while length(profitable_type_options) > 0
                sort!(profitable_type_options, by = x -> -get_expected_utility(get_finance_data(x))[iteration_year])
                most_profitable_option = deepcopy(profitable_type_options[1])
                zone = get_zone(get_tech(most_profitable_option))
                set_name!(most_profitable_option, "$(investor_name)_$(get_type(get_tech(most_profitable_option)))_$(zone)_year_$(iteration_year)_$(Int(counter_by_zone[zone]))")
                set_investment_cost!(most_profitable_option, fill(get_effective_investment_cost(get_finance_data(most_profitable_option)), inv_cost_array_length))
                push!(profitable_options, most_profitable_option)
                counter_by_zone[zone] += 1

                for option in profitable_type_options
                    set_effective_investment_cost!(option, base_inv_cost * (1 + capital_cost_multiplier * total_counter))
                end

                for (idx, option) in enumerate(profitable_type_options)
                    annual_profit, project_utility = calculate_expected_utility(option,
                                                                        scenario_data,
                                                                        market_prices,
                                                                        risk_preference,
                                                                        rep_hour_weight,
                                                                        iteration_year,
                                                                        yearly_horizon,
                                                                        queue_cost,
                                                                        capacity_forward_years,
                                                                        solver)
                    if project_utility < 0 || annual_profit < -1 || counter_by_zone[get_zone(get_tech(option))] > max_new_options[get_name(option)]
                        deleteat!(profitable_type_options, idx)
                    end
                end

                total_counter += 1

            end

        end
    end

    return profitable_options
end
"""
This function makes the investment decisions each year by sending a top N profitable options to the queue,
where N is the maximum number of new projects an investor can invest in within a year.
Returns nothing.
"""
function make_investments!(investor::Investor,
                           max_new_options::Dict{String, Int64},
                           iteration_year::Int64,
                           yearly_horizon::Int64,
                           simulation_years::Int64,
                           capacity_forward_years::Int64,
                           solver::JuMP.MOI.OptimizerWithAttributes)

    scenario_data = get_scenario_data(get_forecast(investor))
    market_prices = get_market_prices(investor)
    risk_preference = get_risk_preference(investor)
    projects = get_projects(investor)

    option_projects = get_options(investor)

    types = unique(get_type.(get_tech.(option_projects)))
    projects_by_type = Dict(type => Project[] for type in types)

    for type in types
        for project in option_projects
            if get_type(get_tech(project)) == type
                push!(projects_by_type[type], project)
            end
        end
    end

    profitable_options = Project[]
    for type in types
        profitable_options = add_profitable_option(projects_by_type[type],
                                                    max_new_options,
                                                    profitable_options,
                                                    scenario_data,
                                                    market_prices,
                                                    risk_preference,
                                                    get_cap_cost_multiplier(investor),
                                                    get_name(investor),
                                                    iteration_year,
                                                    yearly_horizon,
                                                    simulation_years,
                                                    get_rep_hour_weight(investor),
                                                    capacity_forward_years,
                                                    solver)
    end

        sort!(profitable_options, by = x -> -get_expected_utility(get_finance_data(x))[iteration_year])  # sort in descending order of project utility
        for i in 1:min(length(profitable_options), get_max_annual_projects(investor))
            push!(projects, convert(Project{Queue}, profitable_options[i]))
        end

    return
end


text:
n = rand(1:10)

A = rand(n, n)

r = polarfact(A, alg = :halley);

# Test unitary matrix U

U = r.U
H = r.H

@test_approx_eq_eps U'*U eye(n) 1e-7

# Test Hermitian positive semifefinite matrix H

@test issym(H)

for i in eigvals(H)
    @test i >= 0.
end

@test_approx_eq A U*H

println("Halley method passed test...")

##########################################################

m = rand(1:10)

B = rand(m, m)

r2 = polarfact(B, alg = :qdwh);

# Test unitary matrix U

U2 = r2.U
H2 = r2.H

@test_approx_eq_eps U2'*U2 eye(m) 1e-7

# Test Hermitian positive semifefinite matrix H

@test issym(H2)

for i in eigvals(H2)
    @test i >= 0.
end

@test_approx_eq B U2*H2

println("QDWH method passed test...")


text:
## Exercise 6-1
## Write a compare function takes two values, x and y, and returns 1 if x > y, 0 if x == y, and -1 if x < y.
println("Ans: ")

function compare(x, y)
    if x > y
        return 1
    elseif x < y
        return -1
    else
        return 0
    end
end

println(compare(4, 6))
println(compare(90, 56))

println("End.")


text:
# Rounding
function Base.round(x::Decimal, dpts::Int=0; normal::Bool=false)
    shift = BigInt(dpts) + x.q
    if shift > BigInt(0) || shift < x.q
        (normal) ? x : norm(x, rounded=true)
    else
    c = Base.round(x.c / BigInt(10)^(-shift))
    d = Decimal(x.s, BigInt(c), x.q - shift)
    (normal) ? d : norm(d, rounded=true)
    end
end


text:
using Plots

# Model
include_model("planar_push_block_v3")

# Horizon
T = 26

# Time step
h = 0.1

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[model.idx_u[1:2]] .= 5.0
_uu[model.idx_λ[1:4]] .= μ_surface * gravity * h * 0.25
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[model.idx_u[1:2]] .= -5.0
_ul[model.idx_λ[1:4]] .= μ_surface * gravity * h * 0.25

ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.0, 0.0, -r - 1.0e-8, 0.0] # straight
x_goal = 1.0
y_goal = 0.0
θ_goal = 0.0 * π
qT = [x_goal, y_goal, θ_goal, x_goal-r, y_goal-r]

q1 = [0.0, 0.0, 0.0, -r - 1.0e-8, 0.0] # straight
x_goal = 1.0
y_goal = 0.0
θ_goal = 0.0 * π
qT = [x_goal, y_goal, θ_goal, x_goal-r, y_goal-r]

x1 = [q1; q1]
xT = [qT; qT]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)
xl[T][4:5] .= -Inf
xl[T][9:10] .= -Inf
xu[T][4:5] .= Inf
xu[T][9:10] .= Inf

# Objective
include_objective("velocity")
obj_velocity = velocity_objective(
    [t > T / 2 ? Diagonal(1.0 * ones(model.nq)) : Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3]))

_Q_track = [1.0, 1.0, 1.0, 0.1, 0.1]
Q_track = 1.0 * Diagonal([_Q_track; _Q_track])
obj_tracking = quadratic_tracking_objective(
    [Q_track for t = 1:T],
    # [Diagonal(0.1 * ones(model.m)) for t = 1:T-1],
	[Diagonal([0.1 * ones(model.nu);
		zeros(model.nc);
		ones(model.nb);
		zeros(model.m - model.nu - model.nc - model.nb)]) for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj_penalty = PenaltyObjective(1.0e4, model.m)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity])

# Constraints
include_constraints(["contact", "stage"])
t_idx = vcat([t for t = 1:T-1])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])#, con_ctrl_comp, con_ctrl_lim])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.1 * randn(model.m) for t = 1:T-1] # random controls
for t = 1:T-1
    if t <= 5
        u0[t][1:2] = [1.0; 0.0]
    else
        u0[t][1:2] .= 0.0
    end
end

# Pack trajectoriesff into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄, info = solve(prob, copy(z0),
	tol = 1.0e-3, c_tol = 1.0e-3, max_iter = 2500)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
h̄ = h

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
#open(vis)
visualize!(vis, model,
    q, u,
    Δt = h,
	r = model.block_dim[1] + model.block_rnd)

plot(hcat(q̄...)', color = :red, width = 1.0, labels = "")
plot(hcat([ū..., ū[end]]...)[model.idx_u[1:2], :]', linetype = :steppost, color = :black, width = 1.0, labels = "")

# q_array = hcat(q...)
# u_array = hcat(u...)

# traj = Dict("q" => q_array, "u" => u_array)

# using NPZ
# i = 2
# file_path = joinpath(pwd(), "examples/contact_implicit/manipulation/", "traj$i.npz")
# npzwrite(file_path, traj)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
@add_benchmark function sdp_nuclear_norm_atom(handle_problem)
    y = Semidefinite(3)
    p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2)
    handle_problem(p)
end

@add_benchmark function sdp_operator_norm_atom(handle_problem)
    y = Variable((3,3))
    p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
    handle_problem(p)
end

@add_benchmark function sdp_sigma_max_atom(handle_problem)
    y = Variable((3,3))
    p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12)
    handle_problem(p)
end

@add_benchmark function sdp_lambda_max_atom(handle_problem)
    y = Semidefinite(3)
    p = minimize(lambdamax(y), y[1,1]>=4)
    handle_problem(p)
end

@add_benchmark function sdp_lambda_min_atom(handle_problem)
    y = Semidefinite(3)
    p = maximize(lambdamin(y), tr(y)<=6)
    handle_problem(p)
end

@add_benchmark function sdp_matrix_frac_atom(handle_problem)
    probs = [
        let
            x = [1, 2, 3]
            P = Variable(3, 3)
            p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3))
        end,
        let
            x = Variable(3)
            P = Variable(3, 3)
            p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1)
        end
    ]

    handle_problem.(probs)
end

@add_benchmark function sdp_sum_squares_atom(handle_problem)
    n = 10
    A = rand(n,n) + im*rand(n,n)
    A = A + A' # now A is hermitian
    x = ComplexVariable(n,n)
    objective = sumsquares(A - x)
    c1 = x in :SDP
    p = minimize(objective, c1)
    handle_problem(p)
end


text:
#PS: feels like an  NLP course # be more #creative [hint: meta-Idea (explaining the idea ) ]

abstract type abstractIdea end
struct idea <: abstractIdea
    title::string

    content::Any
    Context::Any # ;)
end # what makes an idea=  {Title, content} #Context <--- how about it , too? ;)

struct metaIdea <: idea #TODO:Elaborate ...
"""
an abstract construct(ion)
Rose from an idea
Could be useful (later...) #Q.what makes a concept, a concept?
"""
struct concept <: idea

end

"""

this is a seed, help us by growing it with your support- thank you!

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using DifferentialEquations
using LinearAlgebra
using Random
using Test
include("../src/twoscales_pure_ab.jl")
include("../src/henon_heiles.jl")

function testtwoscales_pure_ab()

    u0=[big"0.12345678", big"0.13787878", big"0.120099999001", big"0.12715124"]

    epsilon = big"0.01"
    t_max = big"1.0"
    println("launching of julia solver")
    prob = ODEProblem(henon_heiles_julia, u0, (big"0.0", t_max), epsilon)
    @time sol = solve(prob , abstol=1e-25, reltol=1e-25)
    sol_ref = sol.u[end]
    nb=10
    parphi = PreparePhi(
        epsilon,
        128,
        [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0],
        henon_heiles
    )
    for order=2:10
        nb = 100
        par_u0 = PrepareU0(parphi, order+2, u0)
        res_err = zeros(BigFloat,5)
        ind = 1
        while nb <= 1000
            pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
            result = twoscales_pure_ab(pargen)
            res_err[ind] = norm(result[:,end]-sol_ref,Inf)
            println("\nnb=$nb order=$order err=$(res_err[ind])")
            nb *= 10
            ind += 1
        end
        coef = Float64(log(res_err[1]) - log(res_err[2]))/log(10)
        println("order=$order coef log = $coef")
    end
end

# function testtwoscales_pure_ab2()

#     u0=[big"0.12345678", big"0.13787878", big"0.120099999001", big"0.12715124"]

#     B = [ big"0.12984599677" big"0.9277" big"0.32984110099677" big"0.142984599677"
#     big"0.4294599677" big"0.127337" big"0.4298411009977" big"0.99484599677"
#     big"0.2298499677" big"0.327667" big"0.1298410099677" big"0.342984599677"
#     big"0.7298459677" big"0.027887" big"0.7294110099677" big"0.66294599677"
#     ]
#     A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]

#     epsilon = big"0.1"

#     C = 1/epsilon * A + B

#     t_max = big"1.0"
#     fct_jul = (u,p,t) -> C*u
#     println("launching of julia solver")
#     prob = ODEProblem(fct_jul, u0, (big"0.0", t_max), epsilon)
#     @time sol = solve(prob , abstol=1e-20, reltol=1e-20)
#     sol_ref = sol.u[end]

#     sol_cal = exp(t_max*C)*u0
#     println("norm = $(norm(sol_ref-sol_cal))")
# end
function testtwoscales_pure_ab3()

    u0=[big"0.12345678", big"0.13787878", big"0.120099999001", big"0.12715124"]

    Random.seed!(9988)

    # B = [ big"-0.12984599677" big"-0.9277" big"0.32984110099677" big"0.142984599677"
    # big"-0.4294599677" big"0.127337" big"0.4298411009977" big"0.99484599677"
    # big"0.2298499677" big"0.327667" big"0.1298410099677" big"-0.342984599677"
    # big"0.7298459677" big"-0.027887" big"0.7294110099677" big"-0.66294599677"
    # ]

    B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)

    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]

    epsilon = big"0.000001"

    C = 1/epsilon * A + B

    t_max = big"1.0"

    sol_ref = exp(t_max*C)*u0

    fct = (u,p,t) -> B*u

    parphi = PreparePhi(
        epsilon,
        32,
        [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0],
        fct,
        B
    )
    for order=2:10
        nb = 1000
        par_u0 = PrepareU0(parphi, order+2, u0)
        res_err = zeros(BigFloat,5)
        ind = 1
        while nb <= 10000
            pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
            result = twoscales_pure_ab(pargen)
            res_err[ind] = norm(result[:,end]-sol_ref,Inf)
            println("\nnb=$nb order=$order err=$(res_err[ind])")
            nb *= 10
            ind += 1
        end
        coef = Float64(log(res_err[1]) - log(res_err[2]))/log(10)
        println("order=$order coef log = $coef")
    end
end
function testtwoscales_pure_ab_epsilon()
    Random.seed!(19988)
    u0=rand(BigFloat,4)
    println("u0=$u0")
    B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    fct = (u,p,t) -> B*u
    order = 5
    for i=1:5
        epsilon = big"1.0"/big"10"^(i*3)
        eps_v = convert(Float32, epsilon)
        t_max = big"1.0"
        parphi = PreparePhi(
            epsilon,
            32,
            A,
            fct,
            B
        )
        sol_ref = getexactsol(parphi, u0, t_max)
        println("epsilon=$eps_v sol_ref=$sol_ref")
        nb = 100
        par_u0 = PrepareU0(parphi, order+2, u0)
        res_err = zeros(BigFloat,5)
        ind = 1
        while nb <= 1000
            pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
            result = twoscales_pure_ab(pargen)
            res_err[ind] = norm(result[:,end]-sol_ref,Inf)
            println("\nnb=$nb epsilon=$eps_v err=$(res_err[ind])")
            nb *= 10
            ind += 1
        end
        coef = Float64(log(res_err[1]) - log(res_err[2]))/log(10)
        eps_v
        println("epsilon=$eps_v coef log = $coef")
    end
end
function testtwoscales_interpolate()
    seed=9818871
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 6
    nb = 100
    t_max = big"1.0"
    for i=1:1
        seed += 10
        Random.seed!(seed)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        fct = (u,p,t) -> B*u
        epsilon = big"0.4"/2.0^i
        eps_v = convert(Float32, epsilon)
        t_max = big"1.0"
        parphi = PreparePhi(epsilon, 32, A, fct, B)
        par_u0 = PrepareU0(parphi, order+2, u0)
        pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
        @time result, tfft, tabu = twoscales_pure_ab(pargen,
    only_end=false, res_fft=true)
        reftol=norm(getexactsol(parphi, u0, t_max)-result[:,end], Inf)*10
        for j=1:10
            t=rand(BigFloat)
            res_ex=getexactsol(parphi, u0, t)
            res_ap=_getresult(tabu, t, parphi, 0, t_max, order)
            println("t=$t norm=$(norm(res_ex-res_ap, Inf))")
            @test isapprox(res_ex, res_ap, atol=reftol, rtol=reftol*10)
        end
    end
end
function testtwoscales_short()
    seed=981885
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 7
    nb = 5
     for i=1:1
        seed += 10
        Random.seed!(seed)
        u0=rand(BigFloat,4)
        println("u0=$u0")
        B = 2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        fct = (u,p,t) -> B*u
        epsilon = big"0.00004"/2.0^i
        eps_v = convert(Float32, epsilon)
        t_max = big"0.01"
        parphi = PreparePhi(epsilon, 32, A, fct, B)
        par_u0 = PrepareU0(parphi, order+2, u0)
        pargen = PrepareTwoScalesPureAB(nb, t_max, order, par_u0)
        @time sol = twoscales_pure_ab(pargen, only_end=true)
        resnorm=norm(getexactsol(parphi, u0, t_max)-sol, Inf)
        println("resnorm=$resnorm")
        @test resnorm < 1e-10
    end
end
function tts_time(t_begin, t_end)
    @time @testset "test twoscales from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = rand(BigFloat, 4)
        fct = (u,p,t) -> B*u
        epsilon=big"0.00345"
        nb = 100
        order = 4
        ordprep= order+2
        parphi = PreparePhi(epsilon, 32, A, fct, B, t_0=t_begin)
        par_u0 = PrepareU0(parphi, order+2, u0)
        pargen = PrepareTwoScalesPureAB(nb, t_end, order, par_u0)
        sol = twoscales_pure_ab(pargen, only_end=true)
        solref = getexactsol(parphi, u0, t_end)
        println("sol=$sol solref=$solref norm=$(norm(sol-solref,Inf))")
        @test isapprox(sol, solref,atol=1e-7, rtol=1e-6)
        result, tfft, tabu = twoscales_pure_ab(pargen, only_end=false, res_fft=true)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=getexactsol(parphi, u0, t)
            res_ap=_getresult(tabu, t, parphi, t_begin, t_end, order)
            println("i=$i")
            println("t=$t")
            println("res_ex=$res_ex")
            println("res_ap=$res_ap")
            @test isapprox(res_ex, res_ap, atol=1e-6, rtol=1e-5)
        end
    end
end
function tts_time_time(t_begin, t_end)
    @time @testset "test twoscales time time from $t_begin to $t_end" begin
        A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
        B =  2rand(BigFloat,4,4)-ones(BigFloat,4,4)
        u0 = 2rand(BigFloat, 4)-ones(BigFloat,4)
        fct = (u,p,t) -> B*u +t*p[1] +p[2]
        tuple_p = (2rand(BigFloat,4)-ones(BigFloat,4),
    2rand(BigFloat,4)-ones(BigFloat,4))
        epsilon=big"0.0003467"
        nb = 100
        order = 4
        ordprep= order+2
        parphi = PreparePhi(epsilon, 32, A, fct, B, paramfct=tuple_p, t_0=t_begin)
        par_u0 = PrepareU0(parphi, order+2, u0)
        pargen = PrepareTwoScalesPureAB(nb, t_end, order, par_u0)
        sol = twoscales_pure_ab(pargen, only_end=true)
        solref = getexactsol(parphi, u0, t_end)
        println("sol=$sol solref=$solref norm=$(norm(sol-solref,Inf))")
        @test isapprox(sol, solref,atol=1e-8, rtol=1e-7)
        result, tfft, tabu = twoscales_pure_ab(pargen, only_end=false, res_fft=true)
        for i=1:10
            t = rand(BigFloat)*(t_end-t_begin) + t_begin
            res_ex=getexactsol(parphi, u0, t)
            res_ap=_getresult(tabu, t, parphi, t_begin, t_end, order)
            @test isapprox(res_ex, res_ap, atol=1e-7, rtol=1e-6)
        end
    end
end
function testtwoscales_time()
    seed=7887
    Random.seed!(seed)
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 5
    nb = 10
    t_max = big"1.0"
    tts_time(big"0.0",big"1.0")
    tts_time(big"0.0",big"1.4565656565")
    tts_time(big"2.234566",big"3.45766790123")
    tts_time(-big"0.8457676",big"0.56716")
    tts_time(-big"1.8111457676",-big"0.345456716")
    tts_time(big"1.8457676",big"0.56716")
    tts_time(-big"0.8457676",-big"1.56716")
    tts_time(-big"0.845722676",-big"0.56716")
end
function testtwoscales_time_time()
    seed=788227
    Random.seed!(seed)
    A =  [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    order = 5
    nb = 10
    t_max = big"1.0"
    tts_time_time(big"0.0",big"1.0")
    tts_time_time(big"0.0",big"1.4565656565")
    tts_time_time(big"11.234566",big"12.45766790123")
    tts_time_time(-big"0.8457676",big"0.56716")
    tts_time_time(-big"1.8111457676",-big"0.345456716")
    tts_time_time(big"1.8457676",big"0.56716")
    tts_time_time(-big"0.8457676",-big"1.56716")
    tts_time_time(-big"0.845722676",-big"0.56716")
end
function test_exact()
    A = [0 0 1 0; 0 0 0 0; -1 0 0 0; 0 0 0 0]

    epsilon = big"0.1"
    n_tau = 32
    parphi = PreparePhi(
        epsilon,
        n_tau,
        A,
        henon_heiles
    )
    @test !isexactsol(parphi)

    u0=[big"0.12345678", big"0.13787878", big"0.120099999001", big"0.12715124"]

    B = [ big"-0.12984599677" big"-0.9277" big"0.32984110099677" big"0.142984599677"
    big"-0.4294599677" big"0.127337" big"0.4298411009977" big"0.99484599677"
    big"0.2298499677" big"0.327667" big"0.1298410099677" big"-0.342984599677"
    big"0.7298459677" big"-0.027887" big"0.7294110099677" big"-0.66294599677"
    ]

    fct = u -> B*u

    parphi = PreparePhi(
        epsilon,
        n_tau,
        A,
        fct,
        B
    )
    @test isexactsol(parphi)

    C = 1/epsilon * A + B

    t = big"0.12347171717"

    @test  isapprox(exp(t*C)*u0, getexactsol(parphi, u0, t), atol=1e-77,rtol=1e-77)

end
test_exact()
testtwoscales_time()

testtwoscales_time_time()
testtwoscales_short()
testtwoscales_interpolate()
# testtwoscales_pure_ab()
# testtwoscales_pure_ab2()
#testtwoscales_pure_ab3()
testtwoscales_pure_ab_epsilon()


text:
###########################################################
# This is my module that contains some wrapper functions.
# Author: Tokuma Suzuki
###########################################################

using LinearAlgebra, SparseArrays

function speye(size, dtype=Float64)
    """
    Return sparse identity matrix.

    This function needs LinearAlgebra and SparseArrays.
    Arguments:
    size: Int64
    dtype: keywords, Float64 or Int64
    """

    return sparse(Matrix{dtype}(I, size, size))
end

function linspace(start,stop,size=50)
    """

    Return linealy spaced 1d array

    Arguments:
    start: Real
    end: Real
    size: Int
    """

    return collect(LinRange(start,stop,size))
end

function interp1(xpt, ypt, x; method="linear", extrapvalue=nothing)

    """
    This function needs Interpolations and works like interp1 in MATLAB.
    If you want to input array{Float64,1} as xpt, choose "linear"

    Input
        xpt: grid
        ypt: vector that we want to interpolate
        x: evaluation point
        method: "linear" or "cubic"
        extravalue: default is nothing

    Return
        y: interpolated function
    """

    if extrapvalue == nothing
        y = zeros(x)
        idx = trues(x)
    else
        y = extrapvalue*ones(x)
        idx = (x .>= xpt[1]) .& (x .<= xpt[end])
    end

    if method == "linear"
        intf = interpolate((xpt,), ypt, Gridded(Linear()))
        y[idx] = intf[x[idx]]

    elseif method == "cubic"
        itp = interpolate(ypt, BSpline(Cubic(Natural())), OnGrid())
        intf = Interpolations.scale(itp, xpt)
        y[idx] = [intf[xi] for xi in x[idx]]
    end

    return y
end


text:
## Simple hello world showing
## * window
## * layout
## * button
## * callback
## * dialog

## This is done twice: first with PyCall syntax, then with qtextras

using PySide

w = Qt[:QWidget]()
w[:setWindowTitle]("Hello world example")
lyt = Qt[:QVBoxLayout](w)
w[:setLayout](lyt)

btn = Qt[:QPushButton]("Click me", w)
lyt[:addWidget](btn)

qconnect(btn, :clicked) do
  msg = Qt[:QMessageBox](btn)
  msg[:setWindowTitle]("A message for you...")
  msg[:setText]("Hello world!")
  msg[:setInformativeText]("Thanks for clicking.")
  msg[:setIcon](Qt[:QMessageBox]()[:Information]) # how to pick out Qt::QMessageBox::Informatino enumeration
  convert(Function, msg[:exec])()       # Must convert to a function
end

raise(w)

## Or in the qtextras style (not so different here)
w = Widget()
setWindowTitle(w, "Hello world example (redux)")
lyt = VBoxLayout(w)
setLayout(w, lyt)

btn = Button(w)
setText(btn, "Click me")
push!(lyt, btn)

qconnect(btn, :clicked) do
    println("Hi htere")
  MessageBox(btn, "Hi there", icon=:Information)
end

raise(w)


text:
### Abstract variable NStates

abstract type VariableNState{F<:VariateForm} end

add_dimension(n::Number) = eltype(n)[n]
add_dimension(a::Array, sa::Tuple=size(a)) = reshape(a, sa..., 1)

### Basic variable NState subtypes

## BasicUnvVariableNState

mutable struct BasicUnvVariableNState{N<:Number} <: VariableNState{Univariate}
  value::Vector{N}
  n::Integer
end

BasicUnvVariableNState(value::Vector{N}) where {N<:Number} = BasicUnvVariableNState{N}(value, length(value))

BasicUnvVariableNState(n::Integer, ::Type{N}=Float64) where {N<:Number} = BasicUnvVariableNState{N}(Array{N}(n), n)

eltype{N<:Number}(::Type{BasicUnvVariableNState{N}}) = N
eltype{N<:Number}(::BasicUnvVariableNState{N}) = N

copy!(nstate::BasicUnvVariableNState, state::BasicUnvVariableState, i::Integer) = (nstate.value[i] = state.value)

function show(io::IO, nstate::BasicUnvVariableNState{N}) where N<:Number
  indentation = "  "

  println(io, "BasicUnvVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  print(io, indentation*"number of states = $(nstate.n)")
end

## BasicMuvVariableNState

mutable struct BasicMuvVariableNState{N<:Number} <: VariableNState{Multivariate}
  value::Matrix{N}
  size::Integer
  n::Integer
end

BasicMuvVariableNState(value::Matrix{N}) where {N<:Number} = BasicMuvVariableNState{N}(value, size(value)...)

BasicMuvVariableNState(size::Integer, n::Integer, ::Type{N}=Float64) where {N<:Number} =
  BasicMuvVariableNState{N}(Array{N}(size, n), size, n)

eltype{N<:Number}(::Type{BasicMuvVariableNState{N}}) = N
eltype{N<:Number}(::BasicMuvVariableNState{N}) = N

copy!(nstate::BasicMuvVariableNState, state::BasicMuvVariableState, i::Integer) =
  (nstate.value[1+(i-1)*state.size:i*state.size] = state.value)

function show(io::IO, nstate::BasicMuvVariableNState{N}) where N<:Number
  indentation = "  "

  println(io, "BasicMuvVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"state size = $(nstate.size)")
  print(io, indentation*"number of states = $(nstate.n)")
end

## BasicMavVariableNState

mutable struct BasicMavVariableNState{N<:Number} <: VariableNState{Matrixvariate}
  value::Array{N, 3}
  size::Tuple{Integer, Integer}
  n::Integer
end

BasicMavVariableNState(value::Array{N, 3}) where {N<:Number} =
  BasicMavVariableNState{N}(value, (size(value, 1), size(value, 2)), size(value, 3))

BasicMavVariableNState(size::Tuple, n::Integer, ::Type{N}=Float64) where {N<:Number} =
  BasicMavVariableNState{N}(Array{N}(size..., n), size, n)

eltype{N<:Number}(::Type{BasicMavVariableNState{N}}) = N
eltype{N<:Number}(::BasicMavVariableNState{N}) = N

copy!(nstate::BasicMavVariableNState, state::BasicMavVariableState, i::Integer, statelen::Integer=prod(state.size)) =
  (nstate.value[1+(i-1)*statelen:i*statelen] = state.value)

function show(io::IO, nstate::BasicMavVariableNState{N}) where N<:Number
  indentation = "  "

  println(io, "BasicMavVariableNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"state size = $(nstate.size)")
  print(io, indentation*"number of states = $(nstate.n)")
end


text:
import ImpulseResponse

using HDF5
using Test
@everywhere using DeviceUnderTest
@everywhere using Soundcard

let (f,h,d,t) = ImpulseResponse.expsinesweep_simulate()
end

function asio_test()
    ms = zeros(1,8)
    ms[1,1] = 1
    mm = zeros(8,1)
    mm[1,1] = 1
    ImpulseResponse.expsinesweep_asio(ms, mm, 48000, 22, 22000, 10, 3, [([1.0],[1.0])], -6)
end
let (f,h,d,t) = asio_test()
end

function fileio_test()
    ms = zeros(1,2)
    ms[1,1] = 1
    mm = zeros(8,8)
    for i = 1:8
        mm[i,i] = 1
    end
    r = DeviceUnderTest.register()
    expsinesweep_fileio(r[0x9fefe994b7e95bf1], ms, mm, 47999.6, 22, 22000, 10, 3, [([1.0],[1.0])], -6, -10, 800, 2000, 0.5, 2, 3)
end

function fileio_asio_test()
    ms = zeros(1,2)
    ms[1,2] = 1
    mm = zeros(8,1)
    mm[1,1] = 1
    r = DeviceUnderTest.register()
    ImpulseResponse.expsinesweep_fileio_asio(r[0x9fefe994b7e95bf1], ms, mm, 48000, 47999.6, 22, 22000, 10, 3, [([1.0],[1.0])], -6, -10, 800, 2000, 0.5, 2, 3)
end

function asio_fileio_test()
    ms = zeros(1,8)
    ms[1,1] = 1
    mm = zeros(8,8)
    for i = 1:8
        mm[i,i] = 1
    end
    r = DeviceUnderTest.register()
    ImpulseResponse.expsinesweep_asio_fileio(r[0x9fefe994b7e95bf1], ms, mm, 48000, 47999.6, 22, 22000, 10, 3, [([1.0],[1.0])], -6, -10, 800, 2000, 0.5, 2, 3)
end

function clockdrift_test()
    ms = zeros(1,2)
    ms[1,1] = 1
    mm = zeros(8,1)
    mm[1,1] = 1
    r = DeviceUnderTest.register()
    rept = 1
    ImpulseResponse.measureclockdrift(r[0x9fefe994b7e95bf1], ms, mm, rept, 48000, 800, 2000, 0.5, -6)
end

function minimalphase_eq_test()

    ms = [0 0 0 0 1 0 0 0]
    mm = [1,0]
    mm = mm[:,:]
    fs = 48000
    atten = -20
    n = 1024

    # transfer function without any equalization
    f0 = 22
    f1 = 22000
    f, m, ϕ, xds = ImpulseResponse.minimalphase_fir_original(ms, mm, fs, f0, f1, atten, n);

    # tuning process
    f0 = 142
    f1 = 16000
    fx = 214
    h = ImpulseResponse.minimalphase_fir_design(f, m, ϕ, xds, fs, f0, f1, fx, n);
    # ...
    f1 = 13564
    h = ImpulseResponse.minimalphase_fir_design(f, m, ϕ, xds, fs, f0, f1, fx, n);

    # plot the verification
    f0 = 22
    f1 = 22000
    ImpulseResponse.minimalphase_fir_verification(h, ms, mm, fs, f0, f1, atten, n);

    h5write(joinpath("equalizationfilters-20181123.h5"), "mth1", convert(Matrix{Float64},h[:,:]))

end


text:
# # [Converting Models to Computation Graphs](@id bilayernetwork_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/covid/bilayerconversion.ipynb)

# using Pkg
# Pkg.status()
# Pkg.activate(".")

# Pkg.instantiate()

# Pkg.develop(path="../../")
# Pkg.instantiate()

using AlgebraicPetri
using AlgebraicPetri.Epidemiology
using AlgebraicPetri.BilayerNetworks

using Catlab
using Catlab.CategoricalAlgebra
using Catlab.Graphics
import Catlab.CategoricalAlgebra: migrate!
using Catlab.WiringDiagrams
using Catlab.Programs.RelationalPrograms
display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

using PrettyTables
function printsoln(bn::AbstractLabelledBilayerNetwork, soln::Vector)
    pretty_table(soln)
end

# #### SIR Model:

# define model
sir = @relation (s,i,r) begin
    infection(s,i)
    recovery(i,r)
end
display_uwd(sir)
#-
#

# Extract the Petri network form of the model.
psir = apex(oapply_epi(sir))
psir
Graph(psir)
#-

# Convert the Petri network into a bilayer network and draw it.
# This model uses a computation graph to express the computation of the vector field of the Petri net.
bnsir = BilayerNetwork()
migrate!(bnsir, psir)
bnsir
to_graphviz(bnsir)
#-

# We can hand code a Bilayer netowork using the `@acset` macro provided by Catlab. As you can see from the code,
# there is a lot of typing to specify the incidence of all these wires. The Petri Net is more compact.
# This notion of Bilayer network comes from the definition of mass action kinetics for reaction networks.

# hand coded Bilayer network
bnsir_test = @acset BilayerNetwork begin
    Qin = 3
    Qout = 3
    Win = 3
    Wa = 3
    Wn = 3
    Box = 2
    arg = [1,2,2]
    call = [1,1,2]
    efflux = [1,1,2]
    effusion = [1,2,2]
    influx = [1,1,2]
    infusion = [2,2,3]
end

@assert bnsir == bnsir_test

function roundtrip(pn::AbstractPetriNet, bn::AbstractBilayerNetwork)
    roundtrippetri = PetriNet()
    migrate!(roundtrippetri, bn)
    pn_structure = PetriNet()
    copy_parts!(pn_structure, pn)
    return roundtrippetri, pn_structure
end
#-

# We can specify the SEIR model with the relation macro and functorial semantics as usual.

seir = @relation (s,e,i,r) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
end

pseir = apex(oapply_epi(seir))

# By converting this to a Bilayer network, we are able to visualize differences in the computational pattern
# of data flow between different reaction network models.

#
bnseir = BilayerNetwork()
migrate!(bnseir, pseir)

bnrt,pnstr = roundtrip(pseir, bnseir)

display_uwd(sir)

Graph(psir)
#-

to_graphviz(bnsir)
#-

to_graphviz(bnseir)
#-

Graph(bnrt)
#-

qm = @relation (s,q) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
    quarexp(e,q)
    quarinf(i,q)
end
display_uwd(qm)
#-

import AlgebraicPetri.Epidemiology: exposure_petri, spontaneous_petri
semantics = Dict(
    :infection => exposure_petri(:S, :I, :I, :inf),
    :exposure  => exposure_petri(:S, :I, :E, :exp),
    :illness   => spontaneous_petri(:E,:I,:ill),
    :recovery  => spontaneous_petri(:I,:R,:rec),
    :death     => spontaneous_petri(:I,:D,:death),
    :quarexp   => spontaneous_petri(:E, :Q, :qe),
    :quarinf   => spontaneous_petri(:I, :Q, :qi),
    :quarrec   => spontaneous_petri(:Q, :R, :qr)
)
pn_quar = oapply(qm, semantics)  |> apex
Graph(pn_quar)
#-

bnquar = LabelledBilayerNetwork()
migrate!(bnquar, pn_quar)
to_graphviz(bnquar)

qm = @relation (s,i,q) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
    quarexp(e,q)
    quarinf(i,q)
    quarrec(q,r)
end
display_uwd(qm)
#-

pn_quar = oapply(qm, semantics)  |> apex
Graph(pn_quar)
#-

bnquar = LabelledBilayerNetwork()
migrate!(bnquar, pn_quar)
to_graphviz(bnquar)
#-

quarrt = LabelledPetriNet()
migrate!(quarrt, bnquar) |> Graph

bnquar
#-

import AlgebraicPetri.BilayerNetworks: evaluate!, evaluate

bnsir = LabelledBilayerNetwork()
migrate!(bnsir, psir)

evaluate(bnsir, [10.0, 1, 0.0], inf=0.1, rec=0.3)

evaluate(bnquar, [10.0, 1, 0, 0,0,0], exp=0.1, rec=0.3, qi=0.2, ill=0.7, qe=0.23, qr=0.3)

function euler(bn::AbstractLabelledBilayerNetwork, state, nsteps::Integer, stepsize::Real; params...)
    #preallocate storage so that each step is nonallocating
    #create a storage space for steps of euler's method store intermediate
    #states as named tuples so that you can integrate with julia Tables.jl data  analysis ecosystem
    du = zeros(length(state))
    ϕ = ones(nparts(bn, :Box))
    u = tuple(state...)
    results = Vector{NamedTuple{tuple(bn[:,:variable]...)}}()
    for i in 1:nsteps
        u = u .+ stepsize.*evaluate!(du, ϕ, bn, u; params...)
        push!(results, NamedTuple{tuple(bn[:,:variable]...)}(u))
    end
    return results
end

soln = euler(bnquar, (S=10.0, I=1, E=0, R=0, Q=0), 30, 0.15, exp=0.1, rec=0.03, qi=0.37, ill=0.7, qe=0.23, qr=0.03)
printsoln(bnquar, soln)
#-

import AlgebraicPetri.BilayerNetworks: compile
compile(bnquar, :du, :ϕ, :u, :p)
#-
#
compile(bnquar, :du, :ϕ, :u, exp=0.1, rec=0.03, qi=0.37, ill=0.7, qe=0.23, qr=0.03)
#-

function eulers(bn::AbstractLabelledBilayerNetwork, funcname::Symbol; params...)
    f = compile(bn, :du, :ϕ, :u; params...)
    varnames = tuple(bn[:,:variable]...)
    nϕ = nparts(bn, :Box)
    quote
    function $funcname(state, nsteps::Integer, stepsize::Real)
        $f
        #preallocate storage so that each step is nonallocating
        du = zeros(length(state))
        ϕ = ones($nϕ)
        u = tuple(state...)
        #create a storage space for steps of euler's method
        #store intermediate states as named tuples so that you can integrate
        #with julia Tables.jl data analysis ecosystem
        results = Vector{NamedTuple{$varnames}}()
        for i in 1:nsteps
            Δ = f!(du, ϕ, u, 0)
            u = u .+ stepsize.*Δ
            push!(results, NamedTuple{$varnames}(u))
        end
        return results
    end
    end
end

eulers(bnsir, :eulsir, inf=0.3, rec=0.2)
#-

eulseirqexp = eulers(bnquar, :eulseirq, exp=0.1, rec=0.03, qi=0.37, ill=0.7, qe=0.23, qr=0.03)
#-

eulsirexp = eulers(bnsir, :eulsir, inf=0.3, rec=0.2)
#-

eval(eulsirexp)
soln_codegen = eulsir((S=10.0, I=1, R=0), 30, 0.15)
pretty_table(soln_codegen)
#-

eval(eulseirqexp)
soln_codegen = eulseirq((S=10.0, I=1, E=0, R=0, Q=0), 30, 0.15)
pretty_table(soln_codegen)
#-

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -65.76627584345755 <= q <= 59.26265199263469)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.4079491356153899 + -0.18134844314194787*1/(1 + exp(-(0.2817697329190594 + 0.989386077401103*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.6068435500303653*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.8530798012021474*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.8376755183348754*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.6594644406874322*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.2976334456396539*1/(1 + exp(-(-0.4984715860832294 + -0.6460783499544172*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.0952691245496573*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.5434919689873383*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.2400778955734424*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.904210133814042*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.11940205422195982*1/(1 + exp(-(-0.9402956128299329 + 0.19522910845028774*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.5236489579886725*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.09802794343243137*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.5359271533353085*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.5362141752531211*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.6350558501595733*1/(1 + exp(-(-0.842982128122256 + 0.6833459994096582*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.15560340828874075*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.5368387122163134*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.6086259049043772*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.5403730914454505*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.4683577877900529*1/(1 + exp(-(0.6067063734003204 + 0.8661978259504042*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.06088102390312189*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.9836506901320718*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.7952476421961925*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.6703006410599177*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))))))))) + 1/(1 + exp(-(0.9379679817528435 + -0.4660824895766815*1/(1 + exp(-(0.2817697329190594 + 0.989386077401103*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.6068435500303653*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.8530798012021474*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.8376755183348754*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.6594644406874322*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.32761179552569075*1/(1 + exp(-(-0.4984715860832294 + -0.6460783499544172*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.0952691245496573*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.5434919689873383*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.2400778955734424*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.904210133814042*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.6880683972685175*1/(1 + exp(-(-0.9402956128299329 + 0.19522910845028774*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.5236489579886725*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.09802794343243137*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.5359271533353085*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.5362141752531211*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.5888723964207281*1/(1 + exp(-(-0.842982128122256 + 0.6833459994096582*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.15560340828874075*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.5368387122163134*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.6086259049043772*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.5403730914454505*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.05595383404471299*1/(1 + exp(-(0.6067063734003204 + 0.8661978259504042*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.06088102390312189*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.9836506901320718*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.7952476421961925*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.6703006410599177*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))))))))) + 1/(1 + exp(-(-0.3048668135843893 + -0.3680078347752045*1/(1 + exp(-(0.2817697329190594 + 0.989386077401103*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.6068435500303653*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.8530798012021474*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.8376755183348754*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.6594644406874322*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.09989706667818732*1/(1 + exp(-(-0.4984715860832294 + -0.6460783499544172*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.0952691245496573*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.5434919689873383*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.2400778955734424*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.904210133814042*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.6733365176188579*1/(1 + exp(-(-0.9402956128299329 + 0.19522910845028774*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.5236489579886725*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.09802794343243137*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.5359271533353085*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.5362141752531211*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.7967060490111804*1/(1 + exp(-(-0.842982128122256 + 0.6833459994096582*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.15560340828874075*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.5368387122163134*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.6086259049043772*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.5403730914454505*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.8483184805894783*1/(1 + exp(-(0.6067063734003204 + 0.8661978259504042*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.06088102390312189*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.9836506901320718*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.7952476421961925*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.6703006410599177*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))))))))) + 1/(1 + exp(-(-0.8260325513973243 + 0.9282935832138137*1/(1 + exp(-(0.2817697329190594 + 0.989386077401103*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.6068435500303653*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.8530798012021474*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.8376755183348754*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.6594644406874322*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.5968467249691596*1/(1 + exp(-(-0.4984715860832294 + -0.6460783499544172*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.0952691245496573*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.5434919689873383*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.2400778955734424*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.904210133814042*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.003449906971975558*1/(1 + exp(-(-0.9402956128299329 + 0.19522910845028774*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.5236489579886725*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.09802794343243137*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.5359271533353085*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.5362141752531211*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.5427745947324003*1/(1 + exp(-(-0.842982128122256 + 0.6833459994096582*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.15560340828874075*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.5368387122163134*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.6086259049043772*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.5403730914454505*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.4794292841893828*1/(1 + exp(-(0.6067063734003204 + 0.8661978259504042*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.06088102390312189*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.9836506901320718*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.7952476421961925*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.6703006410599177*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))))))))) + 1/(1 + exp(-(-0.3442450542855955 + 0.3679171232252214*1/(1 + exp(-(0.2817697329190594 + 0.989386077401103*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.6068435500303653*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.8530798012021474*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.8376755183348754*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.6594644406874322*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.4969759600474797*1/(1 + exp(-(-0.4984715860832294 + -0.6460783499544172*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.0952691245496573*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.5434919689873383*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.2400778955734424*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.904210133814042*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.17345395101703964*1/(1 + exp(-(-0.9402956128299329 + 0.19522910845028774*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.5236489579886725*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + -0.09802794343243137*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.5359271533353085*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + 0.5362141752531211*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + -0.24637355767613478*1/(1 + exp(-(-0.842982128122256 + 0.6833459994096582*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + 0.15560340828874075*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.5368387122163134*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + -0.6086259049043772*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.5403730914454505*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3]))))))) + 0.4183638306902213*1/(1 + exp(-(0.6067063734003204 + 0.8661978259504042*1/(1 + exp(-(-0.7838782818766057 + 0.6083304932652234*$(x[1]) + 0.4982278286707449*$(x[2]) + -0.2123905230276213*$(x[3])))) + -0.06088102390312189*1/(1 + exp(-(-0.6220948510874038 + -0.22981743777707297*$(x[1]) + 0.33188126395442463*$(x[2]) + -0.8264122617752747*$(x[3])))) + 0.9836506901320718*1/(1 + exp(-(0.5453301275165048 + -0.2798145416678528*$(x[1]) + -0.9787415760841154*$(x[2]) + 0.5329097417839521*$(x[3])))) + 0.7952476421961925*1/(1 + exp(-(-0.31934053056836653 + 0.3102265855539219*$(x[1]) + -0.9236015463563882*$(x[2]) + -0.5629164596540628*$(x[3])))) + -0.6703006410599177*1/(1 + exp(-(0.17795541107351465 + 0.6245951080916892*$(x[1]) + -0.14381738742917527*$(x[2]) + 0.44898943772247746*$(x[3])))))))))) - $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