text:
"""
    LBP([mode], block_size)

Construct the classic local binary pattern methods. This includes from O-LBP (original LBP)
[1-3], mLBP (modified LBP), MB-LBP (multi-block LBP)[5-6].

# Arguments

- `mode`: a function that defines how block mode value is computed. By default it uses the
center pixel value as the block mode value, a.k.a., the original LBP [1]. Currently
supported modes are: [`center_mode`](@ref), [`average_mode`](@ref).
- `block_size::Dims`: positive integers that used to specify the size of each block [5-6].
When `all(block_size .== 1)`, it degenerates to the pixel-version [1].

# Examples

```jldoctest; setup=:(using LocalBinaryPatterns)
julia> using LocalBinaryPatterns: center_mode, average_mode

julia> X = [6 7 9; 5 6 2; 2 1 6]
3×3 $(Matrix{Int}):
 6  7  9
 5  6  2
 2  1  6

julia> local_binary_pattern(LBP(center_mode, (1, 1)), X)[2, 2] # mode value: 6
0xa9

julia> local_binary_pattern(LBP(average_mode, (1, 1)), X)[2, 2] # mode value: 44//9 ≈ 4.8
0xab
```

# References

- [1] T. Ojala, M. Pietikäinen, and D. Harwood, “A comparative study of texture measures with classification based on featured distributions,” _Pattern Recognition_, vol. 29, no. 1, pp. 51–59, Jan. 1996, doi: 10.1016/0031-3203(95)00067-4.
- [2] T. Ojala, M. Pietikäinen, and T. Mäenpää, “A Generalized Local Binary Pattern Operator for Multiresolution Gray Scale and Rotation Invariant Texture Classification,” in _Advances in Pattern Recognition — ICAPR 2001, vol. 2013, S. Singh, N. Murshed, and W. Kropatsch, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg_, 2001, pp. 399–408. doi: 10.1007/3-540-44732-6_41.
- [3] Pietikäinen, Matti, Timo Ojala, and Zelin Xu. "Rotation-invariant texture classification using feature distributions." _Pattern recognition_ 33.1 (2000): 43-52.
- [4] T. Ojala, M. Pietikainen, and T. Maenpaa, “Multiresolution gray-scale and rotation invariant texture classification with local binary patterns,” _IEEE Trans. Pattern Anal. Machine Intell._, vol. 24, no. 7, pp. 971–987, Jul. 2002, doi: 10.1109/TPAMI.2002.1017623.
- [5] Zhang, Lun, et al. "Face detection based on multi-block lbp representation." _International conference on biometrics_. Springer, Berlin, Heidelberg, 2007.
- [6] Yan, Shengye, et al. "Locally assembled binary (LAB) feature with feature-centric cascade for fast and accurate face detection." _2008 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2008.
"""
struct LBP{N,F} <: LocalBinaryPattern
    mode::F # mode(X::AbstractArray{T,N}, I::CartesianIndex{N}, offsets::Tuple)
    block_size::Dims{N}
    function LBP(mode::F, block_size::Dims{N}) where {N,F}
        all(block_size .> 0) || throw(ArgumentError("block size `block_size=$(block_size)` should be positive integers."))
        new{N, F}(mode, block_size)
    end
end
LBP(block_size::Dims) = LBP(center_mode, block_size)

Base.show(io::IO, lbp::LBP) = print(io, "LBP(", lbp.mode, ", ", lbp.block_size, ")")

Base.@propagate_inbounds function _lbp_encode(lbp::LBP, X::AbstractMatrix, I::CartesianIndex{2}, offsets::Tuple)
    rst = 0
    gc = lbp.mode(X, I, offsets)
    for i in 1:length(offsets)
        p = I.I .+ offsets[i]
        @boundscheck checkbounds(Bool, X, p...) || continue
        rst += ifelse(gc <= _inbounds_getindex(X, p), 1, 0) << (i-1)
    end
    return rst
end

Base.@propagate_inbounds function _mlbp_encode_inbounds(
        lbp::LBP,
        iX::IntegralArray,
        I::CartesianIndex{2},
        bo::CartesianIndex{2},
        offsets::Tuple)
    rst = 0
    gc = lbp.mode(iX, I, bo, offsets)
    n = prod(bo.I .+ 1)
    for i in 1:length(offsets)
        p = CartesianIndex(I.I .+ offsets[i])
        rst += ifelse(gc <= iX[p..p+bo]/n, 1, 0) << (i-1)
    end
    return rst
end
function _mlbp_encode(
        lbp::LBP,
        iX::IntegralArray,
        I::CartesianIndex{2},
        bo::CartesianIndex{2},
        offsets::Tuple)
    rst = 0
    Rlast = last(CartesianIndices(iX))
    gc = lbp.mode(iX, I, bo, offsets)
    for i in 1:length(offsets)
        p = CartesianIndex(I.I .+ offsets[i])
        checkbounds(Bool, iX, p.I...) || continue
        Rp = p:min(p+bo, Rlast)
        gp = iX[first(Rp)..last(Rp)]/length(Rp)
        rst += ifelse(gc <= gp, 1, 0) << (i-1)
    end
    return rst
end


text:
abstract type BandGapModel end

"""
Experimental Data for the Silicon Band Gap on eV and Kelvin
"""
TexpBG = [3.902699812884201,19.24084469393206,75.93156909917133,88.55564465829099,110.23790430366208,168.89779916243427,192.92167869553595,248.05845139445782,289.0171968279426,360.29047491758,411.88274080014264]
BgexpD = [1.1557025624729165,1.1555143933922498,1.1530489952245262,1.1520632962588473,1.1492423449784364,1.1405983597348612,1.1351549126101401,1.1235500821966806,1.113266923299148,1.093717856420627,1.0793094811010167].+0.013587437527083468

struct VarshniModel <: BandGapModel
    p::Vector{Number}
    f::Function
end

struct GengModel <: BandGapModel
    p::Vector{Number}
    f::Function
end

"""
# Arguments
    T:: Number Latice Temperature [K]
    f:: GengModel
Calculate the Thermal Capacity at constant pressure for the GengBandGap Model
"""
function Cp(T::Number,m::GengModel)
    Θ = m.p[5]
    if(T <= 0)
        return 4*pi^4/15
    else
        G(x) = x^4*exp(x)/(exp(x)-1)^2
        return quadgk(G,0,Θ/T)[1]
    end
end

"""
# Arguments
    p:: [Band Gap at T0[eV], T0 [K], Band Gap at T1[eV], T1[K]]
Calculation of Geng Band Gap Model for low doping Semicondurs
"""
function GengModel(p=[BgexpD[1],TexpBG[1],BgexpD[end],TexpBG[end],480])
    return GengModel(p,x-> x)
end

"""
# Arguments
    p:: [α [eV/K²], β [K] , Eg(0) [eV]]
Calculation of Band Gap Temperature Dependence Based on Varshin Model Eg(T) = Eg(0) - αT²/(T+β)
"""
VarshniModel(;α=4.73e-4*Eunit/Tunit/Tunit,β=636*Tunit,Eg0 =1.17*Eunit) = VarshniModel([Eg0,α,β],(T,p)-> p[1] - p[2]*T^2/(T+p[3]))
#VarshniModel(p = [1.16929,4.9e-4,650]) = VarshniModel(p,(T,p)-> p[1] - p[2]*T^2/(T+p[3]))
Eg(T::Number,m::BandGapModel) = m.f(T,m.p)

function IntGen(T::Number,m::GengModel)
    Θ = m.p[5]
    T0 = m.p[2]
    if T == 0
        return 0.0
    else
        X = collect(LinRange(1.0,T,450))
        Y = (X/Θ).^3 .* Cp.(X,Ref(m))
        return NumericalIntegration.integrate(X,Y)
    end
end

function Eg(T::Number,m::GengModel)
    T = ustrip(T)
    Eg0 = m.p[1]
    T0 = m.p[2]
    Eg1 = m.p[3]
    T1 = m.p[4]
    N = IntGen(T,m)
    D = IntGen(T1,m)
    Egn = (Eg0 - (Eg0-Eg1)*N/D)*Eunit
    return Egn
end

ext: jl


text:
export tJSite,
       tJSites

struct tJSite <: Site
  s::Index
  tJSite(I::Index) = new(I)
end
tJSite(n::Int) = tJSite(Index(3,"Site,tJ,n=$n"))

function tJSites(N::Int;kwargs...)::SiteSet
  sites = SiteSet(N)
  for n=1:N
    set(sites,n,tJSite(n))
  end
  return sites
end

