text:
push!(LOAD_PATH, pwd())

import PyPlot
import HDF5
import .InterpolationFunctions
import .BaselineFunctions
import .PeakshapeFunctions

function baselineAndPeakshape(
  filepath;
  outputfilename="results/_result.hdf5",
  peakshapeRegions=8,
  peakshapeRegionStretch=0.5,
  peakshapeQuantileValue = 0.1,
  peakfindingNoiseThresholdValue = 25,
  peakWindowWidth = 200,
  peakfindingSignalLimit = 0.2
  )

  file = joinpath(filepath, outputfilename)
  massAxis = HDF5.h5read(file, "MassAxis")
  avgSpectrum = HDF5.h5read(file, "AvgSpectrum")

  baselinePoints, baselineValues, baselineNoise = BaselineFunctions.calculateBaseline(massAxis,avgSpectrum,baselinePointWidth = 0.8, threshold=0.2)
  baselineNoiseInterpolated = InterpolationFunctions.interpolate(massAxis, baselinePoints, baselineNoise)
  baselineInterpolated = InterpolationFunctions.interpolate(massAxis, baselinePoints, baselineValues)
  baselineCorrectedAvgSpec = avgSpectrum[:,1] - baselineInterpolated;

  peakIndices = PeakshapeFunctions.findPeakIndices(massAxis, avgSpectrum, baselineInterpolated, baselineNoiseInterpolated, noiseThreshold = peakfindingNoiseThresholdValue, signalLimit = peakfindingSignalLimit)
  peakShapesCenterMass, peakShapesY = PeakshapeFunctions.calculatePeakshapes(massAxis, baselineCorrectedAvgSpec, peakIndices, nbrMassRegions = peakshapeRegions,regionStretch=peakshapeRegionStretch, peakWindowWidth = peakWindowWidth, quantileValue = peakshapeQuantileValue)

  PyPlot.figure()
  PyPlot.title("Baseline Correction")
  PyPlot.semilogy(baselinePoints,baselineValues,".-")
   PyPlot.semilogy(baselinePoints,baselineValues + baselineNoise,".")
  #semilogy(peakMasses, peakValues, "x")
   PyPlot.semilogy(massAxis,avgSpectrum)

  ############ delete h5 data that will be overwritten ###########
  fh = HDF5.h5open(file,"r+")
  if Base.haskey(fh, "AvgBaseline")
    HDF5.o_delete(fh,"AvgBaseline")
  end
  if Base.haskey(fh, "MassDepPeakshape")
  HDF5.o_delete(fh,"MassDepPeakshape")
  end
  if Base.haskey(fh, "MassDepPeakshapeCenterMasses")
  HDF5.o_delete(fh,"MassDepPeakshapeCenterMasses")
  end
  if Base.haskey(fh, "BaseLines")
  HDF5.o_delete(fh,"BaseLines")
  end
  HDF5.close(fh)

  HDF5.h5write(file, "AvgBaseline", baselineInterpolated)
  HDF5.h5write(file, "MassDepPeakshape", peakShapesY)
  HDF5.h5write(file, "MassDepPeakshapeCenterMasses", peakShapesCenterMass)

  fh = HDF5.h5open(file,"r")
  ds = fh["SumSpecs"]
  spectra = size(ds)[2]
  HDF5.close(fh)
  if (spectra > 0)
    SumSpecs = HDF5.h5read(file, "SumSpecs")
    fill!(SumSpecs,0)
    HDF5.h5write(file, "BaseLines", SumSpecs)

  else
    HDF5.h5write(file, "SumSpecs", avgSpectrum)
    fill!(avgSpectrum,0)
    HDF5.h5write(file, "BaseLines", avgSpectrum)
  end

  println("DONE")
end


text:
#-------------------------------------------------------------------------------
# File: sat_encoding.jl
# Description: This file contains all the fuctions needed to encode
# an SGP instance into an SAT model.
# Date: December 10, 2019
# Author: Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
# TBH :  This was mostly inspired by "DOI 10.1007/s10479-015-1914-5"
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# Genenarl Idea of the model :
# One boolean variable denoted g[q', p', g', w'] states if :
# player q' is the p'-th player of the group number g' of week w' or not

#-------------------------------------------------------------------------------
const g, p, w = 2,2,3

"""
    sat_modeling(p::Int, g::Int, w::Int)

Generate a SAT model (in CNF) of the SGP using the given parameter :
    _ p : number of players per group
    _ g : number of groups per week
    _ w : number of weeks considered

Output :
    A file containing the model in order to use a sat solver :
        _ Each line is a clause and is ended by a "0"
        _ Each boolean variable is denoted by a Integer in [ 1 , (q*p*g*w) ] (with q the total number of players)
        _ positives values "i" states for "i"
        _ negatives values "-i" states for "NOT(i)"

Example : (x1 v NOT(x5) v x4 v x1) AND (NOT(x1) v x5 v x3 v x4)
    1 -5 4 0
    -1 5 3 4 0
"""
function sat_modeling()
    q = g * p # q denote the total number of players
    println(" Encoding Model ", w ," weeks ", g, " groups ", p, " players per groups ")
    println(" Total number of player : ", q)
    println(" Total number of variable : ", q*p*g*w)
    nb_clause = 0
    # Writting in the file
    open("sgp.cnf","w") do f

        # Each golfer plays once a week
        for q_var in 1:q, w_var in 1:w
            clause = ""
            for p_var in 1:p, g_var in 1:g
                clause = clause*string( var_to_int(q_var, p_var, g_var, w_var)," ")
            end
            clause = clause*"0\n"
            write(f,clause)
            nb_clause += 1
        end

        for q_var in 1:q, w_var in 1:w, p_var in 1:p, g_var in 1:g
            for p_second in p_var+1:p
                clause = string( not(var_to_int(q_var, p_var, g_var, w_var))," ",not(var_to_int(q_var, p_second, g_var, w_var))," 0\n")
                write(f,clause)
                nb_clause += 1
            end
        end

        for q_var in 1:q, w_var in 1:w, p_var in 1:p, g_var in 1:g
            for g_second in g_var+1:g, p_second in p_var+1:p
                clause = string(not(var_to_int(q_var, p_var, g_var, w_var))," ",not(var_to_int(q_var, p_second, g_second, w_var))," 0\n")
                write(f,clause)
                nb_clause += 1
            end
        end

        # Groups are correct
        for w_var in 1:w, p_var in 1:p, g_var in 1:g
            clause = ""
            for player in 1:q
                clause = clause*string(var_to_int(player, p_var, g_var, w_var)," ")
            end
            clause = clause*"0\n"
            write(f,clause)
            nb_clause += 1
        end

        for w_var in 1:w, p_var in 1:p, g_var in 1:g, player in 1:q
            for player_second in player+1:q
                clause = string(not(var_to_int(player, p_var, g_var, w_var))," ",not(var_to_int(player_second, p_var, g_var, w_var))," 0\n")
                write(f,clause)
                nb_clause += 1
            end
        end

        # Socialization constraint
        for w_var in 1:w, g_var in 1:g
            for w_second in w_var+1:w, g_second in 1:g, player in 1:q, p1_var in 1:p, p1_second in 1:p
                for player_second in player+1:q, p2_var in 1:p, p2_second in 1:p
                    clause = string(not(var_to_int(player, p1_var, g_var, w_var)), " ",not(var_to_int(player_second, p2_var, g_var, w_var)), " ",not(var_to_int(player, p1_second, g_second, w_second)), " ",not(var_to_int(player_second, p2_second, g_second, w_second)), " 0\n")
                    write(f,clause)
                    nb_clause += 1
                end
            end
        end
    end
end

function not(variable::Int)
    return -variable
end

function var_to_int(q_var::Int, p_var::Int, g_var::Int, w_var::Int)
    q = g * p
    return q_var + (p_var-1)*(q) +  (g_var-1)*(q)*(p) + (w_var-1)*(q)*(p)*(g)
end

sat_modeling()

ext: jl

max_stars_repo_name: remi-garcia/ProjetCSP

max_issues_repo_name: remi-garcia/ProjetCSP

max_forks_repo_name: remi-garcia/ProjetCSP


text:
#!jl # # EFTfitter.jl - Empty Template
#jl # EFTfitter.jl - Empty Template

#nb using EFTfitter
#nb using BAT            # for sampling
#nb using IntervalSets   # for specifying the prior
#nb using Distributions  # for specifying the prior
#nb using Plots          # for plotting

#!jl include("lit_inputs.jl")

#!nb using EFTfitter
#!nb using BAT            # for sampling
#!nb using IntervalSets   # for specifying the prior
#!nb using Distributions  # for specifying the prior
#!nb using Plots # for plotting
#jl
#jl #~include definitions of observables, measurements,
#jl #~uncertainties and correlations from the inputs file:
#jl include("inputs.jl")

#~create an `EFTfitterModel` object:
model = EFTfitterModel(parameters, measurements, correlations)
#~create posterior distribution:
posterior = PosteriorDensity(model);

#~sample the posterior distribution with BAT.jl:
algorithm = MCMCSampling(mcalg = MetropolisHastings(), nsteps = 10^5, nchains = 4)
samples = bat_sample(posterior, algorithm).result;

#~create and display a `SampledDensity` object for a quick overview of results:
sd = SampledDensity(posterior, samples)
display(sd)

#~ plot the posterior distribution:
p = plot(samples)
#!nb savefig(p, "plot.pdf")


text:
function _build_complement_constraint(
    _error::Function,
    constraint,
)
    set = JuMP.moi_set(constraint)

    jump_fct = JuMP.jump_function(constraint)
    moi_fct = JuMP.moi_function(constraint)

    return JuMP.VectorConstraint(
        [jump_fct...], ComplementSet{typeof(moi_fct)}(set)
    )
end

function JuMP.parse_one_operator_constraint(_error::Function, vectorized::Bool, ::Val{:!}, constraint)
    _error1 = deepcopy(_error)
    if vectorized
        _error("`$(constraint)` should be non vectorized. There is currently no vectorized support for `!` constraints. Please open an issue at ConstraintSolver.jl")
    end

    vectorized, inner_parsecode, inner_buildcall = JuMP.parse_constraint_expr(_error1, constraint)

    buildcall = :($(esc(:(CS._build_complement_constraint)))(
        $_error,
        $inner_buildcall,
    ))

    return inner_parsecode,buildcall
end


text:
function high_degree_adaptive(g)
	attack_nodes = Int[]
	deg = degree(g)
	degmax = maximum(deg)
	H = [Set{Int}() for i=1:degmax]

	for i in vertices(g)
		deg[i] > 0 && push!(H[deg[i]], i)
	end

	while degmax > 0
		if !isempty(H[degmax])
			i = rand(H[degmax])
			delete!(H[degmax], i)
			deg[i] = 0
			push!(attack_nodes, i)
			for j in neighbors(g,i)
				if deg[j] > 0
					delete!(H[deg[j]], j)
					deg[j] -= 1
					deg[j] > 0 && push!(H[deg[j]], j)
				end
			end
		else
			degmax -= 1
		end
	end
	return attack_nodes
end


text:
import Pkg

"""
The version of the `OMOPCommonDataModel.jl` Julia package.
"""
function version()::VersionNumber
    package_directory = dirname(dirname(@__FILE__))
    project_file = joinpath(package_directory, "Project.toml")
    version_string = Pkg.TOML.parsefile(project_file)["version"]
    version_number = VersionNumber(version_string)
    return version_number
end

"""
The version of the OMOP Common Data Model (CDM) being implemented.
"""
@inline function cdm_version()::VersionNumber
    return CDM_VERSION
end


text:
@testset "1009.complement-of-base-10-integer.jl" begin
    @test bitwise_complement(7) == 0
    @test bitwise_complement(5) == 2
    @test bitwise_complement(10) == 5
