text:
"""
ULMFiT - FINE-TUNING
"""

using Flux

cd(@__DIR__)
include("pretrain_lm.jl")    # importing LanguageModel and useful functions
include("custom_layers.jl")      # importing AWD_LSTM, VarDrop and DroppedEmbeddings
include("utils.jl")         # importing utilities

function discriminative_step!(layers, ηL::Float64, l, gradient_clip::Float64, opts::Vector)
    # Applying gradient clipping
    l = Tracker.hook(x -> grad_clipping(x, gradient_clip), l)

    # Gradient calculation
    grads = Tracker.gradient(() -> l, get_trainable_params(layers))

    # discriminative step
    ηl = ηL/(2.6^(length(layers)-1))
    for (layer, opt) in zip(layers, opts)
        opt.eta = ηl
        for ps in get_trainable_params([layer])
            Tracker.update!(opt, ps, grads[ps])
        end
        ηl *= 2.6
    end
    return
end

# Fine-Tuning Language Model
function fine_tune_lm!(lm::LanguageModel; batchsize::Integer=64, bptt::Integer=70, data_loader::Channel=imdb_fine_tune_data, gradient_clip::Float64=0.25,
        ηL::Float64=4e-3, stlr_cut_frac::Float64=0.1, stlr_ratio::Float32=32, stlr_η_max::Float64=0.01, epochs::Integer=1, checkpoint_itvl::Integer=5000)

    opts = [ADAM(0.001, (0.7, 0.99)) for i=1:4]
    cut = num_of_iters * epochs * stlr_cut_frac
    gpu!.(lm.layers)

    # Fine-Tuning loops
    for epoch=1:epochs
        gen = data_loader()
        num_of_iters = take!(gen)
        T = num_of_iters-Int(floor((num_of_iters*2)/100))
        set_trigger!.(T, lm.layers)
        for i=1:num_of_iters

            # FORWARD
            l = loss(lm, gen)

            # Slanted triangular learning rate step
            t = i + (epoch-1)*num_of_iters
            p_frac = (i < cut) ? i/cut : (1 - ((i-cut)/(cut*(1/stlr_cut_frac-1))))
            ηL = stlr_η_max*((1+p_frac*(stlr_ratio-1))/stlr_ratio)

            # Backprop with discriminative fine-tuning step
            discriminative_step!(lm.layers[[1, 3, 5, 7]], ηL, l, gradient_clip, opts)

            # ASGD Step, after Triggering
            asgd_step!.(i, lm.layers)

            # Resetting dropout masks for all the layers with DropOut or DropConnect
            reset_masks!.(lm.layers)

            println("loss: $l", " iteration completed: $i")

            # Saving checkpoints
            if i == checkpoint_itvl save_model!(lm) end
        end
        println("\nEpoch: $epoch")
    end
end


text:
n_actions(p::LaserTagPOMDP) = p.diag_actions ? 9 : 5
POMDPs.actions(p::LaserTagPOMDP) = 1:n_actions(p)
POMDPs.actionindex(p::LaserTagPOMDP, a::Int) = a

const ACTION_NAMES = SVector("north",
                             "east",
                             "south",
                             "west",
                             "tag",
                             "northeast",
                             "southeast",
                             "southwest",
                             "northwest")

const TAG_ACTION = 5

const ACTION_DIRS = SVector(Coord(0,1),
                            Coord(1,0),
                            Coord(0,-1),
                            Coord(-1,0),
                            Coord(0,0),
                            Coord(1,1),
                            Coord(1,-1),
                            Coord(-1,-1),
                            Coord(-1,1))

const CARDINALS = SVector(Coord(0,1),
                            Coord(1,0),
                            Coord(0,-1),
                            Coord(-1,0))

const DIAGONALS = SVector(Coord(1,1),
                         Coord(1,-1),
                         Coord(-1,-1),
                         Coord(-1,1))

const DIR_TO_ACTION = Dict{Coord, Int}(c=>i for (i,c) in enumerate(ACTION_DIRS))

ext: jl


text:
import ..MadNLP

function MadNLP.Optimizer(m::Model)
    _obj,_grad!,_con!,_jac!,jac_sparsity,_hess!,hess_sparsity = nlp_evaluator(m.obj,m.con)
    _jac_fill_sparsity! = @inline (I,J)->fill_sparsity!(I,J,jac_sparsity)
    _hess_fill_sparsity! = @inline (I,J)->fill_sparsity!(I,J,hess_sparsity)
    p = m.p

    prob = MadNLP.NonlinearProgram(
        m.n,m.m,length(hess_sparsity),length(jac_sparsity),
        0.,m.x,m.g,m.l,m.zl,m.zu,m.xl,m.xu,m.gl,m.gu,
        x->_obj(x,p),(g,x)->_grad!(g,x,p),(c,x)->_con!(c,x,p),(jac,x)->_jac!(jac,x,p),
        (hess,x,lag,sig)->_hess!(hess,x,p,lag,sig),
        _hess_fill_sparsity!,_jac_fill_sparsity!,
        MadNLP.INITIAL,Dict())

    prob.ext[:solver] = MadNLP.Solver(prob;m.opt...)

    return prob
end

optimize!(prob::MadNLP.NonlinearProgram) = MadNLP.optimize!(prob.ext[:solver])


text:
module BHAPtfem

# package code goes here

using Reexport, DataFrames, Statistics, Distributed
using SparseArrays, LinearAlgebra

@reexport using PtFEM

# General utilities

include("util/types.jl")

# input routines (called from BHAPtfem.jl)

include("input/creatematerialdict.jl")
include("input/createmediadict.jl")
include("input/input.jl")
#include("input/updatematerialtable.jl")
#include("input/updatemediatable.jl")

# problem input handling

include("problem/problem.jl")
include("problem/createmesh.jl")
include("problem/createsegmentdf.jl")
include("problem/createpropertydf.jl")
include("problem/createnodedf.jl")
include("problem/createelementdf.jl")
include("problem/createcasetable.jl")

# user visible (called from e.g. example projects)

include("solve/solve.jl")

# core ptfem based mp runs

include("threads/runcase.jl")

# These table are created here, a user might update and/or add to these dicts

bha = Dict{Symbol, Any}()
materials = creatematerialdict()
media = createmediadict()

export
  bha,
  materials,
  media,
  #updatedmaterialstable,
  #updatemedatable,
  input!,
  problem,
  solve

end # module


text:
function coords(idx::Int, s::Int)::Tuple
    row = ((idx - 1) ÷ s) + 1
    col = ((idx - 1) % s) + 1
    return (row, col)
end

function coord_index(row::Int, col::Int, s::Int)::Int
    s * (row - 1) + col
end

function rot(idx::Int, s::Int, angle::Int)::Int
    row, col = coords(idx, s)
    if angle == 0
        return idx
    elseif angle == 90
        rot_row = s + 1 - col
        rot_col = row
    elseif angle == 180
        rot_row = s + 1 - row
        rot_col = s + 1 - col
    elseif angle == 270
        rot_row = col
        rot_col = s + 1 - row
    else
        throw(DomainError(angle, "Angle must be 0º, 90º, 180º, or 270º"))
    end
    coord_index(rot_row, rot_col, s)
end

function rot(combo::AbstractArray, s::Int, angle::Int)::AbstractArray
    [rot(v, s, angle) for v in combo]
end

function reflect(idx::Int, s::Int, horizontal::Bool)::Int
    row, col = coords(idx, s)
    if horizontal
        ref_row = row
        ref_col = s + 1 - col
    else
        ref_row = s + 1 - row
        ref_col = col
    end
    coord_index(ref_row, ref_col, s)
end

function reflect(combo::AbstractArray, s::Int, horizontal::Bool)::AbstractArray
    [reflect(v, s, horizontal) for v in combo]
end

function transpose(idx::Int, s::Int)::Int
    row, col = coords(idx, s)
    coord_index(col, row, s)
end

function transpose(combo::AbstractArray, s::Int)::AbstractArray
    [transpose(v, s) for v in combo]
end

"""
Determines if `left` is less than `right`--that is, `left` becomes
`right` when the vote distributions are in sorted order.
"""
function voter_dist_lt(left::AbstractArray, right::AbstractArray)
    left_s = sort(left)
    right_s = sort(left)
    for (l, r) in zip(left_s, right_s)
        if l > r
            return false
        elseif l < r
            return true
        end
    end
    false  # exactly equal
end

function combo_duplicated(combo::AbstractArray, s::Int)
    voter_dist_lt(reflect(combo, s, true), s) ||
    voter_dist_lt(reflect(combo, s, false), s) ||
    voter_dist_lt(rot(combo, s, 90), s) ||
    voter_dist_lt(rot(combo, s, 180), s) ||
    voter_dist_lt(rot(combo, s, 270), s) ||
    voter_dist_lt(transpose(combo, s), s) ||
    voter_dist_lt(rot(transpose(combo, s), s, 180), s)
end

"""
Determines the number of unique combinations resulting from
transforming a combination (reflection, rotation, transpose).
"""
function combo_count(combo::AbstractArray)
    combos = Set([
        sort(combos),
        sort(reflect(combo, s, true)),
        sort(reflect(combo, s, false)),
        sort(rot(combo, s, 90)),
        sort(rot(combo, s, 180)),
        sort(rot(combo, s, 270)),
        sort(transpose(combo, s)),
        sort(rot(transpose(combo, s), s, 180))
    ])
    length(combos)
end


text:
using Plots

"""
    Strategy(decisions, scores, energies, activeMons)

Struct for information about an entire match and the strategy used throughout.
This includes the decisions made, the scores at each point along the battle,
the energies of the mons (to determine fast and charged move applications), and
the active mons used by each team.
"""
mutable struct IndividualStrategy
    decisions::Array{Tuple{UInt8,UInt8}}
    scores::Array{Float64}
    hps::Array{Tuple{UInt16, UInt16}}
end

"""
    plot_strategy(strat, static_s)

Given the Strategy and the StaticState, plot the battle using a PvPoke-like
notation.
"""
function plot_strategy(strat::IndividualStrategy,
    static_s::StaticIndividualState)
    gr()
    plt1 = plot(1:length(strat.scores), strat.scores, width = 2,
        label = "possible battle scores", ylims = [0, 1],
        ylabel = "Battle Score", xlabel = "Decisions", size = (950, 400))
    plot!(plt1, 1:length(strat.scores), map(i -> strat.hps[i][1],
        1:length(strat.scores)) ./ static_s[0x01].stats.hitpoints,
        color = :green, label = "")
    plot!(plt1, 1:length(strat.scores), map(i -> strat.hps[i][2],
        1:length(strat.scores)) ./ static_s[0x02].stats.hitpoints,
        color = :purple, label = "")
    hline!(plt1, [0.5], label = "win/loss")
    plt2 = plot(xlims = [0, length(strat.scores)], ylims = [-3, 0],
        legend = false, size = (950, 200), axis = nothing)
    for i = 1:length(strat.scores), j = 0x01:0x02
        if strat.decisions[i][j] == 0x01
            scatter!(plt2, [i], [-Int64(j)], markershape = :hexagon,
                markersize = 12, alpha = 0.5, color = shieldColor)
        elseif strat.decisions[i][j] == 0x05 || strat.decisions[i][j] == 0x06
            color = colors[strat.decisions[i][j] == 0x05 ?
                static_s[j].charged_move_1.moveType :
                static_s[j].charged_move_2.moveType]
            scatter!(plt2, [i], [-Int64(j)], markershape = :circle,
                markersize = 10, alpha = 0.5, color = color)
        end
        if i > 1
            if strat.decisions[i][j] != 0x05 &&
                strat.decisions[i][j] != 0x06 &&
                strat.hps[i][get_other_agent(j)] <
                strat.hps[i - 1][get_other_agent(j)]
                color = colors[static_s[j].fast_move.moveType]
                scatter!(plt2, [i], [-Int64(j)], markershape = :square,
                    alpha = 0.5, color = color)
            end
            if strat.hps[i][j] == 0x0000
                scatter!(plt2, [i], [-Int64(j)], markershape = :xcross,
                    alpha = 0.5, markersize = 10, color = :red)
            end
        end
    end

    l = @layout [a; b{0.2h}]
    plot(plt1, plt2, layout = l, size = (950, 600))