function operator(site::tJSite,
                  opname::AbstractString)
  s = site.s
  sP = prime(site.s)
  Emp = s(1)
  EmpP = sP(1)
  Up = s(2)
  UpP = sP(2)
  Dn = s(3)
  DnP = sP(3)

  Op = ITensor(dag(s), s')

  if opname == "Nup"
    Op[Up, UpP] = 1.
  elseif opname == "Ndn"
    Op[Dn, DnP] = 1.
  elseif opname == "Ntot"
    Op[Up, UpP] = 1.
    Op[Dn, DnP] = 1.
  elseif opname == "Cup" || opname == "Aup"
    Op[Up, EmpP] = 1.
  elseif opname == "Cdagup" || opname == "Adagup"
    Op[Emp, UpP] = 1.
  elseif opname == "Cdn" || opname == "Adn"
    Op[Dn, EmpP] = 1.
  elseif opname == "Cdagdn" || opname == "Adagdn"
    Op[Emp, DnP] = 1.
  elseif opname == "FermiPhase" || opname == "FP"
    Op[Up, UpP] = -1.
    Op[Emp, EmpP] = 1.
    Op[Dn, DnP] = -1.
  elseif opname == "Fup"
    Op[Up, UpP] = -1.
    Op[Emp, EmpP] = 1.
    Op[Dn, DnP] = 1.
  elseif opname == "Fdn"
    Op[Up, UpP] = 1.
    Op[Emp, EmpP] = 1.
    Op[Dn, DnP] = -1.
  elseif opname == "Sᶻ" || opname == "Sz"
    Op[Up, UpP] = 0.5
    Op[Dn, DnP] = -0.5
  elseif opname == "Sˣ" || opname == "Sx"
    Op[Up, DnP] = 1.0
    Op[Dn, UpP] = 1.0
  elseif opname == "S⁺" || opname == "Splus"
    Op[Dn, UpP] = 1.
  elseif opname == "S⁻" || opname == "Sminus"
    Op[Up, DnP] = 1.
  elseif opname == "Emp" || opname == "0"
    pEmp = ITensor(s)
    pEmp[Emp] = 1.
    return pEmp
  elseif opname == "Up" || opname == "↑"
    pU = ITensor(s)
    pU[Up] = 1.
    return pU
  elseif opname == "Dn" || opname == "↓"
    pD = ITensor(s)
    pD[Dn] = 1.
    return pD
  else
      throw(ArgumentError("Operator name '$opname' not recognized for SpinOneSite"))
  end
  return Op
end


text:
# Everytime I encounter something that raises some error and I fix it, I should add that
# specific problem to this list to make sure, that this specific run then works without
# bugs.
using ODEFilters
using Test
using LinearAlgebra
using UnPack
using ParameterizedFunctions

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_fitzhughnagumo, prob_ode_vanstiff, prob_ode_mm_linear

@testset "Smoothing with small constant steps" begin
    prob = ODEFilters.remake_prob_with_jac(prob_ode_fitzhughnagumo)
    @test solve(prob, EK0(order=4, diffusionmodel=:fixed, smooth=true),
                adaptive=false, dt=1e-3) isa ODEFilters.ProbODESolution
    @test solve(prob, EK1(order=4, diffusionmodel=:fixed, smooth=true),
                adaptive=false, dt=1e-3) isa ODEFilters.ProbODESolution
end

@testset "Problem with analytic solution" begin
    prob = ODEFilters.remake_prob_with_jac(prob_ode_mm_linear)
    @test solve(prob, EK0(order=4)) isa ODEFilters.ProbODESolution
    @test solve(prob, EK1(order=4)) isa ODEFilters.ProbODESolution
end

@testset "Stiff Vanderpol" begin
    prob = ODEFilters.remake_prob_with_jac(prob_ode_vanstiff)
    @test solve(prob, EK1(order=3)) isa ODEFilters.ProbODESolution
end

@testset "Big Float" begin
    prob = prob_ode_fitzhughnagumo
    prob = remake(prob, u0=big.(prob.u0))
    sol = solve(prob, EK0(order=3))
    @test eltype(eltype(sol.u)) == BigFloat
    @test eltype(eltype(sol.pu.μ)) == BigFloat
    @test eltype(eltype(sol.pu.Σ)) == BigFloat
    @test sol isa ODEFilters.ProbODESolution
end

@testset "OOP problem definition" begin
    prob = ODEProblem((u, p, t) -> ([p[1] * u[1] .* (1 .- u[1])]), [1e-1], (0.0, 5), [3.0])
    @test solve(prob, EK0(order=4)) isa ODEFilters.ProbODESolution
    prob = ODEFilters.remake_prob_with_jac(prob)
    @test solve(prob, EK1(order=4)) isa ODEFilters.ProbODESolution
end

@testset "Callback: Harmonic Oscillator with condition on E=2" begin
    u0 = ones(2)
    function harmonic_oscillator(du,u,p,t)
        du[1] = u[2]
        du[2] = -u[1]
    end
    prob = ODEProblem(harmonic_oscillator, u0, (0.0,100.0))

    function Callback()
        function affect!(integ)
            @unpack dt = integ
            @unpack x_filt, Proj, Precond = integ.cache
            E0 = Proj(0)

            PI = inv(Precond(dt))
            x = x_filt

            m, P = x.μ, x.Σ

            m0, P0 = E0*m, ODEFilters.X_A_Xt(P, E0)

            e = m0'm0
            H = 2m0'E0
            S = H*P*H'

            S_inv = inv(S)
            K = P * H' * S_inv

            mnew = m + K * (2 .- e)
            Pnew = ODEFilters.X_A_Xt(P, (I-K*H)) # + X_A_Xt(R, K)

            # @info m P e S K mnew
            copy!(m, mnew)
            copy!(P, Pnew)
        end
        condtion = (t,u,integrator) -> true
        save_positions = (true, true)
        DiscreteCallback(condtion,affect!,save_positions=save_positions)
    end

    @test solve(prob, EK0(order=3)) isa ODEFilters.ProbODESolution
    @test solve(prob, EK0(order=3), callback=Callback()) isa ODEFilters.ProbODESolution
end

@testset "2nd Order ODE" begin
    function vanderpol!(ddu, du, u, p, t)
        μ = p[1]
        @. ddu = μ * ((1-u^2) * du - u)
    end
    du0 = [0.0]
    u0 = [2.0]
    tspan = (0.0, 6.3)
    p = [1e1]
    prob = SecondOrderODEProblem(vanderpol!, du0, u0, tspan, p)
    @test solve(prob, EK0(order=3)) isa ODEFilters.ProbODESolution
end

@testset "Problem definition with ParameterizedFunctions.jl" begin
    f = @ode_def LotkaVolterra begin
        dx = a*x - b*x*y
        dy = -c*y + d*x*y
    end a b c d
    p = [1.5,1,3,1]
    tspan = (0.0,10.0)
    u0 = [1.0,1.0]
    prob = ODEProblem(f,u0,tspan,p)
    @test solve(prob, EK1(order=3)) isa ODEFilters.ProbODESolution
end


text:
"""
Description of TSP, and the setup used to solve it.
"""
function tsp()
  # Compute the absolute location of the TSP example file.
  file_name = "$(dirname(@__FILE__))/tsp/berlin52.tsp"

  def = algorithm.genetic() do alg
    alg.population = population.simple() do pop
      pop.size = 100
      pop.species = species.simple() do sp
        sp.fitness = fitness.scalar(Float, false)
        sp.representation = representation.permutation(1:52)
      end
      pop.breeder = breeder.flat() do br
        br.selection = selection.tournament(2)
        br.mutation = mutation.single_swap()
        br.crossover = crossover.null()
      end
    end
    alg.evaluator = evaluator.tsp(file_name; threads = 8)
    alg.termination = [criterion.generations(1000)]
    alg.loggers = [logger.fitness()]
  end

  # Build and return the algorithm; somewhat redundant.
  algorithm.compose!(def)
end


text:
"""
    checkMethodsDefined(GenJ)

documentation
"""
function checkParametersDefined(GenJ::GeneticJuliaStructure, index::Integer = 0)

    if index == 0
        if !isdefined(GenJ._experimentInfo, :_individualType)
            error("The experiment individualType must be set")
        end

        if !isdefined(GenJ._experimentInfo, :_algorithm)
            error("The experiment algorithm must be set")
        end

        if !isdefined(GenJ, :_evaluator)
            error("The evaluator must be set")
        end

        if !isdefined(GenJ, :_generator)
            error("The generator must be set")
        end

        if !isdefined(GenJ, :_selector)
            error("The selector must be set")
        end

        if !isdefined(GenJ, :_crossoverOp)
            error("The crossover operator must be set")
        end

        if !isdefined(GenJ, :_mutationOp)
            error("The mutation operator must be set")
        end

        if !isdefined(GenJ, :_replacementOp)
            error("The replacement operator must be set")
        end
    else
        if !isdefined(GenJ._experimentInfo, :_individualType)
            error("The experiment individualType must be set " * string(index))
        end

        if !isdefined(GenJ._experimentInfo, :_algorithm)
            error("The experiment algorithm must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_evaluator)
            error("The evaluator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_generator)
            error("The generator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_selector)
            error("The selector must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_crossoverOp)
            error("The crossover operator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_mutationOp)
            error("The mutation operator must be set in experiment " * string(index))
        end

        if !isdefined(GenJ, :_replacementOp)
            error("The replacement operator must be set in experiment " * string(index))
        end
    end
end # function

function checkParametersCorrectGP(GenJ::GeneticJuliaStructure, index::Integer = 0)

    rng = Random.MersenneTwister(1234)

    ####### GA AND GP DISTINCTION #######

    if index == 0
        if !isdefined(GenJ._experimentInfo, :_GPExperimentInfo)
            error("Information about Genetic Programming problem has not been given for a GP problem ($(GenJ._experimentInfo._individualType)), this information must be specified in setGPInfo")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: CGPInfo && GenJ._experimentInfo._individualType != CGPGenotype
            error("Information about Canonical Genetic Programming problem has been given for a non CGP problem ($(GenJ._experimentInfo._individualType))")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: STGPInfo && GenJ._experimentInfo._individualType != STGPGenotype
            error("Information about Strongly Typed Genetic Programming problem has been given for a non STGP problem ($(GenJ._experimentInfo._individualType))")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: GEPInfo && GenJ._experimentInfo._individualType != GEPGenotype
            error("Information about Gene Expression Programming problem has been given for a non GEP problem ($(GenJ._experimentInfo._individualType))")
        end
    else
        if !isdefined(GenJ._experimentInfo, :_GPExperimentInfo)
            error("Information about Genetic Programming problem has not been given for a GP problem ($(GenJ._experimentInfo._individualType)), this information must be specified in setGPInfo in experiment $index")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: CGPInfo && GenJ._experimentInfo._individualType != CGPGenotype
            error("Information about Canonical Genetic Programming problem has been given for a non CGP problem ($(GenJ._experimentInfo._individualType)) in experiment $index")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: STGPInfo && GenJ._experimentInfo._individualType != STGPGenotype
            error("Information about Strongly Typed Genetic Programming problem has been given for a non STGP problem ($(GenJ._experimentInfo._individualType)) in experiment $index")
        elseif typeof(GenJ._experimentInfo._GPExperimentInfo) <: GEPInfo && GenJ._experimentInfo._individualType != GEPGenotype
            error("Information about Gene Expression Programming problem has been given for a non GEP problem ($(GenJ._experimentInfo._individualType)) in experiment $index")
        end
    end

    if GenJ._experimentInfo._individualType == GEPGenotype || GenJ._experimentInfo._individualType == CGPGenotype
        functionSet = GenJ._experimentInfo._GPExperimentInfo._functionSet
        terminalSet = GenJ._experimentInfo._GPExperimentInfo._terminalSet

        if index == 0
            type = getType(terminalSet[1])
            nTerminals = length(terminalSet)

            for i=2:nTerminals
                if !(getType(terminalSet[i]) <: type)
                    error("A terminal has type different from $type. All terminals must have the same type ($(getType(terminalSet[i])))")
                end
            end

            for func in functionSet
                if !(getType(func) <: type)
                    error("A function returns a type different from $type ($(getType(func)))")
                end

                for argType in func._argTypes
                    if !(argType <: type)
                        error("A function receives an argument of type different from $type ($argType)")
                    end
                end
            end
        else
            type = getType(terminalSet[1])
            nTerminals = length(terminalSet)

            for i=2:nTerminals
                if !(getType(terminalSet[i]) <: type)
                    error("A terminal has type different from $type. All terminals must have the same type in experiment $index ($(getType(func)))")
                end
            end

            for func in functionSet
                if !(getType(func) <: type)
                    println(func)
                    error("A function returns a type different from $type in experiment $index ($(getType(func)))")
                end

                for argType in func._argTypes
                    if !(argType <: type)
                        error("A function receives an argument of type different from $type in experiment $index ($(getType(func)))")
                    end
                end
            end
        end
    end

    OneByOne = true
    All = true
    ind = 0
#GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._generator._varArgs...)
    try
        ind = GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            OneByOne = false
        end
    end

    try
        ind = GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, 1, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            All = false
        end
    end

    if index == 0
        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified")
        end

        ind = All ? GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, 1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals")
        end
    else

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified in experiment $index")
        end

        ind = All ? GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, 1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true in experiment $index")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual in experiment $index")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false in experiment $index")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals in experiment $index")
        end
    end

    margin = 1e-6
    nParents = GenJ._crossoverOp._nParents
    nChildren = GenJ._crossoverOp._nChildren
    individuals = Array{Genotype}(undef, nParents)
    popSize = GenJ._generator._popSize
    k = 0
    fitness = 0

    nSelectedParents = getNSelectedParents(GenJ._selector)
    if !(typeof(nSelectedParents) <: Integer)
        nSelectedParents = convert(Integer, round(nSelectedParents * popSize))
        remainder = nSelectedParents % nParents

        if remainder != 0
            nSelectedParents = nSelectedParents + nParents - remainder
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            nSelectedParents -= nParents
        end
    end

    ####### GENERATION #######

    if GenJ._generator._generateOneByOne
        for i=1:nParents
            individuals[i] = GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._generator._varArgs...)
        end
    else
        individuals = GenJ._generator._method(GenJ._experimentInfo._GPExperimentInfo, nParents, rng, GenJ._generator._varArgs...)
    end

    ####### CROSSOVER, MUTATION AND EVALUATION #######

    if index == 0

        if GenJ._replacementOp._eliteSize > GenJ._generator._popSize
            error("Size of elite set is greater than the size of the population")
        end

        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType))")
        end
#fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[1], GenJ._experimentInfo._GPExperimentInfo, individuals[1])
        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], GenJ._experimentInfo._GPExperimentInfo, individuals[1])

                if !(typeof(fitness) <: Real)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e

            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents)")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement")
        end