end


text:
"""
Calculates the SVD deomposition of a matrix by using the
Jacobi method. Returns a `SVD` factorization object.
"""
gesvj(x) = gesvj!(copy(x))
export gesvj
"""
Same as `gesvj` but saves space by overwriting the input matrix.
"""
function gesvj!(x)
  JacobiSVD.jsvd!(x)
end
export gesvj!


text:
using PkgName
using Base.Test

# write your own tests here
@test "Test for A*V[:,1:m] is almost equal to V[:,1:m+1]*H" begin include("approximate.jl") end
@test "Test for orthogonality" begin include("ortho.jl") end


text:
function triangle(t, a = 1)
    t = mod(abs(2 * t/(2*pi)), 2)
    abs(2 * a - 2 * t) - a
end

function squarewave(t, a = 1)
    t = 2*t
    mt = mod(t/(2*pi),2)

    if 0 <= mt <= 1
        return -a
    elseif 1 < mt <= 2
        return a
    end

end


text:
"""
    CoulombInteractionMultipole(k, g)

Represents the `k`th multipole of the multipole expansion of the
Coulomb interaction `g`.
"""
struct CoulombInteractionMultipole{G<:CoulombInteraction} <: TwoBodyOperator
    k::Int
    g::G
end

"""
    CoulombPotentialMultipole

Type alias for contraction of the
[`CoulombInteractionMultipole`](@ref) over one coordinate, thereby
forming a potential in the other coordinate.
"""
const CoulombPotentialMultipole{A,B} = ContractedOperator{1,2,1,A,<:CoulombInteractionMultipole,B}

Base.show(io::IO, ci::CoulombInteractionMultipole) = write(io, "ĝ", to_superscript(ci.k))

function Base.show(io::IO, me::OrbitalMatrixElement{2,A,<:CoulombInteractionMultipole,B}) where {A,B}
    if me.a == me.b # Direct interaction
        write(io, "F",to_superscript(me.o.k),"($(me.a[1]),$(me.a[2]))")
    elseif me.a[1] == me.b[2] && me.a[2] == me.b[1] # Exchange interaction
        write(io, "G",to_superscript(me.o.k),"($(me.a[1]),$(me.a[2]))")
    else # General case
        write(io, "R",to_superscript(me.o.k),"(", join(string.(me.a), ","), ";", join(string.(me.b), ","), ")")
    end
end

Base.show(io::IO, me::CoulombPotentialMultipole{A,B}) where {A,B}=
write(io, "r⁻¹×Y",to_superscript(me.o.k),"($(me.a[1]),$(me.b[1]))")

radial_integral(a, (k,g)::Tuple{<:Integer,<:CoulombInteraction}, b) =
    OrbitalMatrixElement(a, CoulombInteractionMultipole(k,g), b)

"""
    CoulombTensor(k)

Construct a Coulomb interaction tensor of rank `k`.
"""
struct CoulombTensor{k} <: Tensor{k,'K'} end

"""
    system(::Type{CoulombTensor})

A Coulomb tensor only acts on the coordinates ``r``, ``\\theta`` and
``\\phi``.
"""
system(::Type{<:CoulombTensor}) = SpatialSubSystem()

@doc raw"""
    RadialCoulombMatrixElement

This represents the matrix element of the radial component of the
Coulomb tensor operator:

```math
\tensor{K}^{(k)}(i) \defd
\left\{[1-\Heaviside(r_j-r_i)]r_i^k +
\frac{\Heaviside(r_j-r_i)}{r_i^{k+1}}\right\}
\tensor{C}^{(k)}(i),
\quad
i = 1,2,
\quad
j = 3-i,
```
"""
struct RadialCoulombMatrixElement{k} <: OneBodyOperator end

RadialCoulombMatrixElement(k) = RadialCoulombMatrixElement{k}()

function Base.:(*)(a::RadialCoulombMatrixElement{k₁}, b::RadialCoulombMatrixElement{k₂}) where {k₁,k₂}
    @assert k₁ == k₂
    CoulombInteractionMultipole(k₁)
end

@tensor(CoulombTensor{k} where k) do
    begin
        n′ ~ n # The Coulomb interaction couples orbitals of different
               # n, but there is no selection rule.
        ℓ′ ∈ abs(ℓ - k):2:(ℓ+k)
    end

    raw"""
    rme((n′,ℓ′), ::CoulombTensor{k}, (n,ℓ))

Computes the reduced matrix element of `𝐊` in terms of
[`RadialCoulombMatrixElement`](@ref).
"""
    rme(ℓ′, SphericalTensor(k), ℓ) # *RadialCoulombMatrixElement(k)
end

"""
    ranks(a, ::Type{CoulombTensor}, b)

Return which tensor ranks for Coulomb tensors that fulfill the
triangle condition between spin-orbitals `a` and `b`.
"""
ranks(a::SpinOrbital, ::Type{CoulombTensor}, b::SpinOrbital) =
    triangle_range(a.orb.ℓ, b.orb.ℓ)

export CoulombInteractionMultipole, CoulombPotentialMultipole, CoulombTensor