end


text:
# ========================Adjoint=========================================
# output windowed sources, input is shotv
function MultiStepAdjoint(oz::Int64, ox::Int64, wnz::Int64, wnx::Int64, wnt::Int64,
                          shot::ShotV, fidMtx::FidMtx)
    dt = shot.dt; nt = shot.nt;
    nz = shot.nz; nx = shot.nx;
    ext= shot.ext; iflag = shot.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       zupper = 0
    end
    Nx = nx + 2*ext
    wsc = InitWSrcCube(oz, ox, wnz, wnx, ext, iflag, 0.0, dt, wnt)
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, nt)
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, nt)
    tmp  = zeros(length(spt1.Txxx))
    tmp1 = zeros(tmp);
    AddShotV2Spt!(spt1, shot)
    for it = nt-1 : -1 : 1
        OneStepAdjoint!(spt2, spt1, fidMtx, tmp, tmp1)
        AddShotV2Spt!(spt2, shot)
        CopySnapShot!(spt1, spt2)
        if it <= wnt
           spt2Wsc!(wsc, spt1)
        end
    end
    return wsc
end

# output wave field, input shotv
function MultiStepAdjoint(path::String, shotv::ShotV, fidMtx::FidMtx)
    nz=shotv.nz; nx=shotv.nx; ext=shotv.ext; iflag=shotv.iflag
    ot = shotv.ot; dt = shotv.dt; nt = shotv.nt
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, nt  )
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, nt-1)
    AddShotV2Spt!(spt1, shotv)
    fid = WriteWfd(path, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx, tmp, tmp1)
        AddShotV2Spt!(spt2, shotv)
        CopySnapShot!(spt1, spt2)
        WriteWfd(fid, spt1)
    end
    close(fid)
    ReverseOrderWfd(path)
    return nothing
end

# output two images, input shotv
function MultiStepAdjoint(shotv::ShotV, fidMtx::FidMtx, path::String)
    (nz, nx, ext, iflag, dt, ns) = InfoPv(path);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dm = zeros(Nz*Nx); du = zeros(Nz*Nx);
    ot = shotv.ot; dt = shotv.dt; nt = shotv.nt
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, nt  )
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, nt-1)
    AddShotV2Spt!(spt1, shotv)
    if nt <= ns
       imaging!(dm, du, spt1, path)
    end
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx, tmp, tmp1)
        AddShotV2Spt!(spt2, shotv)
        CopySnapShot!(spt1, spt2)
        if it <= ns
           imaging!(dm, du, spt1, path)
        end
    end
    return dm, du
end

# function MultiStepAdjoint(shotv::ShotV, fidMtx::FidMtx, path::String)
#     (nz, nx, ext, iflag, dt, ns) = InfoPv(path);
#     dm = zeros(nz, nx); du = zeros(nz, nx);
#     ot = shotv.ot; dt = shotv.dt; nt = shotv.nt
#     spt1 = InitSnapShot(nz, nx, ext, iflag, dt, nt  )
#     spt2 = InitSnapShot(nz, nx, ext, iflag, dt, nt-1)
#     AddShotV2Spt!(spt1, shotv)
#     if nt <= ns
#        (dm, du) = imaging(dm, du, spt1, path)
#     end
#     tmp = zeros(length(spt1.Vxx))
#     tmp1= zeros(tmp)
#     for it = nt-1: -1: 1
#         OneStepAdjoint!(spt2, spt1, fidMtx, tmp, tmp1)
#         AddShotV2Spt!(spt2, shotv)
#         CopySnapShot!(spt1, spt2)
#         if it <= ns
#            (dm, du) = imaging(dm, du, spt1, path)
#         end
#     end
#     return dm, du
# end

# output SnapShots, input shotv
function MultiStepAdjoint_spt(path::String, shotv::ShotV, fidMtx::FidMtx)
    nz=shotv.nz; nx=shotv.nx; ext=shotv.ext; iflag=shotv.iflag
    ot = shotv.ot; dt = shotv.dt; nt = shotv.nt
    spt1 = InitSnapShot(nz, nx, ext, iflag, dt, nt  )
    spt2 = InitSnapShot(nz, nx, ext, iflag, dt, nt-1)
    AddShotV2Spt!(spt1, shotv)
    fid = WriteSnapShot(path, spt1)
    tmp = zeros(length(spt1.Vxx))
    tmp1= zeros(tmp)
    for it = nt-1: -1: 1
        OneStepAdjoint!(spt2, spt1, fidMtx, tmp, tmp1)
        AddShotV2Spt!(spt2, shotv)
        CopySnapShot!(spt1, spt2)
        WriteSnapShot(fid, spt1)
    end
    close(fid)
    ReverseOrderSnapShots(path)
    return nothing
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
function get_biotic_phenotype(species, epistasisMat, pleiotropyMat, expressionArray, model::ABM)
  abp = model.biotic_phenotypes[species]
  x = pleiotropyMat[abp, abp] * (epistasisMat[abp, abp] * expressionArray[abp])  # phenotypic values
  d = model.E[species]
  z = x .+ rand(d)
  return z
end

get_biotic_phenotype(ag::Ind, model::ABM) = get_biotic_phenotype(ag.species, ag.epistasisMat, ag.pleiotropyMat, ag.q, model)

function get_abiotic_phenotype(species, epistasisMat, pleiotropyMat, expressionArray, model::ABM)
  abp = model.abiotic_phenotypes[species]
  x = pleiotropyMat[abp, abp] * (epistasisMat[abp, abp] * expressionArray[abp])  # phenotypic values
  d = model.E[species]
  if d.σ == 0 && d.μ == 0
    randd = 0.0
  else
    randd = rand(d)
  end
  z = x .+ randd
  return z
end

get_abiotic_phenotype(ag::Ind, model::ABM) = get_abiotic_phenotype(ag.species, ag.epistasisMat, ag.pleiotropyMat, ag.q, model)

"""
Calculates the average distance between two sets of phenotypes.

Distance is a probability and is calculated by the number of standard deviation that phenotype 2 is different from phenotype 1.

Variance can stand for selection coefficient, with a larger variance, less sever the difference.
"""
function phenotypic_distance(phenotype1, phenotype2, interaction_coeff, variance::Real)
  distance = 0.0
  counter = 0.0
  for ph1 in phenotype1
    for ph2 in phenotype2
      d= abs(0.5 - cdf(Normal(ph1,  variance), ph2)) # 0.5 is the value when ph1 and ph2 are identical
      distance += d
      counter += 1.0
    end
  end
  distance += distance  # Double the distance so that it is in range 0 to 1.
  distance /= counter  # average distance
  if interaction_coeff < 0
    distance = 1.0 - distance
  end
  return distance
end

phenotypic_distance(ag1::Ind, ag2::Ind, model) = phenotypic_distance(ag1.biotic_phenotype, ag2.biotic_phenotype, model.interactions[ag1.species, ag2.species], variance)

"""
Abiotic distance to the optimal values
"""
function abiotic_distance(phenotype, optimal, variance)
  distance = 0.0
  counter = 0.0
  for i in eachindex(phenotype)
    d = abs(0.5 - cdf(Normal(optimal[i],  variance), phenotype[i]))
    distance += d
    counter += 1.0
  end
  distance += distance
  distance /= counter
  return distance
end

function abiotic_fitness(abiotic_phenotype, species::Int, pos, model::ABM)
  rawW = 1.0 - abiotic_distance(abiotic_phenotype, return_opt_phenotype(species, model.step[1], pos, model), variance)
  # W = adjust_abiotic_fitness(rawW, model.selectionCoeffs[species])
  return rawW
end

abiotic_fitness(ag::Ind, model::ABM) = abiotic_fitness(ag.abiotic_phenotype, ag.species, ag.pos, model)

# adjust_abiotic_fitness(rawfitness::AbstractFloat, selection_coeff::AbstractFloat) = 1.0 - ((1.0 - rawfitness) * selection_coeff)

function interaction_power(ag1, ag2, d, model)
  prob = (1.0 - d) * abs(model.interactions[ag1.species, ag2.species])
  return prob
end

"""
    interact!(ag1::Ind, ag2::Ind, model::ABM)

Changes the two interacting individuals' survivability or decide their number of children.

Each individual has a "raw" survival prob. that is determined by its abiot phenotypes.
Interactions increase or decrease that probablity.
"""
function interact!(ag1::Ind, ag2::Ind, model::ABM)
  sp1 = ag1.species
  sp2 = ag2.species

  # record interaction history
  ag1.interaction_history[sp2] = model.step[1]
  ag2.interaction_history[sp1] = model.step[1]

  if sp1 != sp2 && model.food_sources[sp1, sp2] > 0  # predation
    d = phenotypic_distance(ag1, ag2, model)
    if rand() < d
      eat!(ag1, ag2, model)
    end
  elseif sp1 != sp2 && model.food_sources[sp2, sp1] > 0 # predation
    d = phenotypic_distance(ag1, ag2, model)
    if rand() < d
      eat!(ag2, ag1, model)
    end
  else # interaction
    if sp1 == sp2 && ag1.sex != ag2.sex && model.ploidy[sp1] == 2 && in_reproduction_age(ag1, model) && in_reproduction_age(ag2, model)  # reproduce
      reproduce!(ag1::Ind, ag2::Ind, model::ABM)
    else
      ix_value1 = model.interactions[sp1, sp2]
      ix_value2 = model.interactions[sp2, sp1]
      if ix_value1 != 0 || ix_value2 != 0
        d = phenotypic_distance(ag1, ag2, model)
        inx_prob1 = interaction_power(ag1, ag2, d, model)
        inx_prob2 = interaction_power(ag2, ag1, d, model)
        abiotic1 = abiotic_fitness(ag1, model)
        abiotic2 = abiotic_fitness(ag2, model)
        ix_dir1 = sign(ix_value1)
        ix_dir2 = sign(ix_value2)

        # update agents' fitness
        ag1.W = abiotic1 + (inx_prob1 * ix_dir1)
        ag2.W = abiotic2 + (inx_prob2 * ix_dir2)
        adjust_fitness!(ag1, model)
        adjust_fitness!(ag2, model)
      end
    end
  end
end

"""
`ag1` eats `ag2`.
"""
function eat!(ag1, ag2, model)
  ag1.energy += model.food_sources[ag1.species, ag2.species]
  remove_agent!(ag2, model)
end

"""
Returns as many random ids from the agent site as the number of species.
"""
function target_species_ids(agent, model::ABM)
  nspecies = model.nspecies
  allids = collect(nearby_ids(agent, model, 0))
  foundlen = length(allids)
  if foundlen == 0
    return Int[]
  elseif foundlen < nspecies
    return allids
  else
    species_ids = sample(allids, nspecies, replace=false)
    return species_ids
  end
end

"""
    interact!(agent::Ind, model::ABM)

Includes interaction with other species such as competition, cooperation, and predation. It also includes sexual reproduction (if meeting an individual from other sex) and within species competition/cooperation (if meeting an individual from the same sex).
"""
function interact!(agent::Ind, model::ABM)
  sp = agent.species
  target_ids = target_species_ids(agent, model)
  for id in target_ids
    target = model[id]
    target_sp = target.species
    if agent.interaction_history[target_sp] != model.step[1] && target.interaction_history[sp] != model.step[1] # if agent and target have not interacted with such species before
      interact!(agent, target, model)
      # if agent was a prey, check whether it is still alive
      if model.food_sources[target_sp, sp] > 0
        if !agent.isalive
          return
        end
      end
    end
  end