#childCross = GenJ._crossoverOp._method(individuals..., GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._crossoverOp._varArgs...)
        try
            childCross = GenJ._crossoverOp._method(individuals..., GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren))")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType))")
            end
        end

#childMut = GenJ._mutationOp._method(individuals[1], GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._mutationOp._varArgs...)
        try

            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType))")
            end
        end

    else

        if GenJ._replacementOp._eliteSize > GenJ._generator._popSize
            error("Size of elite set is greater than the size of the population in experiment $index")
        end

        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0 in experiment $index")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType)) in experiment $index")
        end
#fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[2], GenJ._experimentInfo._GPExperimentInfo, individuals[1])
        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], GenJ._experimentInfo._GPExperimentInfo, individuals[1])

                if !(typeof(fitness) <: Real)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical in experiment $index")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents) in experiment $index")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement in experiment $index")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren)) in experiment $index")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], GenJ._experimentInfo._GPExperimentInfo, rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end
    end
end

function checkParametersCorrectGA(GenJ::GeneticJuliaStructure, index::Integer = 0)
    rng = Random.MersenneTwister(1234)
    OneByOne = true
    All = true
    ind = 0

    ####### GA AND GP DISTINCTION #######

    if index == 0
        if isdefined( GenJ._experimentInfo, :_GPExperimentInfo)
            @warn "Information about Genetic Programming problem has been given for a non-GP problem ($(GenJ._experimentInfo._individualType)), this information will be dismissed"
        end
    else
        if isdefined( GenJ._experimentInfo, :_GPExperimentInfo)
            @warn "Information about Genetic Programming problem has been given for a non-GP problem ($(GenJ._experimentInfo._individualType)), this information will be dismissed in experiment $index"
        end
    end

    try
        ind = GenJ._generator._method(rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            OneByOne = false
        end
    end

    try
        ind = GenJ._generator._method(1, rng, GenJ._generator._varArgs...)
    catch e
        if isa(e, MethodError)
            All = false
        end
    end

    if index == 0

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified")
        end

        ind = All ? GenJ._generator._method(1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals")
        end
    else

        if !(All || OneByOne)
            error("Generation method ($(GenJ._generator._method)) introduced is not compatible with arguments specified in experiment $index")
        end

        ind = All ? GenJ._generator._method(1, rng, GenJ._generator._varArgs...) : GenJ._generator._method(rng, GenJ._generator._varArgs...)

        if !OneByOne && GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) receives the population size, however, generateOneByOne is set to true in experiment $index")
        end

        if typeof(ind) <: Array && OneByOne
            error("Generation method ($(GenJ._generator._method)) returns an array of individuals while generateOneByOne is set to true, it should return a single individual in experiment $index")
        end

        if !All && !GenJ._generator._generateOneByOne
            error("Generation method ($(GenJ._generator._method)) doesn't receive the population size, however, generateOneByOne is set to false in experiment $index")
        end

        if !(typeof(ind) <: Array) && All
            error("Generation method ($(GenJ._generator._method)) returns a single individual while generateOneByOne is set to false, it should return an array of individuals in experiment $index")
        end
    end

    margin = 1e-6
    nParents = GenJ._crossoverOp._nParents
    nChildren = GenJ._crossoverOp._nChildren
    individuals = Array{Genotype}(undef, nParents)
    popSize = GenJ._generator._popSize
    fitness = 0
    k = 0

    nSelectedParents = getNSelectedParents(GenJ._selector)
    if !(typeof(nSelectedParents) <: Integer)
        nSelectedParents = convert(Integer, round(nSelectedParents * popSize))
        remainder = nSelectedParents % nParents

        if remainder != 0
            nSelectedParents = nSelectedParents + nParents - remainder
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            nSelectedParents -= nParents
        end
    end

    ####### GENERATION #######

    if GenJ._generator._generateOneByOne
        for i=1:nParents
            individuals[i] = GenJ._generator._method(rng, GenJ._generator._varArgs...)
        end

    else
        individuals = GenJ._generator._method(nParents, rng, GenJ._generator._varArgs...)
    end

    ####### CROSSOVER, MUTATION AND EVALUATION #######

    if index == 0

        if GenJ._replacementOp._eliteSize > GenJ._generator._popSize
            error("Size of elite set is greater than the size of the population")
        end

        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType))")
        end

        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], individuals[1])

                if !(typeof(fitness) <: Real)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents)")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren))")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType))")
            end
        end

        try

            childMut = GenJ._mutationOp._method(individuals[1], rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType))")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType))")
            end
        end
    else

        if GenJ._replacementOp._eliteSize > GenJ._generator._popSize
            error("Size of elite set is greater than the size of the population in experiment $index")
        end

        if GenJ._crossoverOp._probability < margin && GenJ._mutationOp._probability < margin
            error("Either crossover or mutation or both must have a probability above 0 in experiment $index")
        end

        if typeof(individuals[1]) != GenJ._experimentInfo._individualType
            error("Generation method ($(GenJ._generator._method)) should generate individuals (you generate $(typeof(individuals[1]))) of the same type as the one specified in individualType ($(GenJ._experimentInfo._individualType)) in experiment $index")
        end

        try
            len = length(GenJ._evaluator._fitnessFunctions)
            k = 1

            for i=1:len
                fitness = fitnessEval(GenJ._evaluator._fitnessFunctions[i], individuals[1])

                if !(typeof(fitness) <: Real)
                    throw(TypeError(fitness, ""))
                end
                k = k + 1
            end
        catch e
            if isa(e, MethodError)
                error("Fitness function ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) is not compatible with the individual type ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The fitness ($fitness) of an individual generated by ($(GenJ._evaluator._fitnessFunctions[k]._fitnessFunction)) must be numerical in experiment $index")
            end
        end

        if (nSelectedParents % nParents) != 0
            error("The number of selected parents ($nSelectedParents) is not divisible by number of parents per cross ($nParents) in experiment $index")
        end

        if !GenJ._selector._samplingWithRep && nSelectedParents > popSize
            error("Number of selected parents ($nSelectedParents) should not be greater than population size ($popSize) when selection method is without replacement in experiment $index")
        end

        try
            childCross = GenJ._crossoverOp._method(individuals..., rng, GenJ._crossoverOp._varArgs...)
            len = length(childCross)
            if len != GenJ._crossoverOp._nChildren
                error("Children generated by crossover ($len) must match the number of children given as argument ($(GenJ._crossoverOp._nChildren)) in experiment $index")
            end
            typeassert(childCross[1], typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The crossover method ($(GenJ._crossoverOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childCross[1]))) generated by crossover does not match the type of its parents ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end

        try
            childMut = GenJ._mutationOp._method(individuals[1], rng, GenJ._mutationOp._varArgs...)
            typeassert(childMut, typeof(individuals[1]))
        catch e
            if isa(e, MethodError)
                error("The mutation method ($(GenJ._mutationOp._method)) is not compatible with the individual type generated ($(GenJ._experimentInfo._individualType)) in experiment $index")
            elseif isa(e, TypeError)
                error("The type of the child ($(typeof(childMut))) generated by mutation does not match the type of its parent ($(GenJ._experimentInfo._individualType)) in experiment $index")
            end
        end
    end