text:
"""
    make_optimization_model(
        model::MetabolicModel,
        optimizer;
        sense = MOI.MAX_SENSE,
    )

Convert `MetabolicModel`s to a JuMP model, place objectives and the equality
constraint.
"""
function make_optimization_model(model::MetabolicModel, optimizer; sense = MOI.MAX_SENSE)

    precache!(model)

    m, n = size(stoichiometry(model))
    xl, xu = bounds(model)

    optimization_model = COBREXA.JuMP.Model(optimizer)
    @variable(optimization_model, x[i = 1:n])
    @objective(optimization_model, sense, objective(model)' * x)
    @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
    @constraint(optimization_model, lbs, xl .<= x) # lower bounds
    @constraint(optimization_model, ubs, x .<= xu) # upper bounds

    C = coupling(model) # empty if no coupling
    cl, cu = coupling_bounds(model)
    isempty(C) || @constraint(optimization_model, c_lbs, cl .<= coupling(model) * x) # coupling lower bounds
    isempty(C) || @constraint(optimization_model, c_ubs, coupling(model) * x .<= cu) # coupling upper bounds

    return optimization_model
end

"""
    optimize_model(
        model::MetabolicModel,
        optimizer;
        sense = MOI.MIN_SENSE,
    )

Use JuMP to solve an instance of CoreModel
"""
function optimize_model(model::MetabolicModel, optimizer; sense = MOI.MIN_SENSE)
    optimization_model = make_optimization_model(model, optimizer; sense = sense)
    COBREXA.JuMP.optimize!(optimization_model)
    return optimization_model
end

"""
    is_solved(optmodel)

Return `true` if `optmodel` solved successfully (solution is optimal or locally
optimal).  Return `false` if any other termination status is reached.
Termination status is defined in the documentation of `JuMP`.
"""
function is_solved(optmodel)
    COBREXA.JuMP.termination_status(optmodel) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ? true :
    false
end

"""
    get_bound_vectors(opt_model)

Returns vectors of the lower and upper bounds of `opt_model` constraints, where
`opt_model` is a JuMP model constructed by e.g.
[`make_optimization_model`](@ref) or [`flux_balance_analysis`](@ref).
"""
get_bound_vectors(opt_model) = (
    [-normalized_rhs(lb) for lb in opt_model[:lbs]],
    [normalized_rhs(ub) for ub in opt_model[:ubs]],
)

"""
    set_bound(index, optimization_model;
        ub=_constants.default_reaction_rate,
        lb=-_constants.default_reaction_rate)

Helper function to set the bounds of variables.
The JuMP `set_normalized_rhs` function is a little confusing,
so this function simplifies setting constraints. In short, JuMP
uses a normalized right hand side representation of constraints,
which means that lower bounds have their sign flipped. This function
does this for you, so you don't have to remember to do this whenever you
change the constraints.

Just supply the constraint `index` and the JuMP model (`opt_model`) that
will be solved, and the variable's bounds will be set to `ub` and `lb`.
"""
function set_bound(
    vind,
    opt_model;
    ub = _constants.default_reaction_rate,
    lb = -_constants.default_reaction_rate,
)
    set_normalized_rhs(opt_model[:lbs][vind], -lb)
    set_normalized_rhs(opt_model[:ubs][vind], ub)
end


text:
using Test

if !@isdefined TESTING
    TESTING = false
end

macro testing(ex)
    @assert ex.head === :call
    @assert length(ex.args) == 2
    quote
        ans = $(esc(ex.args[2]))
        @test $(esc(ex.args[1])) == ans
        ans
    end
end


text:
function myloss(
    par,
    model,
    mult
    )
    #imp_det = det(model.imply.imp_cov)
    if !isposdef(model.imply.imp_cov)
        F = Inf
    else
        a = cholesky!(model.imply.imp_cov)
        ld = logdet(a)
        model.imply.imp_cov .= LinearAlgebra.inv!(a)
        mul!(mult, model.observed.obs_cov, model.imply.imp_cov)
        #mul!()
        F = ld +
            tr(mult)# -
            #semml.logdet_obs -
            #model.observed.n_man
    end
    return F
end

mult = zeros(11,11)
inverse = zeros(11,11)

model.imply(solution.minimizer)

mat = copy(model.imply.imp_cov)

model.imply.imp_cov .= copy(mat)

@btime myloss(solution.minimizer, model, mult, inverse)

a = 11 + model.loss.functions[1].logdet_obs

mat2 = copy(mat)

mat3 = cholesky!(mat2)

function myfun(mat3, res)
    res .= LinearAlgebra.inv!(mat3)
end

res = zeros(11,11)

@time myfun(mat3, res)

c = rand(11,11); d = rand(11,11)

mul!(c, d, mat3)

mat3


text:
"""
    eigenvaluesInGrid(
        m::AIMSteppingMethod,
        p::NumericAIMProblem{N,T},
        c::AIMCache{N,T},
        grid::Tuple{T,T,Int64,Int64};
        xtol::Real = 1.0e-10,
        ftol::Real = 1.0e-10,
        iterations::Int = 1000
        ) where {N <: Unsigned, T <: Complex}

Attempts to find eigenvalues using a grid of complex plane data points as initial guesses passed to nlsolve.

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::NumericAIMProblem{N,T}`: The previously defined problem data.
- `c::AIMCache{N,T}`: The cache constructed from p.
- `grid::Tuple{T,T,Int64,Int64}`: A tuple consisting of (start point, end point, num. of real pts., num. of imag. pots.).
- `xtol::Real`: Norm difference in x between two successive iterates under which convergence is declared.
- `ftol::Real`: Infinite norm of residuals under which convergence is declared.
- `iterations::Int`: Maximum number of iterations performed by NLsolve.

# Output
An object of type `Array{T,1}` containing the modes found within the grid.
"""
function eigenvaluesInGrid(
    m::AIMSteppingMethod,
    p::NumericAIMProblem{N,T},
    c::AIMCache{N,T},
    grid::Tuple{T,T,Int64,Int64};
    xtol::Real=1.0e-10,
    ftol::Real=1.0e-10,
    iterations::Int=1000
    ) where {N <: Unsigned,T <: Complex}

    if grid[3] < 0 || grid[4] < 0
        error("The number of points in the grid must be positive")
    end

    # Extract grid quantities
    re_start = real(grid[1])
    re_end = real(grid[2])
    re_size = grid[3]
    re_step = (re_end - re_start) / (re_size - 1)
    re_range = re_start:re_step:re_end

    im_start = imag(grid[1])
    im_end = imag(grid[2])
    im_size = grid[4]
    im_step = (im_end - im_start) / (im_size - 1)
    im_range = im_start:im_step:im_end

    # Check that the grid is  valid
    if (re_end <= re_start || re_step < 0 ) || (im_end <= im_start || im_step < 0 )
        error("The grid specified must be tuple in the form (start, end, real points, imag. points) where end > start.")
    end

    # Build the output data array with it's expected size
    eigenvalues::Array{T,1} = []

    for realPart in re_range
        for imagPart in im_range
            # Compute the solution, if any
            sol = computeEigenvalues(
                m,
                p,
                c,
                T(realPart, imagPart),
                nls_xtol=xtol,
                nls_ftol=ftol,
                nls_iterations=iterations
                )

            # If nlsolve converge the result is stored. If not take the appropriate action
            if sol.f_converged || sol.x_converged
                push!(eigenvalues, T(sol.zero[1], sol.zero[2]))
            end
        end
    end

    return eigenvalues
end

"""
    eigenvaluesInGrid(
        m::AIMSteppingMethod,
        p::NumericAIMProblem{N,T},
        c::AIMCache{N,T},
        grid::Tuple{T,T};
        roots_atol::Real = 1.0e-10,
        roots_rtol::Real = 1.0e-10,
        roots_xatol::Real = 1.0e-10,
        roots_xrtol::Real = 1.0e-10,
        roots_maxevals::Int = 100
        ) where {N <: Unsigned, T <: Real}

Attempts to find eigenvalues using a range of real data points as a search region to find_zeros.
For details on convergence settings see [Roots.jl](https://juliahub.com/docs/Roots/o0Xsi/1.0.7/reference/#Convergence).

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::NumericAIMProblem{N,T}`: The previously defined problem data.
- `c::AIMCache{N,T}`: The cache constructed from p.
- `grid::Tuple{T,T}`: A tuple consisting of (start point, end point).
- `roots_atol::Real`: Absolute tolerance.
- `roots_rtol::Real`: Relative tolerance.
- `roots_xatol::Real`: Floating point comparison absolute tolerance.
- `roots_xrtol::Real`: Floating point comparison relative tolerance.
- `roots_maxevals::Int`: Number of algorithm iterations performed.

# Output
An object of type `Array{T,1}` containing the eigenvalues found within the grid.
"""
function eigenvaluesInGrid(
    m::AIMSteppingMethod,
    p::NumericAIMProblem{N,T},
    c::AIMCache{N,T},
    grid::Tuple{T,T};
    roots_atol::Real=1.0e-10,
    roots_rtol::Real=1.0e-10,
    roots_xatol::Real=1.0e-10,
    roots_xrtol::Real=1.0e-10,
    roots_maxevals::Int=100,
    ) where {N <: Unsigned,T <: Real}

    if grid[1] > grid[2]
        error("The interval for searching eigenvalues must be ordered.")
    end

    try
        find_zeros(
            x -> computeDelta!(m, p, c, x),
            grid[1],
            grid[2],
            atol=roots_atol,
            rtol=roots_rtol,
            xatol=roots_xatol,
            xrtol=roots_xrtol,
            maxevals=roots_maxevals
            )
    catch
        println("find_zeros was unable to converge to an eigenvalue")
    end
end


text:
const ColIndex = Union{Int, AbstractString}
const KeyType = Union{Symbol, String}
const TransitionType = Union{String, NTuple{2,String}, NTuple{3,String}}
# const DbIndexType = Union{String, Integer}
# const DbIndexArrayType = Union{Array{String}, Array{Int}}


text:
# This definitely needs cleaning
using OrdinaryDiffEq, DiffEqDevTools, DiffEqBase, DiffEqProblemLibrary, Base.Test
probArr = Vector{ODETestProblem}(2)
probArr[1] = prob_ode_linear

probArr[2] = prob_ode_2Dlinear
srand(100)
## Convergence Testing
println("Convergence Test on Linear")
dts = 1.//2.^(8:-1:4)
testTol = 0.2

for i = 1:2
  prob = probArr[i]
  println("Special RKs")
  sim = test_convergence(dts,prob,Euler())
  @test abs(sim.𝒪est[:final]-1) < testTol
  sim2 = test_convergence(dts,prob,Midpoint())
  @test abs(sim2.𝒪est[:l∞]-2) < testTol
  sim3 = test_convergence(dts,prob,RK4())
  @test abs(sim3.𝒪est[:l∞]-4) < testTol
  sim4 = test_convergence(dts,prob,BS3())
  @test abs(sim4.𝒪est[:l2]-3) < testTol

  ### Stiff Solvers

  println("Convergence Test on Stiff")
  dts = 1.//2.^(8:-1:4)

  sim12 = test_convergence(dts,prob,ImplicitEuler(autodiff=true))
  @test abs(sim12.𝒪est[:final]-1) < testTol
  sim122 = test_convergence(dts,prob,ImplicitEuler(autodiff=false))
  @test abs(sim122.𝒪est[:final]-1) < testTol
  sim13 = test_convergence(dts,prob,Trapezoid(autodiff=true))
  @test abs(sim13.𝒪est[:final]-2) < testTol
  sim132 = test_convergence(dts,prob,Trapezoid(autodiff=false))
  @test abs(sim132.𝒪est[:final]-2) < testTol
end


text:
"""
Utility functions for seeding databases.
"""
module DatabaseSeeding

using SearchLight

export random_seeder

"""
    random_seeder(m::Module, quantity = 10, save = false)

Generic random database seeder. `m` must expose a `random()` function which returns a SearchLight instance.
If `save` the data will be persisted to the database, as configured for the current environment.
"""
function random_seeder(m::Module, quantity::Int = 10, save::Bool = true; skip_validation::Bool = true)
  Core.eval(parse("using $(m)"))

  seeds = []
  for i in 1:quantity
    item = m.random()
    push!(seeds, item)

    save && SearchLight.save!(item, skip_validation = skip_validation)
  end

  seeds
end

end


text:
# TODO Import Indicators.jl to get rolling average

module Forecasts

using Base.Iterators: take, cycle
using MultivariateStats: llsq
using ShiftedArrays: lag
using Missings: skipmissing

# package code goes here
export forecast, naivef, snaivef, meanf, rwf, fit!, update!, take, fit, ar

abstract type Forecast end

Base.eltype(::Type{<:Forecast}) = Float64
Base.iteratorsize(::Type{<:Forecast}) = Base.IsInfinite()
Base.done(::Forecast, state) = false

mutable struct NaiveForecast <: Forecast
    value::Float64
end

Base.start(fc::NaiveForecast) = nothing
Base.next(fc::NaiveForecast, state) = (fc.value, nothing)

function naivef(x)
    NaiveForecast(last(x))
end

function forecast(m::NaiveForecast, h::Int)
    repeat([m.value], inner=h)
end

function fit!(m::NaiveForecast, x::AbstractArray{Number})
    m.value = last(x)
end

function update!(m::NaiveForecast, x::Number)
    m.value = last(x)
end

### Seasonal Naive Forecast ###

struct SeasonalNaiveForecast <: Forecast
    # TODO Enforce that length(values) == s
    values::Array{Number,1}
    s::Int  # length of season
end

function snaivef(y, s)
    SeasonalNaiveForecast(y[(end-s+1):end], s)
end

Base.start(fc::SeasonalNaiveForecast) = start(cycle(fc.values))
# TODO Figure out if there's a better way to dispatch this `next`
Base.next(fc::SeasonalNaiveForecast, state) = next(cycle(fc.values), state)

### Mean Forecast ###

mutable struct MeanForecast <: Forecast
    sum::Number
    nobs::Number
end

function meanf(y)
    MeanForecast(sum(y), length(y))
end

Base.start(fc::MeanForecast) = (fc.sum / fc.nobs)
Base.next(fc::MeanForecast, state) = (state, state)

function update!(fc::MeanForecast, y)
    m.sum += sum(y)
    m.mobs += length(y)
end

### RwDriftForecast ###

mutable struct RwDriftForecast <: Forecast
    value::Number
    drift::Number
end

function rwf(y)
    RwDriftForecast(last(y), mean(diff(y)))
end

Base.start(fc::RwDriftForecast) = fc.value
Base.done(fc::RwDriftForecast, state) = false

function Base.next(fc::RwDriftForecast, prev)
    (prev + fc.drift, prev + fc.drift)
end

## AR Model ###
mutable struct ARForecast <: Forecast
    coefs::Array{Float64,2}
end

function ar(y, p=1)
    fit(ARForecast, y, order=p)
end

function fit(::Type{ARForecast}, y::Array{Float64,2}; order=1)
    trunc = x -> x[(order+1):endof(x),:]
    Ly = trunc(lag(y, order))
    Ly = oftype(y, Ly)
    coefs = llsq(Ly, trunc(y))
    ARForecast(coefs)
end

fit(T::Type{ARForecast}, y::Array{Float64,1}; order=1) = fit(T, reshape(y, length(y), 1), order=order)

function forecast(fc::Forecast, h::Int)
    collect(take(fc, h))
end

end # module


text:
#####
## Utilities for profiling
#####

# Chrome tracing reference:
# - https://www.gamedeveloper.com/programming/in-depth-using-chrome-tracing-to-view-your-inline-profiling-data
# - https://docs.google.com/document/d/1CvAClvFfyA5R-PhYUmn5OOQtYMH4h6I0nSsKchNAySU/preview

"""
Export some profiling info to be visualized with chrome://tracing.
"""
module ProfUtils

import JSON3
using LoggingExtras: FormatLogger, EarlyFilteredLogger, TeeLogger, global_logger

const LOG_GROUP = :chrome_tracing

function event(; name::String, cat::String, typ::String, pid::Int, tid::Int)
  @assert typ in ["B", "E", "i"]
  ts = Int(time_ns()) ÷ 1000
  msg = name
  typ == "B" && (msg *= " (started)")
  typ == "E" && (msg *= " (ended)")
  args = (; name, cat, ph=typ, ts, pid, tid)
  @debug msg _group=LOG_GROUP args...
end

event_start(; name::String, cat::String, pid::Int, tid::Int) =
  event(; name, cat, typ="B", pid, tid)

event_end(; name::String, cat::String, pid::Int, tid::Int) =
  event(; name, cat, typ="E", pid, tid)

instant_event(; name::String, cat::String, pid::Int, tid::Int) =
  event(; name, cat, typ="i", pid, tid)

function log_event(f; name::String, cat::String, pid::Int, tid::Int)
  event_start(; name, cat, pid, tid)
  res = f()
  event_end(; name, cat, pid, tid)
  return res
end

function chrome_tracing_logger(file::String; always_flush=false)
  stream = open(file, "w")
  println(stream, "[")
  formatter = FormatLogger(stream; always_flush) do stream, log
    args = values(log.kwargs)
    dict = (; args.name, args.cat, args.ph, args.ts, args.pid, args.tid)
    JSON3.write(stream, dict)
    println(stream, ",")
  end
  return EarlyFilteredLogger(formatter) do log
    return log.group == LOG_GROUP
  end
end

function set_chrome_tracing(file::String; always_flush=false)
  tracing_logger = chrome_tracing_logger(file; always_flush)
  logger = TeeLogger(global_logger(), tracing_logger)
  global_logger(logger)
  return
end

end


text:
module GPXFile

import ..GPX
import DataFrames: DataFrame
import Geodesy: LLA
import TimeZones: astimezone, localzone
import PooledArrays: PooledArray

function read(filename; tz=localzone())
  if isdir(filename)
    filenames = String[]
    for f ∈ readdir(filename; join=true)
      endswith(lowercase(f), ".gpx") && push!(filenames, f)
    end
    return read(filenames)
  end
  gpx = GPX.read_gpx_file(filename)
  length(gpx.tracks) > 1 && throw(ErrorException("More than one tracks in GPX file"))
  length(gpx.tracks[1].segments) > 1 && throw(ErrorException("More than one segments in GPX file"))
  p = gpx.tracks[1].segments[1].points
  df = DataFrame(time=[astimezone(p1.time, tz) for p1 ∈ p], location=[LLA(p1.lat, p1.lon) for p1 ∈ p])
  df.filename = repeat(PooledArray([filename]), size(df, 1))
  df
end

function read(filenames::AbstractVector)
  df = vcat(read.(filenames)...)
  df.filename = PooledArray(df.filename)
  sort!(df, :time)
end

end # module


text:
@doc_str """
    edit_distance(G₁::AbstractGraph, G₂::AbstractGraph)

Compute the edit distance between graphs `G₁` and `G₂`. Return the minimum
edit cost and edit path to transform graph `G₁` into graph `G₂``.
An edit path consists of a sequence of pairs of vertices ``(u,v) ∈ [0,|G₁|] × [0,|G₂|]``
representing vertex operations:

- ``(0,v)``: insertion of vertex ``v ∈ G₂``
- ``(u,0)``: deletion of vertex ``u ∈ G₁``
- ``(u>0,v>0)``: substitution of vertex ``u ∈ G₁`` by vertex ``v ∈ G₂``

### Optional Arguments
- `insert_cost::Function=v->1.0`
- `delete_cost::Function=u->1.0`
- `subst_cost::Function=(u,v)->0.5`

By default, the algorithm uses constant operation costs. The
user can provide classical Minkowski costs computed from vertex
labels μ₁ (for G₁) and μ₂ (for G₂) in order to further guide the
search, for example:

```
edit_distance(G₁, G₂, subst_cost=MinkowskiCost(μ₁, μ₂))
```
- `heuristic::Function=DefaultEditHeuristic`: a custom heuristic provided to the A*
search in case the default heuristic is not satisfactory.

### Performance
- Given two graphs ``|G₁| < |G₂|``, `edit_distance(G₁, G₂)` is faster to
compute than `edit_distance(G₂, G₁)`. Consider swapping the arguments
if involved costs are equivalent.
- The use of simple Minkowski costs can improve performance considerably.
- Exploit vertex attributes when designing operation costs.

### References
- RIESEN, K., 2015. Structural Pattern Recognition with Graph Edit Distance: Approximation Algorithms and Applications. (Chapter 2)

### Author
- Júlio Hoffimann Mendes (juliohm@stanford.edu)
"""
function edit_distance(G₁::AbstractGraph, G₂::AbstractGraph;
                       insert_cost::Function=v->1.0,
                       delete_cost::Function=u->1.0,
                       subst_cost::Function=(u,v)->0.5,
                       heuristic::Function=DefaultEditHeuristic)

  # A* search heuristic
  h(λ) = heuristic(λ, G₁, G₂)

  # initialize open set
  OPEN = PriorityQueue(Vector{Tuple}, Float64)
  for v in 1:nv(G₂)
    enqueue!(OPEN, [(1,v)], subst_cost(1,v) + h([(1,v)]))
  end
  enqueue!(OPEN, [(1,0)], delete_cost(1) + h([(1,0)]))

  while true
    # minimum (partial) edit path
    λ, cost = DataStructures.peek(OPEN)
    dequeue!(OPEN)

    if is_complete_path(λ, G₁, G₂)
      return cost, λ
    else
      k, _ = λ[end]
      vs = setdiff(1:nv(G₂), [v for (u,v) in λ])

      if k < nv(G₁) # there are still vertices to process in G₁?
        for v in vs
          λ⁺ = [λ; (k+1,v)]
          enqueue!(OPEN, λ⁺, cost + subst_cost(k+1,v) + h(λ⁺) - h(λ))
        end
        λ⁺ = [λ; (k+1,0)]
        enqueue!(OPEN, λ⁺, cost + delete_cost(k+1) + h(λ⁺) - h(λ))
      else
        # add remaining vertices of G₂ to the path
        λ⁺ = [λ; [(0,v) for v in vs]]
        total_insert_cost = sum(insert_cost, vs)
        enqueue!(OPEN, λ⁺, cost + total_insert_cost + h(λ⁺) - h(λ))
      end
    end
  end
end

function is_complete_path(λ, G₁, G₂)
  us = Set(); vs = Set()
  for (u,v) in λ
    push!(us, u)
    push!(vs, v)
  end
  delete!(us, 0)
  delete!(vs, 0)

  return length(us) == nv(G₁) && length(vs) == nv(G₂)
end

function DefaultEditHeuristic(λ, G₁::AbstractGraph, G₂::AbstractGraph)
  vs = Set([v for (u,v) in λ])
  delete!(vs, 0)

  return nv(G₂) - length(vs)
end

#-------------------------
# Edit path cost functions
#-------------------------

"""
    MinkowskiCost(μ₁, μ₂; p::Real=1)

For labels μ₁ on the vertices of graph G₁ and labels μ₂ on the vertices
of graph G₂, compute the p-norm cost of substituting vertex u ∈ G₁ by
vertex v ∈ G₂.

### Optional Arguments
`p=1`: the p value for p-norm calculation.
"""
function MinkowskiCost(μ₁::AbstractVector, μ₂::AbstractVector; p::Real=1)
  (u,v) -> norm(μ₁[u] - μ₂[v], p)
end

"""
    BoundedMinkowskiCost(μ₁, μ₂)

Return value similar to `MinkowskiCost`, but ensure costs smaller than 2τ.

### Optional Arguments
`p=1`: the p value for p-norm calculation.
`τ=1`: value specifying half of the upper limit of the Minkowski cost.
"""
function BoundedMinkowskiCost(μ₁::AbstractVector, μ₂::AbstractVector; p::Real=1, τ::Real=1)
  (u,v) -> 1 / (1/(2τ) + exp(-norm(μ₁[u] - μ₂[v], p)))
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
function weights(p::Array{<:Particle{T}}; logscale = true, normalscale = true) where {T<:Any}
    logw = vec(getfield.(p, :ℓw))
    if normalscale
        logw = normalise(logw)
    end
    if logscale
        return logw
    else
        return exp.(logw)
    end
end

function location(p::Array{Particle{T}}; flatten = false) where {T<:Any}
    if(flatten)
        loc = vec.(getfield.(p, :x))
    else
        loc = getfield.(p, :x)
    end
    return loc
end

function logdensity(p::Array{DynamicParticle{T}}) where {T<:Any}
    logdensityval = vec(getfield.(p, :cacheℓdens))
    return logdensityval
end


text:
using DataDeps

register(DataDep(
  "cpu-small",
  "https://www.openml.org/d/227",
  "https://www.openml.org/data/get_csv/3664/dataset_2213_cpu_small.arff",
  "5183289a479c5f2effbdfb1a221964bab6510ccd1d8e4264b48a15b3807611ae",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=13,
        feature_indices=1:12,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using EntropicToy
using Documenter

DocMeta.setdocmeta!(EntropicToy, :DocTestSetup, :(using EntropicToy); recursive=true)

makedocs(;
    modules=[EntropicToy],
    authors="Rethna Pulikkoonattu et al.",
    repo="https://github.com/geekymode/EntropicToy.jl/blob/{commit}{path}#{line}",
    sitename="EntropicToy.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "nothing") == "true",
        canonical="https://geekymode.github.io/EntropicToy.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/geekymode/EntropicToy.jl",
    devbranch="main",
)


text:
using PEGParser

@grammar mlgrammar begin
  start = sexpr{ _1 }
  sexpr  = "(" + list(expr, white_space) + ")"
  expr = r"[0-9]+"
  white_space = r"[ \t\n]+"
end


text:
const DEFAULT_PRIORITY = 1000
const DEFAULT_TEMPLATE_DIR = Ref{String}(joinpath(dirname(dirname(pathof(PkgTemplates))), "templates"))

"""
    @plugin struct ... end

Define a plugin subtype with keyword constructors and default values.

For details on the general syntax, see
[Parameters.jl](https://mauro3.github.io/Parameters.jl/stable/manual/#Types-with-default-values-and-keyword-constructors-1).

There are a few extra restrictions:

- Before using this macro, you must have imported `@with_kw_noshow`
  via `using PkgTemplates: @with_kw_noshow`
- The type must be a subtype of [`Plugin`](@ref) (or one of its abstract subtypes)
- The type cannot be parametric
- All fields must have default values

## Example

```julia
using PkgTemplates: @plugin, @with_kw_noshow, Plugin
@plugin struct MyPlugin <: Plugin
    x::String = "hello!"
    y::Union{Int, Nothing} = nothing
end
```

## Implementing `@plugin` Manually

If for whatever reason, you are unable to meet the criteria outlined above,
you can manually implement the methods that `@plugin` would have created for you.
This is only mandatory if you want to use your plugin in interactive mode.

### Keyword Constructors

If possible, use `@with_kw_noshow` to create a keyword constructor for your type.
Your type must be capable of being instantiated with no arguments.

### Default Values

If your type's fields have sensible default values, implement `defaultkw` like so:

```julia
using PkgTemplates: PkgTemplates, Plugin
struct MyPlugin <: Plugin
    x::String
end
PkgTemplates.defaultkw(::Type{MyPlugin}, ::Val{:x}) = "my default"
```

Remember to add a method to the function belonging to PkgTemplates,
rather than creating your own function that PkgTemplates won't see.

If your plugin's fields have no sane defaults, then you'll need to implement
[`prompt`](@ref) appropriately instead.
"""
macro plugin(ex::Expr)
    @assert ex.head === :struct "Expression must be a struct definition"
    @assert ex.args[2] isa Expr && ex.args[2].head === :<: "Type must have a supertype"
    T = ex.args[2].args[1]
    @assert T isa Symbol "@plugin does not work for parametric types"

    msg = "Run `using PkgTemplates: @with_kw_noshow` before using this macro"
    @assert isdefined(__module__, Symbol("@with_kw_noshow")) msg
    block = :(begin @with_kw_noshow $ex end)

    foreach(filter(arg -> arg isa Expr, ex.args[3].args)) do field
        @assert field.head === :(=) "Field must have a default value"
        name = QuoteNode(field.args[1].args[1])
        default = field.args[2]
        def = :(PkgTemplates.defaultkw(::Type{$T}, ::Val{$name}) = $default)
        push!(block.args, def)
    end

    return esc(block)
end

function Base.:(==)(a::T, b::T) where T <: Plugin
    return all(n -> getfield(a, n) == getfield(b, n), fieldnames(T))
end

struct Disabled{P<:Plugin} end
Base.:(!)(P::Type{<:Plugin}) = Disabled{P}()

"""
    Secret(name::AbstractString)

Represents a GitHub repository secret.
When converted to a string, yields `\${{ secrets.<name> }}`.
"""
struct Secret
    name::String
end

Base.print(io::IO, s::Secret) = print(io, "\${{ secrets.$(s.name) }}")

"""
A simple plugin that, in general, creates a single file.
"""
abstract type FilePlugin <: Plugin end

"""
    default_file(paths::AbstractString...) -> String

Return a path relative to the default template file directory
(`PkgTemplates/templates`).
"""
default_file(paths::AbstractString...) = joinpath(DEFAULT_TEMPLATE_DIR[], paths...)

"""
    view(::Plugin, ::Template, pkg::AbstractString) -> Dict{String, Any}

Return the view to be passed to the text templating engine for this plugin.
`pkg` is the name of the package being generated.

For [`FilePlugin`](@ref)s, this is used for both the plugin badges
(see [`badges`](@ref)) and the template file (see [`source`](@ref)).
For other [`Plugin`](@ref)s, it is used only for badges,
but you can always call it yourself as part of your [`hook`](@ref) implementation.

By default, an empty `Dict` is returned.
"""
view(::Plugin, ::Template, ::AbstractString) = Dict{String, Any}()

"""
    user_view(::Plugin, ::Template, pkg::AbstractString) -> Dict{String, Any}

The same as [`view`](@ref), but for use by package *users* for extension.

Values returned by this function will override those from [`view`](@ref)
when the keys are the same.
"""
user_view(::Plugin, ::Template, ::AbstractString) = Dict{String, Any}()

"""
    combined_view(::Plugin, ::Template, pkg::AbstractString) -> Dict{String, Any}

This function combines [`view`](@ref) and [`user_view`](@ref) for use in text templating.
If you're doing manual file creation or text templating (i.e. writing [`Plugin`](@ref)s
that are not [`FilePlugin`](@ref)s), then you should use this function
rather than either of the former two.

!!! note
    Do not implement this function yourself!
    If you're implementing a plugin, you should implement [`view`](@ref).
    If you're customizing a plugin as a user, you should implement [`user_view`](@ref).
"""
function combined_view(p::Plugin, t::Template, pkg::AbstractString)
    return merge(view(p, t, pkg), user_view(p, t, pkg))
end

"""
    tags(::Plugin) -> Tuple{String, String}

Return the delimiters used for text templating.
See the [`Citation`](@ref) plugin for a rare case where changing the tags is necessary.

By default, the tags are `"{{"` and `"}}"`.
"""
tags(::Plugin) = "{{", "}}"

"""
    priority(::Plugin, ::Union{typeof(prehook), typeof(hook), typeof(posthook)}) -> Int

Determines the order in which plugins are processed (higher goes first).
The default priority (`DEFAULT_PRIORITY`), is `$DEFAULT_PRIORITY`.

You can implement this function per-stage (by using `::typeof(hook)`, for example),
or for all stages by simply using `::Function`.
"""
priority(::Plugin, ::Function) = DEFAULT_PRIORITY

"""
    gitignore(::Plugin) -> Vector{String}

Return patterns that should be added to `.gitignore`.
These are used by the [`Git`](@ref) plugin.

By default, an empty list is returned.
"""
gitignore(::Plugin) = String[]

"""
    badges(::Plugin) -> Union{Badge, Vector{Badge}}

Return a list of [`Badge`](@ref)s, or just one, to be added to `README.md`.
These are used by the [`Readme`](@ref) plugin to add badges to the README.

By default, an empty list is returned.
"""
badges(::Plugin) = Badge[]

"""
    source(::FilePlugin) -> Union{String, Nothing}

Return the path to a plugin's template file, or `nothing` to indicate no file.

By default, `nothing` is returned.
"""
source(::FilePlugin) = nothing

"""
    destination(::FilePlugin) -> String

Return the destination, relative to the package root, of a plugin's configuration file.

This function **must** be implemented.
"""
function destination end

"""
    Badge(hover::AbstractString, image::AbstractString, link::AbstractString)

Container for Markdown badge data.
Each argument can contain placeholders,
which will be filled in with values from [`combined_view`](@ref).

## Arguments
- `hover::AbstractString`: Text to appear when the mouse is hovered over the badge.
- `image::AbstractString`: URL to the image to display.
- `link::AbstractString`: URL to go to upon clicking the badge.
"""
struct Badge
    hover::String
    image::String
    link::String
end

Base.string(b::Badge) = "[![$(b.hover)]($(b.image))]($(b.link))"

# Format a plugin's badges as a list of strings, with all substitutions applied.
function badges(p::Plugin, t::Template, pkg::AbstractString)
    bs = badges(p)
    bs isa Vector || (bs = [bs])
    return map(b -> render_text(string(b), combined_view(p, t, pkg)), bs)
end

"""
    validate(::Plugin, ::Template)

Perform any required validation for a [`Plugin`](@ref).

It is preferred to do validation here instead of in [`prehook`](@ref),
because this function is called at [`Template`](@ref) construction time,
whereas the prehook is only run at package generation time.
"""
validate(::Plugin, ::Template) = nothing

"""
    prehook(::Plugin, ::Template, pkg_dir::AbstractString)

Stage 1 of the package generation process (the "before" stage, in general).
At this point, `pkg_dir` is an empty directory that will eventually contain the package,
and neither the [`hook`](@ref)s nor the [`posthook`](@ref)s have run.

!!! note
    `pkg_dir` only stays empty until the first plugin chooses to create a file.
    See also: [`priority`](@ref).
"""
prehook(::Plugin, ::Template, ::AbstractString) = nothing

"""
    hook(::Plugin, ::Template, pkg_dir::AbstractString)

Stage 2 of the package generation pipeline (the "main" stage, in general).
At this point, the [`prehook`](@ref)s have run, but not the [`posthook`](@ref)s.

`pkg_dir` is the directory in which the package is being generated
(so `basename(pkg_dir)` is the package name).

!!! note
    You usually shouldn't implement this function for [`FilePlugin`](@ref)s.
    If you do, it should probably `invoke` the generic method
    (otherwise, there's not much reason to subtype `FilePlugin`).
"""
hook(::Plugin, ::Template, ::AbstractString) = nothing

"""
    posthook(::Plugin, ::Template, pkg_dir::AbstractString)

Stage 3 of the package generation pipeline (the "after" stage, in general).
At this point, both the [`prehook`](@ref)s and [`hook`](@ref)s have run.
"""
posthook(::Plugin, ::Template, ::AbstractString) = nothing

function validate(p::T, ::Template) where T <: FilePlugin
    src = source(p)
    src === nothing && return
    isfile(src) || throw(ArgumentError("$(nameof(T)): The file $src does not exist"))
end

function hook(p::FilePlugin, t::Template, pkg_dir::AbstractString)
    source(p) === nothing && return
    pkg = basename(pkg_dir)
    path = joinpath(pkg_dir, destination(p))
    text = render_plugin(p, t, pkg)
    gen_file(path, text)
end

function render_plugin(p::FilePlugin, t::Template, pkg::AbstractString)
    return render_file(source(p), combined_view(p, t, pkg), tags(p))
end

"""
    gen_file(file::AbstractString, text::AbstractString)

Create a new file containing some given text.
Trailing whitespace is removed, and the file will end with a newline.
"""
function gen_file(file::AbstractString, text::AbstractString)
    mkpath(dirname(file))
    text = strip(join(map(rstrip, split(text, "\n")), "\n")) * "\n"
    write(file, text)
end

"""
    render_file(file::AbstractString view::Dict{<:AbstractString}, tags=nothing) -> String

Render a template file with the data in `view`.
`tags` should be a tuple of two strings, which are the opening and closing delimiters,
or `nothing` to use the default delimiters.
"""
function render_file(file::AbstractString, view::Dict{<:AbstractString}, tags=nothing)
    return render_text(read(file, String), view, tags)
end

"""
    render_text(text::AbstractString, view::Dict{<:AbstractString}, tags=nothing) -> String

Render some text with the data in `view`.
`tags` should be a tuple of two strings, which are the opening and closing delimiters,
or `nothing` to use the default delimiters.
"""
function render_text(text::AbstractString, view::Dict{<:AbstractString}, tags=nothing)
    return tags === nothing ? render(text, view) : render(text, view; tags=tags)
end

"""
    needs_username(::Plugin) -> Bool

Determine whether or not a plugin needs a Git hosting service username to function correctly.
If you are implementing a plugin that uses the `user` field of a [`Template`](@ref),
you should implement this function and return `true`.
"""
needs_username(::Plugin) = false

include(joinpath("plugins", "project_file.jl"))
include(joinpath("plugins", "src_dir.jl"))
include(joinpath("plugins", "tests.jl"))
include(joinpath("plugins", "readme.jl"))
include(joinpath("plugins", "license.jl"))
include(joinpath("plugins", "git.jl"))
include(joinpath("plugins", "tagbot.jl"))
include(joinpath("plugins", "develop.jl"))
include(joinpath("plugins", "coverage.jl"))
include(joinpath("plugins", "ci.jl"))
include(joinpath("plugins", "compat_helper.jl"))
include(joinpath("plugins", "citation.jl"))
include(joinpath("plugins", "documenter.jl"))
include(joinpath("plugins", "badges.jl"))
include(joinpath("plugins", "benchmarks.jl"))


text:
@cache mutable struct KenCarp3ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = KenCarp3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getoopnlsolvefields

  KenCarp3ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct KenCarp3Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F,kType} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  k1::kType
  k2::kType
  k3::kType
  k4::kType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = KenCarp3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getiipnlsolvefields

  if typeof(f) <: SplitFunction
    k1 = similar(u); k2 = similar(u)
    k3 = similar(u); k4 = similar(u)
  else
    k1 = nothing; k2 = nothing
    k3 = nothing; k4 = nothing
    uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  end

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = z
  atmp = similar(u,uEltypeNoUnits)

  KenCarp3Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,k1,k2,k3,k4,dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end

@cache mutable struct Kvaerno4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Kvaerno4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getoopnlsolvefields
  Kvaerno4ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct Kvaerno4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Kvaerno4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getiipnlsolvefields

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u); z₅ = z
  atmp = similar(u,uEltypeNoUnits)

  Kvaerno4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end

@cache mutable struct KenCarp4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = KenCarp4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getoopnlsolvefields
  KenCarp4ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct KenCarp4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F,kType} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  z₆::uType
  k1::kType
  k2::kType
  k3::kType
  k4::kType
  k5::kType
  k6::kType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = KenCarp4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getiipnlsolvefields

  if typeof(f) <: SplitFunction
    k1 = similar(u); k2 = similar(u)
    k3 = similar(u); k4 = similar(u)
    k5 = similar(u); k6 = similar(u)
  else
    k1 = nothing; k2 = nothing
    k3 = nothing; k4 = nothing
    k5 = nothing; k6 = nothing
    uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  end

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u); z₅ = zero(u)
  z₆ = z
  atmp = similar(u,uEltypeNoUnits)

  KenCarp4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,z₆,k1,k2,k3,k4,k5,k6,
                dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end