end


text:
  cd(@__DIR__)
  import Pkg
  Pkg.activate(".")

  function main()
    include(joinpath("src", "MyGenieExample.jl"))
  end; main()


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control
# Author: Martin Otter, DLR-SR
#
# Provide information about the structure of the ODE or DAE equations.

import DataFrames
import OrderedCollections
using  LinearAlgebra
import TimerOutputs
import RecursiveFactorization

export StateCategory, XD, XALG, XLAMBDA, XMUE
export ResidualCategory, FD, FC_ALG, FC_LOW_HIGH, FC_MUE

export LinearEquations

export EquationInfo, StateElementInfo, update_equationInfo!, get_stateNames, get_x_table

export EquationInfoStatus, MANUAL, CODE_GENERATION, SOLVER_MODEL

"""
    @enum StateCategory

Type of a state.

| value     | description                                                                       |
|:----------|:----------------------------------------------------------------------------------|
| `XD`      | Differential variable (derivative of variable appears in equations)               |
| `XALG`    | Algebraic variable that is included in x                                          |
| `XLAMBDA` | Algebraic variable that is included in der_x                                      |
| `XMUE`    | Algebraic variable used for stabilization and included in der_x (exact value = 0) |
"""
@enum StateCategory XD=1 XALG XLAMBDA XMUE

"""
    @enum ResidualCategory

Type of a residual (see also [`StateCategory`](@ref) for the type of a state).

| value         | description                                                              |
|:--------------|:-------------------------------------------------------------------------|
| `FD`          | Differential equation (depends on XD, XLAMBDA, XMUE variables)           |
| `FC_ALG`      | Constraint equation that is purely algebraic (depends on XALG variables) |
| `FC_LOW_HIGH` | Constraint equation on lowest and on highest derivative level            |
| `FC_MUE`      | Constraint equation that is associated with a stabilizer XMUE            |
"""
@enum ResidualCategory FD=1 FC_ALG FC_LOW_HIGH FC_MUE

const niter_max = 20  # Maximum number of fixed-point iterations to solve A*x = b

"""
    leq = LinearEquations{FloatType}(vTear_names::Vector{String},  vTear_lengths::Vector{Int},
                                     nResiduals::Int, A_is_constant::Bool;
                                     useRecursiveFactorizationUptoSize = 0)

Define linear equation system "A*x=b" with `length(x) = sum(vTear_lengths)`.
If `A_is_constant = true` then `A` is a matrix that is constant after
initialization.

If length(x) <= useRecursiveFactorizationUptoSize, then linear equation systems will be solved with
`RecursiveFactorization.jl` instead of the default `lu!(..)` and `ldiv!(..)`.

For details of its usage for code generation see [`LinearEquationsIteration`](@ref).
"""
mutable struct LinearEquations{FloatType <: Real}
    A_is_constant::Bool             # = true, if A-matrix is constant
    vTear_names::Vector{String}     # The names of the tearing variables (for error message)
    vTear_lengths::Vector{Int}      # Lengths of the vTear_names elements (length(vTear_value) = sum(vTear_lengths))
    vTear_value::Vector{FloatType}  # Values of vTear variables (needed during execution of generated code)
    nResiduals::Int                 # Number of residuals

    A::Matrix{FloatType}
    b::Vector{FloatType}
    pivots::Vector{Int}             # Pivot vector if recursiveFactorization = true
    residuals::Vector{FloatType}    # Values of the residuals FloatType vector; length(residuals) = sum(residuals_length) = sum(vTear_lengths)
    residual_value::AbstractVector  # Values of the residual variables ::Vector{Any}, length(residual_values) = nResiduals

    # Iteration status of for-loop
    mode::Int       # Operational mode (see function iterate below)
    niter::Int      # Current number of iterations in the fix point iteration scheme
    niter_max::Int  # Maximum number of iterations
    success::Bool   # = true, if solution of A*x = b is successfully computed
                    # = false, if solution is not computed; continue with fix point iteration

    # For warning message if niter > niter_max
    inconsistentPositive::Vector{String}
    inconsistentNegative::Vector{String}

    # Constructed during initialization
    residual_unitRanges::Vector{UnitRange{Int}} # residuals[residual_unitRanges[i]] = residual_value[i], if residual is a vector
    residual_indices::Vector{Int}               # residuals[residual_indices[i]] = residual_value[i], if residual is a scalar
    useRecursiveFactorization::Bool             # = true, if RecursiveFactorization.jl shall be used to solve the linear equation system

    luA::LU{FloatType,Array{FloatType,2}}       # lu-Decomposition of A

    function LinearEquations{FloatType}(vTear_names::Vector{String}, vTear_lengths::Vector{Int},
                                        nResiduals::Int, A_is_constant::Bool;
                                        useRecursiveFactorizationUptoSize::Int = 0) where {FloatType <: Real}
        @assert(length(vTear_names) > 0)
        @assert(length(vTear_names) == length(vTear_lengths))
        nx = sum(vTear_lengths)
        @assert(nx > 0)
        state = -1
        useRecursiveFactorization = nx <= useRecursiveFactorizationUptoSize

        new(A_is_constant, vTear_names, vTear_lengths, zeros(FloatType,nx), nResiduals,
            zeros(FloatType,nx,nx), zeros(FloatType,nx), fill(0,nx), zeros(FloatType,nx),
            Vector{Any}(undef,nResiduals), -4, -1, niter_max, false, String[], String[],
            fill(0:0,nResiduals), fill(0,nResiduals), useRecursiveFactorization)
    end
end
LinearEquations(args...) = LinearEquations{Float64}(args...)