end

"""
    checkParameters(GenJ)

Checks that everything has been set correctly before running an experiment.
"""
function checkParametersCorrect(GenJ::GeneticJuliaStructure, index::Integer = 0)

    if GenJ._experimentInfo._algorithm == SPEA && GenJ._evaluator._globalFitnessFunction._fitnessFunction == noFunc
        error("SPEA algorithm needs a global fitness function for pareto frontier setting")
    end

    if GenJ._experimentInfo._algorithm == SPEA && length(GenJ._evaluator._fitnessFunctions) == 1
        @warn "Using SPEA with only one fitness function ($(GenJ._evaluator._fitnessFunctions[1]._fitnessFunction))"
    end

    if GenJ._experimentInfo._individualType <: GPGenotype
        checkParametersCorrectGP(GenJ, index)
    else
        checkParametersCorrectGA(GenJ, index)
    end
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# This implements the contractions time estimator for the "platonic"
# operator L, the original small matrix method in Galatolo-Nisoli-Saussol
# An elementary way to rigorously estimate convergence to equilibrium and escape rates

abstract type TrinormWeakCompare end

struct WeakStrongerTrinorm <: TrinormWeakCompare end
struct WeakTrinormUncomparable <: TrinormWeakCompare end

comparetrinormweak(Bas::Basis) = WeakTrinormUncomparable
comparetrinormweak(Bas::Ulam) = WeakStrongerTrinorm

#convergencerateabstract(Bas::Basis, D::Dynamic, norms) = _convergenceratesabstract(Bas, D, norms, comparetrinormweak(Bas))

#function _convergenceratesabstract(Bas::Basis, D::Dynamic, norms, ::WeakStrongerTrinorm)
#    boundL = bound_weak_norm_abstract(Bas)
#    A, B = dfly(strong_norm(Bas), aux_norm(Bas), D)
#    m = length(coarse_norms)
#
#    Kh =  BasisDefinition.weak_projection_error(coarse_basis)
#
#    strong_norms = fill(NaN, m+1)
#    weak_norms = fill(NaN, m+1)

#    strong_norms[1] = 1.
#    weak_norms[1] = 1.

#    for i in 1:m

using LinearAlgebra
function eig_costants_small_matrix(A)
    Aint = Interval.(A)
    #@info "1, 1 entry" A[1,1]
    #@info "2, 2 entry" A[1,2]
    ρ = Aint[1, 1]+Aint[2, 2]+sqrt((Aint[1,1]-Aint[2,2])^2- 4*Aint[2,1]*Aint[1, 2])

    B = Aint'-ρ*I

    v = [B[2, 1] ; ρ-B[1,1]]

    v = v/(v[1]+v[2])

    return ρ, v
end

function convergencerateabstract(Bas::Ulam, D::Dynamic, norms)
    boundL = BasisDefinition.bound_weak_norm_abstract(Bas)
    A, B = dfly(strong_norm(Bas), aux_norm(Bas), D)
    #@info A, B
    m = length(norms)
    Kh =  BasisDefinition.weak_projection_error(Bas)

    C = (A⊕₊ 1.0)/(1.0 ⊖₋A)
    D = B ⊗₊ (A⊕₊2.0)

    small_matrices = [[A^n B; Kh*C Kh*n*D+norms[n]] for n in 1:length(norms)]

    weak_norms = ones(length(norms)) ### This is specific to Ulam
    strong_norms = [A^i+B for i in 1:length(norms)]

    for i in 1:length(small_matrices)
        ρ, v = eig_costants_small_matrix(small_matrices[i])
        if ρ<1 && ρ != ∅
            #@info ρ
            strong_norms_here = [strong_norms[1:i-1];[((1/v[1]+B/v[2])*ρ^floor(j÷i)).hi for j in i:length(norms)]]
            strong_norms = min.(strong_norms_here, strong_norms)
            weak_norms_here = [weak_norms[1: i-1]; [((B/v[2])*ρ^floor(j÷i)).hi for j in i:length(norms)]]
            weak_norms = min.(weak_norms_here, weak_norms)
        end
    end
    weak_norms = refine_norms_of_powers(weak_norms, length(weak_norms))
    strong_norms = refine_norms_of_powers(strong_norms, length(strong_norms))
    return weak_norms , strong_norms
end


text:
module TestShow

using Phylo
using DataFrames

using Compat.Test
using Compat: @info

@info """
These tests only check that the show() commands
do not give warnings or errors, not that they are correct!
"""

ntips = 10
a = IOBuffer()
@testset "Nonultrametric()" begin
    nt = rand(Nonultrametric(ntips))
    @test_nowarn show(a, nt)
    @test_nowarn show(a, first(nodeiter(nt)))
    @test_nowarn show(a, first(branchiter(nt)))
    @test_nowarn show(a, first(nodenameiter(nt)) => first(nodeiter(nt)))
    @test_nowarn show(a, first(branchnameiter(nt)) => first(nodeiter(nt)))
    bt = rand(Nonultrametric{BinaryTree{DataFrame, Vector{String}}}(ntips))
    @test_nowarn show(a, bt)
    @test_nowarn show(a, first(nodeiter(bt)))
    @test_nowarn show(a, first(branchiter(bt)))
    @test_nowarn show(a, first(nodenameiter(bt)) => first(nodeiter(bt)))
    @test_nowarn show(a, first(branchnameiter(bt)) => first(nodeiter(bt)))
    ts = rand(Nonultrametric(10), 10)
    @test_nowarn show(a, ts)
    @test_nowarn show(a, parsenewick("((,),(,,));", NamedPolytomousTree))
end
end


text:
function rewrite_backend_specific{F}(m::GLMethod, f::F, types, expr)
    if f == (erf) && length(types) == 1 && all(t-> t <: AbstractFloat, types)
        expr.args[1] = LazyMethod(m, gl_erf, types)
        return expr
    elseif f == (erfc) && length(types) == 1 && all(t-> t <: AbstractFloat, types)
        expr.args[1] = LazyMethod(m, gl_erfc, types)
        return expr
    # StaticArrays.norm is in OpenGL length
    elseif f == norm && length(types) == 1 && all(t-> is_fixedsize_array(t), types)
        expr.args[1] = LazyMethod(m, length, types)
        return expr
    elseif f in (muladd, fma) && length(types) == 3 && all(t-> t <: cli.Floats, types)
        expr.args[1] = LazyMethod(m, gl_fma, types)
        return expr
    # Constructors
    else
        false, expr
    end
end


text:
"""
    Order0()

The `Order0` method is engineered to be a more robust, though possibly
slower, alternative to the other derivative-free root-finding
methods. The implementation roughly follows the algorithm described in
*Personal Calculator Has Key to Solve Any Equation f(x) = 0*, the
SOLVE button from the
[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).
The basic idea is to use a secant step. If along the way a bracket is
found, switch to bisection, using `AlefeldPotraShi`.  If the secant
step fails to decrease the function value, a quadratic step is used up
to 4 times.

This is not really 0-order: the secant method has order
1.6...[Wikipedia](https://en.wikipedia.org/wiki/Secant_method#Comparison_with_other_root-finding_methods_
and the the bracketing method has order
1.6180...[Wikipedia](http://www.ams.org/journals/mcom/1993-61-204/S0025-5718-1993-1192965-2/S0025-5718-1993-1192965-2.pdf)
so for reasonable starting points and functions, this algorithm should be
superlinear, and relatively robust to non-reasonable starting points.

"""
struct Order0 <: AbstractSecant end