@cache mutable struct Kvaerno5ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Kvaerno5Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getoopnlsolvefields

  Kvaerno5ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct Kvaerno5Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  z₆::uType
  z₇::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Kvaerno5Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getiipnlsolvefields

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u); z₅ = zero(u)
  z₆ = zero(u); z₇ = z
  atmp = similar(u,uEltypeNoUnits)

  Kvaerno5Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,z₆,z₇,dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end

@cache mutable struct KenCarp5ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = KenCarp5Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getoopnlsolvefields

  KenCarp5ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct KenCarp5Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F,kType} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  z₆::uType
  z₇::uType
  z₈::uType
  k1::kType
  k2::kType
  k3::kType
  k4::kType
  k5::kType
  k6::kType
  k7::kType
  k8::kType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::KenCarp5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = KenCarp5Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.c3
  J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
  nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
  @getiipnlsolvefields

  if typeof(f) <: SplitFunction
    k1 = similar(u); k2 = similar(u)
    k3 = similar(u); k4 = similar(u)
    k5 = similar(u); k6 = similar(u)
    k7 = similar(u); k8 = similar(u)
  else
    k1 = nothing; k2 = nothing
    k3 = nothing; k4 = nothing
    k5 = nothing; k6 = nothing
    k7 = nothing; k8 = nothing
    uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  end

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u)
  z₅ = zero(u); z₆ = zero(u); z₇ = zero(u); z₈ = z
  atmp = similar(u,uEltypeNoUnits)

  KenCarp5Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,
                k1,k2,k3,k4,k5,k6,k7,k8,
                dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end