"""
     iterating = LinearEquationsIteration(leq::LinearEquations{FloatType}, isInitial::Bool, time, timer)

This function solves a linear equation system in residual form "residual = A*x - b"
by iterating with a while loop over this system:

```julia
function getDerivatives!(_der_x, _x, _m, _time)::Nothing
    _leq::Union{Nothing,LinearEquations{FloatType}} = nothing
    ...
    _leq      = _m.linearEquations[<nr>]   # leq::LinearEquations{FloatType}
    _leq.mode = -2  # initializes the iteration
    while LinearEquationsIteration(_leq, _m.isInitial, _m.time)
        v_tear1 = _leq.vTear_value[1:3]
        v_tear2 = _leq.vTear_value[4]
        ...
        v_solved1 = f(v_tear1, v_tear2, ..., positive(v_tear1,.., _leq))
        v_solved2 = f(v_tear1, v_tear2, ..., v_solved1)
        ...
        _leq.residuals[..] = < residual equation(v_tear1,...) >
    end
    ...
end
```
"""
function LinearEquationsIteration(leq::LinearEquations{FloatType}, isInitial::Bool, time, timer)::Bool where {FloatType}
    #=
    while LinearEquationsIteration(leq, isInitial,time)
        # body of while-loop
    end

    Solve "residuals = A*x - b" where
    A and b can be functions of positive(c_i'*x - d_i). In this case this system
    is solved by a fixed point iteration scheme, that is "A*x = b" is solved
    until potentially present positive(c_i'*x - d_i) calls are consistent to x.
    Iteration is terminated after 20 iterations if no convergence is found.
    In this case a warning message is triggered and simulation is continued
    (it might be that simulation is still successful, even if "x" and
    positive(..) are temporarily not consistent to each other).

    Note, the current values of A,x,b,residuals are stored in leq.
    If A is fixed after initialization (leq.A_is_constant = true), then A is computed
    only once at initialization, the LU decomposition of A is stored in leq
    and used in subsequent calls to solve the equation system.

        On input (mode = leq.mode, x = leq.vTear_value):

        mode = -2: # Initialize fixed point iteration
                   leq.niter   = 0
                   leq.success = false
                   leq.mode    = 0
                   x .= 0
                   return true

        mode = -1: if leq.success
                      return false   # Terminate while-loop
                   elseif leq.niter > 20
                      <warning>
                      return false   # Terminate while-loop
                   end
                   leq.mode = 0      # Re-initialize next solver-iteration
                   x .= 0
                   return true

        mode =  0: b    = -residuals  # Compute b
                   mode = 1           # Compute A[:,1] after the next iteration
                   x[1] = 1.0
                   return true

        mode >  1: j      = mode
                   A[:,j] = residuals + b  # Compute A[:,j].
                   x[j]   = 0.0
                   if j < nx
                       x[j+1]   = 1.0
                       leq.mode = j+1
                       return true         # Compute A[:,j+1] after the next iteration
                   elseif j == nx
                       x = A\b            # Solve linear equation system A*x = b for x.
                       leq.mode    = -1   # Compute all variables IN the next iteration as function of x
                       leq.niter  += 1    # Increment number of iterations to solve A*x = b
                       leq.success = true # Terminate for-loop at the beginning of the next iteration,
                                          # provided no positive(..) call changes its value.
                                          # (leq.success is set to false in positive(..), if the return value changes).
                       return true
                   end
    =#
    mode = leq.mode
    x    = leq.vTear_value
    nx   = length(x)

    if mode == -2
        # Initialize fixed point iteration
        leq.niter   = 0
        leq.success = false
        empty!(leq.inconsistentPositive)
        empty!(leq.inconsistentNegative)
        leq.mode = 0
        x .= 0
        return true

    elseif mode == -1
        # Either terminate fixed point iteration, or start a new iteration
        if leq.success
            return false
        elseif leq.niter > niter_max
            str = ""
            if length(leq.inconsistentPositive) > 0
                str = str * "positive(expr) is inconsistent for expr = $(leq.inconsistentPositive)."
            end
            if length(leq.inconsistentNegative) > 0
                if length(str) > 0
                    str = str * "\n"
                end
                str = str * "negative(expr) is inconsistent for expr = $(leq.inconsistentNegative)."
            end
            @warn "At time = $time, no consistent solution found for mixed linear equation system.\n" *
                  "Simulation is continued although some variables might not be correct at this time instant.\n$str"
            return false
        end
        leq.mode = 0
        x .= 0
        return true

    elseif mode < -2 || mode > nx
        @goto ERROR
    end

    A   = leq.A
    b   = leq.b
    nResiduals          = leq.nResiduals
    residuals           = leq.residuals
    residual_value      = leq.residual_value
    residual_unitRanges = leq.residual_unitRanges
    residual_indices    = leq.residual_indices

    if isInitial && mode == 0
        # Construct unit ranges for the residual variables vector to copy values into the residuals vector
        j = 1
        for i = 1:nResiduals
            res_value = residual_value[i]
            if typeof(res_value) <: Number
                residual_indices[i] = j
                j = j+1
            else
                len = length(res_value)
                if length(res_value) == 0
                    error("Residual variable value $i has length zero")
                end
                k = j+len-1
                residual_unitRanges[i] = j:k
                j = k+1
            end
        end

        if j-1 != nx
            len_res = j-1
            error("The length of the residuals vector (= $len_res) is not equal to the length of the vTear_value vector (= $nx)")
        end
    end

    # Copy residual variable values to residuals vector
    index = 0
    for i = 1:nResiduals
        res_value = residual_value[i]
        index     = residual_indices[i]
        if index > 0
            residuals[index] = res_value
        else
            residuals[residual_unitRanges[i]] = res_value
        end
    end

    if !leq.A_is_constant || isInitial  # A is not constant or A is constant and isInitial = true
        if mode == 0
            # residuals = A*x - b -> b = -residuals)
            for i = 1:nx
                b[i] = -residuals[i]
            end
            leq.mode = 1
            x[1] = convert(FloatType, 1)
            return true
        end

        # residuals = A*x - b -> A[:,j] = residuals + b)
        j = mode
        for i = 1:nx
            A[i,j] = residuals[i] + b[i]
        end
        x[j] = 0

        if j < nx
            leq.mode += 1
            x[leq.mode] = convert(FloatType, 1)
            return true
        end

        # Solve linear equation system
        if nx == 1
            x[1] = b[1]/A[1,1]
            if !isfinite(x[1])
                error("Linear scalar equation system is singular resulting in: ", leq.vTear_names[1], " = ", x[1])
            end
        else
            x .= b
            if leq.useRecursiveFactorization
                ModiaBase.TimerOutputs.@timeit timer "solve A*x=b (RecursiveFactorization)" begin
                    leq.luA = RecursiveFactorization.lu!(A, leq.pivots)
                    ldiv!(leq.luA, x)
                end
            else
                ModiaBase.TimerOutputs.@timeit timer "solve A*x=b (LinearAlgebra)" begin
                    leq.luA = lu!(A)
                    ldiv!(leq.luA, x)
                end
            end
        end

    elseif leq.A_is_constant && !isInitial  # isInitial=false, LU decomposition of A is available in leq.luA
        for i = 1:nx
            x[i] = -residuals[i]
        end

        # Solve linear equation system
        if nx == 1
            x[1] = x[1]/A[1,1]
            if !isfinite(x[1])
                error("Linear scalar equation system is singular resulting in: ", leq.vTear_names[1], " = ", x[1])
            end
        else
            ldiv!(leq.luA, x)
        end

    else
        @goto ERROR
    end

    # Linear equation system solved
    leq.mode    = -1   # Compute all variables IN the next iteration as function of x
    leq.niter  += 1    # Increment number of iterations to solve A*x = b
    leq.success = true # Terminate for-loop at the beginning of the next iteration,
                       # provided no positive(..) call changes its value.
                       # (leq.success is set to false in positive(..), if the return value changes).
    return true

    @label ERROR
    @error "Should not occur (Bug in file ModiaBase/src/EquationAndStateInfo.jl,\nleq.mode = $(leq.mode); vTear = $(leq.vTear_names), time=$time, isInitial=$isInitial)."
end

"""
    @enum EquationInfoStatus

Status of an EquationInfo instance:

- `MANUAL`: Is defined manually. The following variables in x_info::Vector{StateElementInfo}
   are **not** defined:
   x_names_julia, der_x_names_julia, length, unit, startIndex.
   Also variables nx and x_infoByIndex are not defined.
   With function [`update_equationInfo!`](@ref), the variables length, unit, startIndex
   in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex.

- `CODE_GENERATION`: Is constructed during code generation with getSortedAndSolvedAST(..).
  The following variables in x_info::Vector{StateElementInfo}
  are **not** defined: startIndex.
  Also variables nx and x_infoByIndex are not defined.
  With function [`update_equationInfo!`](@ref), the variables startIndex
  in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex.

- `SOLVER_MODEL`: Is used during simulation in a SolverModel. With function
  [`update_equationInfo!`](@ref), missing variables are constructed depending
  on the information given with `MANUAL` or `CODE_GENERATION` and the actual
  modelValues instance (here unit and length information is available).
"""
@enum EquationInfoStatus MANUAL=1 CODE_GENERATION SOLVER_MODEL

"""
    xe_info = StateElementInfo(...)

Return an instance of the mutable struct `StateElementInfo` that defines the information
for one element of the state vector. There are three constructors:

- Default constructor (all variables under section Arguments are given;
  used to write/read the complete information).

- All variables under section Arguments are given with exception of startIndex
  (used for `EquationInfoStatus = CODE_GENERATION`).

- StateElementInfo(x_name, der_x_name, stateCategory; fixed=nothing, nominal=NaN, unbounded=false)
  (used for `EquationInfoStatus = MANUAL`):

# Arguments
- x_name: Name of x-element or "" if no name (if stateCatebory = XLAMBDA or XMUE)
- x_name_julia: Julia name of x-element in getDerivatives!/getResiduals! function
  or "" if not needed (since no code generation).
- der_x_name: Name of der_x-element or "" if either `der(x_name)` or if no name
  (if stateCategory = XALG).
- der_x_name_julia: Julia name of der_x-element in getDerivatives!/getResiduals! function
  or "" if not needed (since no code generation).
- stateCategory::StateCategory: Category of the state
- length: length of x-element (or -1 if not yet known)
- unit: unit of XD, XALG (x_name) or XLAMBDA, XMUE (der_x_name) or "" if not yet known)
- fixed: false (= guess value) or true (= not changed by initialization).
         Only relevant for ode=false, otherwise ignored.
- nominal: Nominal value (NaN if determined via start value)
- unbounded: false or true
- startIndex: start index of x-element with respect to x-vector
              or -1 if not yet known.
"""
mutable struct StateElementInfo
    x_name::String                # Modia name of x-element or "" if no name (for )
    x_name_julia                  # Julia name of x-element in getDerivatives! function
                                  # or "", if not needed (since no code generation).
    der_x_name::String            # Modia name of der_x-element or "" if either "der(x_name)" or if no name,
    der_x_name_julia              # Julia name of der_x-element in getDerivatives! function
                                  # or not needed (since no code generation)
    stateCategory::StateCategory  # category of the state
    length::Int                   # length of x-element (or -1 if not yet known)
    unit::String                  # unit of x-element as string (or "" if not yet known)
    fixed::Bool                   # false or true
    nominal::Float64              # nominal value (NaN if determined via start value)
    unbounded::Bool               # false or true
    startIndex::Int               # start index of x-element with respect to x-vector
                                  # or -1 if not yet known.
end

# Constructor for code-generation
StateElementInfo(x_name, x_name_julia, der_x_name, der_x_name_julia,
                 stateCategory, length, unit, fixed, nominal, unbounded) = StateElementInfo(
                 x_name, x_name_julia, der_x_name, der_x_name_julia,
                 stateCategory, length, unit, fixed, nominal, unbounded, -1)

# Constructor for reading StateElementInfo after code-generation (x_name_julia and der_x_name_julia are not included
StateElementInfo(x_name, der_x_name,
                 stateCategory, length, unit, fixed, nominal, unbounded) = StateElementInfo(
                 x_name, :(), der_x_name, :(),
                 stateCategory, length, unit, fixed, nominal, unbounded, -1)

# Constructor for manually generated equation info.
StateElementInfo(x_name, der_x_name, stateCategory=XD; fixed=false, nominal=NaN, unbounded=false) =
    StateElementInfo(x_name, :(), der_x_name, :(), stateCategory, -1, "", fixed, nominal, unbounded, -1)

function Base.show(io::IO, xe_info::StateElementInfo)
    print(io, "ModiaBase.StateElementInfo(")
    show( io, xe_info.x_name)
	#print(io, ",")
    #show( io, xe_info.x_name_julia)
	print(io, ",")
    show( io, xe_info.der_x_name)
	#print(io, ",")
    #show( io, xe_info.der_x_name_julia)
    print(io, ",ModiaBase.", xe_info.stateCategory)
    print(io, ",", xe_info.length)
	print(io, ",")
    show( io, xe_info.unit)
    print(io, ",", xe_info.fixed)
    print(io, ",", xe_info.nominal)
    print(io, ",", xe_info.unbounded)
    #if xe_info.startIndex > 0
    #    print(io, ",", xe_info.startIndex)
    #end
    print(io, ")")
    return nothing
end

"""
    x_table = get_x_table(x_info::Vector{StateElementInfo}

Return the state element info as DataFrames.DataFrame table, for example
to print it as:

```
show(x_table, allrows=true, allcols=true, summary=false, eltypes=false)
```
"""
function get_x_table(x_info::Vector{StateElementInfo})
    x_table = DataFrames.DataFrame(name=String[], length=Int[], unit=String[], fixed=Bool[], nominal=Float64[], unbounded=Bool[])

    for xe_info in x_info
        push!(x_table, (xe_info.x_name, xe_info.length, xe_info.unit, xe_info.fixed, xe_info.nominal, xe_info.unbounded))
    end

    return x_table
end

"""
    eqInfo = EquationInfo(;
                status               = MANUAL,
                ode                  = true,
                nz                   = 0,
                x_info               = StateElementInfo[],
                residualCategories   = ResidualCategory[],
                linearEquations      = Tuple{Vector{String},Bool}[],
                vSolvedWithFixedTrue = String[],
                defaultParameterAndStartValues = nothing,
                ResultType = nothing,
                ResultTypeHasFloatType = false)

Return instance `eqInfo` that defines the information for the equation system.

# Arguments
- status::EquationInfoStatus: Defines the variables that have a value.
- ode: = true if ODE-interface (`getDerivatives!`),
       = false if DAE-interface (`getResiduals!`).
- nz: Number of zero crossing functions.
- x_info: Vector of StateElementInfo elements provding info for every x-element
- residualCategories: If ode=true, length(residualCategories) = 0.
             If ode=false: residualCategories[i] is the `ResidualCategory`](@ref) of x-element "i".
- `linearEquations::Vector{Tuple{Vector{String},Vector{Int},Int,Bool}}`:
               linearEquations[i] defines a
               ModiaBase.LinearEquations system, where the first tuple value
               is a vector of the names of the unknowns, the second tuple value
               is a vector with the lengths of the unknowns, the third tuple value is the number
               of residuals and the fourth tuple value
               defines whether the coefficient matrix A
               has only constant entries (=true) or not (=false).
- `vSolvedWithFixedTrue::Vector{String}`: Vector of variables that are computed
               from other variables and have `fixed=true`. During initialization
               it is checked whether the calculated values and the start values of
               these variables agree. If this is not the case, an error is triggered.
- `defaultParameterAndStartValues`: Dictionary of default paramter and default start values.
- `ResultType::AbstractModelValues`: ModelValues type that shall be used for the result generation
               (= ModelValues struct without parameters).
- `ResultTypeHasFloatType::Bool`: = false, if `ResultType` is not parameterized.
                                  = true, if `ResultType` has parameter `FloatType`.
"""
mutable struct EquationInfo
    status::EquationInfoStatus
    ode::Bool                                      # = true if ODE model interface, otherwise DAE model interface
    nz::Int                                        # Number of crossing functions
    x_info::Vector{StateElementInfo}
    residualCategories::Vector{ResidualCategory}   # If ode=true, length(residualCategories) = 0
                                                   # If ode=false, residualCategories[j] is the ResidualCategory of residual[j]
    linearEquations::Vector{Tuple{Vector{String},Vector{Int},Int,Bool}}
    vSolvedWithFixedTrue::Vector{String}
    nx::Int                                        # = length(x) or -1 if not yet known
    x_infoByIndex::Vector{Int}                     # i = x_infoByIndex[j] -> x_info[i]
                                                   # or empty vector, if not yet known.
    x_dict::OrderedCollections.OrderedDict{String,Int} # x_dict[x_name] returns the index of x_name with respect to x_info
    defaultParameterAndStartValues::Union{AbstractDict,Nothing}
    ResultType
    ResultTypeHasFloatType::Bool
end

EquationInfo(; status                = MANUAL,
               ode                   = true,
               nz                    = 0,
               x_info                = StateElementInfo[],
               residualCategories    = ResidualCategory[],
               linearEquations       = Tuple{Vector{String},Vector{Int},Int,Bool}[],
               vSolvedWithFixedTrue  = String[],
               nx                    = -1,
               x_infoByIndex         = Int[],
               defaultParameterAndStartValues::Union{AbstractDict,Nothing} = nothing,
               ResultType = nothing,
               ResultTypeHasFloatType = false) =
               EquationInfo(status, ode, nz, x_info,
                            residualCategories, linearEquations,
                            vSolvedWithFixedTrue, nx, x_infoByIndex,
                            OrderedCollections.OrderedDict{String,Int}(),
                            defaultParameterAndStartValues,
                            ResultType, ResultTypeHasFloatType)

"""
    updateEquationInfo!(eqInfo::EquationInfo)

Set eqInfo.x_dict, eqInfo.nx and eqInfo.x_info[:].startIndex
"""
function updateEquationInfo!(eqInfo::EquationInfo)::Nothing
    x_dict = eqInfo.x_dict
    startIndex = 1
    for (i, xi_info) in enumerate(eqInfo.x_info)
        x_dict[xi_info.x_name] = i
        xi_info.startIndex = startIndex
        startIndex += xi_info.length
    end
    eqInfo.nx = startIndex - 1
    return nothing
end

get_x_startIndexAndLength(eqInfo::EquationInfo, name::String) = begin
    xe_info = eqInfo.x_info[ eqInfo.x_dict[name] ]
    (xe_info.startIndex, xe_info.length)
end

function Base.show(io::IO, eqInfo::EquationInfo; indent=4)
    indentation  = repeat(" ", indent)
    indentation2 = repeat(" ", 2*indent)
    indentation3 = repeat(" ", 3*indent)
    println(io, "ModiaBase.EquationInfo(")
    println(io, indentation2, "status = ModiaBase.", eqInfo.status, ",")
    println(io, indentation2, "ode = ", eqInfo.ode, ",")
    if eqInfo.nz > 0
        println(io, indentation2, "nz = ", eqInfo.nz, ",")
    end
    print(io, indentation2, "x_info = ModiaBase.StateElementInfo[")
    for (i, xe_info) in enumerate(eqInfo.x_info)
        if i == 1
            print(io, "\n", indentation3)
        else
            print(io, ",\n", indentation3)
        end
        show(io, xe_info)
    end
    print(io,"]")

    if length(eqInfo.residualCategories) > 0
        print(io, ",\n", indentation2, "residualCategories = [")
        for (i, rcat) in enumerate(eqInfo.residualCategories)
            if i > 1
                print(io, ",")
            end
            show(io, rcat)
        end
        print(io,"]")
    end

    leqs = eqInfo.linearEquations
    if length(leqs) > 0
        println(io, ",\n", indentation2, "linearEquations = [")
        for (i, leq) in enumerate(leqs)
            print(io, indentation3, "(", leq[1], ",\n",
                      indentation3, " ", leq[2], ",\n",
                      indentation3, " ", leq[3], ", ", leq[4], ")")
            if i < length(leqs)
                print(io, ",\n")
            end
        end
        print(io, "]")
    end

    if length(eqInfo.vSolvedWithFixedTrue) > 0
        print(io, ",\n", indentation2, "vSolvedWithFixedTrue = ")
        show(io, eqInfo.vSolvedWithFixedTrue)
    end

    if eqInfo.nx > 0
        print(io, ",\n", indentation2, "nx = ", eqInfo.nx)
    end

    if length(eqInfo.x_infoByIndex) > 0
        print(io, ",\n", indentation2, "x_infoByIndex = ")
        show(io, eqInfo.x_infoByIndex)
    end

    if !isnothing(eqInfo.defaultParameterAndStartValues)
        print(io, ",\n", indentation2, "defaultParameterAndStartValues = ")
        show(io, eqInfo.defaultParameterAndStartValues, indent=12, finalLineBreak=false)
    end

    if !isnothing(eqInfo.ResultType)
        print(io, ",\n", indentation2, "ResultType = ")
        show(io, eqInfo.ResultType)

        if eqInfo.ResultTypeHasFloatType
            print(io, ",\n", indentation2, "ResultTypeHasFloatType = ")
            show(io, eqInfo.ResultTypeHasFloatType)
        end
    end

    println(io, "\n", indentation, ")")
    return nothing
end

"""
    names = get_stateNames(equationInfo::EquationInfo)

Return the names of the states defined in `equationInfo` as a Vector of strings.
"""
get_stateNames(eqInfo::EquationInfo) = String[xi_info.x_name for xi_info in eqInfo.x_info]

"""
    names = get_xNames(eqInfo::EquationInfo)

Return the names of all elements of the x-vector as a vector of strings.
"""
function get_xNames(eqInfo::EquationInfo)::Vector{String}
    xNames = Vector{String}(undef, eqInfo.nx)
    for xe_info in eqInfo.x_info
        if xe_info.length == 1
            xNames[xe_info.startIndex] = xe_info.x_name
        else
            for i = 1:xe_info.length
                xNames[xe_info.startIndex+i-1] = xe_info.x_name*"["*string(i)*"]"
            end
        end
    end
    return xNames
end

#=
"""
    update_equationInfo!(eqInfo::EquationInfo)

Update equationInfo as needed for [`EquationInfoStatus`](@ref)` = SOLVER_MODEL`.
"""
function update_equationInfo!(eqInfo::EquationInfo, modelValues::ModiaBase.AbstractModelValues)::Nothing
    x_info = eqInfo.x_info

    if eqInfo.status == MANUAL
        # Determine length and unit
		count = 0
        for xi_info in x_info
            if xi_info.x_name != "" # XD or XALG
                (component,name) = ModiaBase.get_modelValuesAndName(modelValues, xi_info.x_name)
            elseif xi_info.der_x_name != ""  # XLAMBDA or XMUE
                (component,name) = ModiaBase.get_modelValuesAndName(modelValues, xi_info.der_x_name)
            else
				# Pure algebraic with dummy equation der(x) = -x, x(0) = 0
				xi_info.length = 1
				count = count+1
				continue
            end

	        xi_type        = fieldtype(typeof(component), name)
			xi_numberType  = ModiaBase.numberType(xi_type)
            xi_info.unit   = replace(string(unit(xi_numberType)), " " => "*")
            xi_info.length = Base.length(xi_type)  # Base.length(::Type{<:Number})=1 is defined in ModiaBase, since not defined in Base
        end
		if count == 1 && length(x_info) > 1 || count >  1
			error("Error in update_equationInfo!: x_info is wrong.\n",
				  "x_info = $x_info")
		end
    end

    # Set startIndex in x_info and compute nx
    startIndex = 1
    for xi_info in x_info
        xi_info.startIndex = startIndex
        startIndex += xi_info.length
    end
    nx = startIndex - 1
    eqInfo.nx = nx

    # Set x_infoByIndex::Vector{Int}
    eqInfo.x_infoByIndex = zeros(Int,nx)
    i1 = 0
    i2 = 0
    for (i,xi_info) in enumerate(x_info)
        i1 = xi_info.startIndex
        i2 = i1 + xi_info.length - 1
        for j = i1:i2
            eqInfo.x_infoByIndex[j] = i
        end
    end

    eqInfo.status = SOLVER_MODEL

    return nothing
end
=#

max_stars_repo_name: ModiaSim/ModiaBase

max_issues_repo_name: ModiaSim/ModiaBase

max_forks_repo_name: ModiaSim/ModiaBase


text:
struct BundleMethod <: UnconstrainedNonSmoothSolver
  solver_factory # Must support QPs with linear constraints.
  µ::Float64 # Used in the proximal step. Considered to be constant for now. TODO: schedules
  K::Float64
  n_iter::Int
end

function solve(p::UnconstrainedNonSmoothProblem, params::BundleMethod, x0::Vector{Float64}; info_callback::Union{Nothing, Function}=nothing)
  @assert p.sense == Minimise # TODO: current limitation.
  @assert length(x0) == p.dimension

  # Build the bundle model.
  m = JuMP.Model(params.solver_factory)
  f = @JuMP.variable(m)
  bx = @JuMP.variable(m, [1:p.dimension])

  # Create the first cutting plane with the starting iterate (i.e. the program will never be unbounded).
  @JuMP.constraint(m, f >= p.f(x0) + dot(p.g(x0), bx - x0))

  # Start iterating.
  x = copy(x0)
  for k in 1:params.n_iter
    t0 = time_ns()

    # Get the new test point (proximal step).
    @JuMP.objective(m, Min, f + params.µ * sum((bx[i] - x[i])^2 for i in 1:3))
    JuMP.optimize!(m)
    @assert JuMP.termination_status(m) == MOI.OPTIMAL

    y = value.(bx)
    @JuMP.constraint(m, f >= p.f(y) + dot(p.g(y), bx - y))

    # Is this point optimum?
    @JuMP.objective(m, Min, f)
    c = @JuMP.constraint(m, bx .== y)
    JuMP.optimize!(m)
    @assert JuMP.termination_status(m) == MOI.OPTIMAL

    f_bundle_y = JuMP.value(f)
    JuMP.delete(m, c)

    f = p.f(x)
    v = f_bundle_y - f
    if v >= 0
      break
    end

    # Perform a step:
    if p.f(y) <= f + params.K * v
      # Descent step.
      x = y
    else
      # Short step.
      nothing
    end

    if info_callback !== nothing
      info_callback(k, f, x, g, f_best, x_best, time_ns() - t0)
    end
  end

  return x
end


text:
################################################################################
# Particle
################################################################################
dir = joinpath(@__DIR__, "particle")
model = deepcopy(particle)

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Particle 2D
################################################################################
dir = joinpath(@__DIR__, "particle_2D")
model = deepcopy(particle_2D)

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Hopper (2D)
################################################################################
dir = joinpath(@__DIR__, "hopper_2D")
model = deepcopy(hopper_2D)

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=true)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Hopper (3D)
################################################################################
dir = joinpath(@__DIR__, "hopper_3D")
model = deepcopy(hopper_3D)

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=true)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Quadruped
################################################################################
dir = joinpath(@__DIR__, "quadruped")
model = deepcopy(quadruped)
path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=false)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Quadruped Payload
################################################################################
dir = joinpath(@__DIR__, "quadruped")
model = deepcopy(quadruped_payload)
path_base = joinpath(dir, "dynamics_payload/base.jld2")
path_dyn = joinpath(dir, "dynamics_payload/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=false)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# Flamingo
################################################################################
dir = joinpath(@__DIR__, "flamingo")
model = deepcopy(flamingo)
path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=false)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)