function find_zero(fs, x0, method::Order0;
                   p=nothing,
                   tracks::AbstractTracks=NullTracks(),
                   verbose=false,
                   kwargs...)
    M = Order1()
    N = AlefeldPotraShi()
    F = Callable_Function(M, fs, p)
    find_zero(F, x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
end

function Base.copy(state::UnivariateZeroState{T,S}) where {T, S}
    UnivariateZeroState(state.xn1, state.xn0, state.xstar, copy(state.m),
                        state.fxn1, state.fxn0, state.fxstar, copy(state.fm),
                        state.steps, state.fnevals,
                        state.stopped, state.x_converged,
                        state.f_converged, state.convergence_failed,
                        state.message)
end

function Base.copy!(dst::UnivariateZeroState{T,S}, src::UnivariateZeroState{T,S}) where {T,S}
    dst.xn1 = src.xn1
    dst.xn0 = src.xn0
    empty!(dst.m); append!(dst.m, src.m)
    dst.fxn1 = src.fxn1
    dst.fxn0 = src.fxn0
    empty!(dst.fm); append!(dst.fm, src.fm)
    dst.steps = src.steps
    dst.fnevals = src.fnevals
    dst.stopped = src.stopped
    dst.x_converged = src.x_converged
    dst.f_converged = src.f_converged
    dst.convergence_failed = src.convergence_failed
    dst.message = src.message
    nothing
end

# todo: consolidate this with find_zero(M,N,f, x0)...
function find_zero(fs, x0,
                   M::AbstractUnivariateZeroMethod,
                   N::AbstractBracketing;
                   p = nothing,
                   tracks::AbstractTracks=NullTracks(),
                   verbose=false,
                   kwargs...)

    F = Callable_Function(M, fs, p)
    #    state = init_state(M, F, x0)
    state = init_state(M, F, x0)
    options = init_options(M, state; kwargs...)
    l = Tracks(verbose, tracks, state)

    xstar = find_zero(M, N, F, state, options, l)
    verbose &&  show_trace(M, N, state, l)

    isnan(xstar) && throw(ConvergenceFailed("Stopped at: xn = $(state.xn1). $(state.message)"))

    return xstar

end

# Robust version using some tricks: idea from algorithm described in
# [The SOLVE button from the
# HP-34]C(http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).
# * use bracketing method if one identifed
# * limit steps so as not too far or too near the previous one
# * if not decreasing, use a quad step upto 4 times to bounce out of trap, if possible
# First uses M, then N if bracket is identified
function find_zero(M::AbstractUnivariateZeroMethod,
                   N::AbstractBracketing,
                   F,
                   state::AbstractUnivariateZeroState,
                   options::UnivariateZeroOptions,
                   l::AbstractTracks=NullTracks()
                   )

    log_step(l, M, state, :init)
    state0 = copy(state)
    quad_ctr = 0

    while true

        if assess_convergence(M, state, options)
            break
        end

        copy!(state0, state)
        update_state(M, F, state0, options) # state0 is proposed step

        adj = false
        ## did we find a zero or a bracketing interval?
        if iszero(state0.fxn1)
            copy!(state, state0)
            state.xstar, state.fxstar = state.xn1, state.fxn1
            state.f_converged = true
            break
        elseif sign(state0.fxn0) * sign(state0.fxn1) < 0
            copy!(state, state0)
            a, b = state0.xn0, state0.xn1 # could save some fn calls here
            run_bisection(N, F, (a, b), state)
            break
        end

        ## did we move too far?
        r, a, b = state0.xn1, state.xn0, state.xn1
        Δr = abs(r - b)
        Δx = abs(b - a)
        ts,TB = 1e-3, 1e2 # too small, too big
        if  Δr >= TB * Δx
            adj = true
            r = b + sign(r-b) * TB * Δx  ## too big
            state0.xn1 = r
            state0.fxn1 = F(r)
            incfn(state)
        elseif Δr  <= ts *  Δx
            adj = true
            r = b + sign(r - b) * ts * Δx
            fr = first(F(r))
            incfn(state)
            state0.xn1 = r
            state0.fxn1 = fr
        end

        # a sign change after shortening?
        if sign(state.fxn1) * sign(state0.fxn1) < 0
            state.xn0, state.fxn0 = state.xn1, state.fxn1
            state.xn1, state.fxn1 = state0.xn1, state0.fxn1
            a, b = state.xn0, state.xn1
            run_bisection(N, F, (a, b), state)
            break
        end

        ## did we improve?
        if adj || abs(state0.fxn1) < abs(state.fxn1)
            if isnan(state0.xn1) || isnan(state0.fxn1) || isinf(state0.xn1) || isinf(state0.fxn1)
                break
            end
            copy!(state, state0)
            log_step(l, M, state)
            incsteps(state)
            quad_ctr = 0
            continue
        end

        ## try quad_vertex, unless that has gotten old
        if quad_ctr > 4
            copy!(state, state0)
            state.stopped = true
            break
        else
            quad_ctr += 1
            r = quad_vertex(state0.xn1, state0.fxn1, state.xn1, state.fxn1, state.xn0, state.fxn0)

            if isnan(r) || isinf(r)
                copy!(state, state0)
            else

                fr = F(r)
                incfn(state)

                state0.xn1 = r
                state0.fxn1 = fr
                incfn(state)
                copy!(state, state0)
            end
        end

        log_step(l, M, state)
        incsteps(state)
    end

    decide_convergence(M, F, state, options)
end

# Switch to bracketing method
function run_bisection(N::AbstractBracketing, f, ab, state)
    steps, fnevals = state.steps, state.fnevals
    f = Callable_Function(N, f)
    init_state!(state, N, f; clear=true)
    find_zero(N, f, state, init_options(N, state))
    a, b = _extrema(ab)
    u,v = a > b ? (b, a) : (a, b)
    state.message *= "Bracketing used over ($u, $v), those steps not shown. "
    return nothing
end


text:
module Model

include("./grammar.jl")
include("./CausalModels.jl")
using .Grammar, .CausalModels, Random, Distributions, .SEMLang
export Node, NonTerminalNode, TerminalNode, TaggedParseTree, generateTree, getPriorLogProb, getPriorProb, getConditionalLogProb, proposeTree, getExpr, isValid, getLikelihood, markovChain

""" ----- STRUCTS ----- """

abstract type Node end

""" Non-Terminal Node Struct """
mutable struct NonTerminalNode <: Node
    tag:: Tuple  # binary tuple (i, k) from non_terminal_id (N_i) and production_rule_id (R_ik)
    symbol  # string representing grammar symbol
    node_position  # tuple representing position of node in parse tree (a_E)
    prob:: Float64  # probability associated with production_rule_id k
    parent  # parent node
    children:: Array{Node} # ordered list of child nodes
end

""" Terminal Node Struct """
mutable struct TerminalNode <: Node
    tag:: Tuple  # binary tuple (i, k) from non_terminal_id (N_i) and production_rule_id (R_ik)
    symbol  # string representing grammar symbol
    node_position # tuple representing position of node in parse tree (a_E)
    prob:: Float64 # probability associated with terminal distribution
    parent # parent node
    value # value of terminal node
end

""" Empty Node Initializers """
NonTerminalNode() = NonTerminalNode((), (), 0, 1, nothing, [])
TerminalNode() = TerminalNode((), (), 1, 0, nothing, :())

""" Tagged Parse Tree Struct """
mutable struct TaggedParseTree
    root_node::NonTerminalNode # root node of tagged parse tree
    node_positions # list of position tuples of nodes (A)
end

""" ----- METHODS ----- """
# Set all variables to nothing
function setVariablesToNothing()
    strs = []
    for index=1:5
        for prefix in ["bool_var_", "float_var_", "int_var_"]
            push!(strs, prefix)
            push!(strs, string(index))
            push!(strs, " = nothing;")
        end
    end
    expr = Meta.parse(join(strs,""))
    eval(expr)
end

# Carries out the calculations for a markov chain with unknown
# variables and bayesian-synthesis.
# X = observed data, n = number of iterations, variables = how many
# float, int, and bool variables to use
function markovChain(observed_data::NamedTuple, num_iterations::Int64=1000, variables=nothing)
    setVariablesToNothing()
    if variables === nothing
        E = generateTree()
    else
        E = generateTree(variables)
    end

    for i=1:num_iterations
        #println(join(["iteration: ", string(i)],""))
        E = generateNewExpression(E, observed_data)
    end
    if getLikelihood(E, observed_data) > 0
        E
    else
        nothing
    end
end

# Creates a new tree and evaluates the likelihood then accepts
# or returns the old tree
function generateNewExpression(E, observed_data, variables=nothing)
    if variables === nothing
        newE = proposeTree(E)
    else
        newE = proposeTree(E, variables)
    end
    lik = getLikelihood(E, observed_data)
    setVariablesToNothing()
    newLik = getLikelihood(newE, observed_data)
    setVariablesToNothing()
    println(join(["lik: ", string(lik), " newLik: ", string(newLik), " "]," "))
    if lik == 0
        p = 1
    else
        p = min(1, length(E.node_positions)*newLik/(length(newE.node_positions)*lik))
    end
    r = rand(Distributions.Uniform(0,1))
    if r<=p
        return newE
    else
        return E
    end
end

""" Compute likelihood of data given TaggedParseTree
    CURRENTLY SHOULD ONLY BE USED WITH OBSERVATION DATA THAT ARE
    DISCRETE DISTRIBUTIONS, E.G. BOOLEAN VARIABLES AND CATEGORICAL
    VARIABLES (need to add support for cateogorical variables to grammar)
"""
function getLikelihood(tree::TaggedParseTree, observed_data::NamedTuple)
    if (!isValid(tree, observed_data))
        0
    else
        treeExpr = getExpr(tree)
        eval(SEM(treeExpr))

        # compute probability of data given tree
        # deconstruct data into array of variables and corresponding array of values
        variables = []
        values = []
        for (var_name, val) in zip(keys(observed_data), observed_data)
            push!(variables, eval(var_name))
            push!(values, val)
        end
        prob(variables, values)
    end
end

""" Check if TaggedParseTree is valid w.r.t variable order/contains data variables """
function isValid(tree::TaggedParseTree, observed_data::NamedTuple)
    try
        expr = getExpr(tree)
        eval(SEM(expr))

        # extract variables from expr; check if any extracted variable
        # is nothing; if so, return false
        str_expr = repr(expr)
        for index=1:5
            for prefix in ["bool_var_", "int_var_", "float_var_"]
                var_name = join([prefix, string(index)], "")
                if (occursin(var_name, str_expr))
                    if (eval(Meta.parse(var_name)) === nothing)
                        return false
                    end
                end
            end
        end

        for name in keys(observed_data)
            #println(name)
            if !(occursin(string(name), str_expr))
                return false
            end
        end
        #println("here!")
        true
    catch
        false
    end
end

""" Recursively construct random TaggedParseTree """
function generateTree(rng::Int64)
    Random.seed!(rng)
    # initialize tree object and root node object
    node_positions = []
    root = generateTreeHelper(1, "expr", nothing, node_positions)
    tree = TaggedParseTree(root, node_positions)
    tree
end

generateTree() = generateTree(Random.GLOBAL_RNG)

# generateTree(rng::Int64) = generateTree(Random.AbstractRNG(rng))

function generateTree(variables::NamedTuple)
    # do some stuff
end

""" Recursive helper function for constructing random TaggedParseTree """
function generateTreeHelper(node_index, symbol, parent_node, node_positions)
    node = isTerminalIndex(node_index) ? TerminalNode() : NonTerminalNode()
    node.parent = parent_node

    # sample production rule
    production_rule_index, prob = getProductionRuleIndexAndProb(node_index)
    node.tag = (node_index, production_rule_index)
    node.symbol = symbol
    node.prob = prob

    if isTerminalIndex(node_index)
        # set terminal value
        node.value = getTerminalValue(node_index, production_rule_index)

        # construct new node_position
        node_position = tuplejoin(parent_node.node_position, node.tag)
        node.node_position = node_position
        push!(node_positions, node_position)
    else
        # construct new node_position
        if (parent_node === nothing)
            node_position = node.tag
        else
            node_position = tuplejoin(parent_node.node_position, node.tag)
        end
        node.node_position = node_position
        push!(node_positions, node_position)

        # construct child nodes
        child_node_symbols = getSymbolsOfProductionRule(node_index, production_rule_index)
        for child_index in 1:length(child_node_symbols)
            child_node_symbol = child_node_symbols[child_index]
            child_node_index = getNodeIndexFromSymbol(child_node_symbol)
            child_node = generateTreeHelper(child_node_index, child_node_symbol, node, node_positions)
            push!(node.children, child_node)
        end
    end
    node
end

""" Compute log probability of sampling TaggedParseTree """
function getPriorLogProb(tree::TaggedParseTree)
    getPriorLogProbHelper(tree.root_node)
end

function getPriorLogProbHelper(node::Node)
    node_index = getNodeIndexFromSymbol(node.symbol)
    log_prob = log(node.prob)
    if (!isTerminalIndex(node_index))
        log_prob = log(node.prob)
        for child in node.children
            log_prob += getPriorLogProbHelper(child)
        end
    end
    log_prob
end

""" Compute probability of sampling TaggedParseTree """
function getPriorProb(tree::TaggedParseTree)
    getPriorProbHelper(tree.root_node)
end

function getPriorProbHelper(node::Node)
    node_index = getNodeIndexFromSymbol(node.symbol)
    prob = node.prob
    if (!isTerminalIndex(node_index))
        for child in node.children
            prob *= getPriorProbHelper(child)
        end
    end
    prob
end

"""Compute conditional probability of transitioning between TaggedParseTrees """
function getConditionalLogProb(tree1::TaggedParseTree, tree2::TaggedParseTree)
    node_positions_1 = tree1.node_positions
    node_positions_2 = tree2.node_positions

    #find the node_positions only in tree2
    tree2_only = setdiff(node_positions_2, node_positions_1)

    sum = 0.0
    #Loop through the node positions in tree2 and calculate the probability of
    #getting the resulting tree from the current node if it is an ancestor of
    #the node that is changed
    for pos in node_positions_2
        if length(pos)<=length(tree2_only[1])
            if pos == tree2_only[1][1:length(pos)]
                node = findNodeWithPosition(tree2.root_node, pos)
                sum += getPriorProbHelper(node)
            end
        end
    end
    return log(sum/length(node_positions_1))
end

function proposeTree(tree::TaggedParseTree)
    copied_tree = deepcopy(tree)
    node_positions = copied_tree.node_positions
    node_position = rand(node_positions)

    if (length(node_position) == 2)
        generateTree()
    else
        # find node with uniformly randomly chosen node position
        node = findNodeWithPosition(copied_tree.root_node, node_position)
        #find the current position of the node so the new node can be inserted at the same location
        node_index = findfirst(map(x->x == node, node.parent.children))

        # remove that node from node.parent.children
        if (node.parent !== nothing)
            filter!(child_node -> child_node.node_position != node_position, node.parent.children)
        end

        filter!(pos -> !(length(pos) >= length(node_position) && pos[1:length(node_position)] == node_position), node_positions)
        # generate new node at node_position
        new_node = generateTreeHelper(node.tag[1], node.symbol, node.parent, node_positions)

        # add new_node to new_node.parent.children
        if (node.parent !== nothing)
            if node_index>1
                store = node.parent.children[1:node_index-1]
            else
                store = []
            end
            append!(append!(store, [new_node]), node.parent.children[node_index:length(node.parent.children)])
            node.parent.children = store

        end
        copied_tree
    end
end

function findNodeWithPosition(node, pos)
    if (node.node_position == pos)
        return node
    else
        if typeof(node) == NonTerminalNode
            for child_node in node.children
                ret_val = findNodeWithPosition(child_node, pos)
                if ret_val !== nothing
                    return ret_val
                end
            end
            return nothing
        else
            return nothing
        end
    end
end

""" Extract native Julia Expr object from TaggedParseTree """
function getExpr(tree::TaggedParseTree)
    expr_tokens = []
    getExprHelper(tree.root_node, expr_tokens)
    lines = split(join(expr_tokens), "\n")
    exprs = map(Meta.parse, lines)
    Expr(:block, exprs...)
end

function getExprHelper(node::TerminalNode, arr)
    str_value = repr(node.value)
    if (str_value[1] == ':')
        push!(arr, str_value[2:length(str_value)])
    else
        push!(arr, str_value)
    end
end

function getExprHelper(node::NonTerminalNode, arr)
    node_index = getNodeIndexFromSymbol(node.symbol)
    rule_index = node.tag[2]
    children = node.children

    if (node_index == 1) # line \n expr | line
        if (rule_index == 1)
            getExprHelper(children[1], arr)
            push!(arr, "\n")
            getExprHelper(children[2], arr)
        elseif (rule_index == 2)
            getExprHelper(children[1], arr)
        end
    elseif (node_index == 2) # exo_line | endo_line
        getExprHelper(children[1], arr)
    elseif (node_index == 3) # bool_var '~' dist | float_var '~' float_dist
        getExprHelper(children[1], arr)
        push!(arr, " ~ ")
        getExprHelper(children[2], arr)
    elseif (node_index == 4) # unary_endo_expr | binary_endo_expr | ternary_endo_expr
        getExprHelper(children[1], arr)
    elseif (node_index == 5) # '!' bool_var | var
        if (rule_index == 1)
            getExprHelper(children[1], arr)
            push!(arr, " = !")
            getExprHelper(children[2], arr)
        elseif (rule_index == 2 || rule_index == 3 || rule_index == 4)
            getExprHelper(children[1], arr)
            push!(arr, " = ")
            getExprHelper(children[2], arr)
        end
    elseif (node_index == 6) # bool_binary_op bool_var bool_var_list | num_binary_op int_var int_var_list | num_binary_op float_var float_var_list
        getExprHelper(children[1], arr)
        push!(arr, " = ")
        getExprHelper(children[2], arr)
        push!(arr, "(")
        getExprHelper(children[3], arr)
        push!(arr, ", ")
        getExprHelper(children[4], arr)
        push!(arr, ")")
    elseif (node_index == 7) # bool_var '?' var ':' var
        getExprHelper(children[1], arr)
        push!(arr, " = ")
        getExprHelper(children[2], arr)
        push!(arr, " ? ")
        getExprHelper(children[3], arr)
        push!(arr, " : ")
        getExprHelper(children[4], arr)
    elseif (node_index == 8 || node_index == 9 || node_index == 10) # bool_var ',' bool_var_list | bool_var
        if (rule_index == 1)
            getExprHelper(children[1], arr)
            push!(arr, ", ")
            getExprHelper(children[2], arr)
        elseif (rule_index == 2)
            getExprHelper(children[1], arr)
        end
    elseif (node_index == 11) # bool_var | int_var | float_var
        getExprHelper(children[1], arr)
    elseif (node_index == 12) # Bernoulli( bernoulli_params )
        push!(arr, "Bernoulli(")
        getExprHelper(children[1], arr)
        push!(arr, ")")
    elseif (node_index == 21)  # Normal ( normal_params ) | Uniform( uniform_params )
        if (rule_index == 1)
            push!(arr, "Normal")
            getExprHelper(children[1], arr)
        elseif (rule_index == 2)
            push!(arr, "Uniform")
            getExprHelper(children[1], arr)
        end
    else
        throw(ArgumentError(""))
    end
end

end


text:
using Turing.VarReplay
using Turing.VarReplay: set_flag!, unset_flag!
using Distributions
using Test

vi = VarInfo()
vn_x = VarName(gensym(), :x, "", 1)
dist = Normal(0, 1)
r = rand(dist)
gid = 0

push!(vi, vn_x, r, dist, gid)

# del is set by default
@test is_flagged(vi, vn_x, "del") == false

set_flag!(vi, vn_x, "del")
@test is_flagged(vi, vn_x, "del") == true

unset_flag!(vi, vn_x, "del")
@test is_flagged(vi, vn_x, "del") == false


text:
# Autogenerated wrapper script for AMReX_jll for x86_64-linux-gnu-cxx11
export libamrex

using CompilerSupportLibraries_jll
using MPICH_jll
JLLWrappers.@generate_wrapper_header("AMReX")
JLLWrappers.@declare_library_product(libamrex, "libamrex.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, MPICH_jll)
    JLLWrappers.@init_library_product(
        libamrex,
        "lib/libamrex.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text: include("histograms.jl")


text:
export SVMClassifier, SVMRegressor,
       SVMNuClassifier, SVMNuRegressor,
       SVMLClassifier, SVMLRegressor

"""
    SVMClassifier(; kwargs...)
C-Support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html#sklearn.svm.SVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMNuClassifier, SVMLClassifier, SVMRegressor
"""
mutable struct SVMClassifier <: MMI.Deterministic
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    decision_function_shape::String
    random_state
end

# constructor:
#> all arguments are kwargs with a default value
function SVMClassifier(
    ;C          = 1.0
    ,kernel     = "rbf"
    ,degree     = 3
    ,gamma      = "auto"
    ,coef0      = 0.0
    ,shrinking  = true
    ,tol        = 1e-3
    ,cache_size = 200
    ,max_iter   = -1
    ,decision_function_shape="ovr"
    ,random_state = nothing)

    model = SVMClassifier(
        C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , decision_function_shape
        , random_state
        )

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMClassifier)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
        model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
        model.gamma="auto"
    end
    if(!(model.decision_function_shape in ["ovo","ovr"]))
            warning *="decision_function_shape parameter is not valid, setting to default=\"ovr\" \n"
        model.decision_function_shape="ovr"
    end
    return warning
end

"""
    SVMNuClassifier(; kwargs...)
NU-Support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVC.html#sklearn.svm.NuSVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVC.html#sklearn.svm.NuSVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMLClassifier, SVMNuRegressor
"""
mutable struct SVMNuClassifier <: MMI.Deterministic
    nu::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    decision_function_shape::String
    random_state
end

# constructor:
#> all arguments are kwargs with a default value
function SVMNuClassifier(
    ;nu         = 0.5
    ,kernel     = "rbf"
    ,degree     = 3
    ,gamma      = "auto"
    ,coef0      = 0.0
    ,shrinking  = true
    ,tol        = 1e-3
    ,cache_size = 200
    ,max_iter   = -1
    ,decision_function_shape = "ovr"
    ,random_state = nothing)

    model = SVMNuClassifier(
        nu
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , decision_function_shape
        , random_state
        )

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMNuClassifier)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
        model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
        model.gamma="auto"
    end
    if(!(model.decision_function_shape in ["ovo","ovr"]))
            warning *="decision_function_shape parameter is not valid, setting to default=\"ovr\" \n"
        model.decision_function_shape="ovr"
    end
    return warning
end

"""
    SVMLClassifier(; kwargs...)
Linear-support Vector classifier from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVC.html#sklearn.svm.LinearSVC). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMNuClassifier, SVMLRegressor
"""

mutable struct SVMLClassifier <: MMI.Deterministic
    C::Float64
    loss::String
    dual::Bool
    penalty::String
    tol::Float64
    max_iter::Int
    intercept_scaling::Float64
    random_state
end

# constructor:
#> all arguments are kwargs with a default value

function SVMLClassifier(
    ;C        = 1.0
    ,loss     = "squared_hinge"
    ,dual     = true
    ,penalty  = "l2"
    ,tol      = 1e-3
    ,max_iter = -1
    ,intercept_scaling = 1.0
    ,random_state = nothing)

    model = SVMLClassifier(
        C
        , loss
        , dual
        , penalty
        , tol
        , max_iter
        , intercept_scaling
        , random_state
        )

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMLClassifier)
    warning = ""
    if(!(model.loss in ["hinge","squared_hinge"]))
            warning *="loss parameter is not valid, setting to default=\"squared_hinge\" \n"
        model.loss="squared_hinge"
    end
    if(!(model.penalty in ["l1","l2"]))
            warning *="penalty parameter is not valid, setting to default=\"l2\" \n"
        model.penalty="l2"
    end
    return warning