text:
@testset "191.number-of-1-bits.jl" begin
    @test hamming_weight(0b00000000000000000000000000001011) == 3
    @test hamming_weight(0b00000000000000000000000010000000) == 1
    @test hamming_weight(0b11111111111111111111111111111101) == 31
end


text:
let sys = DenseSystem(test_params...),
    beta = 12.34,
    P = 256,
    sampling_sys = DiagonalSystem(test_params_sampling...)

    @test isapprox(ground_pes(sys, (-1.0, 2.0, -3.0, 4.0); lengths=(2, 3)), [11.73057536 14.61381238; -1.50969268 1.59347336; 22.29332112 25.61509171])

    @test isapprox(path_mean_std(sampling_sys, beta, P, 1, 1), 0.201370345)
    @test isapprox(path_mean_std(sampling_sys, beta, P, test_S_sampling, test_M), 0.108103579)
end


text:
using Parameters: @unpack
using ..Continuations: find_zero!
using ..BifurcationsBase: AbstractSpecialPoint, special_intervals,
    special_points, push_special_point!,
    contkind, FixedPointCont, ContinuationKind

function default_resolve_exception_handler(err, interval, warn_exceptions)
    if any(isa.(Ref(err), warn_exceptions))
        @warn(
            """
            Failed to find bifurcation point within:
            $(interval)
            """,
            exception = err)
        return true
    end
    return false