################################################################################
# PushBot
################################################################################
dir = joinpath(@__DIR__, "pushbot")
model = deepcopy(pushbot)

path_base = joinpath(dir, "dynamics/base.jld2")
path_dyn = joinpath(dir, "dynamics/dynamics.jld2")

expr_base = generate_base_expressions(model, M_analytical=true, C_analytical=false)
save_expressions(expr_base, path_base, overwrite=true)
instantiate_base!(model, path_base)

expr_dyn = generate_dynamics_expressions(model)
save_expressions(expr_dyn, path_dyn, overwrite=true)
instantiate_dynamics!(model, path_dyn)


text:
"""
    StaticLineSearch

`StaticLineSearch` directly uses the initial guess as the new step length.
"""
struct StaticLineSearch end

(sls::StaticLineSearch)(f, state, p, α₀) = α₀


text:
function debugging end
enable_debug() = @eval debugging() = true
disable_debug() = @eval debugging() = false
enable_debug()
# disable_debug()

macro _assert(args...)
    call = Expr(
        :macrocall,
        getfield(ArgCheck, Symbol("@check")),  # injecting the callable for esc
        __source__,
        args...,
    )
    ex = Expr(:block, __source__, call)
    quote
        if $debugging()
            $ex
        end
        nothing
    end |> esc
end

function ceillog2(n::Integer)
    n > 0 || throw(DomainError(n))
    i = trailing_zeros(n)
    j = 8 * sizeof(n) - leading_zeros(n) - 1
    if i == j
        return i
    else
        return j + 1
    end
end

abstract type Cell{T} end
mutable struct MutableCell{T} <: Cell{T}
    value::T
    # TODO: pad
    MutableCell{T}() where {T} = new{T}()
end
struct SingletonCell{T} <: Cell{T} end

@inline Base.getindex(cell::MutableCell{T}) where {T} = cell.value
@inline Base.getindex(::SingletonCell{T}) where {T} = T.instance
@inline Base.setindex!(cell::MutableCell{T}, value::T) where {T} = cell.value = value
@inline Base.setindex!(::SingletonCell{T}, value::T) where {T} = value