end

"""
    SVMRegressor(; kwargs...)
Epsilon-Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html#sklearn.svm.SVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMClassifier, SVMNuRegressor, SVMLRegressor
"""
mutable struct SVMRegressor <: MMI.Deterministic
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
    epsilon::Float64
end

# constructor:
#> all arguments are kwargs with a default value
function SVMRegressor(
    ;C          = 1.0
    ,kernel     = "rbf"
    ,degree     = 3
    ,gamma      = "auto"
    ,coef0      = 0.0
    ,shrinking  = true
    ,tol        = 1e-3
    ,cache_size = 200
    ,max_iter   = -1
    ,epsilon    = 0.1)

    model = SVMRegressor(
        C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter
        , epsilon)

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMRegressor)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
        model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
        model.gamma="auto"
    end
    return warning
end

"""
    SVMNuRegressor(; kwargs...)
Nu Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVR.html#sklearn.svm.NuSVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.NuSVR.html#sklearn.svm.NuSVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMNuClassifier, SVMRegressor, SVMLRegressor
"""

mutable struct SVMNuRegressor <: MMI.Deterministic
    nu::Float64
    C::Float64
    kernel::Union{String,Function}
    degree::Int
    gamma::Union{Float64,String}
    coef0::Float64
    shrinking::Bool
    tol::Float64
    cache_size::Float64
    max_iter::Int