end

default_resolve_exception_handler(warn_exceptions::Tuple) =
    (args...) -> default_resolve_exception_handler(args..., warn_exceptions)

function resolving_points(
        solver::BifurcationSolver, args...;
        warn_exceptions = (SingularException,),
        exception_handler = default_resolve_exception_handler(warn_exceptions),
        )
    Iterators.filter(
        point -> point !== nothing,
        try
            resolve_point(interval, solver)
        catch err
            # [[../continuations/branching.jl::SingularException]]
            if ! exception_handler(err, interval)
                rethrow()
            end
            nothing
        end for interval in special_intervals(solver, args...))
end

resolved_points(solver::BifurcationSolver, args...; kwargs...) =
    collect(resolving_points(solver, args...; kwargs...))

function resolve_point(point::AbstractSpecialPoint, solver::BifurcationSolver)
    super = as(solver, ContinuationSolver)
    return resolve_point(point, super.cache, super.opts)
end

"""
    resolve_points!(solver; kwargs...)

Resolve special points if not.  It is a no-op if all points are
already resolved.
"""
function resolve_points!(
        solver::BifurcationSolver;
        warn_exceptions = (SingularException,),
        exception_handler = default_resolve_exception_handler(warn_exceptions),
        )
    for sweep in solver.sol.sweeps
        ofs = length(sweep.special_points) + 1
        for interval in @view sweep.special_intervals[ofs:end]
            point = try
                resolve_point!(interval,
                               as(solver.cache, ContinuationCache),
                               solver.opts)
            catch err
                if !exception_handler(err, interval)
                    rethrow()
                end
                continue
            end
            push_special_point!(sweep, point)
        end
    end
end

function special_points!(solver, args...; kwargs...)
    resolve_points!(solver; kwargs...)
    return special_points(solver, args...)
end

resolve_point(point::AbstractSpecialPoint,
              cache::ContinuationCache,
              args...) =
    resolve_point!(point, deepcopy(cache), args...)

resolve_point!(point::SpecialPoint, _...) = point  # already resolved

function resolve_point!(point::SpecialPointInterval,
                        cache::ContinuationCache,
                        opts::ContinuationOptions,
                        ) :: SpecialPoint
    @unpack point_type, point_index = point
    f = testfn_for(point_type, timekind(cache), contkind(cache),
                   cache.prob_cache)
    direction = as(point.sweep.value, ContinuationSweep).direction
    u, tJ, L, Q, J =
        find_zero!(cache, opts, f, point.u0, point.u1, direction)
    return SpecialPoint(timekind(point),
                        point_type, point_index, u, J,
                        WeakRef(point.sweep.value))
end

function testfn_for(point_type::Enum, tkind::TimeKind, ckind::ContinuationKind,
                    prob_cache)
    ptype = Val{point_type}()
    # For call signature of `f`, see:
    # [[../continuations/zero_point.jl::f(]]
    return (u, J, L, Q) -> testfn(ptype, tkind, ckind, prob_cache, u, J, L, Q)
end

testfn(::Val{PointTypes.saddle_node}, ::Continuous, ::FixedPointCont,
       prob_cache, u, J, L, Q) =
    det(ds_jacobian(J))

function testfn(::Val{PointTypes.hopf}, ::Continuous, ::FixedPointCont,
                prob_cache, u, J, L, Q)
    vals = _eigvals(ds_jacobian(J))
    _, i = findmin(abs.(real.(vals)))
    return real(vals[i])
end

testfn(::Val{PointTypes.saddle_node}, ::Discrete, ::FixedPointCont,
       prob_cache, u, J, L, Q) =
    det(ds_jacobian(J) - I)

testfn(::Val{PointTypes.period_doubling}, ::Discrete, ::FixedPointCont,
       prob_cache, u, J, L, Q) =
    det(ds_jacobian(J) + I)


text:
function get_fourbar(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    model=:fourbar,
    spring=0.0,
    damper=0.0,
    parse_damper=true,
    T=Float64)

    path = joinpath(@__DIR__, "../deps/$(String(model)).urdf")
    mech = Mechanism(path; floating=false, T,
        gravity,
        timestep,
        parse_damper)

    # Adding springs and dampers
    set_springs!(mech.joints, spring)
    set_dampers!(mech.joints, damper)

    return mech
end

function initialize_fourbar!(mechanism::Mechanism;
    angle=0.0,
    angular_velocity=szeros(2)) where T

    zero_velocity!(mechanism)
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :jointb1);
        xmin=[ -angle, angular_velocity[1]])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :joint12);
        xmin=[+2angle, 0])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :jointb3);
        xmin=[ +angle, angular_velocity[2]])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :joint34);
        xmin=[-2angle, 0])
end


text:
module sparsematrix

const SparseMatrix = SparseMatrixCSC{Complex128, Int}

function gemm_sp_dense_small(alpha::Complex128, M::SparseMatrix, B::Matrix{Complex128}, result::Matrix{Complex128})
    if alpha == Complex128(1.)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = alpha*M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    end
end

function gemm_sp_dense_big(alpha::Complex128, M::SparseMatrix, B::Matrix{Complex128}, result::Matrix{Complex128})
    if alpha == Complex128(1.)
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    else
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = alpha*B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    end
end