USE_PADDED_ATOMICS = true
# USE_PADDED_ATOMICS = false
if USE_PADDED_ATOMICS
    # TODO: use RecordArrays to allocate atomics
    # TODO: load cache line size
    const ATOMICS_NPADS = 7

    # The allocator is very likely to hand over objects allocated at consecutive locations
    function cache_aligned_atomic(x::T) where {T}
        a1 = Threads.Atomic{T}(x)
        a2 = Threads.Atomic{T}(x)
        a3 = Threads.Atomic{T}(x)
        a4 = Threads.Atomic{T}(x)
        a5 = Threads.Atomic{T}(x)
        a6 = Threads.Atomic{T}(x)
        a7 = Threads.Atomic{T}(x)
        a8 = Threads.Atomic{T}(x)
        if mod(UInt(pointer_from_objref(a1)), 64) == 0
            return (a1, (a2, a3, a4, a5, a6, a7, a8))
        elseif mod(UInt(pointer_from_objref(a2)), 64) == 0
            return (a2, (a1, a3, a4, a5, a6, a7, a8))
        elseif mod(UInt(pointer_from_objref(a3)), 64) == 0
            return (a3, (a1, a2, a4, a5, a6, a7, a8))
        else
            return (a4, (a1, a2, a3, a5, a6, a7, a8))
        end
    end
else
    const ATOMICS_NPADS = 0
    cache_aligned_atomic(x::T) where {T} = (Threads.Atomic{T}(x), ())
end

""" Something like Transducers.Reduced """
struct Break{T}
    value::T
end

right(_, x) = x

mutable struct ExperimentalMutableNTuple{N,T}
    values::NTuple{N,T}
    ExperimentalMutableNTuple{N,T}(values::NTuple{N,T}) where {N,T} = new{N,T}(values)
    ExperimentalMutableNTuple{N,T}() where {N,T} = new{N,T}()
end

@inline Base.getindex(mut::ExperimentalMutableNTuple, i) = mut.values[i]
# TODO: need to use unsafe_store! to make it data race-free
#=
@inline function Base.setindex!(mut::ExperimentalMutableNTuple{N,T}, v, i) where {T}
    @boundscheck 1 <= i <= N || throw(BoundsError(mut, i))
    if Base.issingletontype(T)
        return
    elseif Base.isimmutable(T)
    else
    end
    mut.values = Base.setindex(mut.values, v, i)
end
=#

struct FallbackMutableNTuple{N,T}
    values::Vector{T}
    FallbackMutableNTuple{N,T}(values::NTuple{N,T}) where {N,T} = new{N,T}(collect(values))
    FallbackMutableNTuple{N,T}() where {N,T} = new{N,T}(Vector{T}(undef, N))
end

Base.@propagate_inbounds Base.getindex(mut::FallbackMutableNTuple, i) = mut.values[i]
Base.@propagate_inbounds Base.setindex!(mut::FallbackMutableNTuple, v, i) =
    mut.values[i] = v

USE_EXPERIMENTAL_MUTABLE_NTUPLE = false
if USE_EXPERIMENTAL_MUTABLE_NTUPLE
    const MutableNTuple = ExperimentalMutableNTuple
else
    const MutableNTuple = FallbackMutableNTuple
end

function define_docstrings()
    docstrings = [:SyncBarriers => joinpath(dirname(@__DIR__), "README.md")]
    docsdir = joinpath(@__DIR__, "docs")
    for filename in readdir(docsdir)
        stem, ext = splitext(filename)
        ext == ".md" || continue
        name = Symbol(stem)
        name in names(SyncBarriers, all=true) || continue
        push!(docstrings, name => joinpath(docsdir, filename))
    end
    for (name, path) in docstrings
        include_dependency(path)
        doc = read(path, String)
        doc = replace(doc, r"^```julia"m => "```jldoctest $name")
        doc = replace(doc, "<kbd>TAB</kbd>" => "_TAB_")
        @eval SyncBarriers $Base.@doc $doc $name
    end
end


text:
import LightGraphs

function prufer_code_to_tree(p)
    n = length(p)

    deg = [sum(p .== i) + 1 for i=1:(n+2)]
    G = LG.SimpleGraph(n+2)
    for i in p
        for j=1:(n+2)
            if deg[j] == 1
                LG.add_edge!(G, i,j)
                deg[i] -= 1
                deg[j] -= 1
                break
            end
        end
    end

    u = 0
    v = 0
    for i=1:(n+2)
        if deg[i] == 1
            if u==0
                u = i
            else
                v = i
                break
            end
        end
    end
    LG.add_edge!(G,u,v)

    G
end

function get_all_undirected_trees(num_verts)
    possible_codes = collect(Base.Iterators.product([1:num_verts for _ in 1:(num_verts-2)]...))
    possible_graphs = prufer_code_to_tree.(possible_codes)[:]
end

export prufer_code_to_tree, get_all_undirected_trees


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