end

# constructor:
#> all arguments are kwargs with a default value
function SVMNuRegressor(
    ;nu         = 0.5
    ,C          = 1.0
    ,kernel     = "rbf"
    ,degree     = 3
    ,gamma      = "auto"
    ,coef0      = 0.0
    ,shrinking  = true
    ,tol        = 1e-3
    ,cache_size = 200
    ,max_iter   = -1)

    model = SVMNuRegressor(
        nu
        , C
        , kernel
        , degree
        , gamma
        , coef0
        , shrinking
        , tol
        , cache_size
        , max_iter)

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMNuRegressor)
    warning = ""
    if(typeof(model.kernel)==String && (!(model.kernel  in
            ["linear","poly","rbf","sigmoid","precomputed"])))
            warning *="kernel parameter is not valid, setting to default=\"rbf\" \n"
        model.kernel="rbf"
    end
    if(typeof(model.gamma)==String && (!(model.gamma  in
            ["auto","scale"])))
            warning *="gamma parameter is not valid, setting to default=\"auto\" \n"
        model.gamma="auto"
    end
    return warning
end

"""
    SVMLRegressor(; kwargs...)
Linear Support Vector Regression from
[https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR](https://scikit-learn.org/stable/modules/generated/sklearn.svm.LinearSVR.html#sklearn.svm.LinearSVR). Implemented hyperparameters as per
package documentation cited above.
See also, SVMRegressor, SVMNuRegressor, SVMLClassifier
"""

mutable struct SVMLRegressor <: MMI.Deterministic
    C::Float64
    loss::String
    fit_intercept::Bool
    dual::Bool
    tol::Float64
    max_iter::Int
    epsilon::Float64
end

# constructor:
#> all arguments are kwargs with a default value
function SVMLRegressor(
    ;C        = 1.0
    ,loss     = "epsilon_insensitive"
    ,fit_intercept = true
    ,dual     = true
    ,tol      = 1e-3
    ,max_iter = -1
    ,epsilon  = 0.1)

    model = SVMLRegressor(
        C
        , loss
        , fit_intercept
        , dual
        , tol
        , max_iter
        , epsilon)

    message = MMI.clean!(model)       #> future proof by including these
    isempty(message) || @warn message #> two lines even if no clean! defined below

    return model