function gemm!(alpha::Complex128, M::SparseMatrix, B::Matrix{Complex128}, beta::Complex128, result::Matrix{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(result, beta)
    end
    if nnz(M) > 1000
        gemm_sp_dense_big(alpha, M, B, result)
    else
        gemm_sp_dense_small(alpha, M, B, result)
    end
end

function gemm!(alpha::Complex128, B::Matrix{Complex128}, M::SparseMatrix, beta::Complex128, result::Matrix{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(result, beta)
    end
    dimB = size(result,1)
    if alpha == Complex128(1.)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]*alpha
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    end
end

function gemv!(alpha::Complex128, M::SparseMatrix, v::Vector{Complex128}, beta::Complex128, result::Vector{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(result, beta)
    end
    if alpha == Complex128(1.)
        @inbounds for colindex = 1:M.n
            vj = v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    else
        @inbounds for colindex = 1:M.n
            vj = alpha*v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    end
end

function gemv!(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrix, beta::Complex128, result::Vector{Complex128})
    if beta == Complex128(0.)
        fill!(result, beta)
    elseif beta != Complex128(1.)
        scale!(result, beta)
    end
    if alpha == Complex128(1.)
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*v[M.rowval[i]]
            end
        end
    else
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
            end
        end
    end
end

function sub2sub(shape1::NTuple{N, Int}, shape2::NTuple{M, Int}, I::CartesianIndex{N}) where {N, M}
    linearindex = sub2ind(shape1, I.I...)
    CartesianIndex(ind2sub(shape2, linearindex)...)
end

function ptrace(x, shape_nd::Vector{Int}, indices::Vector{Int})
    shape_nd = (shape_nd...)
    N = div(length(shape_nd), 2)
    shape_2d = (x.m, x.n)
    shape_nd_after = ([i ∈ indices || i-N ∈ indices ? 1 : shape_nd[i] for i=1:2*N]...)
    shape_2d_after = (prod(shape_nd_after[1:N]), prod(shape_nd_after[N+1:end]))
    I_nd_after_max = CartesianIndex(shape_nd_after...)
    y = spzeros(Complex128, shape_2d_after...)
    for I in eachindex(x)
        I_nd = sub2sub(shape_2d, shape_nd, I)
        if I_nd.I[indices] != I_nd.I[indices + N]
            continue
        end
        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max))
        y[I_after] += x[I]
    end
    y
end

function permutedims(x, shape, perm)
    shape = (shape...)
    shape_perm = ([shape[i] for i in perm]...)
    y = spzeros(Complex128, x.m, x.n)
    for I in eachindex(x)
        linear_index = sub2ind((x.m, x.n), I.I...)
        cartesian_index = ind2sub(shape, linear_index)
        cartesian_index_perm = [cartesian_index[i] for i=perm]
        linear_index_perm = sub2ind(shape_perm, cartesian_index_perm...)
        J = ind2sub((x.m, x.n), linear_index_perm)
        y[J...] = x[I.I...]
    end
    y
end

end # module


text:
@testset "Function fields" begin
  include("FunField/DegreeLocalization.jl")
  include("FunField/Round2.jl")
end


text:
using SparseArrays, LightXML

function osm2lar(
        filename::String,
        doCrop::Bool = false,
        doFrame::Bool = false
    )::Lar.Model
    xdoc  = parse_file(filename)
    xroot = root(xdoc)

    cropBox   = get_elements_by_tagname(xroot, "bounds")[1]
    node_list = get_elements_by_tagname(xroot, "node"  )
    way_list  = get_elements_by_tagname(xroot, "way"   )

    #=== EXTRACT VERTEX SET AND VERTEX POSITION MAP ===#

    Vidxs = [parse(Int, attribute(v, "id")) for v in node_list]
    Vmax  = max(Vidxs...)
    Vmap  = SparseArrays.spzeros(Int, Vmax)
    Vnum  = length(Vidxs)
    for i = 1 : Vnum   Vmap[Vidxs[i]] = i   end

    #=== EXTRACT (X,Y) AND RESCALE POINTS IN [0, 1] ===#
    lons  = [parse(Float64, attribute(v, "lon")) for v in node_list]
    lmin  = parse(Float64, attribute(cropBox, "minlon"))
    lmax  = parse(Float64, attribute(cropBox, "maxlon"))
    lons  = (lons .- lmin) ./ (lmax - lmin)

    lats  = [parse(Float64, attribute(v, "lat")) for v in node_list]
    lmin  = parse(Float64, attribute(cropBox, "minlat"))
    lmax  = parse(Float64, attribute(cropBox, "maxlat"))
    lats  = (lats .- lmin) ./ (lmax - lmin)

    #=== BUILD LAR MODEL ===#

    model = Lar.Model([lons'; lats'])

    for street in way_list
        nodes = [
                 Vmap[parse(Int, attribute(n, "ref"))]
                 for n in get_elements_by_tagname(street, "nd")
                ]
        for i = 1 : length(nodes) - 1
            newcell = SparseArrays.spzeros(Int8, Vnum)
            newcell[nodes[  i]] = 1
            newcell[nodes[i+1]] = 1
            Lar.addModelCells!(model, 1, convert(Lar.ChainOp, newcell'))
        end
    end

    if doCrop
        Lar.mapBoxCrop!(model, [0.0 1.0;0.0 1.0], doFrame)
    end
    return model
end

function mapBoxCrop!(
        model::Lar.Model,
        cropBox::Lar.Points,
        doFrame::Bool
    )::Nothing

    toDelV = convert(Array{Int64,1},
        [i  for i = 1 : size(model, 0, 2)
            if model.G[1, i] <= cropBox[1, 1]
            || model.G[1, i] >= cropBox[1, 2]
            || model.G[2, i] <= cropBox[2, 1]
            || model.G[2, i] >= cropBox[2, 2]
        ]
    )
    Lar.deleteModelVertices!(model, toDelV)

    if doFrame
        idx = size(model, 0, 2) + 1
        Lar.addModelVertices!(model, [
            cropBox[1, 1] cropBox[1, 1] cropBox[1, 2] cropBox[1, 2]
            cropBox[2, 1] cropBox[2, 2] cropBox[2, 2] cropBox[2, 1]
        ])
        Lar.addModelCells!(model, 1, [
            [idx, idx+1], [idx+1, idx+2], [idx+2, idx+3], [idx+3, idx]
        ])
    end
    return
end


text:
using Test
using CartesianGeneticProgramming
import Cambrian
import Random

cfg = get_config("test.yaml")

@testset "Mutation" begin

    parent = CGPInd(cfg)

    # Uniform mutation
    child = uniform_mutate(cfg, parent)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Goldman mutation : ensure structural difference
    child = goldman_mutate(cfg, parent)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    # Profiling mutation: ensure output different for provided inputs
    inputs = rand(cfg.n_in, 1)
    child = profiling_mutate(cfg, parent, inputs)
    @test any(parent.chromosome .!= child.chromosome)
    @test any(parent.genes .!= child.genes)

    out_parent = process(parent, inputs[:, 1])
    out_child = process(child, inputs[:, 1])
    @test any(out_parent .!= out_child)
end


text:
# julia --project=@. test/runtests.jl
#
# julia> include("test/runtests.jl")
using Test

# see https://discourse.julialang.org/t/writing-tests-in-vs-code-workflow-autocomplete-and-tooltips/57488
# see https://github.com/julia-vscode/julia-vscode/issues/800
if isdefined(@__MODULE__, :LanguageServer)
    # invoked by VS Code
    include("../src/Classic.jl")
    using .Classic
else
    # invoked during test
    using Classic
end

include("test_ch2.jl")
include("test_ch4.jl")
include("test_gene.jl")
include("test_kmeans.jl")


text:
using SpecialFunctions, IterTools

"Normalises a vector"
function normalise(x::AbstractArray)
    if length(x) == 0
        error("cannot normalise a vector of length 0")
    end
    return x/sum(x)
end

lgamma_local(x) = SpecialFunctions.logabsgamma(x)[1]
const precomputed_lfact = SpecialFunctions.logfactorial.(1:10000)
function lgamma_local(x::Integer)
    if x < 10000 && x > 2
        return precomputed_lfact[x-1]
    else
        # return SpecialFunctions.logfactorial(x)
        return SpecialFunctions.logabsgamma(x)[1]
    end
end

function kmax(x::AbstractArray{T, 1}, k::Integer) where T <: Number
    if length(x) < k
        error("length(x) < $k, cannot take the $k largest elements of a vector of size $(length(x))")
    else
        res = x[(end-k+1):end]
        return kmax_rec(x, k, findmin(res), res)
    end
end

# import Base.length
#
# function length(x::Union{IterTools.Distinct})
#     l = 0
#     for k in x
#         l +=1
#     end
#     return l
# end

function kmax_rec(x::AbstractArray{T, 1}, k::Integer, smallest::Tuple{T,U}, res::AbstractArray{T, 1}) where {T <: Number, U <: Integer}
    # println("x = $x")
    # println("smallest = $smallest")
    # println("res = $res")
    if length(x) == k
        return res
    else
        if x[1] > smallest[1]
            res[smallest[2]] = x[1]
        end
        return kmax_rec(x[2:end], k, findmin(res), res)
    end
end

# function precompute_lgamma_α(α, data)
#     return lgamma_local.(i + α for i in 0:sum(sum.(values(data))))
# end
#
# function precompute_lfactorial(data)
#     return  SpecialFunctions.logfactorial.(1:maximum(maximum.(values(data))))
# end

"Converts a dictionary to a new dictionary with log values"
function convert_weights_to_logweights(weights_dict)
    return Dict(k => log.(v) for (k,v) in weights_dict)
end

"Converts a dictionary to a new dictionary with exponential values"
function convert_logweights_to_weights(weights_dict)
    return Dict(k => exp.(v) for (k,v) in weights_dict)
end

using Memoization
@memoize function log_pochammer_rec(x::Real, n::Integer)
    # @info "n=$n"
    if n==0
        return 0.
    elseif n==1
        return log(x)
    else
        return log(x + n - 1)  + log_pochammer_rec(x, n-1)
    end
end

function test_equal_spacing_of_observations(data; override = false, digits_after_comma_for_time_precision = 4)
    if !override&&(data |> keys |> collect |> sort |> diff |> x -> round.(x; digits = digits_after_comma_for_time_precision) |> unique |> length > 1)
        println(data |> keys |> collect |> sort |> diff |> x -> truncate_float.(x,digits_after_comma_for_time_precision) |> unique)
        error("Think twice about precomputing all terms, as the time intervals are not equal. You can go ahead using the option 'override = true.'")
    end
end

function log_binomial_safe_but_slow(n::Int64, k::Int64)
    @assert n >= 0
    @assert k >= 0
    @assert k <= n
    if k == 0 || k == n
        return 0
    elseif k == 1 || k == n-1
        return log(n)
    else
        return sum(log(i) for i in (n-k+1):n) - sum(log(i) for i in 2:k)
    end
end

function log_descending_fact_no0(x::Real, n::Int64)
    return sum(log(x-i) for i in 0:(n-1))
end

function log_pochammer(x::Real, n::Integer)
    if n==0
        return 0.
    elseif n==1
        return log(x)
    else
        return sum(log(x + i) for i in 0:(n-1))
    end
end

function loghypergeom_pdf_using_precomputed(i, m, si::Integer, sm::Integer, log_binomial_coeff_ar_offset::Array{Float64,2})
    return sum(log_binomial_coeff_ar_offset[m[k]+1,i[k]+1] for k in eachindex(m)) - log_binomial_coeff_ar_offset[sm+1, si+1]
end

function assert_constant_time_step_and_compute_it(data)
    Δts = keys(data) |> collect |> sort |> diff |> unique
    if length(Δts) > 1
        test_equal_spacing_of_observations(data; override = false)
    end
    Δt = mean(Δts)
    return Δt
end

function sample_from_Gamma_mixture(δ, θ, Λ, wms)
    #use 1/θ because of the way the Gamma distribution is parameterised in Julia Distributions.jl

    latent_mixture_idx = rand(Categorical(wms))
    return rand(Gamma(δ/2 + Λ[latent_mixture_idx], 1/θ))
end
function create_gamma_mixture_parameters(δ, θ, Λ)
    α = [δ/2 + m for m in Λ]
    β = [θ for m in Λ]
    return α, β
end
function create_gamma_mixture_pdf(δ, θ, Λ, wms)
    #use 1/θ because of the way the Gamma distribution is parameterised in Julia Distributions.jl
    return x -> sum(wms.*Float64[pdf(Gamma(δ/2 + m, 1/θ),x) for m in Λ])
end

function create_dirichlet_mixture(α::Array{T, 1}, Λ::Array{Array{U,1},1}) where {T <: Real, U <:Integer}
    α_mixt = Array{Array{T,1},1}(undef, length(Λ))
    for i in eachindex(Λ)
        α_mixt[i] = α .+ Λ[i]
    end
    return α_mixt
end

create_dirichlet_mixture_parameters = create_dirichlet_mixture #alis for consistency with CIR

function get_quantiles_from_mass(mass)
    qinf = 0.5*(1-mass)
    return (qinf, 1-qinf)
end

function create_Gamma_mixture_density(δ, θ, Λ, wms)
    #use 1/θ because of the way the Gamma distribution is parameterised in Julia Distributions.jl
    return x -> sum(wms.*Float64[pdf(Gamma(δ/2 + m, 1/θ),x) for m in Λ])
end

function compute_quantile_mixture_hpi(δ, θ, Λ, wms, q::Float64)
    #use 1/θ because of the way the Gamma distribution is parameterised in Julia Distributions.jl
    f = x -> sum(wms.*Float64[cdf(Gamma(δ/2 + m, 1/θ),x) for m in Λ])
    return fzero(x -> f(x)-q, 0, 10^9)
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Test, Random, EndowmentsLH

function quantile_test()
    rng = MersenneTwister(45);
    @testset "Quantiles" begin
        n = 50;
        ed = EndowmentsLH.make_test_endowment_draws(n);
        # pctV = LinRange(0.1, 0.9, n);
        for dName in names(ed)
            drawV = get_draws(ed, dName);
            if eltype(drawV) <: Real
                qV = endow_quantiles(ed, dName, [0.3, 0.7]);
                @test eltype(qV) == eltype(drawV)
                @test size(qV) == (2,)
                @test all(diff(qV) .> 0.0)

                q = endow_quantiles(ed, dName, 0.7);
                @test isapprox(qV[2], q)
                nLess = sum(get_draws(ed, dName) .<= q);
                @test abs(nLess / n - 0.7) < 1.5 / n
            end
        end
    end
end

@testset "Quantiles" begin
    quantile_test();
end

# -----------


text:
#these functions are copied, and slightly modified, from the code included with the GMD paper

###############################################################################
#
# As the land model is modularized, there could be multiple places that conatin
#     exact the same information. Thus, these information needs to be synced.
#     This type of information includes leaf area, hydraulic traits, and
#     photosynthetic traits.
#
# The following function is meant for leaf area
#
###############################################################################
"""
    update_LAI!(node::SPACMono{FT}, lai::FT) where {FT<:AbstractFloat}
Update SPAC leaf area index
"""
function update_LAI!(node::SPACMono{FT}, lai::FT) where {FT<:AbstractFloat}
    # 1. Update the general information of the SPAC
    node.la = node.ga * lai;

    # 2. Update the LAI for canopy radiation model
    node.canopy_rt.LAI = lai;

    # 3. Update the LA and LAI for stomatal and photosynthesis models
    for _iPS in node.plant_ps
        _iPS.LA   = node.la / node.n_canopy;
        _iPS.LAI  = lai / node.n_canopy;
        _iPS.tLAI = lai;
    end

    # 4. Update the LA for hydraulic model
    #    TODO abstract this to avoid memory allocation
    for _iHS in node.plant_hs.leaves
        _iHS.area = node.la / node.n_canopy;
    end

    return nothing
end

"""
    update_Kmax!(node::SPACMono{FT}, kmax::FT) where {FT<:AbstractFloat}
Update the maximal hydraulic conductance for SPAC
TODO need to abstractize this to PlantHydraulics.jl
"""
function update_Kmax0!(node::SPACMono{FT}, kmax::FT) where {FT<:AbstractFloat}
    # Root:Stem:Leaf conductance ratio is set as 1:2:2.
    #    For trees, roots: 2X, trunk: 8X, branch: 8X, leaves: 4X
    #    For palms, roots: 2X, trunk: 4X, leaves: 4X
    #    For grasses, roots: 2X, leaves: 2X

    # 1. update the hydraulic conductance in Roots
    for root in node.plant_hs.roots
        root.k_max      = 2 * kmax / node.plant_hs.n_root;
        root.k_element .= root.k_max * root.N;
    end

    # 2. If is a tree
    if typeof(node.plant_hs) <: TreeLikeOrganism
        node.plant_hs.trunk.k_max      = 8 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for stem in node.plant_hs.branch
            stem.k_max      = 8 * kmax / node.plant_hs.n_canopy;
            stem.k_element .= stem.k_max * stem.N;
        end
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 3. If is a palm
    if typeof(node.plant_hs) <: PalmLikeOrganism
        node.plant_hs.trunk.k_max      = 4 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 4. If is a grass
    if typeof(node.plant_hs) <: GrassLikeOrganism
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 2 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    return nothing
end

function update_Kmax_ratio!(node::SPACMono{FT}, kmax::FT, ratio::Array{FT}) where {FT<:AbstractFloat}
    # Root:Stem:Leaf conductance ratio is set as 1:2:2.
    #    For trees, roots: 2X, trunk: 8X, branch: 8X, leaves: 4X
    #    For palms, roots: 2X, trunk: 4X, leaves: 4X
    #    For grasses, roots: 2X, leaves: 2X

	resist_sum = sum(1 ./ ratio);
	ratio_norm = ratio * resist_sum;

    # 1. update the hydraulic conductance in Roots
    for root in node.plant_hs.roots
        root.k_max      = ratio_norm[1] * kmax / node.plant_hs.n_root;
        root.k_element .= root.k_max * root.N;
    end

    # 2. If is a tree
    if typeof(node.plant_hs) <: TreeLikeOrganism
        node.plant_hs.trunk.k_max      = ratio_norm[2] * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for stem in node.plant_hs.branch
            stem.k_max      = ratio_norm[3] * kmax / node.plant_hs.n_canopy;
            stem.k_element .= stem.k_max * stem.N;
        end
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = ratio_norm[4] * kmax / node.plant_hs.n_canopy / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    return nothing
end

"""
    update_VJRWW!(node::SPACMono{FT}, vcmax::FT) where {FT<:AbstractFloat}
Update Vcmax25(WW), Jmax25(WW), and Rd25(WW) from a given Vcmax25
"""
function update_VJRWW!(node::SPACMono{FT}, vcmax::FT) where {FT<:AbstractFloat}
    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)
    for _iPS in node.plant_ps
        _iPS.ps.Vcmax     = vcmax;
        _iPS.ps.Vcmax25   = vcmax;
        _iPS.ps.Vcmax25WW = vcmax;
        _iPS.ps.Jmax      = vcmax * 1.67;
        _iPS.ps.Jmax25    = vcmax * 1.67;
        _iPS.ps.Jmax25WW  = vcmax * 1.67;
        _iPS.ps.Rd        = vcmax * 0.015;
        _iPS.ps.Rd25      = vcmax * 0.015;
        _iPS.ps.Rd25WW    = vcmax * 0.015;
    end

    return nothing
end

"""
    update_VJR!(node::SPACMono{FT}, ratio::FT) where {FT<:AbstractFloat}
Update effective Vcmax25, Jmax25, and Rd25 from a given tune factor
"""
function update_VJR!(node::SPACMono{FT}, ratio::FT) where {FT<:AbstractFloat}
    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)
    for _iPS in node.plant_ps
        _iPS.ps.Vcmax   = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Vcmax25 = _iPS.ps.Vcmax25WW * ratio;
        _iPS.ps.Jmax    = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Jmax25  = _iPS.ps.Jmax25WW * ratio;
        _iPS.ps.Rd      = _iPS.ps.Rd25WW * ratio;
        _iPS.ps.Rd25    = _iPS.ps.Rd25WW * ratio;
    end

    return nothing
end

function update_VJR_leaf!(iPS::CanopyLayer{FT}, ratio::FT) where {FT<:AbstractFloat}
    # TODO change the ratio accordingly to photo_set
    # TODO add another ratio V2J in photo_set
    # Update Vcmax25, Jmax25 (1.67 Vcmax), and Rd25 (0.015 Vcmax)

	iPS.ps.Vcmax   = iPS.ps.Vcmax25WW * ratio;
	iPS.ps.Vcmax25 = iPS.ps.Vcmax25WW * ratio;
	iPS.ps.Jmax    = iPS.ps.Jmax25WW * ratio;
	iPS.ps.Jmax25  = iPS.ps.Jmax25WW * ratio;
	iPS.ps.Rd      = iPS.ps.Rd25WW * ratio;
	iPS.ps.Rd25    = iPS.ps.Rd25WW * ratio;

    return nothing
end

"""
    update_Cab!(node::SPACMono{FT}, cab::FT) where {FT<:AbstractFloat}
Update leaf chlorophyll content, and then rerun `fluspect!`
"""
function update_Cab!(node::SPACMono{FT}, cab::FT) where {FT<:AbstractFloat}
    for _leaf in node.leaves_rt
        _leaf.Cab = cab;
        fluspect!(_leaf, node.wl_set);
    end

    return nothing
end

"""
    update_Weibull!(node::SPACMono{FT}, b::FT, c::FT) where {FT<:AbstractFloat}
Update Weibull B and C for all components
"""
function update_Weibull!(node::SPACMono{FT}, b::FT) where {FT<:AbstractFloat}
    # 1. update B and C for roots
    for _root in node.plant_hs.roots
        _root.vc.b = b;
    end

    # 2. update B and C for trunk (if exists)
    if !(typeof(node.plant_hs) <: GrassLikeOrganism)
        node.plant_hs.trunk.vc.b = b;
    end

    # 3. update B and C for branches (if exist)
    if typeof(node.plant_hs) <: TreeLikeOrganism
        for _stem in node.plant_hs.branch
            _stem.vc.b = b;
        end
    end

    # 4. update B and C for leaves
    for _leaf in node.plant_hs.leaves
        _leaf.vc.b = b;
    end

    return nothing
end

function update_Weibull!(node::SPACMono{FT}, b::FT, c::FT) where {FT<:AbstractFloat}
    # 1. update B and C for roots
    for _root in node.plant_hs.roots
        _root.vc.b = b;
        _root.vc.c = c;
    end

    # 2. update B and C for trunk (if exists)
    if !(typeof(node.plant_hs) <: GrassLikeOrganism)
        node.plant_hs.trunk.vc.b = b;
        node.plant_hs.trunk.vc.c = c;
    end

    # 3. update B and C for branches (if exist)
    if typeof(node.plant_hs) <: TreeLikeOrganism
        for _stem in node.plant_hs.branch
            _stem.vc.b = b;
            _stem.vc.c = c;
        end
    end

    # 4. update B and C for leaves
    for _leaf in node.plant_hs.leaves
        _leaf.vc.b = b;
        _leaf.vc.c = c;
    end

    return nothing
end

function update_Kmax!(node::SPACMono{FT}, kmax::FT) where {FT<:AbstractFloat}
    # Root:Stem:Leaf conductance ratio is set as 1:2:2.
    #    For trees, roots: 2X, trunk: 8X, branch: 8X, leaves: 4X
    #    For palms, roots: 2X, trunk: 4X, leaves: 4X
    #    For grasses, roots: 2X, leaves: 2X

    # 1. update the hydraulic conductance in Roots
    for root in node.plant_hs.roots
        root.k_max      = 2 * kmax / node.plant_hs.n_root;
        root.k_element .= root.k_max * root.N;
    end

    # 2. If is a tree
    if typeof(node.plant_hs) <: TreeLikeOrganism
        node.plant_hs.trunk.k_max      = 8 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for stem in node.plant_hs.branch
            stem.k_max      = 8 * kmax / node.plant_hs.n_canopy;
            stem.k_element .= stem.k_max * stem.N;
        end
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 3. If is a palm
    if typeof(node.plant_hs) <: PalmLikeOrganism
        node.plant_hs.trunk.k_max      = 4 * kmax;
        node.plant_hs.trunk.k_element .= node.plant_hs.trunk.k_max *
                                         node.plant_hs.trunk.N;
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 4 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    # 4. If is a grass
    if typeof(node.plant_hs) <: GrassLikeOrganism
        for leaf in node.plant_hs.leaves
            leaf.k_sla      = 2 * kmax / node.la;
            leaf.k_element .= leaf.k_sla * leaf.N;
        end
    end

    return nothing
end


text:
using ASTInterpreter2
using Test
using REPL

using DebuggerFramework

@testset "Main tests" begin
    include("utils.jl")
    include("evaling.jl")
    include("stepping.jl")
    include("interpret.jl")
    include("misc.jl")
    include("toplevel.jl")
end
include("ui.jl")


text:
### Default filtering method selects all MCMC samples

select(c::MCChain, s::MCSampler) = trues(size(c.samples, 1))

### Sampler-specific selectors
# select(c::MCChain, s::ARS) = c.diagnotics["accept"]


text:
# Examples presented in class - Lecture 1
# Example 3 - Robust Knapsack Problem

using Distributions, ECOS, LaTeXStrings
# Distributions includes probability distributions,
# ECOS is a solver that we need because of the stucture of the problem
# LatexStrings is to write LaTeX in our plots

N = 18

#Input data
value = [50, 40, 70, 55, 80, 35, 65, 50, 60, 85, 20, 45, 55, 25, 80, 45, 45, 65]
weight_average = [3.5, 4, 5.5, 5, 6, 4.5, 6, 4, 5.5, 7, 5, 4.5, 7, 3.5, 5.5, 3.5, 4, 6.5]
weight_stdev = weight_average*0.3

#Randomly generating initial data using normal distribution
weight_data = Array{Float64}(100,N)
for j= 1:N
   weight_data[:,j] = rand(Normal(weight_average[j], weight_stdev[j]), 100)
end

#Input data: weight limit (capacity)
capacity = 20

# Protection elipsoid: adds 50% of weight as a protection level
P = 0.5*weight_average

# Declaring the model and solving it as a function.
function solve_robust_model(N, value, weight_average, P, Γ)
    m = Model(solver = ECOSSolver()) #creates the model, select the solver

    @variable(m, 0 <= x[1:N] <= 1 ) # creates the binary variables x, one for each item

    @constraint(m, sum(weight_average[j]*x[j] for j = 1:N) + Γ*norm(P'*x) <= capacity) # declare the knapsack constraint

    @objective(m, Max, sum(value[j]*x[j] for j = 1:N)) # declare the objective function

    solve(m)

    return getvalue(x)
end

# This function simulates the item selection againts feasibility.
function feasibility_estimate(solution, repetitions)
    feasible_count = 0
    actual_weight = zeros(18)
    for n= 1:repetitions
        for j= 1:N
           # generate random weights according to distribution
           actual_weight[j] = rand(Normal(weight_average[j], weight_stdev[j]))
        end
        #if total weight more than capacity => problem infeasible.
        if actual_weight'*solution <= capacity
            feasible_count += 1
        end
    end
    return feasible_count/repetitions
end

# generate a range of Gammas to try
Γ_range = linspace(0,1,10)

# storing the results of each run
feas = []
total_value = []
Γ_used = []

# for each Gamma, solve te model and store the results
for Γ in Γ_range
    println(Γ)
    x = solve_robust_model(N, value, weight_average, P, Γ)
    push!(Γ_used, Γ)
    push!(feas, feasibility_estimate(x,5000))
    push!(total_value, value'x)
end

# plotting the results from simulation
p1 = plot(Γ_used, feas, xlabel = L"\Gamma", ylabel = "feas. prob.", legend=false)
p2 = plot(Γ_used, total_value, xlabel = L"\Gamma", ylabel = "total value", legend=false, color=:orange)
plot(p1,p2, layout = (2,1))


text:
const BASE_PORT = 2011 # nasa_daq base_port
const MONITOR_PORT = BASE_PORT+2 # nasa_daq monitor.py port
const ZMQ_PORT = BASE_PORT+4 # pope zmq port port


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Dex"
version = v"2.23.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/dexidp/dex.git", "d820fd45d80cef74d4c65f5fcc5766ddb1fa514e"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir -p "${bindir}" "${prefix}/share"
cd dex/
atomic_patch -p1 ../patches/01-allow-optional-github-gitlab-scopes.patch
install_license LICENSE
make
mv bin/dex "$bindir/dex${exeext}"
mv bin/example-app "$bindir/example-app${exeext}"
mv bin/grpc-client "$bindir/grpc-client${exeext}"
tar -czvf $prefix/share/webtemplates.tar.gz -C ./web static templates themes
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:musl),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("example-app", :exampleapp),
    ExecutableProduct("grpc-client", :grpcclient),
    ExecutableProduct("dex", :dex),
    FileProduct("share/webtemplates.tar.gz", :webtemplates),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers = [:go, :c])