mutable struct IoK8sApiAppsV1Deployment <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiAppsV1DeploymentSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiAppsV1DeploymentStatus } # spec name: status

    function IoK8sApiAppsV1Deployment(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1Deployment, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAppsV1Deployment, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAppsV1Deployment, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiAppsV1Deployment, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiAppsV1Deployment, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiAppsV1Deployment

const _property_map_IoK8sApiAppsV1Deployment = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiAppsV1Deployment = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiAppsV1DeploymentSpec", Symbol("status")=>"IoK8sApiAppsV1DeploymentStatus")
Base.propertynames(::Type{ IoK8sApiAppsV1Deployment }) = collect(keys(_property_map_IoK8sApiAppsV1Deployment))
Swagger.property_type(::Type{ IoK8sApiAppsV1Deployment }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1Deployment[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1Deployment }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1Deployment[property_name]

function check_required(o::IoK8sApiAppsV1Deployment)
    true
end

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


text:
#!/usr/bin/env julia

using ComplexSystems
using Plots
using StatsBase

avg_p_1 = Array{Tuple{Float64,Float64}}(undef,11)
for p in 0:0.1:1
    avg = Array{Float64}(undef,30)
    for j=1:30
        space = initialize_population(200,200,p)
        for i=1:100
            sznajd_one_site(space,p)
        end
        avg[j]=sum(space)/length(space)
    end
    avg_p_1[Int(p*10+1)]=mean_and_std(avg)
end

avg_p_1_stupid = Array{Tuple{Float64,Float64}}(undef,11)
for p in 0:0.1:1
    avg = Array{Float64}(undef,30)
    for j=1:30
        space = initialize_population(200,200,p)
        for i=1:100
            if i in [50,60,70,80,90]
                stupid_bots(space,0.1)
            end
            sznajd_one_site(space,p)
        end
        avg[j]=sum(space)/length(space)
    end
    avg_p_1_stupid[Int(p*10+1)]=mean_and_std(avg)
end

avg_p_1_smart = Array{Tuple{Float64,Float64}}(undef,11)
for p in 0:0.1:1
    avg = Array{Float64}(undef,30)
    for j=1:30
        space = initialize_population(200,200,p)
        for i=1:100
            if i in [50,60,70,80,90]
                smart_bots(space,0.1)
            end
            sznajd_one_site(space,p)
        end
        avg[j]=sum(space)/length(space)
    end
    avg_p_1_smart[Int(p*10+1)]=mean_and_std(avg)
end

avg_p_1_cambridge = Array{Tuple{Float64,Float64}}(undef,11)
for p in 0:0.1:1
    avg = Array{Float64}(undef,30)
    for j=1:30
        space = initialize_population(200,200,p)
        for i=1:100
            if i == 75
                cambridge_analytica(space,0.1)
            end
            sznajd_one_site(space,p)
        end
        avg[j]=sum(space)/length(space)
    end
    avg_p_1_cambridge[Int(p*10+1)]=mean_and_std(avg)
end

theavg=map(x -> x[1],avg_p_1) |> collect
thetime=0:0.1:1|>collect
σ=map(x->x[2],avg_p_1)|> collect

theavg2=map(x -> x[1],avg_p_2) |> collect
σ2=map(x->x[2],avg_p_2)|> collect

theavg3=map(x -> x[1],avg_p_3) |> collect
σ3=map(x->x[2],avg_p_3)|> collect

plot(thetime,theavg,ribbon=σ,fillalpha=.5,lab="1 site")
plot!(thetime,theavg2,ribbon=σ2,fillalpha=.5,lab="2 site")
plot!(thetime,theavg3,ribbon=σ3,fillalpha=.5,lab="3 site")


text:
# This is for comparison with Leo Breiman paper "RANDOM FORESTS"
using UrlDownload
using DataFrames
using CSV
using StableRNGs
using StatsBase

datasets = [
   # "glass" dataset
   (url = "http://archive.ics.uci.edu/ml/machine-learning-databases/glass/glass.data",
    kwargs = [:header => ["Id", "RI", "Na","Mg", "Al", "Si", "K", "Ca", "Ba", "Fe", "glass_type"]],
    transform = x -> select!(x, Not(1)),
    out = "uci_glass",
    split = true),

   # "sonar" dataset
   (url = "https://archive.ics.uci.edu/ml/machine-learning-databases/undocumented/connectionist-bench/sonar/sonar.all-data",
    kwargs = [:header => false],
    transform = identity,
    out = "uci_sonar",
    split = true),

   # "image" dataset
   (url = "http://archive.ics.uci.edu/ml/machine-learning-databases/image/segmentation.data",
      kwargs = [:skipto => 6, :header => replace.(lowercase.(["target", "REGION-CENTROID-COL", "REGION-CENTROID-ROW", "REGION-PIXEL-COUNT", "SHORT-LINE-DENSITY-5", "SHORT-LINE-DENSITY-2", "VEDGE-MEAN", "VEDGE-SD", "HEDGE-MEAN", "HEDGE-SD", "INTENSITY-MEAN", "RAWRED-MEAN", "RAWBLUE-MEAN", "RAWGREEN-MEAN", "EXRED-MEAN", "EXBLUE-MEAN", "EXGREEN-MEAN", "VALUE-MEAN", "SATURATION-MEAN", "HUE-MEAN"]), "-" => "_")],
      transform = x -> select!(x, [names(x)[2:end]..., names(x)[1]]),
      out = "uci_images_train.csv",
      split = false
   ),

   (url = "http://archive.ics.uci.edu/ml/machine-learning-databases/image/segmentation.test",
      kwargs = [:skipto => 6, :header => replace.(lowercase.(["target", "REGION-CENTROID-COL", "REGION-CENTROID-ROW", "REGION-PIXEL-COUNT", "SHORT-LINE-DENSITY-5", "SHORT-LINE-DENSITY-2", "VEDGE-MEAN", "VEDGE-SD", "HEDGE-MEAN", "HEDGE-SD", "INTENSITY-MEAN", "RAWRED-MEAN", "RAWBLUE-MEAN", "RAWGREEN-MEAN", "EXRED-MEAN", "EXBLUE-MEAN", "EXGREEN-MEAN", "VALUE-MEAN", "SATURATION-MEAN", "HUE-MEAN"]), "-" => "_")],
      transform = x -> select!(x, [names(x)[2:end]..., names(x)[1]]),
      out = "uci_images_test.csv",
      split = false
   ),

   # "sat-images" dataset
   (url = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/satimage/sat.trn",
      kwargs = [:header => false, :delim => ' '],
      transform = identity,
      out = "uci_satimage_train.csv",
      split = false),

   (url = "https://archive.ics.uci.edu/ml/machine-learning-databases/statlog/satimage/sat.tst",
      kwargs = [:header => false, :delim => ' '],
      transform = identity,
      out = "uci_satimage_test.csv",
      split = false),

   # "letters" dataset
   (url = "https://archive.ics.uci.edu/ml/machine-learning-databases/letter-recognition/letter-recognition.data",
      kwargs = [:header => false],
      transform = x -> select!(x, [names(x)[2:end]..., names(x)[1]]),
      out = "uci_letters",
      split = true),

   # "iris" dataset
   (url = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data",
      kwargs = [:header => false],
      transform = identity,
      out = "uci_iris_all.csv",
      split = false
   )
]

function load_dataset(dataset)
   df = urldownload(dataset.url,
      format = :CSV; dataset.kwargs...) |> DataFrame
   dataset.transform(df)
   if dataset.split
      rng = StableRNG(2020)
      test_ids = sample(rng, 1:size(df, 1), Int(floor(0.1*size(df, 1))), replace = false)
      train_ids = setdiff(1:size(df, 1), test_ids)
      CSV.write(joinpath(@__DIR__, "data", dataset.out*"_test.csv"), df[test_ids, :])
      CSV.write(joinpath(@__DIR__, "data", dataset.out*"_train.csv"), df[train_ids, :])
   else
      CSV.write(joinpath(@__DIR__, "data", dataset.out), df)
   end
end

for dataset in datasets
   load_dataset(dataset)
end

####
# letter recognition
load_dataset(datasets[end])


text:
function parseInput(filename)
    values = readlines(filename*".txt")
    polymer = values[1]
    replacement = Dict{String, String}()
    for i in 3:length(values)
        st = values[i]
        st = split(st, "->")
        first = strip(st[1])
        scd = strip(st[2])
        replacement[first] = scd
    end
    return polymer, replacement
end

function applyReplacement(st, repval)

end

function solve(filename)
    polymer, dict = parseInput(filename)
    steps = 40
    for _ in 1:steps
        newpol = ""
        for p in 1:length(polymer)-1
            pair = polymer[p:p+1]
            rep = dict[pair]
            newpol = newpol *pair[1]*rep
            if p == length(polymer) - 1
                newpol = newpol * pair[2]
            end
        end
        polymer = newpol
    end
    #println(polymer)
    countdict = Dict{String,Int}()
    for char in polymer
        c = count(x -> x == char, polymer)
        countdict[string(char)] = c
    end
    mincount = minimum(values(countdict))
    maxcount = maximum(values(countdict))
    res = maxcount - mincount
    println(res)
end

function solve2(filename)
    polymer, dict = parseInput(filename)
    steps = 40
    k = collect(keys(dict))
    pairs = Dict{String, Int}(i => 0 for i in k)
    for i in 1:length(polymer)-1
        pair = polymer[i:i+1]
        pairs[pair] += 1
    end
    for _ in 1:steps
        change = Dict{String,Int}(i => 0 for i in k)
        for (p,c) in pairs
            oldpair = p
            rep = dict[p]
            new1 = p[1]*rep
            new2 = rep*p[2]
            change[oldpair] -= c
            change[new1] += c
            change[new2] += c
        end
        for (k, v) in pairs
            pairs[k] += change[k]
        end
    end
    letters = unique(collect(values(dict)))
    countdict = Dict{String, Int}(l => 0 for l in letters)
    for p in pairs
        pstring = p[1]
        countdict[string(pstring[1])] += p[2]
    end
    last = string(polymer[length(polymer)])
    countdict[last] += 1
    mincount = minimum(values(countdict))
    maxcount = maximum(values(countdict))
    res = maxcount - mincount
    println(res)
end

#solve("example")
#solve("input")
solve2("example")
solve2("input")


text:
"""
   DifferentialExtinction
   Original Python code (solar_diffex.py) by Andrea Lin
   Ported to Julia by Eric Ford
   Based on Collier-Cameron et al. (2019) MNRAS 487, 1082.  doi:10.1093/mnras/stz1215
   See also Astronomical Algorithms (Meeus, 1991) p151 & p178
"""
__precompile__() # this module is safe to precompile
module DifferentialExtinction
using PyCall
using Dates
using DataFrames, CSV

const JplHorizons = PyNULL()

function __init__()
	copy!(JplHorizons , pyimport("astroquery.jplhorizons") )
end

import ..NeidSolarScripts
import ..NeidSolarScripts.SolarRotation

AU = 149597870.7  # km
R_sun = 695700.0  # km # should we use better, but non-standard value, e.g, 695660?
valid_obs = [:WIYN, :HARPSN]

""" get_obs_loc(obs::Symbol)
 Returns a Dict with long & lat (degrees) and elevation (km)
 Warning: Currently only has info for :WIYN and :HARPSN.
"""
function get_obs_loc(obs::Symbol)
	@assert obs ∈ valid_obs
	if obs == :HARPSN
		loc = Dict("lon"=> -17.88905, "lat"=> 28.754, "elevation"=> 2.3872)
	elseif obs == :WIYN
		loc = Dict("lon"=> -111.600562, "lat"=> 31.958092, "elevation"=> 2.091)
	else
		@error("Don't have coordinates for obs = " * string(obs))
	end
	return loc
end

f_z(alt) = π/2-alt

function f_airmass(z)
	secz = sec(z)
	return secz * (1.0 - 0.0012*(secz^2 - 1.0))
end

"""fractional brightness diff (from sun center)"""
function f_epsilon(z_center, dist_sun) #
        z_upper = z_center - asin(R_sun/dist_sun)
        z_lower = z_center + asin(R_sun/dist_sun)
        delta_airmass = f_airmass(z_upper) - f_airmass(z_lower)
        k_ex = 0.157482  # WARNING: Hardwired for WIYN.  TODO: Should recalc for each obs?
		#=
		 X = [ones(length(eph.columns["airmass"])) eph.columns["airmass"]]
		 y =  eph.columns["magextinct"]
		 coeff = X \y
		 k_ex = coeff[2]
		 =#
		return k_ex * delta_airmass/2.0 * log(10.)/2.5
end

"parallactic angle"
function f_q(H_sun, z_center, lat)
	#lat = deg2rad(wiyn["lat"]) #* u.deg # observer's latitude # no longer hard-wired
	return asin(sin(H_sun) * cos(lat)/sin(z_center))
end

f_phi(pole_angle, q) = (pole_angle - q)

function f_pole_angle_ecl(alpha_sun)
	e_tilt = deg2rad(23.4) # ecliptic tilt
    asin(sin(e_tilt) * cos(alpha_sun))
end

f_psi(pole_angle_ecl, q) = (pole_angle_ecl - q)

 """
 f_B0(JD)
  see Astronomical Algorithms (Meeus, 1991) p151 & p178
 """
 function f_B0(JD)

    T = (JD - 2451545.) / 36525.
    L0 = deg2rad(280.46645 + 36000.76983*T + 0.0003032*T^2)
    MA = deg2rad(357.52910 + 35999.05030*T - 0.0001559*T^2 - 0.00000048*T^3)
    sun_center = (1.914600 - 0.004817*T - 0.000014*T^2) * sin(MA) +
    			 (0.019993 - 0.000101*T) * sin(2*MA) +
    			 0.000290 * sin(3*MA)
	sun_center = deg2rad(sun_center)
    sun_true_lon = L0 + sun_center
    lam = sun_true_lon - deg2rad(0.00569) - deg2rad(0.00478 * sin(deg2rad(125.04 - 1934.136*T)))
    kappa = deg2rad(73.6667 + 1.3958333 * (JD - 2396758.) / 36525.)
    return asin(sin(lam - kappa) * sin(deg2rad(7.25)))
end

function diff_extinct_corr(JD, alt, HA, RA, pole_angle, dist_sun, lat)
    z = f_z(deg2rad(alt))
    epsilon = f_epsilon(z, dist_sun)
    q = f_q(deg2rad(HA*15), z, deg2rad(lat))

    phi = (f_phi(deg2rad(pole_angle), q))
    psi = f_psi(f_pole_angle_ecl(deg2rad(RA)), q)

    inc_sol_axis = f_B0(JD)

	ld_coeff = 0.6 # limb darkening coeff (linear)

	# sidereal diff rotation (Snodgrass & Ulrich 1990)
	sr_A = 2.972e-6 #* u.rad / u.s
	sr_B = -0.484e-6 #* u.rad / u.s
	sr_C = -0.361e-6 #* u.rad / u.s

	omega_orb = (2*π) / (365.25 * 24. * 3600.) #* u.rad / u.s # omega of Earth around Sun

    diff_rot_fac = sr_A*(7*ld_coeff - 15.) / (20*(3-ld_coeff)) +
      sr_B*(19 *ld_coeff - 35 + cos(inc_sol_axis)^2 * (3*ld_coeff - 35)) / (280*(3-ld_coeff)) +
      sr_C*(187*ld_coeff - 315 - cos(inc_sol_axis)^2 * (630 - 118*ld_coeff) -
      105 *cos(inc_sol_axis)^4 * (ld_coeff-1.)) / (6720*(3-ld_coeff))

    delta_vr = (epsilon * sin(phi) * R_sun * sin(inc_sol_axis) * diff_rot_fac -
       (7*ld_coeff - 15) * epsilon * sin(psi) / (20*(3-ld_coeff)) * R_sun * omega_orb)
	delta_vr *= 1e3 # km/s -> m/s.  WARN: Andrea's python code returns cm/s

    return delta_vr
end

"""
   calc_Δv_diff_extinction(bjd; obs, ...)
   calc_Δv_diff_extinction(start, stop; obs, dt, dt_unit, ...)
Estimate apparent Δv (m/s) in solar observations due to differential extinction.
"""
function calc_Δv_diff_extinction end

function query_horizons(;epochs::Union{Vector{Float64},Dict}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	if typeof(epochs) <: Dict
		@assert haskey(epochs,"start")
		@assert haskey(epochs,"stop")
		@assert haskey(epochs,"step")
	elseif typeof(epochs) <: Vector{Float64}
		@assert 1 <= length(epochs) <= 10000
	end
	@assert 1 < max_airmass <= 10.0
	loc = get_obs_loc(obs)
	# Sometimes this fails.  Why?
	obj = JplHorizons.Horizons(id="10", location=loc, epochs=epochs, id_type="majorbody")
	eph = obj.ephemerides(quantities="1,4,8,17,20,42", airmass_lessthan=max_airmass, refraction=true)
	if verbose  println(eph)  end
	return eph
end

function calc_Δv_diff_extinction_horizons(;epochs::Union{Vector{Float64},Dict}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	loc = get_obs_loc(obs)
	eph = query_horizons(epochs=epochs,obs=obs,max_airmass=max_airmass, verbose=verbose)
	delta_vr = diff_extinct_corr.(eph.columns["datetime_jd"], eph.columns["EL"],
	 	eph.columns["hour_angle"], eph.columns["RA"],
		eph.columns["NPole_ang"], eph.columns["delta"].*AU, loc["lat"])
	return (jd=eph.columns["datetime_jd"], Δv=delta_vr)
end

function cache_horizons_output(year::Int64, obs::Symbol)
	fn = "horizons_" * string(obs) * "_" * string(year) * ".csv"
	filename = joinpath(dirname(pathof(NeidSolarScripts)),"..","data","horizons",fn)
	if isfile(filename)
		return
	end
	df = DataFrame(:jd=>Float64[], :sol_npole_ang=>Float64[], :sol_dist_au=>Float64[])
	samples_per_day = 4
	epoch_start = Dates.datetime2julian(DateTime(year,1,1))
	for i in 1:samples_per_day
		epoch_stop = epoch_start + ceil(Int64,366/samples_per_day)
		epochs = collect(range(epoch_start,stop=epoch_stop,step=1.0/samples_per_day) )
		eph = query_horizons(epochs=epochs, obs=obs, max_airmass=10.0)
		df_tmp = DataFrame(:jd=>eph.columns["datetime_jd"], :sol_npole_ang=>eph.columns["NPole_ang"], :sol_dist_au=>eph.columns["delta"])
		append!(df, df_tmp)
		epoch_start = epoch_stop
	end
	CSV.write(filename,df)
	return df
end

function read_horizons_output(jd::Union{Float64,Vector{Float64}}, obs::Symbol)
	@assert obs ∈ valid_obs
	years = unique(year.(Dates.julian2datetime.(jd)))
	@assert length(years) >= 1
	cache = Dict{Int64,Any}()
	for year in years
		fn = "horizons_" * string(obs) * "_" * string(year) * ".csv"
		filename = joinpath(dirname(pathof(NeidSolarScripts)),"..","data","horizons",fn)
		if isfile(filename)
			cache[year] = CSV.read(filename, DataFrame)
		else
			cache[year] = cache_horizons_output(year, obs)
		end
		cache[year].sol_npole_ang[cache[year].sol_npole_ang .> 180] .-= 360
	end
	if length(years) == 1
		return cache[first(keys(cache))]
	else
		sorted_keys = sort(collect(keys(cache)))
		output = cache[first(sorted_keys)]
		for k in sorted_keys[2:end]
			append!(output,cache[k])
		end
		return output
	end
end

function get_horizons_output(jd::Float64, obs::Symbol)
  df = read_horizons_output(jd, obs)
  idxlo = searchsortedlast(df.jd, jd)
  @assert 1 <= idxlo < size(df,1)
  idxhi = idxlo + 1
  result = Dict{Symbol,Float64}()
  num = jd-df.jd[idxlo]
  denom = df.jd[idxhi]-df.jd[idxlo]
  whi = num/denom
  wlo = 1-whi
  for k in Symbol.(names(df))
	if k == :jd continue end
  	result[k] = df[idxlo,k] * wlo +  df[idxhi,k] * whi
  end
  result[:jd] = jd
  return result
end

function get_horizons_output(jd::Vector{Float64}, obs::Symbol)
  df = read_horizons_output(jd, obs)
  idxlo = map(t->searchsortedlast(df.jd, t),jd)
  @assert all(1 .<= idxlo .< size(df,1))
  result = DataFrame(:jd=>jd)
  whi = (jd.-df.jd[idxlo]) ./ (df.jd[idxlo.+1].-df.jd[idxlo])
  #wlo = 1.-whi
  for k in  Symbol.(names(df))
	if k == :jd continue end
  	result[!,k] = df[idxlo,k] .* (1.0.-whi) .+  df[idxlo.+1,k] .* whi
  end
  return result
end

function calc_Δv_diff_extinction_astropy(;epochs::Vector{Float64}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	#=if typeof(epochs) <: Dict
		@assert haskey(epochs,"start")
		@assert haskey(epochs,"stop")
		@assert haskey(epochs,"step")
	else=#
	if typeof(epochs) <: Vector{Float64}
		@assert 1 <= length(epochs) <= 10000
	end
	@assert 1 < max_airmass <= 10.0
	sol = SolarRotation.get_solar_info(epochs, obs=obs)
	ha = sol.lst .- sol.ra./15
	idx_visible = 0 .< sol.airmass .<= max_airmass
	horizons = get_horizons_output(epochs[idx_visible],obs)
	#delta = SolarRotation.get_earth_sun_dist.(epochs[idx_visible], obs=obs)
	delta = sun.sol_dist_au
	loc = get_obs_loc(obs)
	#println(" jd = ", epochs, " delta_h = ", horizons.sol_dist_au, " delta_ap = ", delta)
	#println(" diff = ", horizons.sol_dist_au.-delta)
	delta_vr = diff_extinct_corr.(epochs[idx_visible], sol.alt[idx_visible], ha[idx_visible], sol.ra[idx_visible], horizons.sol_npole_ang, horizons.sol_dist_au.*AU, loc["lat"]) #eph["NPole_ang"], eph["delta"])
	return delta_vr
end

function calc_Δv_diff_extinction_astropy(epoch::Float64; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	@assert 1 < max_airmass <= 10.0
	sol = SolarRotation.get_solar_info(epoch, obs=obs)
	ha = sol[:lst] .- sol[:ra]./15
	idx_visible = 0 < sol[:airmass] <= max_airmass
	if !idx_visible return (jd=epoch, Δv = 0.0) end
	horizons = get_horizons_output(epoch,obs)
	#delta = SolarRotation.get_earth_sun_dist(epoch, obs=obs)
	loc = get_obs_loc(obs)
	delta_vr = diff_extinct_corr(epoch, sol[:alt], ha, sol[:ra], horizons[:sol_npole_ang], sol[:sol_dist_au].*AU, loc["lat"]) #eph["NPole_ang"], eph["delta"])
	return delta_vr
end

function calc_Δv_diff_extinction(epoch::Float64; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	calc_Δv_diff_extinction_astropy(epoch; obs=obs, max_airmass=max_airmass, verbose=verbose)
end

function compare_Horizons_Astropy(;epochs::Vector{Float64}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	eph = query_horizons(epochs=epochs,obs=obs,max_airmass=max_airmass, verbose=verbose)
	delta_vr = diff_extinct_corr.(eph.columns["datetime_jd"], eph.columns["EL"],
	 	eph.columns["hour_angle"], eph.columns["RA"],
		eph.columns["NPole_ang"], eph.columns["delta"].*AU, loc["lat"])

	sol = SolarRotation.get_solar_info(epochs, obs=obs)
	ha = sol.lst .- sol.ra./15
	idx_visible = 0 .< sol.airmass .<= max_airmass
	delta = SolarRotation.get_earth_sun_dist.(epochs[idx_visible], obs=obs)

	println("Horizons airmass =", eph.columns["airmass"])
	println("Astropy  airmass =", sol.airmass[idx_visible])
	println("Δjd = ", eph.columns["datetime_jd"], " vs ", epochs[idx_visible], " = ", eph.columns["datetime_jd"].-epochs[idx_visible])
	println("Δalt = ", eph.columns["EL"].- sol.alt[idx_visible])
	println("Δha = ", eph.columns["hour_angle"].- ha[idx_visible])
	println("ΔRA = ", eph.columns["RA"].- sol.ra[idx_visible])

	horizons = get_horizons_output(epochs[idx_visible],obs)
	delta_vr_new = diff_extinct_corr.(epochs[idx_visible], sol.alt[idx_visible], ha[idx_visible], sol.ra[idx_visible], zeros(length(epochs[idx_visible])), fill(AU,length(epochs[idx_visible])), loc["lat"]) #eph["NPole_ang"], eph["delta"])

	println("delta_vr_old = ", delta_vr, "\ndelta_vr_new = ", delta_vr_new)
	println("delta = ", delta_vr.- delta_vr_new)
	return (eph=eph, sol=sol, delta_vr_old=delta_vr, delta_vr_new=delta_vr_nw)
end

function calc_Δv_diff_extinction(bjd::Real; obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	epochs = [bjd]
	(bjd_out, Δv) = calc_Δv_diff_extinction(epochs=epochs, obs=obs, max_airmass=max_airmass, verbose=verbose)
	return first(Δv)
end

function calc_Δv_diff_extinction(start_date::String, stop_date::String; obs::Symbol, dt::Integer = 1, dt_unit::String = "m", max_airmass::Float64 = 10.0, verbose::Bool = false)
	@assert 1 <= dt <= 24*60
	@assert dt_unit ∈ ["m", "h", "d"]
	@assert 1 < max_airmass <= 10.0
	epochs = Dict("start"=>start_date, "stop"=>stop_date, "step"=>string(dt) * dt_unit)
	(bjd_out, Δv) = calc_Δv_diff_extinction(epochs=epochs, obs=obs, max_airmass=max_airmass, verbose=verbose)
	return (bjd=bjd_out, Δv=Δv)
end

function calc_Δv_diff_extinction(start_bjd::Real, stop_bjd::Real; obs::Symbol, dt::Integer = 1, dt_unit::String = "m", max_airmass::Float64 = 10.0, verbose::Bool=false)
	start_date = string(Date(Dates.julian2datetime(start_bjd)))
	stop_date = string(Date(Dates.julian2datetime(stop_bjd)))
	calc_Δv_diff_extinction(start_date, stop_date, obs=obs, dt=dt, dt_unit=dt_unit, max_airmass=max_airmass, verbose=verbose)
end

function calc_Δv_diff_extinction_demo(; verbose::Bool = false)
	start_date = "2021-01-01"
	stop_date = "2021-02-11"
	calc_Δv_diff_extinction(start_date, stop_date, obs=:WIYN, verbose=verbose)
end

export calc_Δv_diff_extinction

end # module DifferentialExtinction


text:
### routines for fitting a composite sigmoid to piece-wise linear functions.

function evalcompositelogisticprobitcost(t_range, f::Function, p::Vector{T})::T where T

    s = zero(T)
    for (i,t) in enumerate(t_range)
        s += ( f(t) - evalcompositelogisticprobit(t, p[1], p[2]) )^2
    end

    return s
end

"""
    getcompactsigmoidparameters(infos::Vector{MonotoneMaps.Piecewise2DLineType{T}};
        N_fit_positions::Int = 15,
        max_iters = 5000,
        xtol_rel = 1e-5,
        ftol_rel = 1e-5,
        maxtime = Inf,
        p0 = [0.5; 0.0],
        p_lb = [0.1; -5.0],
        p_ub = [0.6; 5.0]

Given a set of single region piece-wise linear function parameter objects, infos, fit the `a` and `b` parameters of the compact sigmoid (see function `evalcompositelogisticprobit`) for each piece-wise linear function.

∀ i ∈ [length(infos)], infos[i].xs must be 4 elements, with the first and last being the boundary, and the second and third being the end points of the single region.

p0, p_lb, p_ub are two element 1-D arrays.
...
# Select arguments
Two-element 1-D array for the following. The first entry relates to the a parameter, and the second relates to the b parameter.
- `p0::Vector{T}`: initial guess to the optimization.
- `p_lb::Vector{T}`: lower bounds to the optimization.
- `p_ub::Vector{T}`: upper bounds to the optimization.
- `optim_algorithm::Symbol` can be :GN_ESCH, :GN_ISRES, :LN_BOBYQA, :GN_DIRECT_L

...

See optim.jl in the examples folder for usage details.
"""
function getcompactsigmoidparameters(infos::Vector{MonotoneMaps.Piecewise2DLineType{T}};
    N_fit_positions::Int = 15,
    max_iters = 5000,
    xtol_rel = 1e-5,
    ftol_rel = 1e-5,
    maxtime = Inf,
    optim_algorithm::Symbol = :GN_ESCH,
    p0 = [0.5; 0.0],
    p_lb = [0.1; -5.0],
    p_ub = [0.6; 5.0],
    evalcostfunc = evalcompositelogisticprobitcost) where T

    fs = collect( xx->evalpiecewise2Dlinearfunc(xx, infos[i]) for i = 1:length(infos) )

    L = length(fs)
    gs = Vector{Function}(undef, L)
    p_star_set = Vector{Vector{T}}(undef, L)
    rets = Vector{Symbol}(undef, L)

    for l = 1:L

        t_range = LinRange(infos[l].xs[2], infos[l].xs[3], N_fit_positions)

        g = pp->evalcostfunc(t_range, fs[l], pp)
        dg = xx->Zygote.gradient(f, xx)

        opt = NLopt.Opt(optim_algorithm, length(p0))

        minf, minx, ret, N_evals = runNLopt!(opt,
            p0,
            g,
            dg,
            p_lb,
            p_ub;
            max_iters = max_iters,
            xtol_rel = xtol_rel,
            ftol_rel = ftol_rel,
            maxtime = maxtime)
        # update.
        p_star_set[l] = minx
        gs[l] = g
        rets[l] = ret
    end

    return gs, p_star_set, rets
end