end

function MMI.clean!(model::SVMLRegressor)
    warning = ""
    if(!(model.loss in ["epsilon_insensitive","squared_epsilon_insensitive"]))
            warning *="loss parameter is not valid, setting to default=\"epsilon_insensitive\" \n"
        model.loss="epsilon_insensitive"
    end
    return warning
end

function MMI.fit(model::SVMClassifier
                     , verbosity::Int   #> must be here (and typed) even if not used (as here)
                     , X
                     , y)

    Xmatrix = MMI.matrix(X)
    y_plain = MMI.int(y)
    decode  = MMI.decoder(y[1]) # for predict method

    cache = SVC(C=model.C,
                kernel=model.kernel,
                degree=model.degree,
                coef0=model.coef0,
                shrinking=model.shrinking,
                gamma=model.gamma,
                tol=model.tol,
                cache_size=model.cache_size,
                max_iter=model.max_iter,
                decision_function_shape=model.decision_function_shape,
                random_state=model.random_state
                )

    result    = ScikitLearn.fit!(cache, Xmatrix, y_plain)
    fitresult = (result, decode)
    report    = NamedTuple()

    return fitresult, nothing, report
end

function MMI.fit(model::SVMNuClassifier
                     , verbosity::Int   #> must be here (and typed) even if not used (as here)
                     , X
                     , y)

    Xmatrix = MMI.matrix(X)
    y_plain = MMI.int(y)
    decode  = MMI.decoder(y[1]) # for predict method

    cache = NuSVC(nu=model.nu,
                  kernel=model.kernel,
                  degree=model.degree,
                  coef0=model.coef0,
                  shrinking=model.shrinking,
                  gamma=model.gamma,
                  tol=model.tol,
                  cache_size=model.cache_size,
                  max_iter=model.max_iter,
                  decision_function_shape=model.decision_function_shape,
                  random_state=model.random_state
                  )

    result    = ScikitLearn.fit!(cache,Xmatrix,y_plain)
    fitresult = (result, decode)
    report    = NamedTuple()

    return fitresult, nothing, report
end

function MMI.fit(model::SVMLClassifier
                     , verbosity::Int   #> must be here (and typed) even if not used (as here)
                     , X
                     , y)

    Xmatrix = MMI.matrix(X)
    y_plain = MMI.int(y)
    decode  = MMI.decoder(y[1]) # for predict method

    cache = LinearSVC(C=model.C,
                      loss = model.loss,
                      dual=model.dual,
                      penalty=model.penalty,
                      intercept_scaling=model.intercept_scaling,
                      tol=model.tol,
                      max_iter=model.max_iter,
                      random_state=model.random_state
                      )

    result    = ScikitLearn.fit!(cache,Xmatrix,y_plain)
    fitresult = (result, decode)
    report    = NamedTuple()

    return fitresult, nothing, report
end

function MMI.fit(model::SVMRegressor
             , verbosity::Int   #> must be here (and typed) even if not used (as here)
             , X
             , y)

    Xmatrix = MMI.matrix(X)

    cache = SVR(C=model.C,
                kernel=model.kernel,
                degree=model.degree,
                coef0=model.coef0,
                shrinking=model.shrinking,
                gamma=model.gamma,
                tol=model.tol,
                cache_size=model.cache_size,
                max_iter=model.max_iter,
                epsilon=model.epsilon)

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report    = NamedTuple()

    return fitresult, nothing, report
end

function MMI.fit(model::SVMNuRegressor
                     , verbosity::Int   #> must be here (and typed) even if not used (as here)
                     , X
                     , y)

    Xmatrix = MMI.matrix(X)

    cache = NuSVR(nu=model.nu,
                  C=model.C,
                  kernel=model.kernel,
                  degree=model.degree,
                  coef0=model.coef0,
                  shrinking=model.shrinking,
                  gamma=model.gamma,
                  tol=model.tol,
                  cache_size=model.cache_size,
                  max_iter=model.max_iter
                  )

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report    = NamedTuple()

    return fitresult, nothing, report
end

function MMI.fit(model::SVMLRegressor
                     , verbosity::Int   #> must be here (and typed) even if not used (as here)
                     , X
                     , y)

    Xmatrix = MMI.matrix(X)

    cache = LinearSVR(C=model.C,
                      loss=model.loss,
                      fit_intercept=model.fit_intercept,
                      dual=model.dual,
                      tol=model.tol,
                      max_iter=model.max_iter,
                      epsilon=model.epsilon)

    fitresult = ScikitLearn.fit!(cache,Xmatrix,y)
    report    = NamedTuple()

    return fitresult, nothing, report
end

# placeholder types for predict dispatching
SVMC = Union{SVMClassifier, SVMNuClassifier, SVMLClassifier}
SVMR = Union{SVMRegressor, SVMNuRegressor, SVMLRegressor}
SVM  = Union{SVMC, SVMR}

function MMI.predict(model::SVMC
                         , (fitresult, decode)
                         , Xnew)

    xnew       = MMI.matrix(Xnew)
    prediction = ScikitLearn.predict(fitresult, xnew)
    return decode(prediction)
end

function MMI.predict(model::SVMR
                         , fitresult
                         , Xnew)

    xnew       = MMI.matrix(Xnew)
    prediction = ScikitLearn.predict(fitresult,xnew)
    return prediction
end

## METADATA

MMI.load_path(::Type{<:SVMClassifier})   = "MLJModels.ScikitLearn_.SVMClassifier"
MMI.load_path(::Type{<:SVMNuClassifier}) = "MLJModels.ScikitLearn_.SVMNuClassifier"
MMI.load_path(::Type{<:SVMLClassifier})  = "MLJModels.ScikitLearn_.SVMLClassifier"
MMI.load_path(::Type{<:SVMRegressor})    = "MLJModels.ScikitLearn_.SVMRegressor"
MMI.load_path(::Type{<:SVMNuRegressor})  = "MLJModels.ScikitLearn_.SVMNuRegressor"
MMI.load_path(::Type{<:SVMRegressor})    = "MLJModels.ScikitLearn_.SVMRegressor"
MMI.load_path(::Type{<:SVMLRegressor})   = "MLJModels.ScikitLearn_.SVMLRegressor"

MMI.package_name(::Type{<:SVM})    = "ScikitLearn"
MMI.package_uuid(::Type{<:SVM})    = "3646fa90-6ef7-5e7e-9f22-8aca16db6324"
MMI.is_pure_julia(::Type{<:SVM})   = false
MMI.package_url(::Type{<:SVM})     = "https://github.com/cstjean/ScikitLearn.jl"
MMI.input_scitype(::Type{<:SVM})   = Table(Continuous)
MMI.target_scitype(::Type{<:SVMC}) = AbstractVector{<:Finite}
MMI.target_scitype(::Type{<:SVMR}) = AbstractVector{Continuous}


text:
# =========================================================================================
"""
    readmsms{T <: AbstractFloat}(
        vertstream::IOStream,
        facestream::IOStream,
                  ::Type{T}=Float64
    )

Reads a surface model from the given MSMS-generated `.face` and `.vert` files.

!!! note
    This file type does not support charge models! Hence, the charge list of the returning
    `Model` object is empty and has to be set separately.

# Return type
`Model{T, Triangle{T}}`

# Alias

    readmsms{T}(fname::String, ::Type{T}=Float64)

Reads the model using a common file name prefix (`fname.{vert,face}`) for both files
rather than `IOStream` objects.

"""
function readmsms(
        vertstream::IOStream,
        facestream::IOStream,
                  ::Type{T}=Float64
    ) where T <: AbstractFloat
    nodes = readmsms_nodes(vertstream, T)
    Model(nodes, readmsms_elements(facestream, nodes))
end

function readmsms(
        fname::String,
             ::Type{T}=Float64
    ) where T <: AbstractFloat
    open(ff -> open(fv -> readmsms(fv, ff, T), "$fname.vert"), "$fname.face")
end

# =========================================================================================
"""
    readmsms_nodes{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads all nodes from the given MSMS-generated `.vert` file.

# Return type
`Vector{Vector{T}}`
"""
function readmsms_nodes(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    nodes = Vector{T}[]
    # skip header lines
    [readline(stream) for i in 1:3]
    for line in eachline(stream)
        push!(nodes, [parse(T, e) for e in split(line)[1:3]])
    end
    nodes
end

# =========================================================================================
"""
    readmsms_elements{T <: AbstractFloat}(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    )

Reads all elements from the given MSMS-generated `.face` file.

# Return type
`Vector{Triangle{T}}`
"""
function readmsms_elements(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    ) where T <: AbstractFloat
    elements = Triangle{T}[]
    # skip header lines
    [readline(stream) for i in 1:3]
    for line in eachline(stream)
        push!(elements, Triangle([nodes[parse(Int, e)] for e in split(line)[1:3]]...))
    end
    elements
end


text:
@symbol_func function f_RP(cur_reactor::AbstractReactor)
  cur_f = Q_E(cur_reactor)

  cur_f += 1

  cur_f ^= -1

  cur_f
end


text:
@diffusion_process Sine begin
    :parameters
    (a, b, c, σ) --> Float64

    :additional
    constdiff --> true
end

b(t, x, P::Sine) = ℝ{1}(P.a + P.b*sin.(P.c * x))

σ(t, x, P::Sine) = ℝ{1}(P.σ)