text:
@testset "_is_strictly_increasing" begin
    @test TimeDag._is_strictly_increasing([])
    @test TimeDag._is_strictly_increasing([1])
    @test TimeDag._is_strictly_increasing([1, 2, 3])
    @test !TimeDag._is_strictly_increasing([1, 2, 3, 3])
    @test !TimeDag._is_strictly_increasing([1, 2, 2, 3])
    @test !TimeDag._is_strictly_increasing([1, 3, 2])
end

@testset "basic" begin
    times = DateTime(2020, 1, 1):Day(1):DateTime(2020, 1, 5)
    values = 1:length(times)
    block = Block(times, values)

    @test value_type(block) == Int64
    @test length(block) == 5
    @test first(block) == (DateTime(2020, 1, 1), 1)
    @test last(block) == (DateTime(2020, 1, 5), 5)
    @test collect(block) == collect(zip(times, values))

    @testset "getindex" begin
        @test block[1] == first(block)
        @test block[2] == (DateTime(2020, 1, 2), 2)
        @test block[5] == last(block)
        @test block[1:1] == Block([times[1]], [values[1]])
        @test block[1:2] == Block(times[1:2], values[1:2])
        @test block[1:2:5] == Block(times[1:2:5], values[1:2:5])
    end

    @testset "slice" begin
        # Slicing follows open-closed interval semantics.
        block_slice = TimeDag._slice(block, DateTime(2020, 1, 1), DateTime(2020, 1, 3))
        @test value_type(block_slice) == Int64
        @test length(block_slice) == 2
        @test first(block_slice) == (DateTime(2020, 1, 1), 1)
        @test last(block_slice) == (DateTime(2020, 1, 2), 2)
        @test block_slice[1] == first(block_slice)
        @test block_slice[2] == last(block_slice)
        @test collect(block_slice) == collect(zip(times[1:2], values[1:2]))
    end
end

@testset "invalid" begin
    # Mismatched lengths.
    @test_throws ArgumentError Block([DateTime(1990)], Int64[])
    @test_throws ArgumentError Block(DateTime[], [1])

    # Mis-ordered timesdtamps.
    @test_throws ArgumentError Block([DateTime(1990), DateTime(1990)], [1, 2])
    @test_throws ArgumentError Block([DateTime(1991), DateTime(1990)], [1, 2])

    # Creating an invalid block when unchecked should not throw an exception.
    block = Block(TimeDag.unchecked, [DateTime(1990)], Int64[])
    @test block.times == [DateTime(1990)]
    @test block.values == Int64[]
end

@testset "tables" begin
    @test Tables.columnnames(b1) == (:time, :value)

    for (i, row) in enumerate(Tables.namedtupleiterator(b1))
        @test b1.times[i] == row.time
        @test b1.values[i] == row.value
    end

    df1 = DataFrame(b1)
    @test Tables.columnnames(df1) == [:time, :value]
    @test df1[!, :time] == b1.times
    @test df1[!, :value] == b1.values
end


text:
module Select

using FieldMetadata

import FieldMetadata: flattenable, @flattenable, selectable, @selectable, default, label

export select, updateselected, selectable, @selectable

nested(T::Type, expr_builder, expr_combiner) =
    expr_combiner(T, [Expr(:..., expr_builder(T, fn)) for fn in fieldnames(T)])

select_combiner(T, expressions) = Expr(:tuple, expressions...)
select_builder(T, fname) = quote
    selection = selectable($T, Val{$(QuoteNode(fname))})
    if selection != Nothing
        ((selection, typeof(t.$fname).name.wrapper, label($T, Val{$(QuoteNode(fname))})), )
    elseif flattenable($T, Val{$(QuoteNode(fname))})
        select(getfield(t, $(QuoteNode(fname))))
    else
        ()
    end
end
select_inner(T) = nested(T, select_builder, select_combiner)

@generated select(t) = select_inner(t)

updateselected_builder(T, fname) = quote
    selection = selectable($T, Val{$(QuoteNode(fname))})
    if selection != Nothing
        typ = data[n]
        n += 1
        # Update only if the default constructor returns a differrent type
        if typeof(typ()) != fieldtype($T, $(QuoteNode(fname)))
            (typ(),)
        else
            (getfield(t, $(QuoteNode(fname))),)
        end
    elseif flattenable($T, Val{$(QuoteNode(fname))})
        val, n = updateselected(getfield(t, $(QuoteNode(fname))), data, n)
        val
    else
        (getfield(t, $(QuoteNode(fname))),)
    end
end

updateselected_combiner(T, expressions) = :(($(Expr(:call, :($T.name.wrapper), expressions...)),), n)
updateselected_combiner(T::Type{<:Tuple}, expressions) = :(($(Expr(:tuple, expressions...)),), n)

updateselected_inner(::Type{T}) where T = nested(T, updateselected_builder, updateselected_combiner)

updateselected(x::Number, data, n) = x, n
updateselected(x::Symbol, data, n) = x, n
updateselected(x::Nothing, data, n) = x, n
" Reconstruct an object from partial Tuple or Vector data and another object"
updateselected(t, data) = updateselected(t, data, 1)[1][1]
# Also increment vector position counter - the returned n + 1 becomes the new n
@generated updateselected(t, data, n) = updateselected_inner(t)

end # module

ext: jl


text:
"""
	MSMOptions

MSMOptions is a mutable struct that contains options related to the optimization.
# Examples
```julia-repl
julia> options = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes, localOptimizer = :LBFGS)
```
"""
mutable struct MSMOptions
	globalOptimizer::Symbol #algorithm for finding a global maximum
	localOptimizer::Symbol 	#algorithm for finding a local maximum
	maxFuncEvals::Int64			#maximum number of evaluations
	saveName::String				#name under which the optimization should be saved
	showDistance::Bool			#show the distance, everytime the objective function is calculated?
	minBox::Bool						#When looking for a local maximum, use Fminbox ?
	populationSize::Int64		#When using BlackBoxOptim, set the population size
	penaltyValue::Float64 	#Objective function's value when the model fails
	gridType::Symbol				#sampling procedure to use (latin hypercube by default)
	saveStartingValues::Bool #whether or not saving the starting values when using local_to_global
	maxTrialsStartingValues::Int64 #maximum number of attempts when searching for valid starting values
	thresholdStartingValue::Float64 #value under which a point is considered as a valid starting value
end

"""
	MSMOptions

Constructor for MSMOptions. MSMOptions is a mutable struct that contains options related to the optimization.

# Examples
```julia-repl
julia> options = MSMOptions(maxFuncEvals=1000, globalOptimizer = :dxnes, localOptimizer = :LBFGS)
```
"""
function MSMOptions( ;
					globalOptimizer::Symbol=:dxnes,
					localOptimizer::Symbol=:LBFGS,
					maxFuncEvals::Int64=1000,
					saveName::String = get_now(),
					showDistance::Bool = false,
					minBox::Bool = false,
					populationSize::Int64 = 50,
					penaltyValue::Float64 = 999999.0,
					gridType::Symbol = :LHC,
					saveStartingValues::Bool = false,
					maxTrialsStartingValues::Int64 = 20,
					thresholdStartingValue::Float64 = penaltyValue/10.0)

	if thresholdStartingValue > penaltyValue
		error("Please set thresholdStartingValue < penaltyValue.")
	end

	#Fminbox does not work with all optimizers in Optim
	if minBox == true
		listValidLocalOptimizers = [:GradientDescent, :BFGS, :LBFGS, :ConjugateGradient]
		if in(localOptimizer, listValidLocalOptimizers) == false
			error("if minBox == true, localOptimizer must be in $(listValidLocalOptimizers)")
		end
	end

	#Check valid grid types
	listValidGridTypes = [:LHC, :Sobol]
	if in(gridType, listValidGridTypes) == false
		error("gridType must be in $(listValidGridTypes)")
	end

	MSMOptions(globalOptimizer,
				localOptimizer,
				maxFuncEvals,
				saveName,
				showDistance,
				minBox,
				populationSize,
				penaltyValue,
				gridType,
				saveStartingValues,
				maxTrialsStartingValues,
				thresholdStartingValue)

end

"""
	MSMProblem

MSMProblem is a mutable struct that caries all the information needed to
perform the optimization and display the results.
"""
mutable struct MSMProblem
	priors::OrderedDict{String,Array{Float64,1}}
	W::Matrix{Float64} #Weight matrix in MSM objective function
	empiricalMoments::OrderedDict{String,Array{Float64,1}}
	simulatedMoments::OrderedDict{String, Float64}
	distanceEmpSimMoments::Float64
	simulate_empirical_moments::Function					#returns an ordered dict
	simulate_empirical_moments_array::Function		#returns an Array
	objective_function::Function
	options::MSMOptions
	bbSetup::BlackBoxOptim.OptController					#set up when using BlackBoxOptim (global minimum)
	bbResults::BlackBoxOptim.OptimizationResults	#results when using BlackBoxOptim (global minimum)
	optimResults::Optim.OptimizationResults				#results when using Optim (local minimum)
	Sigma0::Array{Float64,2}											#distance matrix (in the terminology of Duffie and Singleton (1993))
	Avar::Array{Float64,2}											  #asymptotic variance of the SMM estimate
end

# Constructor for MSMProblem
#------------------------------------------------------------------------------
function MSMProblem(  ; priors::OrderedDict{String,Array{Float64,1}} = OrderedDict{String,Array{Float64,1}}(),
						W::Matrix{Float64}=Matrix(1.0 .* I(1)),
						empiricalMoments::OrderedDict{String,Array{Float64,1}} = OrderedDict{String,Array{Float64,1}}(),
						simulatedMoments::OrderedDict{String, Float64} = OrderedDict{String,Float64}(),
						distanceEmpSimMoments::Float64 = 0.,
						simulate_empirical_moments::Function = default_function,       #returns an ordered dict
						simulate_empirical_moments_array::Function = default_function, #returns an Array
						objective_function::Function = default_function,
						options::MSMOptions = MSMOptions(),
						bbSetup::BlackBoxOptim.OptController = defaultbbOptimOptController,
						bbResults::BlackBoxOptim.OptimizationResults = defaultbbOptimOptimizationResults,
						optimResults::Optim.OptimizationResults = defaultOptimResults,
						Sigma0::Array{Float64,2} = Array{Float64}(undef,0,0),
						Avar::Array{Float64,2} = Array{Float64}(undef,0,0))

	MSMProblem(priors,
				W,
				empiricalMoments,
				simulatedMoments,
				distanceEmpSimMoments,
				simulate_empirical_moments,
				simulate_empirical_moments_array,
				objective_function,
				options,
				bbSetup,
				bbResults,
				optimResults,
				Sigma0,
				Avar)

end

"""
	default_function(x)

Function x->x. Used to initialize functions.
"""
function default_function(x)
	x
end

"""
	rosenbrock2d(x)

Rosenbrock function. Used to initialize BlackBoxOptim.OptController().
"""
function rosenbrock2d(x)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

# It is quite useful to have "default" BlackBoxOptim.OptController and
# BlackBoxOptim.OptimizationResults objects, since
# BlackBoxOptim.OptController() and BlackBoxOptim.OptimizationResults()
# do not work
#-------------------------------------------------------------------------------
defaultbbOptimOptController = bbsetup(x -> rosenbrock2d(x);
											Method=:dxnes,
											SearchRange = (-5.0, 5.0),
											NumDimensions = 2, MaxFuncEvals = 2,
											TraceMode = :silent);

defaultbbOptimOptimizationResults = bboptimize(defaultbbOptimOptController);

defaultOptimResults = optimize(rosenbrock2d, [0.0, 0.0], LBFGS());

"""
	is_global_optimizer(s::Symbol)

function to check that the global optimizer chosen is available.
"""
function is_global_optimizer(s::Symbol)

	# Global Optimizers using BlackBoxOptim
	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	#------------------------------------------------------------------------------
	listValidGlobalOptimizers = [:dxnes, :adaptive_de_rand_1_bin_radiuslimited, :xnes,
								 :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin,
								 :generating_set_search, :de_rand_1_bin,
								 :separable_nes, :resampling_inheritance_memetic_search,
								 :probabilistic_descent, :resampling_memetic_search,
								 :de_rand_2_bin_radiuslimited, :de_rand_2_bin,
								 :random_search, :simultaneous_perturbation_stochastic_approximation]

	in(s, listValidGlobalOptimizers)

end

"""
	is_bb_optimizer(s::Symbol)

function to check whether the global optimizer is using BlackBoxOptim
"""
function is_bb_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listbbOptimizers = [:dxnes, :adaptive_de_rand_1_bin_radiuslimited, :xnes,
					 :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin,
					 :generating_set_search, :de_rand_1_bin,
					 :separable_nes, :resampling_inheritance_memetic_search,
					 :probabilistic_descent, :resampling_memetic_search,
					 :de_rand_2_bin_radiuslimited, :de_rand_2_bin,
					 :random_search, :simultaneous_perturbation_stochastic_approximation]

	in(s, listbbOptimizers)

end

"""
	is_local_optimizer(s::Symbol)

function to check that the local optimizer chosen is available.
"""
function is_local_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listValidLocalOptimizers = [:NelderMead, :SimulatedAnnealing, :ParticleSwarm,
								:BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
								:MomentumGradientDescent, :AcceleratedGradientDescent]

	in(s, listValidLocalOptimizers)

end

"""
	is_optim_optimizer(s::Symbol)

function to check whether the local minimizer uses the package Optim.
"""
function is_optim_optimizer(s::Symbol)

	# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/examples/benchmarking/latest_toplist.csv
	listOptimOptimizers = [:NelderMead, :SimulatedAnnealing, :ParticleSwarm,
							:BFGS, :LBFGS, :ConjugateGradient, :GradientDescent,
							:MomentumGradientDescent, :AcceleratedGradientDescent]

	in(s, listOptimOptimizers)

end

"""
	convert_to_optim_algo(s::Symbol)

function to convert local optimizer (of type Symbol) to an Optim algo.
"""
function convert_to_optim_algo(s::Symbol)
	eval(Meta.parse("$(s)()"))
end

"""
	convert_to_fminbox(s::Symbol)

function to convert local optimizer (of type Symbol) to a Fminbox usable
by Optim.
"""
function convert_to_fminbox(s::Symbol)

	# Old API (before v0.15.0)
	# To be changed when switching to Julia v0.7
	eval(Meta.parse("Fminbox($(String(s))())"))

end


text:
using Test, Random, Distances
using RandomisationModels

macro include_testset(filename)
    @assert filename isa AbstractString
    quote
        @testset $(filename) begin
            include($(filename))
        end
    end
end

@include_testset "models/complete-randomisation.jl"
@include_testset "models/permuted-block.jl"


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

mutable struct DeploymentProperties <: SwaggerModel
    template::Any # spec type: Union{ Nothing, Any } # spec name: template
    templateLink::Any # spec type: Union{ Nothing, TemplateLink } # spec name: templateLink
    parameters::Any # spec type: Union{ Nothing, Any } # spec name: parameters
    parametersLink::Any # spec type: Union{ Nothing, ParametersLink } # spec name: parametersLink
    mode::Any # spec type: Union{ Nothing, String } # spec name: mode
    debugSetting::Any # spec type: Union{ Nothing, DebugSetting } # spec name: debugSetting
    onErrorDeployment::Any # spec type: Union{ Nothing, OnErrorDeployment } # spec name: onErrorDeployment

    function DeploymentProperties(;template=nothing, templateLink=nothing, parameters=nothing, parametersLink=nothing, mode=nothing, debugSetting=nothing, onErrorDeployment=nothing)
        o = new()
        validate_property(DeploymentProperties, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(DeploymentProperties, Symbol("templateLink"), templateLink)
        setfield!(o, Symbol("templateLink"), templateLink)
        validate_property(DeploymentProperties, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(DeploymentProperties, Symbol("parametersLink"), parametersLink)
        setfield!(o, Symbol("parametersLink"), parametersLink)
        validate_property(DeploymentProperties, Symbol("mode"), mode)
        setfield!(o, Symbol("mode"), mode)
        validate_property(DeploymentProperties, Symbol("debugSetting"), debugSetting)
        setfield!(o, Symbol("debugSetting"), debugSetting)
        validate_property(DeploymentProperties, Symbol("onErrorDeployment"), onErrorDeployment)
        setfield!(o, Symbol("onErrorDeployment"), onErrorDeployment)
        o
    end
end # type DeploymentProperties

const _property_map_DeploymentProperties = Dict{Symbol,Symbol}(Symbol("template")=>Symbol("template"), Symbol("templateLink")=>Symbol("templateLink"), Symbol("parameters")=>Symbol("parameters"), Symbol("parametersLink")=>Symbol("parametersLink"), Symbol("mode")=>Symbol("mode"), Symbol("debugSetting")=>Symbol("debugSetting"), Symbol("onErrorDeployment")=>Symbol("onErrorDeployment"))
const _property_types_DeploymentProperties = Dict{Symbol,String}(Symbol("template")=>"Any", Symbol("templateLink")=>"TemplateLink", Symbol("parameters")=>"Any", Symbol("parametersLink")=>"ParametersLink", Symbol("mode")=>"String", Symbol("debugSetting")=>"DebugSetting", Symbol("onErrorDeployment")=>"OnErrorDeployment")
Base.propertynames(::Type{ DeploymentProperties }) = collect(keys(_property_map_DeploymentProperties))
Swagger.property_type(::Type{ DeploymentProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DeploymentProperties[name]))}
Swagger.field_name(::Type{ DeploymentProperties }, property_name::Symbol) =  _property_map_DeploymentProperties[property_name]

const _allowed_DeploymentProperties_mode = ["Incremental", "Complete"]

function check_required(o::DeploymentProperties)
    true
end

function validate_property(::Type{ DeploymentProperties }, name::Symbol, val)
    if name === Symbol("mode")
        Swagger.validate_param(name, "DeploymentProperties", :enum, val, _allowed_DeploymentProperties_mode)
    end
end


text:
# sqfChol.jl

export sqfChol_to_Σ, PDMat_adj;

## Note this follows my unconventional definition of the square-root-free
# Cholesky decomposition with construction from the back of the vector.
function sqfChol_to_Σ(β::Array{Array{T, 1}, 1}, δ::Array{T, 1},
    pdthrowerror::Bool=true) where T <: Real

    sqrtΔ = Diagonal(sqrt.(δ))
    βmat = BayesInference.xpnd_tri(vcat(β...), false)'

    βinvSqrtΔ = Matrix(βmat \ sqrtΔ)
    Σ = βinvSqrtΔ * βinvSqrtΔ'

    PDMat_adj(Σ, throwerror=pdthrowerror)
end

function lNX_sqfChol(X::Union{Array{T, 1}, Array{T, 2}}, μ::Array{T, 1},
    β::Array{Array{T, 1}, 1}, δ::Array{T, 1}, pdthrowerror::Bool=true) where T <: Real

    size(X,1) == length(μ) || throw("In lNX_sqfChol, the columns of X are the observations.")

    Σ = sqfChol_to_Σ(β, δ, pdthrowerror)

    if isnothing(Σ)
        return nothing
    else
        d = MultivariateNormal(μ, Σ)
        if pdthrowerror
            return logpdf(d, X)
        else
            try logpdf(d, X)
            catch excep
                if isa(excep, SingularException)
                    println("Singular exception in lNX_sqfChol.\n")
                    return nothing
                else
                    return logpdf(d, X)
                end
            end
        end
    end
end

# function lNX_seq(x::Array{T, 1}, μ::Array{T, 1},
#     β::Array{Array{T, 1}, 1}, δ::Array{T, 1}) where T <: Real
#
#     K = length(μ)
#     dev = x .- μ
#
#     out = ldnorm(x[K], μ[K], δ[K])
#
#     for k = (K-1):-1:1
#         out += ldnorm(x[k], μ[k] - sum(β[k] .* dev[(k+1):K]) , δ[k])
#     end
#
#     out
# end

### test
#
# using BayesInference
# using Distributions
# using PDMats
# using LinearAlgebra
#
# n = 1000
# K = 5
# X = randn(n, K)
#
# μ_x = randn(K)
# δ_x = exp.(randn(K))
# β_x = [ randn(k) for k = (K-1):-1:1 ]
#
#
# Σ = sqfChol_to_Σ(β_x, δ_x)
# MvNormal(μ_x, Σ)
#
# using BenchmarkTools
#
# @time out1 = lNX(Array(X'), μ_x, β_x, δ_x)
# @btime lNX(Array(X'), μ_x, β_x, δ_x)
# # lNX(X[1,:], μ_x, β_x, δ_x)
#
# @time out2 = [ lNX_seq(X[i,:], μ_x, β_x, δ_x) for i = 1:size(X, 1) ]
# @btime [ lNX_seq(X[i,:], μ_x, β_x, δ_x) for i = 1:size(X, 1) ]
#
# out1 ≈ out2

# naive implemtation (non-sequential) is far superior

function PDMat_adj(A::Matrix{Float64}; maxadd::Float64=1.0e-6,
        epsfact::Float64=100.0, cumadd::Float64=0.0, throwerror::Bool=true)

    try PDMat(A)
    catch excep
        if isa(excep, PosDefException) && cumadd <= maxadd
            a = epsfact * eps(Float64)
            A += a * I
            cumadd += a
            epsfactnext = 10.0 * epsfact
            return PDMat_adj(A, maxadd=maxadd, epsfact=epsfactnext,
                    cumadd=cumadd, throwerror=throwerror)
        else
            println("Failure to adjust to positive definiteness.\n A = ", A, "\n")
            if throwerror
                return PDMat(A) # just trigger original error
            else
                return nothing
            end
        end
    end
end


text:
module JuliaRL

greet() = print("Hello Reinforcement Learning Julia!")

export FeatureCreators
include("FeatureCreators.jl")

export Learning
include("Learning.jl")

export ExperienceReplay, WeightedExperienceReplay, size, getindex, add!
include("Replay.jl")

export
    AbstractEnvironment,
    MountainCar,
    GymEnv,
    start,
    start!,
    step!,
    step,
    get_reward,
    get_state,
    is_terminal,
    get_actions,
    render

include("Environments.jl")

export AbstractState

abstract type AbstractState end

export AbstractPolicy, AbstractQPolicy, EpsilonGreedyQPolicy, get

abstract type AbstractPolicy end

abstract type AbstractQPolicy <: AbstractPolicy end

mutable struct EpsilonGreedyQPolicy <: AbstractQPolicy
    ϵ::Float64
    actions::AbstractArray
end

function get(policy::AbstractQPolicy, values::AbstractArray; rng=Random.GLOBAL_RNG)
    action = findmax(values)[2]
    if rand(rng) < policy.ϵ
        action = rand(rng, policy.actions)
    end
    return action
end

export Agent

include("Agent.jl")

export my_module

module my_module
export func

"""
    func(x)
Returns double the number `x` plus `1`.
"""
func(x) = 2x + 1
end

end # module

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
@testset "Search" begin

    @testset "Approximate" begin
        seq = dna"ACGTACG"

        @testset "forward" begin
            @test approxsearch(seq, dna"", 0) === 1:0
            @test approxsearch(seq, dna"AC", 0) === 1:2
            @test approxsearch(seq, dna"AC", 0, 2) === 5:6
            @test approxsearch(seq, dna"AC", 0, 2, 5) === 0:-1
            @test approxsearch(seq, dna"AT", 0) === 0:-1
            @test approxsearch(seq, dna"AT", 1) === 1:1
            @test approxsearch(seq, dna"AT", 1, 2) === 3:4
            @test approxsearch(seq, dna"AT", 1, 2, 3) === 0:-1
            @test approxsearch(seq, dna"NG", 0) === 2:3
            @test approxsearch(seq, dna"NG", 1) === 1:1
            @test approxsearch(seq, dna"GN", 0) === 3:4
            @test approxsearch(seq, dna"GN", 1) === 1:1
            @test approxsearch(seq, dna"ACG", 0) === 1:3
            @test approxsearch(seq, dna"ACG", 1) === 1:2
            @test approxsearch(seq, dna"ACG", 2) === 1:1
            @test approxsearch(seq, dna"ACG", 3) === 1:0
            @test approxsearch(seq, dna"ACG", 4) === 1:0

            @test approxsearchindex(seq, dna"", 0) === 1
            @test approxsearchindex(seq, dna"AC", 0) === 1
            @test approxsearchindex(seq, dna"AC", 0, 2) === 5
            @test approxsearchindex(seq, dna"AC", 0, 2, 5) === 0

            query = ApproximateSearchQuery(dna"ACG")
            @test approxsearch(seq, query, 1) === 1:2
            @test approxsearch(seq, query, 1, 2) === 2:3
            @test approxsearch(seq, query, 1, 2, 2) === 0:-1
            @test approxsearchindex(seq, query, 1) === 1
            @test approxsearchindex(seq, query, 1, 2) === 2
            @test approxsearchindex(seq, query, 1, 2, 2) === 0
        end

        @testset "backward" begin
            # TODO: maybe this should return 8:7 like rsearch
            @test approxrsearch(seq, dna"", 0) === 7:6
            @test approxrsearch(seq, dna"AC", 0) === 5:6
            @test approxrsearch(seq, dna"AC", 0, 5) === 1:2
            @test approxrsearch(seq, dna"AC", 0, 5, 2) === 0:-1
            @test approxrsearch(seq, dna"AT", 0) === 0:-1
            @test approxrsearch(seq, dna"AT", 1) === 5:6
            @test approxrsearch(seq, dna"AT", 1, 6) === 5:6
            @test approxrsearch(seq, dna"AT", 1, 5) === 5:5
            @test approxrsearch(seq, dna"AT", 1, 3, 2) === 0:-1
            @test approxrsearch(seq, dna"NG", 0) === 6:7
            @test approxrsearch(seq, dna"NG", 0, 6) === 2:3
            @test approxrsearch(seq, dna"GN", 0) === 3:4
            @test approxrsearch(seq, dna"GN", 1) === 7:7
            @test approxrsearch(seq, dna"ACG", 0) === 5:7
            @test approxrsearch(seq, dna"ACG", 1) === 6:7
            @test approxrsearch(seq, dna"ACG", 2) === 7:7
            @test approxrsearch(seq, dna"ACG", 3) === 7:6
            @test approxrsearch(seq, dna"ACG", 4) === 7:6

            # TODO: maybe this should return 8 like rsearchindex
            @test approxrsearchindex(seq, dna"", 0) === 7
            @test approxrsearchindex(seq, dna"AC", 0) === 5
            @test approxrsearchindex(seq, dna"AC", 0, 5) === 1
            @test approxrsearchindex(seq, dna"AC", 0, 5, 2) === 0

            query = ApproximateSearchQuery(dna"ACG")
            @test approxrsearch(seq, query, 1, 7) === 6:7
            @test approxrsearch(seq, query, 1, 6) === 5:6
            @test approxrsearch(seq, query, 1, 6, 6) === 0:-1
            @test approxrsearchindex(seq, query, 1, 7) === 6
            @test approxrsearchindex(seq, query, 1, 6) === 5
            @test approxrsearchindex(seq, query, 1, 6, 6) === 0
        end
    end

    @testset "Regular Expression" begin
        @test isa(biore"A+"d, Seq.RE.Regex{DNA})
        @test isa(biore"A+"r, Seq.RE.Regex{RNA})
        @test isa(biore"A+"a, Seq.RE.Regex{AminoAcid})
        @test isa(biore"A+"dna, Seq.RE.Regex{DNA})
        @test isa(biore"A+"rna, Seq.RE.Regex{RNA})
        @test isa(biore"A+"aa, Seq.RE.Regex{AminoAcid})
        @test_throws Exception eval(:(biore"A+"))
        @test_throws Exception eval(:(biore"A+"foo))
        @test string(biore"A+"dna) == "biore\"A+\"dna"
        @test string(biore"A+"rna) == "biore\"A+\"rna"
        @test string(biore"A+"aa) == "biore\"A+\"aa"

        @test  ismatch(biore"A"d, dna"A")
        @test !ismatch(biore"A"d, dna"C")
        @test !ismatch(biore"A"d, dna"G")
        @test !ismatch(biore"A"d, dna"T")
        @test  ismatch(biore"N"d, dna"A")
        @test  ismatch(biore"N"d, dna"C")
        @test  ismatch(biore"N"d, dna"G")
        @test  ismatch(biore"N"d, dna"T")
        @test  ismatch(biore"[AT]"d, dna"A")
        @test !ismatch(biore"[AT]"d, dna"C")
        @test !ismatch(biore"[AT]"d, dna"G")
        @test  ismatch(biore"[AT]"d, dna"T")
        @test !ismatch(biore"[^AT]"d, dna"A")
        @test  ismatch(biore"[^AT]"d, dna"C")
        @test  ismatch(biore"[^AT]"d, dna"G")
        @test !ismatch(biore"[^AT]"d, dna"T")

        re = biore"^A(C+G*)(T{2,})N$"d
        @test !ismatch(re, dna"AC")
        @test !ismatch(re, dna"AGTT")
        @test !ismatch(re, dna"CCGTT")
        @test !ismatch(re, dna"ACTT")
        @test !ismatch(re, dna"ACTTGT")
        @test  ismatch(re, dna"ACGTTA")
        @test  ismatch(re, dna"ACGTTT")
        @test  ismatch(re, dna"ACCGGTTT")
        @test  ismatch(re, dna"ACCGGTTT")
        @test  ismatch(re, dna"ACCGGTTTA")
        @test  ismatch(re, dna"ACCGGTTTG")

        @test matched(match(re, dna"ACCGTTTTA")) == dna"ACCGTTTTA"
        @test get(captured(match(re, dna"ACCGTTTTA"))[1]) == dna"CCG"
        @test get(captured(match(re, dna"ACCGTTTTA"))[2]) == dna"TTTT"

        # greedy
        @test matched(match(biore"A*"d, dna"AAA")) == dna"AAA"
        @test matched(match(biore"A+"d, dna"AAA")) == dna"AAA"
        @test matched(match(biore"A?"d, dna"AAA")) == dna"A"
        @test matched(match(biore"A{2}"d, dna"AAA")) == dna"AA"
        @test matched(match(biore"A{2,}"d, dna"AAA")) == dna"AAA"
        @test matched(match(biore"A{2,4}"d, dna"AAA")) == dna"AAA"

        # lazy
        @test matched(match(biore"A*?"d, dna"AAA")) == dna""
        @test matched(match(biore"A+?"d, dna"AAA")) == dna"A"
        @test matched(match(biore"A??"d, dna"AAA")) == dna""
        @test matched(match(biore"A{2}?"d, dna"AAA")) == dna"AA"
        @test matched(match(biore"A{2,}?"d, dna"AAA")) == dna"AA"
        @test matched(match(biore"A{2,4}?"d, dna"AAA")) == dna"AA"

        # search
        @test search(dna"ACGTAAT", biore"A+"d) == 1:1
        @test search(dna"ACGTAAT", biore"A+"d, 1) == 1:1
        @test search(dna"ACGTAAT", biore"A+"d, 2) == 5:6
        @test search(dna"ACGTAAT", biore"A+"d, 7) == 0:-1

        # eachmatch
        matches = [dna"CG", dna"GC", dna"GC", dna"CG"]
        for (i, m) in enumerate(eachmatch(biore"GC|CG"d, dna"ACGTTATGCATGGCG"))
            @test matched(m) == matches[i]
        end
        matches = [dna"CG", dna"GC", dna"GC"]
        for (i, m) in enumerate(collect(eachmatch(biore"GC|CG"d, dna"ACGTTATGCATGGCG", false)))
            @test matched(m) == matches[i]
        end

        # matchall
        @test matchall(biore"A*"d, dna"") == [dna""]
        @test matchall(biore"A*"d, dna"AAA") == [
            dna"AAA", dna"AA", dna"A", dna"",
            dna"AA",  dna"A",  dna"",
            dna"A",   dna""]
        @test matchall(biore"AC*G*T"d, dna"ACCGGGT") == [dna"ACCGGGT"]

        @test matchall(biore"A*"d, dna"", false) == [dna""]
        @test matchall(biore"A*"d, dna"AAA", false) == [dna"AAA"]
        @test matchall(biore"AC*G*T"d, dna"ACCGGGT", false) == [dna"ACCGGGT"]

        # RNA and Amino acid
        @test  ismatch(biore"U(A[AG]|GA)$"r, rna"AUUGUAUGA")
        @test !ismatch(biore"U(A[AG]|GA)$"r, rna"AUUGUAUGG")
        @test  ismatch(biore"T+[NQ]A?P"a, aa"MTTQAPMFTQPL")
        @test  ismatch(biore"T+[NQ]A?P"a, aa"MTTAAPMFTQPL")
        @test !ismatch(biore"T+[NQ]A?P"a, aa"MTTAAPMFSQPL")

        # PROSITE
        @test  ismatch(prosite"[AC]-x-V-x(4)-{ED}", aa"ADVAARRK")
        @test  ismatch(prosite"[AC]-x-V-x(4)-{ED}", aa"CPVAARRK")
        @test !ismatch(prosite"[AC]-x-V-x(4)-{ED}", aa"ADVAARRE")
        @test !ismatch(prosite"[AC]-x-V-x(4)-{ED}", aa"CPVAARK")
        @test  ismatch(prosite"<[AC]-x-V-x(4)-{ED}>", aa"ADVAARRK")
        @test !ismatch(prosite"<[AC]-x-V-x(4)-{ED}>", aa"AADVAARRK")
        @test !ismatch(prosite"<[AC]-x-V-x(4)-{ED}>", aa"ADVAARRKA")
    end
end


text:
using JuMP
using LinearAlgebra
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MathOptInterface.Bridges
const MOIU = MathOptInterface.Utilities

import MathOptInterface.FileFormats

include("sdd_constraint.jl")
include("dd_constraint.jl")
include("change_of_basis_model.jl")


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    let
        rates = [0.057, 0.0755, 0.0837, 0.078, 0.1084, 0.0702, 0.1167]
        tenors = [1,3,5,10,15,20,25]
        Yields.Par(rates,tenors)
        Yields.CMT(rates,tenors)
        Yields.Forward(rates,tenors)
        Yields.OIS(rates,tenors)
        Yields.Zero(NelsonSiegel(), rates,tenors)
        Yields.Zero(NelsonSiegelSvensson(), rates,tenors)
        Yields.Zero(rates,tenors)
        c = Yields.Zero(rates,tenors)
        Yields.zero(c)
        Yields.par(c)
        Yields.forward(c,5,6)
    end
    let
        rates = [0.057, 0.0755, 0.0837, 0.078, 0.1084, 0.0702, 0.1167]
        tenors = [1.0,3,5,10,15,20,25]
        Yields.Par(rates,tenors)
        Yields.CMT(rates,tenors)
        Yields.Forward(rates,tenors)
        Yields.OIS(rates,tenors)
        Yields.Zero(rates,tenors)
        Yields.Zero(NelsonSiegel(), rates,tenors)
        Yields.Zero(NelsonSiegelSvensson(), rates,tenors)

    end

end


text:
using BenchmarkTools

include("helpers.jl")

const SUITE = BenchmarkGroup()

kwargs = Dict(:verbose=>false)

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 60

n = 300
p = 3
covars, counts = bdata(n, p, 10^2)

for d in 10 .^ (1:4)
  SUITE[d] = BenchmarkGroup()
  covars, counts = bdata(n, p, d)
  for f in (DMR, HDMR)
    SUITE[d][string(f)] = BenchmarkGroup()
    for par = ("parallel", "serial")
      parallel = (par == "parallel")
      SUITE[d][string(f)][par] = BenchmarkGroup()
      if parallel
        clustertypes = ["local", "remote"]
      else
        clustertypes = ["local"]
      end
      for clus in clustertypes
        local_cluster = (clus == "local")
        # SUITE[d][f][par][clus] = @benchmarkable $(f)($covars, $counts; parallel=$parallel, local_cluster=$local_cluster, $kwargs...)
        SUITE[d][string(f)][par][clus] = @benchmarkable fit($f, $covars, $counts; parallel=$parallel, local_cluster=$local_cluster, $kwargs...)
      end
    end
  end
end


text:
using Test
using Altro
using RobotDynamics
using RobotZoo
const RD = RobotDynamics

# for T in (Float32,Float64)
@testset "Expansions ($T)" for T in (Float32, Float64)
    n,m = 10,5

    @testset "Dynamics expansion ($(typeof(model))" for model in (RobotZoo.Cartpole(), RobotZoo.Quadrotor())
        N = 11
        dt = 0.1
        model = RobotZoo.Quadrotor()
        dmodel = RD.DiscretizedDynamics{RD.RK4}(model)
        n,m = RD.dims(model)
        e = RD.errstate_dim(model)

        # State Expansion
        E = Altro.StateControlExpansion{T}(e,m)
        @test size(E.xx) == (e,e)
        @test size(E.ux) == (m,e)
        @test size(E.uu) == (m,m)
        @test size(E.x) == (e,)
        @test size(E.u) == (m,)
        @test eltype(E.data) == T

        Eerr = Altro.CostExpansion{T}(e,m,N)
        Efull = Altro.FullStateExpansion(Eerr, dmodel)
        if model isa RobotZoo.Cartpole
            @test Efull === E
        end
        let E = Efull[1]
            @test size(E.xx) == (n,n)
            @test size(E.ux) == (m,n)
            @test size(E.uu) == (m,m)
            @test size(E.x) == (n,)
            @test size(E.u) == (m,)
            @test eltype(E.data) == T
        end

        S = Altro.StateControlExpansion{T}(n)
        @test size(S.xx) == (n,n)
        @test size(S.x) == (n,)
        @test_throws UndefRefError S.u
        @test eltype(S.data) == T

        Z = SampledTrajectory(map(1:N) do k
            x,u = rand(model)
            RD.KnotPoint{n,m}(Vector(x),Vector(u),(k-1)*dt,dt)
        end)
        D = [Altro.DynamicsExpansion{T}(n, e, m) for k = 1:N-1]

        for d in D
            if model isa RobotZoo.Cartpole
                @test d.∇f === d.∇e
            else
                @test d.∇f !== d.∇e
            end
        end

        sig = RD.StaticReturn()
        diffmethod = RD.ForwardAD()
        for k = 1:N-1
            RD.jacobian!(sig, diffmethod, dmodel, D[k], Z[k])
        end
        J = zeros(n,n+m)
        xn = zeros(n)
        for k = 1:N-1
            RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), dmodel, J, xn, Z[k])
            @test J ≈ D[k].∇f
            eltype(D[1].A) == T
        end

        G = [zeros(T,n,e) for k = 1:N]
        models = [dmodel for k = 1:length(Z)-1]
        Altro.errstate_jacobians!(models,G,Z)
        for k = 1:N-1
            Gk = zeros(T,n,e)
            RD.errstate_jacobian!(RD.RotationState(), model, Gk, Z[k])
            G[k] ≈ Gk
        end

        Altro.error_expansion!(models, D, G)
        for k = 1:N-1
            A = D[k].A
            B = D[k].B
            G2 = G[k+1]
            G1 = G[k]
            Aerr = G2'A*G1
            Berr = G2'B
            @test D[k].fx ≈ Aerr
            @test D[k].fu ≈ Berr
        end
    end
end


text:
# [Lith, Boonkkamp, and IJzerman](https://doi.org/10.1016/j.amc.2017.09.003)
# A family of different methods that includes the secant method and Newton's method

"""
    LithBoonkkampIJzerman{S,D} <: AbstractNewtonLikeMethod
    LithBoonkkampIJzerman(S,D)

Specify a linear multistep solver with `S` steps and `D` derivatives following [Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003).

## Examples

```
find_zero(sin, 3, Roots.LithBoonkkampIJzerman(2,0)) # the secant method
find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(1,1)) # Newton's method
find_zero((sin,cos), 3, Roots.LithBoonkkampIJzerman(3,1)) # Faster convergence rate
find_zero((sin,cos, x->-sin(x)), 3, Roots.LithBoonkkampIJzerman(1,2)) # Halley-like method
```

The method can be more robust to the intial condition. This example is from the paper (p13). Newton's method (the `S=1`, `D=1` case) fails if `|x₀| ≥ 1.089` but methods with more memory succeed.

```
fx =  ZeroProblem((tanh,x->sech(x)^2), 1.239) # zero at 0.0
solve(fx, Roots.LithBoonkkampIJzerman(1,1)) # Newton, NaN
p = init(fx, Roots.LithBoonkkampIJzerman(2,1))
solve!(p); p.state.steps # 6
p = init(fx, Roots.LithBoonkkampIJzerman(3,1))
solve!(p); p.state.steps # 7
p = init(fx, Roots.LithBoonkkampIJzerman(4,1), verbose=true)
solve!(p); Roots.tracks(p) # 4 iterations, 16 function evaluations
```

Multiple derivatives can be constructed automatically using automatic differentiation. For example,

```
using ForwardDiff
function δ(f, n::Int=1)
    n <= 0 && return f
    n == 1 && return x -> ForwardDiff.derivative(f,float(x))
    δ(δ(f,1),n-1)
end
fs(f,n) = ntuple(i->δ(f,i-1), Val(n+1))

f(x) = cbrt(x)*exp(-x^2) # cf. Table 6 in paper
fx = ZeroProblem(fs(f,1), 0.1147)
opts = (xatol=2eps(), xrtol=0.0, atol=0.0, rtol=0.0) # converge if |xₙ - xₙ₋₁| <= 2ϵ
solve(fx, Roots.LithBoonkkampIJzerman(1, 1); opts...) # NaN -- no convergence
solve(fx, Roots.LithBoonkkampIJzerman(2, 1); opts...) # converges
fx = ZeroProblem(fs(f,2), 0.06)                       # need better starting point
solve(fx, Roots.LithBoonkkampIJzerman(2, 2); opts...) # converges
```

For the case `D=1`, a bracketing method based on this approach is implemented in [`LithBoonkkampIJzermanBracket`](@ref)

## Reference

In [Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003) an analysis is
given of the convergence rates when using linear multistep methods to
solve `0=f(x)` using `f⁻¹(0) = x` when `f` is a sufficiently smooth
linear function. The reformulation, attributed to Grau-Sanchez, finds
a differential equation for `f⁻¹`: `dx/dy = [f⁻¹]′(y) = 1/f′(x) = F` as
`x(0) = x₀ + ∫⁰_y₀ F(x(y)) dy`.

A linear multi-step method is used to solve this equation
numerically.  Let S be the number of memory steps (S= 1,2,...) and D be
the number of derivatives employed, then, with `F(x) = dx/dy`
`x_{n+S} = ∑_{k=0}^{S-1} aₖ x_{n+k} +∑d=1^D ∑_{k=1}^{S-1} aᵈ_{n+k}F⁽ᵈ⁾(x_{n+k})`.
 The `aₖ`s and `aᵈₖ`s are computed each step.

This table is from Tables 1 and 3 of the paper and gives the
convergence rate for simple roots identified therein:

```
s: number of steps remembered
d: number of derivatives uses
s/d  0    1    2    3    4
1    .    2    3    4    5
2    1.62 2.73 3.79 4.82 5.85
3    1.84 2.91 3.95 4.97 5.98
4    1.92 2.97 3.99 4.99 5.996
5    1.97 .    .    .    .
```

That is, more memory leads to a higher convergence rate; more
derivatives leads to a higher convergence rate. However, the
interval about `α`, the zero, where the convergence rate is guaranteed
may get smaller.

!!! Note:
    For the larger values of `S`, the expressions to compute the next value get quite involved.
    The higher convergence rate is likely only to be of help for finding solutions to high precision.

"""
struct LithBoonkkampIJzerman{S,D} <: AbstractNewtonLikeMethod end
LithBoonkkampIJzerman(s,d) = LithBoonkkampIJzerman{s,d}()

function init_state(L::LithBoonkkampIJzerman{S,D}, fs, x) where {S,D}

    xs, ys = init_lith(L, fs, x) # [x₀,x₁,…,xₛ₋₁], ...
    R = eltype(xs)
    T = eltype(ys[1])

    ys′ = vcat(ys...)
    # skip unit consideration here, as won't fit within storage of ys
    state = UnivariateZeroState(xs[end],    # xₙ
                                S >= 2 ? xs[end-1] : one(R)*NaN, # xₙ₋₁
                                one(R)*NaN, # α, or xtar
                                xs,         # all xs (not all but first 2)
                                ys[1][end], # fₙ
                                one(T)*NaN, # fₙ₋₁
                                one(T)*NaN, # f(α)
                                ys′,        # flattened ys, as a vector
                                0,          # no steps
                                iszero(D) + S*(1+D), # no function calls in initialization
                                false,false,false,false,
                                "")

    state
end

function update_state(L::LithBoonkkampIJzerman{S,D}, fs, o::UnivariateZeroState, options) where {S,D}

    xs, ys′ = o.m, o.fm
    ys = ntuple(i -> view(ys′, 1 + (i-1)*S:(i*S)), Val(1+D)) # unflatten ys′

    xᵢ = lmm(L, xs, ys...)

    for i in 1:S-1
        xs[i] = xs[i+1]
    end
    xs[end] = xᵢ

    for (i,fⁱ) ∈ enumerate(fs)
        for j ∈ 1:S-1
            ys[i][j] = ys[i][j+1]
        end
        ys[i][end] = fⁱ(xᵢ)
    end

    o.xn0, o.xn1 = o.xn1, xᵢ
    o.fxn0, o.fxn1 = o.fxn1, ys[1][end]
    incfn(o, 1+D)

    nothing

end

# manufacture initial xs, ys
# use lower memory terms to boot strap up. Secant uses initial default step
function init_lith(L::LithBoonkkampIJzerman{S,D}, fs, x₀) where {S,D}
    # initialize xs, ys
    x̃₀ = float(x₀)
    R = eltype(x̃₀)
    f = first(fs)
    T = promote_type(R, eltype(f(x̃₀)))

    xs = Vector{T}(undef,S)
    ys = ntuple(i -> Vector{T}(undef,S), D+1)
    xs[1] = x̃₀
    for (i,fⁱ) ∈ enumerate(fs)
        ys[i][1] = fⁱ(x̃₀)
    end

    # build up to get S of them
    N = 1
    if D == 0
        x₀ = _default_secant_step(x̃₀)
        xs[2] = xs[1] # shift to first
        xs[1] = x₀
        ys[1][2] = f(x₀)
        N += 1
    end
    for (i,S′) ∈ enumerate(N:S-1)
        xᵢ = lmm(LithBoonkkampIJzerman(S′,D), xs, ys...)
        xs[N+i] = xᵢ
        for (j,fʲ) ∈ enumerate(fs)
            ys[j][N+i] = fʲ(xᵢ)
        end
    end

    xs, ys
end

"""
    LithBoonkkampIJzermanBracket()

A bracketing method which is a modification of Brent's method due to
[Lith, Boonkkamp, and
IJzerman](https://doi.org/10.1016/j.amc.2017.09.003). The best
possible convergence rate is 2.91.

A function, its derivative, and a bracketing interval need to be specified.

The state includes the 3 points -- a bracket `[a,b]` (`b=xₙ` has
`f(b)` closest to `0`) and `c=xₙ₋₁` -- and the corresponding values
for the function and its derivative at these three points.

The next proposed step is either a S=2 or S=3 selection for the
[`LithBoonkkampIJzerman`](@ref) methods with derivative information
included only if it would be of help. The proposed is modified if it
is dithering. The proposed is compared against a bisection step; the
one in the bracket and with the smaller function value is chosen as
the next step.

"""
struct LithBoonkkampIJzermanBracket <: AbstractBracketing end

function init_state(::LithBoonkkampIJzermanBracket, F::FirstDerivative, xs)
    u, v = promote(float.(xs)...)
    fu,fv = F.f(u), F.f(v)
    isbracket(fu, fv) || throw(ArgumentError(bracketing_error))

    if abs(fu) < abs(fv)
        a,b,fa,fb = v,u,fv,fu
    else
        a,b,fa,fb = u,v,fu,fv
    end
    f′a,f′b = F.fp(a), F.fp(b)
    c,fc,f′c = a,fa,f′a

    # keep bracket, int[a,b] as xn1, xn0, m=[c]
    # store fb,fa, [fc, f′a, f′c,f′b]
    state = UnivariateZeroState(b,a, one(a)*NaN, [c], # b,c=xₙ,xₙ-1; int[a,b] a bracket
                                fb,fa, one(fb)*NaN, [fc, f′a, f′c, f′b],
                                0, 4,
                                false, false, false, false,
                                "")
end

function update_state(M::LithBoonkkampIJzermanBracket, F, state::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

    b::T,c::T,a::T = state.xn1, state.m[1], state.xn0
    fb::S,fc::S,fa::S = state.fxn1, state.fm[1], state.fxn0
    f′a::S, f′c::S, f′b::S = state.fm[2],state.fm[3],state.fm[4]

    # Get next interpolating step
    # decide on S and D;
    # S is 3 if a,b,c are distinct; D=1 unless all derivative info will be of the wrong sign.
    s::Int = ((a == c) || (b == c)) ? 2 : 3

    # which derivatives do we include
    sₘ = sign((fb-fa)/(b-a))
    mc, mb = sign(f′c) == sₘ, sign(f′b) == sₘ

    d₀::S = zero(S)
    if s == 2
        if mc || mb
            # D = 1
            as, bs = lmm_coefficients(LithBoonkkampIJzerman{s,1}(), (c,b), (fc, fb))
            h = -fb

            d₀ = -sum(as .* (c,b))
            mb && (d₀ += h * bs[2]/f′b)
            mc && (d₀ += h * bs[1]/f′c)
        else
            d₀ = lmm(LithBoonkkampIJzerman{s, 0}(), (c,b), (fc, fb))
        end
    else
        ma = sign(f′a) == sₘ
        if mc || mb || ma
            # D = 1
            as, bs = lmm_coefficients(LithBoonkkampIJzerman{3,1}(),(a,c,b), (fa,fc,fb))
            h = -fb

            d₀ = -sum(as .* (a,c,b))
            mb && (d₀ += h * bs[end]/f′b) # only when helpful
            mc && (d₀ += h * bs[end-1]/f′c)
            ma && (d₀ += h * bs[end-2]/f′a)
        else
            d₀ = lmm(LithBoonkkampIJzerman{3, 0}(), (a,c,b), (fa, fc,fb))
        end
    end

    # If the step is smaller than the tolerance, use the tolerance as step size.
    xatol, xrtol = options.xabstol, options.xreltol
    δ = xatol + abs(b) * xrtol
    Δ₀ = b - d₀
    if abs(Δ₀) <= δ
        d₀ = b - sign(Δ₀)*δ
    end

    # compare to bisection step; extra function evalution
    d₁ = a + (b-a)* (0.5) #_middle(a, b)
    f₀, f₁ = F.f(d₀), F.f(d₁)

    # interpolation outside a,b or bisection better use that
    d::T,fd::S,f′d::S = zero(T), zero(S), zero(S)
    if (abs(f₀) < abs(f₁)) && (min(a,b) < d₀ < max(a,b))
        d,fd,f′d = d₀,f₀,F.fp(d₀) # interp
    else
        d,fd,f′d = d₁,f₁,F.fp(d₁)  # bisection
    end

    # either [a,d] a bracket or [d,b]
    # [a < d] < b ...c -- b -> d, c-> b (update?)
    # a < [d < b] ...c -- a -> d (update?)
    if sign(fa) * sign(fd) <= 0
        c, fc, f′c = b, fb, f′b
        b, fb, f′b = d, fd, f′d
    else
        a,fa,f′a = d, fd, f′d
    end

    # a,b bracket; keep |fb| ≤ |fa|
    if abs(fa) < abs(fb)
        c, fc,f′c = b, fb, f′b
        a,b,fa,fb,f′a,f′b = b,a,fb,fa,f′b,f′a
    end

    state.xn1,state.xn0, state.m[1] = b,a,c
    state.fxn1,state.fxn0 = fb, fa
    state.fm .= (fc, f′a, f′c, f′b)

    incfn(state, 3)
    nothing

end

function default_tolerances(::M, ::Type{T}, ::Type{S}) where {M<:LithBoonkkampIJzermanBracket,T,S}
    xatol = 2eps(T)
    xrtol = zero(one(T))
    atol = zero(float(one(S))) * oneunit(S)
    rtol = 2eps(float(one(S))) * one(S)
    maxevals = typemax(Int)
    maxfnevals = typemax(Int)
    strict = true
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

### ------
# Script used to generate expressions
#
# Finds expressions by assuming an interpolating polynomial
# goes through the points. From (20) in the paper
# some cases take a **long** time to run.
# At some point, early on, this gets to be more of an academic exercise than
# a performant solution

#= -------
using SymPy

# see Ansatz (20) on p10
function inverse_polynomial_interpretation(s=2,d=2)

    @vars y
    hs = [Sym("h$i") for i ∈ 0:(1+d)*s-1]
    xs = [Sym("x$i") for i ∈ 0:s-1]
    fs = [Sym("f$i") for i ∈ 0:s-1]
    f′s = [Sym("f′$i") for i ∈ 0:s-1]
    f′′s = [Sym("f′′$i") for i ∈ 0:s-1]
    f′′′s = [Sym("f′′′$i") for i ∈ 0:s-1]
    f′′′′s = [Sym("f′′′′$i") for i ∈ 0:s-1]
    f′′′′′s = [Sym("f′′′′′$i") for i ∈ 0:s-1]
    f′′′′′′s = [Sym("f′′′′′′$i") for i ∈ 0:s-1]
    h0 = first(hs)

    H(y) = sum(hs[i+1]*(y-fs[s])^i for i ∈ 0:(1+d)*s-1)
    Hⁱ = H⁰ = H(y)
    Hⁱs = eltype(H⁰)[]
    for i ∈ 1:d
        Hⁱ = diff(Hⁱ,y)
        push!(Hⁱs, Hⁱ)
    end

    eqs = Sym[subs(H(fs[i]), Dict(h0=>xs[s])) - xs[i] for i ∈ 1:s-1]
    for i ∈ 1:s
        # cf. Liptag
        f1,f2,f3,f4,f5,f6 = f′s[i],f′′s[i],f′′′s[i],f′′′′s[i],f′′′′′s[i],f′′′′′′s[i]
        g′ = 1/f1
        g′′ = -f2/f1^3
        g′′′ = (3*f2^2 - f1*f3)/(f1^5)
        g′′′′ = -(15*f2^2 + 10*f1*f2*f3+f1^2*f4)/f1^7
        g′′′′′ = (105*f2^4 -105*f1*f2^2*f3 + 10*f1^2*f3^2 + 15*f1^2*f2*f4 -f1^3*f5)/f1^9
        g′′′′′′ = (-f1^4*f6 + 21*f1^3*f2*f5 + 35*f1^3*f3*f4 - 210*f1^2*f2^2*f4 - 280*f1^2*f2*f3^2 + 1260*f1*f2^3*f3 - 945*f2^5)/f1^11
        gⁱs = [g′,g′′,g′′′, g′′′′,g′′′′′,g′′′′′′]

        for j ∈ 1:d
            push!(eqs, subs(Hⁱs[j], Dict(y=>fs[i], h0=>xs[s])) - gⁱs[j])
        end

    end

    ϕ = sympy.linsolve(eqs, hs[2:end]...)
    ϕ = first(elements(ϕ))
    ϕ = Sym.(convert(Tuple, ϕ.__pyobject__))
    D = Dict{Any,Any}(h0=>xs[s])
    for i in 1:(d+1)*s-1
        D[hs[i+1]] = ϕ[i]
    end
    subs(H(0), D) |> simplify
end

# For g = f⁻¹ return [g', g'', g''',..., g⁽ⁿ⁾]
# (cf [Liptaj](https://vixra.org/pdf/1703.0295v1.pdf)
function liptag(N)

    @vars x₀ Δₓ
    fs = [Sym("f$i") for i ∈ 1:N]
    gs = [Sym("g$i") for i ∈ 1:N]
    a(i) = fs[i]/factorial(i)
    b(i) = gs[i]/factorial(i)

    gᵏs = [1/fs[1]]
    for n ∈ 2:N
        Δy = sum(a(j) * Δₓ^j for j ∈ 1:n)
        l = x₀ + Δₓ
        r = x₀ + sum(b(i)*Δy^i for i ∈ 1:n)
        ϕ = solve(l-r, gs[n])[1]

        for j ∈ 1:n-1
            ϕ = subs(ϕ, gs[j] => gᵏs[j])
        end
        L = limit(ϕ, Δₓ => 0)
        push!(gᵏs, L)
    end
    gᵏs
end
=#

# have computed these
# S/D 0 1 2 3 4  5 6
# 1   x ✓ ✓ ✓ ✓ ✓ ✓
# 2   ✓ ✓ ✓ ✓ ✓ ✓ ✓
# 3   ✓ ✓ - - x x x
# 4   ✓ - x x x x x
# 5   ✓ x x x x x x
# 6   - x x x x x x

# - can be found with script, but answers are too long for
#   inclusion here

## We have two means to do this:
## Using coefficients as,bs, ... returned by lmm_coefficients
## x = ∑ aᵢxᵢ + ∑ⱼ₊₁ⁿ ∑ᵢ bʲᵢFʲᵢ, where Fʲ is the jth derivative of g⁻¹ (F¹ = 1/f'...)
## Using a polynomial interpolant, H(y), going through (xᵢ,fʲ(xᵢ)), j ∈ 0:N)

# secant
function lmm(::LithBoonkkampIJzerman{2,0}, xs, fs)
    x0,x1 = xs
    f0,f1 = fs

    (f0*x1 - f1*x0)/(f0 - f1)
end

function lmm(::LithBoonkkampIJzerman{3,0}, xs, fs)
    x0,x1,x2 = xs
    f0,f1,f2 = fs

    (f0^2*f1*x2 - f0^2*f2*x1 - f0*f1^2*x2 + f0*f2^2*x1 + f1^2*f2*x0 - f1*f2^2*x0)/(f0^2*f1 - f0^2*f2 - f0*f1^2 + f0*f2^2 + f1^2*f2 - f1*f2^2)

end

function lmm(::LithBoonkkampIJzerman{4,0},xs,fs)
    x0,x1,x2,x3 = xs
    f0,f1,f2,f3 = fs

    (f0^3*f1^2*f2*x3 - f0^3*f1^2*f3*x2 - f0^3*f1*f2^2*x3 + f0^3*f1*f3^2*x2 + f0^3*f2^2*f3*x1 - f0^3*f2*f3^2*x1 - f0^2*f1^3*f2*x3 + f0^2*f1^3*f3*x2 + f0^2*f1*f2^3*x3 - f0^2*f1*f3^3*x2 - f0^2*f2^3*f3*x1 + f0^2*f2*f3^3*x1 + f0*f1^3*f2^2*x3 - f0*f1^3*f3^2*x2 - f0*f1^2*f2^3*x3 + f0*f1^2*f3^3*x2 + f0*f2^3*f3^2*x1 - f0*f2^2*f3^3*x1 - f1^3*f2^2*f3*x0 + f1^3*f2*f3^2*x0 + f1^2*f2^3*f3*x0 - f1^2*f2*f3^3*x0 - f1*f2^3*f3^2*x0 + f1*f2^2*f3^3*x0)/(f0^3*f1^2*f2 - f0^3*f1^2*f3 - f0^3*f1*f2^2 + f0^3*f1*f3^2 + f0^3*f2^2*f3 - f0^3*f2*f3^2 - f0^2*f1^3*f2 + f0^2*f1^3*f3 + f0^2*f1*f2^3 - f0^2*f1*f3^3 - f0^2*f2^3*f3 + f0^2*f2*f3^3 + f0*f1^3*f2^2 - f0*f1^3*f3^2 - f0*f1^2*f2^3 + f0*f1^2*f3^3 + f0*f2^3*f3^2 - f0*f2^2*f3^3 - f1^3*f2^2*f3 + f1^3*f2*f3^2 + f1^2*f2^3*f3 - f1^2*f2*f3^3 - f1*f2^3*f3^2 + f1*f2^2*f3^3)

end

function lmm(::LithBoonkkampIJzerman{5,0},xs,fs)
    x0,x1,x2,x3,x4 = xs
    f0,f1,f2,f3,f4 = fs

    (f0^4*f1^3*f2^2*f3*x4 - f0^4*f1^3*f2^2*f4*x3 - f0^4*f1^3*f2*f3^2*x4 + f0^4*f1^3*f2*f4^2*x3 + f0^4*f1^3*f3^2*f4*x2 - f0^4*f1^3*f3*f4^2*x2 - f0^4*f1^2*f2^3*f3*x4 + f0^4*f1^2*f2^3*f4*x3 + f0^4*f1^2*f2*f3^3*x4 - f0^4*f1^2*f2*f4^3*x3 - f0^4*f1^2*f3^3*f4*x2 + f0^4*f1^2*f3*f4^3*x2 + f0^4*f1*f2^3*f3^2*x4 - f0^4*f1*f2^3*f4^2*x3 - f0^4*f1*f2^2*f3^3*x4 + f0^4*f1*f2^2*f4^3*x3 + f0^4*f1*f3^3*f4^2*x2 - f0^4*f1*f3^2*f4^3*x2 - f0^4*f2^3*f3^2*f4*x1 + f0^4*f2^3*f3*f4^2*x1 + f0^4*f2^2*f3^3*f4*x1 - f0^4*f2^2*f3*f4^3*x1 - f0^4*f2*f3^3*f4^2*x1 + f0^4*f2*f3^2*f4^3*x1 - f0^3*f1^4*f2^2*f3*x4 + f0^3*f1^4*f2^2*f4*x3 + f0^3*f1^4*f2*f3^2*x4 - f0^3*f1^4*f2*f4^2*x3 - f0^3*f1^4*f3^2*f4*x2 + f0^3*f1^4*f3*f4^2*x2 + f0^3*f1^2*f2^4*f3*x4 - f0^3*f1^2*f2^4*f4*x3 - f0^3*f1^2*f2*f3^4*x4 + f0^3*f1^2*f2*f4^4*x3 + f0^3*f1^2*f3^4*f4*x2 - f0^3*f1^2*f3*f4^4*x2 - f0^3*f1*f2^4*f3^2*x4 + f0^3*f1*f2^4*f4^2*x3 + f0^3*f1*f2^2*f3^4*x4 - f0^3*f1*f2^2*f4^4*x3 - f0^3*f1*f3^4*f4^2*x2 + f0^3*f1*f3^2*f4^4*x2 + f0^3*f2^4*f3^2*f4*x1 - f0^3*f2^4*f3*f4^2*x1 - f0^3*f2^2*f3^4*f4*x1 + f0^3*f2^2*f3*f4^4*x1 + f0^3*f2*f3^4*f4^2*x1 - f0^3*f2*f3^2*f4^4*x1 + f0^2*f1^4*f2^3*f3*x4 - f0^2*f1^4*f2^3*f4*x3 - f0^2*f1^4*f2*f3^3*x4 + f0^2*f1^4*f2*f4^3*x3 + f0^2*f1^4*f3^3*f4*x2 - f0^2*f1^4*f3*f4^3*x2 - f0^2*f1^3*f2^4*f3*x4 + f0^2*f1^3*f2^4*f4*x3 + f0^2*f1^3*f2*f3^4*x4 - f0^2*f1^3*f2*f4^4*x3 - f0^2*f1^3*f3^4*f4*x2 + f0^2*f1^3*f3*f4^4*x2 + f0^2*f1*f2^4*f3^3*x4 - f0^2*f1*f2^4*f4^3*x3 - f0^2*f1*f2^3*f3^4*x4 + f0^2*f1*f2^3*f4^4*x3 + f0^2*f1*f3^4*f4^3*x2 - f0^2*f1*f3^3*f4^4*x2 - f0^2*f2^4*f3^3*f4*x1 + f0^2*f2^4*f3*f4^3*x1 + f0^2*f2^3*f3^4*f4*x1 - f0^2*f2^3*f3*f4^4*x1 - f0^2*f2*f3^4*f4^3*x1 + f0^2*f2*f3^3*f4^4*x1 - f0*f1^4*f2^3*f3^2*x4 + f0*f1^4*f2^3*f4^2*x3 + f0*f1^4*f2^2*f3^3*x4 - f0*f1^4*f2^2*f4^3*x3 - f0*f1^4*f3^3*f4^2*x2 + f0*f1^4*f3^2*f4^3*x2 + f0*f1^3*f2^4*f3^2*x4 - f0*f1^3*f2^4*f4^2*x3 - f0*f1^3*f2^2*f3^4*x4 + f0*f1^3*f2^2*f4^4*x3 + f0*f1^3*f3^4*f4^2*x2 - f0*f1^3*f3^2*f4^4*x2 - f0*f1^2*f2^4*f3^3*x4 + f0*f1^2*f2^4*f4^3*x3 + f0*f1^2*f2^3*f3^4*x4 - f0*f1^2*f2^3*f4^4*x3 - f0*f1^2*f3^4*f4^3*x2 + f0*f1^2*f3^3*f4^4*x2 + f0*f2^4*f3^3*f4^2*x1 - f0*f2^4*f3^2*f4^3*x1 - f0*f2^3*f3^4*f4^2*x1 + f0*f2^3*f3^2*f4^4*x1 + f0*f2^2*f3^4*f4^3*x1 - f0*f2^2*f3^3*f4^4*x1 + f1^4*f2^3*f3^2*f4*x0 - f1^4*f2^3*f3*f4^2*x0 - f1^4*f2^2*f3^3*f4*x0 + f1^4*f2^2*f3*f4^3*x0 + f1^4*f2*f3^3*f4^2*x0 - f1^4*f2*f3^2*f4^3*x0 - f1^3*f2^4*f3^2*f4*x0 + f1^3*f2^4*f3*f4^2*x0 + f1^3*f2^2*f3^4*f4*x0 - f1^3*f2^2*f3*f4^4*x0 - f1^3*f2*f3^4*f4^2*x0 + f1^3*f2*f3^2*f4^4*x0 + f1^2*f2^4*f3^3*f4*x0 - f1^2*f2^4*f3*f4^3*x0 - f1^2*f2^3*f3^4*f4*x0 + f1^2*f2^3*f3*f4^4*x0 + f1^2*f2*f3^4*f4^3*x0 - f1^2*f2*f3^3*f4^4*x0 - f1*f2^4*f3^3*f4^2*x0 + f1*f2^4*f3^2*f4^3*x0 + f1*f2^3*f3^4*f4^2*x0 - f1*f2^3*f3^2*f4^4*x0 - f1*f2^2*f3^4*f4^3*x0 + f1*f2^2*f3^3*f4^4*x0)/(f0^4*f1^3*f2^2*f3 - f0^4*f1^3*f2^2*f4 - f0^4*f1^3*f2*f3^2 + f0^4*f1^3*f2*f4^2 + f0^4*f1^3*f3^2*f4 - f0^4*f1^3*f3*f4^2 - f0^4*f1^2*f2^3*f3 + f0^4*f1^2*f2^3*f4 + f0^4*f1^2*f2*f3^3 - f0^4*f1^2*f2*f4^3 - f0^4*f1^2*f3^3*f4 + f0^4*f1^2*f3*f4^3 + f0^4*f1*f2^3*f3^2 - f0^4*f1*f2^3*f4^2 - f0^4*f1*f2^2*f3^3 + f0^4*f1*f2^2*f4^3 + f0^4*f1*f3^3*f4^2 - f0^4*f1*f3^2*f4^3 - f0^4*f2^3*f3^2*f4 + f0^4*f2^3*f3*f4^2 + f0^4*f2^2*f3^3*f4 - f0^4*f2^2*f3*f4^3 - f0^4*f2*f3^3*f4^2 + f0^4*f2*f3^2*f4^3 - f0^3*f1^4*f2^2*f3 + f0^3*f1^4*f2^2*f4 + f0^3*f1^4*f2*f3^2 - f0^3*f1^4*f2*f4^2 - f0^3*f1^4*f3^2*f4 + f0^3*f1^4*f3*f4^2 + f0^3*f1^2*f2^4*f3 - f0^3*f1^2*f2^4*f4 - f0^3*f1^2*f2*f3^4 + f0^3*f1^2*f2*f4^4 + f0^3*f1^2*f3^4*f4 - f0^3*f1^2*f3*f4^4 - f0^3*f1*f2^4*f3^2 + f0^3*f1*f2^4*f4^2 + f0^3*f1*f2^2*f3^4 - f0^3*f1*f2^2*f4^4 - f0^3*f1*f3^4*f4^2 + f0^3*f1*f3^2*f4^4 + f0^3*f2^4*f3^2*f4 - f0^3*f2^4*f3*f4^2 - f0^3*f2^2*f3^4*f4 + f0^3*f2^2*f3*f4^4 + f0^3*f2*f3^4*f4^2 - f0^3*f2*f3^2*f4^4 + f0^2*f1^4*f2^3*f3 - f0^2*f1^4*f2^3*f4 - f0^2*f1^4*f2*f3^3 + f0^2*f1^4*f2*f4^3 + f0^2*f1^4*f3^3*f4 - f0^2*f1^4*f3*f4^3 - f0^2*f1^3*f2^4*f3 + f0^2*f1^3*f2^4*f4 + f0^2*f1^3*f2*f3^4 - f0^2*f1^3*f2*f4^4 - f0^2*f1^3*f3^4*f4 + f0^2*f1^3*f3*f4^4 + f0^2*f1*f2^4*f3^3 - f0^2*f1*f2^4*f4^3 - f0^2*f1*f2^3*f3^4 + f0^2*f1*f2^3*f4^4 + f0^2*f1*f3^4*f4^3 - f0^2*f1*f3^3*f4^4 - f0^2*f2^4*f3^3*f4 + f0^2*f2^4*f3*f4^3 + f0^2*f2^3*f3^4*f4 - f0^2*f2^3*f3*f4^4 - f0^2*f2*f3^4*f4^3 + f0^2*f2*f3^3*f4^4 - f0*f1^4*f2^3*f3^2 + f0*f1^4*f2^3*f4^2 + f0*f1^4*f2^2*f3^3 - f0*f1^4*f2^2*f4^3 - f0*f1^4*f3^3*f4^2 + f0*f1^4*f3^2*f4^3 + f0*f1^3*f2^4*f3^2 - f0*f1^3*f2^4*f4^2 - f0*f1^3*f2^2*f3^4 + f0*f1^3*f2^2*f4^4 + f0*f1^3*f3^4*f4^2 - f0*f1^3*f3^2*f4^4 - f0*f1^2*f2^4*f3^3 + f0*f1^2*f2^4*f4^3 + f0*f1^2*f2^3*f3^4 - f0*f1^2*f2^3*f4^4 - f0*f1^2*f3^4*f4^3 + f0*f1^2*f3^3*f4^4 + f0*f2^4*f3^3*f4^2 - f0*f2^4*f3^2*f4^3 - f0*f2^3*f3^4*f4^2 + f0*f2^3*f3^2*f4^4 + f0*f2^2*f3^4*f4^3 - f0*f2^2*f3^3*f4^4 + f1^4*f2^3*f3^2*f4 - f1^4*f2^3*f3*f4^2 - f1^4*f2^2*f3^3*f4 + f1^4*f2^2*f3*f4^3 + f1^4*f2*f3^3*f4^2 - f1^4*f2*f3^2*f4^3 - f1^3*f2^4*f3^2*f4 + f1^3*f2^4*f3*f4^2 + f1^3*f2^2*f3^4*f4 - f1^3*f2^2*f3*f4^4 - f1^3*f2*f3^4*f4^2 + f1^3*f2*f3^2*f4^4 + f1^2*f2^4*f3^3*f4 - f1^2*f2^4*f3*f4^3 - f1^2*f2^3*f3^4*f4 + f1^2*f2^3*f3*f4^4 + f1^2*f2*f3^4*f4^3 - f1^2*f2*f3^3*f4^4 - f1*f2^4*f3^3*f4^2 + f1*f2^4*f3^2*f4^3 + f1*f2^3*f3^4*f4^2 - f1*f2^3*f3^2*f4^4 - f1*f2^2*f3^4*f4^3 + f1*f2^2*f3^3*f4^4)

end

function lmm(::LithBoonkkampIJzerman{6,0},xs,fs)
    x0,x1,x2,x3,x4,x5 = xs
    f0,f1,f2,f3,f4,f5 = fs

    error("not implemented")

end

## d = 1; Newton-like

# return (as, bs⁰,[bs¹,...,bsⁿ⁻¹])
# where -∑ aᵢ xᵢ + h ⋅ ∑ₙ (∑ bsʲᵢ Fʲᵢ)l
function lmm_coefficients(::LithBoonkkampIJzerman{1,1}, xs, fs)
    a0 = -one(xs[1])
    b0 = one(fs[1])

    return (a0,), (b0,)

end

function lmm_coefficients(::LithBoonkkampIJzerman{2,1}, xs, fs)
    q = fs[1]/fs[2]

    # from the paper
    # x2 + a1 x1 + a0x0 =  h3 * (b1 * 1/fp1 + b0 * 1/fp0)
    a0 = (1-3q)/(q-1)^3
    a1 = -1 - a0
    b0 = q/(q-1)^2
    b1 = q * b0

    return (a0,a1), (b0,b1)

end

function lmm_coefficients(::LithBoonkkampIJzerman{3,1}, xs, fs)

    # from the paper
    q0 = fs[3-2]/fs[3]
    q1 = fs[3-1]/fs[3]

    a0 = (q1^2 * (q0 * (3 + 3q1 - 5q0) - q1)) /  ((q0-1)^3 * (q0-q1)^3)
    a1 = (q0^2 *  (q1 * (5q1 - 3q0 - 3)+q0))  /  ((q1-1)^3 * (q0-q1)^3)
    a2 = (q0^2*q1^2*(3q1 - q0*(q1-3) - 5))   /   ((q0-1)^3 * (q1-1)^3) # minor typo in (27c)

    b0 = (q0*q1^2)   / ((q0-1)^2 * (q0-q1)^2)
    b1 = (q0^2*q1)   / ((q0-q1)^2* (q1-1)^2)
    b2 = (q0^2*q1^2) / ((q0-1)^2 * (q1-1)^2)

    return (a0,a1,a2), (b0,b1,b2)
end

function lmm_coefficients(::LithBoonkkampIJzerman{S,1}, xs, fs) where {S}
    error("not computed")
end

function lmm(L::LithBoonkkampIJzerman{S,1}, xs, fs, f′s) where {S}

    as, bs = lmm_coefficients(L, xs, fs)
    Fs = 1 ./ f′s # F = (g⁻¹)'
    h = -fs[S]

    -sum(as[i] * xs[i] for i ∈ 1:S) + h * sum(bs[i] * Fs[i] for i ∈ 1:S )

end

function lmm(::LithBoonkkampIJzerman{4,1},xs, fs, f′s)
    x0,x1,x2,x3 = xs
    f0,f1,f2,f3 = fs
    f′0,f′1,f′2,f′3 = f′s

    # can get with script, but too long as found
    error("not implemented")

end

## d = 2; Halley-like
function lmm(::LithBoonkkampIJzerman{1,2}, xs, fs, f′s,f′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]

    -f0^2*f′′0/(2*f′0^3) - f0/f′0 + x0

end

function lmm(::LithBoonkkampIJzerman{2,2}, xs, fs, f′s,f′′s)
    x0,x1 = xs
    f0,f1 = fs
    f′0,f′1 = f′s
    f′′0,f′′1 = f′′s

    (-f0^5*f1^2*f′0^3*f′′1/2 - f0^5*f1*f′0^3*f′1^2 + f0^5*f′0^3*f′1^3*x1 + f0^4*f1^3*f′0^3*f′′1 + f0^4*f1^3*f′1^3*f′′0/2 + 5*f0^4*f1^2*f′0^3*f′1^2 - 5*f0^4*f1*f′0^3*f′1^3*x1 - f0^3*f1^4*f′0^3*f′′1/2 - f0^3*f1^4*f′1^3*f′′0 - 4*f0^3*f1^3*f′0^3*f′1^2 + 4*f0^3*f1^3*f′0^2*f′1^3 + 10*f0^3*f1^2*f′0^3*f′1^3*x1 + f0^2*f1^5*f′1^3*f′′0/2 - 5*f0^2*f1^4*f′0^2*f′1^3 - 10*f0^2*f1^3*f′0^3*f′1^3*x0 + f0*f1^5*f′0^2*f′1^3 + 5*f0*f1^4*f′0^3*f′1^3*x0 - f1^5*f′0^3*f′1^3*x0)/(f′0^3*f′1^3*(f0^5 - 5*f0^4*f1 + 10*f0^3*f1^2 - 10*f0^2*f1^3 + 5*f0*f1^4 - f1^5))

end

function lmm(::LithBoonkkampIJzerman{3,2}, xs, fs, f′s,f′′s)
    x0,x1,x2 = xs
    f0,f1,f2 = fs
    f′0,f′1,f′2 = f′s
    f′′0,f′′1,f′′2 = f′′s

    ## can get from script, but too long for inclusion here
    error("not implemented")

end

function lmm(::LithBoonkkampIJzerman{4,2}, xs, fs, f′s,f′′s)
    error("Not computed")
end

## d = 3
function lmm(::LithBoonkkampIJzerman{1,3}, xs,  fs, f′s,f′′s, f′′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]

    (f0^3*(f′0*f′′′0 - 3*f′′0^2)/6 - f0^2*f′0^2*f′′0/2 - f0*f′0^4 + f′0^5*x0)/f′0^5

end

function lmm(::LithBoonkkampIJzerman{2,3}, xs,  fs, f′s, f′′s, f′′′s)
    x0,x1 = xs
    f0,f1 = fs
    f′0,f′1 = f′s
    f′′0,f′′1 = f′′s
    f′′′0,f′′′1 = f′′′s

    (f0^7*f1^3*f′0^5*f′1*f′′′1 - 3*f0^7*f1^3*f′0^5*f′′1^2 - 3*f0^7*f1^2*f′0^5*f′1^2*f′′1 - 6*f0^7*f1*f′0^5*f′1^4 + 6*f0^7*f′0^5*f′1^5*x1 - 3*f0^6*f1^4*f′0^5*f′1*f′′′1 + 9*f0^6*f1^4*f′0^5*f′′1^2 + f0^6*f1^4*f′0*f′1^5*f′′′0 - 3*f0^6*f1^4*f′1^5*f′′0^2 + 21*f0^6*f1^3*f′0^5*f′1^2*f′′1 + 42*f0^6*f1^2*f′0^5*f′1^4 - 42*f0^6*f1*f′0^5*f′1^5*x1 + 3*f0^5*f1^5*f′0^5*f′1*f′′′1 - 9*f0^5*f1^5*f′0^5*f′′1^2 - 3*f0^5*f1^5*f′0*f′1^5*f′′′0 + 9*f0^5*f1^5*f′1^5*f′′0^2 - 33*f0^5*f1^4*f′0^5*f′1^2*f′′1 - 15*f0^5*f1^4*f′0^2*f′1^5*f′′0 - 126*f0^5*f1^3*f′0^5*f′1^4 + 126*f0^5*f1^2*f′0^5*f′1^5*x1 - f0^4*f1^6*f′0^5*f′1*f′′′1 + 3*f0^4*f1^6*f′0^5*f′′1^2 + 3*f0^4*f1^6*f′0*f′1^5*f′′′0 - 9*f0^4*f1^6*f′1^5*f′′0^2 + 15*f0^4*f1^5*f′0^5*f′1^2*f′′1 + 33*f0^4*f1^5*f′0^2*f′1^5*f′′0 + 90*f0^4*f1^4*f′0^5*f′1^4 - 90*f0^4*f1^4*f′0^4*f′1^5 - 210*f0^4*f1^3*f′0^5*f′1^5*x1 - f0^3*f1^7*f′0*f′1^5*f′′′0 + 3*f0^3*f1^7*f′1^5*f′′0^2 - 21*f0^3*f1^6*f′0^2*f′1^5*f′′0 + 126*f0^3*f1^5*f′0^4*f′1^5 + 210*f0^3*f1^4*f′0^5*f′1^5*x0 + 3*f0^2*f1^7*f′0^2*f′1^5*f′′0 - 42*f0^2*f1^6*f′0^4*f′1^5 - 126*f0^2*f1^5*f′0^5*f′1^5*x0 + 6*f0*f1^7*f′0^4*f′1^5 + 42*f0*f1^6*f′0^5*f′1^5*x0 - 6*f1^7*f′0^5*f′1^5*x0)/(6*f′0^5*f′1^5*(f0^7 - 7*f0^6*f1 + 21*f0^5*f1^2 - 35*f0^4*f1^3 + 35*f0^3*f1^4 - 21*f0^2*f1^5 + 7*f0*f1^6 - f1^7))

end

function lmm(::LithBoonkkampIJzerman{3,3} ,xs,  fs, f′s, f′′s, f′′′s)
    x0,x1,x2 = xs
    f0,f1,f2 = fs
    f′0,f′1,f′2 = f′s
    f′′0,f′′1,f′′2 = f′′s
    f′′′0,f′′′1,f′′′2 = f′′′s

    # can get from script, but too long for inclusion here
    error("not implemented")

end

function lmm(::LithBoonkkampIJzerman{4,3}, xs,  fs, f′s, f′′s, f′′′s)
    x0,x1,x2,x3 = xs
    f0,f1,f2,f3 = fs
    f′0,f′1,f′2,f′3 = f′s
    f′′0,f′′1,f′′2,f′′3 = f′′s
    f′′′0,f′′′1,f′′′2,f′′′3 = f′′′s

    error("not computed")
end

## d = 4
function lmm(::LithBoonkkampIJzerman{1,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0 = xs[1]
    f0 = fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]

    (-f0^4*(f′0^2*f′′′′0 + 10*f′0*f′′0*f′′′0 + 15*f′′0^2)/24 + f0^3*f′0^2*(f′0*f′′′0 - 3*f′′0^2)/6 - f0^2*f′0^4*f′′0/2 - f0*f′0^6 + f′0^7*x0)/f′0^7

end

function lmm(::LithBoonkkampIJzerman{2,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0,x1 = xs
    f0,f1 = fs
    f′0,f′1 = f′s
    f′′0,f′′1 = f′′s
    f′′′0,f′′′1 = f′′′s
    f′′′′0,f′′′′1 = f′′′′s

    (-f0^9*f1^4*f′0^7*f′1^2*f′′′′1 - 10*f0^9*f1^4*f′0^7*f′1*f′′1*f′′′1 - 15*f0^9*f1^4*f′0^7*f′′1^2 + 4*f0^9*f1^3*f′0^7*f′1^3*f′′′1 - 12*f0^9*f1^3*f′0^7*f′1^2*f′′1^2 - 12*f0^9*f1^2*f′0^7*f′1^4*f′′1 - 24*f0^9*f1*f′0^7*f′1^6 + 24*f0^9*f′0^7*f′1^7*x1 + 4*f0^8*f1^5*f′0^7*f′1^2*f′′′′1 + 40*f0^8*f1^5*f′0^7*f′1*f′′1*f′′′1 + 60*f0^8*f1^5*f′0^7*f′′1^2 + f0^8*f1^5*f′0^2*f′1^7*f′′′′0 + 10*f0^8*f1^5*f′0*f′1^7*f′′0*f′′′0 + 15*f0^8*f1^5*f′1^7*f′′0^2 - 36*f0^8*f1^4*f′0^7*f′1^3*f′′′1 + 108*f0^8*f1^4*f′0^7*f′1^2*f′′1^2 + 108*f0^8*f1^3*f′0^7*f′1^4*f′′1 + 216*f0^8*f1^2*f′0^7*f′1^6 - 216*f0^8*f1*f′0^7*f′1^7*x1 - 6*f0^7*f1^6*f′0^7*f′1^2*f′′′′1 - 60*f0^7*f1^6*f′0^7*f′1*f′′1*f′′′1 - 90*f0^7*f1^6*f′0^7*f′′1^2 - 4*f0^7*f1^6*f′0^2*f′1^7*f′′′′0 - 40*f0^7*f1^6*f′0*f′1^7*f′′0*f′′′0 - 60*f0^7*f1^6*f′1^7*f′′0^2 + 84*f0^7*f1^5*f′0^7*f′1^3*f′′′1 - 252*f0^7*f1^5*f′0^7*f′1^2*f′′1^2 - 24*f0^7*f1^5*f′0^3*f′1^7*f′′′0 + 72*f0^7*f1^5*f′0^2*f′1^7*f′′0^2 - 432*f0^7*f1^4*f′0^7*f′1^4*f′′1 - 864*f0^7*f1^3*f′0^7*f′1^6 + 864*f0^7*f1^2*f′0^7*f′1^7*x1 + 4*f0^6*f1^7*f′0^7*f′1^2*f′′′′1 + 40*f0^6*f1^7*f′0^7*f′1*f′′1*f′′′1 + 60*f0^6*f1^7*f′0^7*f′′1^2 + 6*f0^6*f1^7*f′0^2*f′1^7*f′′′′0 + 60*f0^6*f1^7*f′0*f′1^7*f′′0*f′′′0 + 90*f0^6*f1^7*f′1^7*f′′0^2 - 76*f0^6*f1^6*f′0^7*f′1^3*f′′′1 + 228*f0^6*f1^6*f′0^7*f′1^2*f′′1^2 + 76*f0^6*f1^6*f′0^3*f′1^7*f′′′0 - 228*f0^6*f1^6*f′0^2*f′1^7*f′′0^2 + 588*f0^6*f1^5*f′0^7*f′1^4*f′′1 + 252*f0^6*f1^5*f′0^4*f′1^7*f′′0 + 2016*f0^6*f1^4*f′0^7*f′1^6 - 2016*f0^6*f1^3*f′0^7*f′1^7*x1 - f0^5*f1^8*f′0^7*f′1^2*f′′′′1 - 10*f0^5*f1^8*f′0^7*f′1*f′′1*f′′′1 - 15*f0^5*f1^8*f′0^7*f′′1^2 - 4*f0^5*f1^8*f′0^2*f′1^7*f′′′′0 - 40*f0^5*f1^8*f′0*f′1^7*f′′0*f′′′0 - 60*f0^5*f1^8*f′1^7*f′′0^2 + 24*f0^5*f1^7*f′0^7*f′1^3*f′′′1 - 72*f0^5*f1^7*f′0^7*f′1^2*f′′1^2 - 84*f0^5*f1^7*f′0^3*f′1^7*f′′′0 + 252*f0^5*f1^7*f′0^2*f′1^7*f′′0^2 - 252*f0^5*f1^6*f′0^7*f′1^4*f′′1 - 588*f0^5*f1^6*f′0^4*f′1^7*f′′0 - 1344*f0^5*f1^5*f′0^7*f′1^6 + 1344*f0^5*f1^5*f′0^6*f′1^7 + 3024*f0^5*f1^4*f′0^7*f′1^7*x1 + f0^4*f1^9*f′0^2*f′1^7*f′′′′0 + 10*f0^4*f1^9*f′0*f′1^7*f′′0*f′′′0 + 15*f0^4*f1^9*f′1^7*f′′0^2 + 36*f0^4*f1^8*f′0^3*f′1^7*f′′′0 - 108*f0^4*f1^8*f′0^2*f′1^7*f′′0^2 + 432*f0^4*f1^7*f′0^4*f′1^7*f′′0 - 2016*f0^4*f1^6*f′0^6*f′1^7 - 3024*f0^4*f1^5*f′0^7*f′1^7*x0 - 4*f0^3*f1^9*f′0^3*f′1^7*f′′′0 + 12*f0^3*f1^9*f′0^2*f′1^7*f′′0^2 - 108*f0^3*f1^8*f′0^4*f′1^7*f′′0 + 864*f0^3*f1^7*f′0^6*f′1^7 + 2016*f0^3*f1^6*f′0^7*f′1^7*x0 + 12*f0^2*f1^9*f′0^4*f′1^7*f′′0 - 216*f0^2*f1^8*f′0^6*f′1^7 - 864*f0^2*f1^7*f′0^7*f′1^7*x0 + 24*f0*f1^9*f′0^6*f′1^7 + 216*f0*f1^8*f′0^7*f′1^7*x0 - 24*f1^9*f′0^7*f′1^7*x0)/(24*f′0^7*f′1^7*(f0^9 - 9*f0^8*f1 + 36*f0^7*f1^2 - 84*f0^6*f1^3 + 126*f0^5*f1^4 - 126*f0^4*f1^5 + 84*f0^3*f1^6 - 36*f0^2*f1^7 + 9*f0*f1^8 - f1^9))

end

function lmm(::LithBoonkkampIJzerman{3,4}, xs, fs, f′s, f′′s, f′′′s, f′′′′s)
    x0,x1,x2 = xs
    f0,f1,f2 = fs
    f′0,f′1,f′2 = f′s
    f′′0,f′′1,f′′2 = f′′s
    f′′′0,f′′′1,f′′′2 = f′′′s
    f′′′′0,f′′′′1,f′′′′2 = f′′′′s

    error("not computed")

end

# n = 5

function lmm(::LithBoonkkampIJzerman{1,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0 = xs[1]
    f0= fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]
    f′′′′′0 = f′′′′′s[1]

    (f0^5*(f′0^3*f′′′′′0 - 15*f′0^2*f′′0*f′′′′0 - 10*f′0^2*f′′′0^2 + 105*f′0*f′′0^2*f′′′0 - 105*f′′0^4) - 5*f0^4*f′0^2*(f′0^2*f′′′′0 + 10*f′0*f′′0*f′′′0 + 15*f′′0^2) + 20*f0^3*f′0^4*(f′0*f′′′0 - 3*f′′0^2) - 60*f0^2*f′0^6*f′′0 - 120*f0*f′0^8 + 120*f′0^9*x0)/(120*f′0^9)

end

function lmm(::LithBoonkkampIJzerman{2,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0,x1 = xs
    f0,f1 = fs
    f′0,f′1 = f′s
    f′′0,f′′1 = f′′s
    f′′′0,f′′′1 = f′′′s
    f′′′′0,f′′′′1 = f′′′′s
    f′′′′′0,f′′′′′1 = f′′′′′s

    (f0^11*f1^5*f′0^9*f′1^3*f′′′′′1 - 15*f0^11*f1^5*f′0^9*f′1^2*f′′1*f′′′′1 - 10*f0^11*f1^5*f′0^9*f′1^2*f′′′1^2 + 105*f0^11*f1^5*f′0^9*f′1*f′′1^2*f′′′1 - 105*f0^11*f1^5*f′0^9*f′′1^4 - 5*f0^11*f1^4*f′0^9*f′1^4*f′′′′1 - 50*f0^11*f1^4*f′0^9*f′1^3*f′′1*f′′′1 - 75*f0^11*f1^4*f′0^9*f′1^2*f′′1^2 + 20*f0^11*f1^3*f′0^9*f′1^5*f′′′1 - 60*f0^11*f1^3*f′0^9*f′1^4*f′′1^2 - 60*f0^11*f1^2*f′0^9*f′1^6*f′′1 - 120*f0^11*f1*f′0^9*f′1^8 + 120*f0^11*f′0^9*f′1^9*x1 - 5*f0^10*f1^6*f′0^9*f′1^3*f′′′′′1 + 75*f0^10*f1^6*f′0^9*f′1^2*f′′1*f′′′′1 + 50*f0^10*f1^6*f′0^9*f′1^2*f′′′1^2 - 525*f0^10*f1^6*f′0^9*f′1*f′′1^2*f′′′1 + 525*f0^10*f1^6*f′0^9*f′′1^4 + f0^10*f1^6*f′0^3*f′1^9*f′′′′′0 - 15*f0^10*f1^6*f′0^2*f′1^9*f′′0*f′′′′0 - 10*f0^10*f1^6*f′0^2*f′1^9*f′′′0^2 + 105*f0^10*f1^6*f′0*f′1^9*f′′0^2*f′′′0 - 105*f0^10*f1^6*f′1^9*f′′0^4 + 55*f0^10*f1^5*f′0^9*f′1^4*f′′′′1 + 550*f0^10*f1^5*f′0^9*f′1^3*f′′1*f′′′1 + 825*f0^10*f1^5*f′0^9*f′1^2*f′′1^2 - 220*f0^10*f1^4*f′0^9*f′1^5*f′′′1 + 660*f0^10*f1^4*f′0^9*f′1^4*f′′1^2 + 660*f0^10*f1^3*f′0^9*f′1^6*f′′1 + 1320*f0^10*f1^2*f′0^9*f′1^8 - 1320*f0^10*f1*f′0^9*f′1^9*x1 + 10*f0^9*f1^7*f′0^9*f′1^3*f′′′′′1 - 150*f0^9*f1^7*f′0^9*f′1^2*f′′1*f′′′′1 - 100*f0^9*f1^7*f′0^9*f′1^2*f′′′1^2 + 1050*f0^9*f1^7*f′0^9*f′1*f′′1^2*f′′′1 - 1050*f0^9*f1^7*f′0^9*f′′1^4 - 5*f0^9*f1^7*f′0^3*f′1^9*f′′′′′0 + 75*f0^9*f1^7*f′0^2*f′1^9*f′′0*f′′′′0 + 50*f0^9*f1^7*f′0^2*f′1^9*f′′′0^2 - 525*f0^9*f1^7*f′0*f′1^9*f′′0^2*f′′′0 + 525*f0^9*f1^7*f′1^9*f′′0^4 - 170*f0^9*f1^6*f′0^9*f′1^4*f′′′′1 - 1700*f0^9*f1^6*f′0^9*f′1^3*f′′1*f′′′1 - 2550*f0^9*f1^6*f′0^9*f′1^2*f′′1^2 - 35*f0^9*f1^6*f′0^4*f′1^9*f′′′′0 - 350*f0^9*f1^6*f′0^3*f′1^9*f′′0*f′′′0 - 525*f0^9*f1^6*f′0^2*f′1^9*f′′0^2 + 1100*f0^9*f1^5*f′0^9*f′1^5*f′′′1 - 3300*f0^9*f1^5*f′0^9*f′1^4*f′′1^2 - 3300*f0^9*f1^4*f′0^9*f′1^6*f′′1 - 6600*f0^9*f1^3*f′0^9*f′1^8 + 6600*f0^9*f1^2*f′0^9*f′1^9*x1 - 10*f0^8*f1^8*f′0^9*f′1^3*f′′′′′1 + 150*f0^8*f1^8*f′0^9*f′1^2*f′′1*f′′′′1 + 100*f0^8*f1^8*f′0^9*f′1^2*f′′′1^2 - 1050*f0^8*f1^8*f′0^9*f′1*f′′1^2*f′′′1 + 1050*f0^8*f1^8*f′0^9*f′′1^4 + 10*f0^8*f1^8*f′0^3*f′1^9*f′′′′′0 - 150*f0^8*f1^8*f′0^2*f′1^9*f′′0*f′′′′0 - 100*f0^8*f1^8*f′0^2*f′1^9*f′′′0^2 + 1050*f0^8*f1^8*f′0*f′1^9*f′′0^2*f′′′0 - 1050*f0^8*f1^8*f′1^9*f′′0^4 + 230*f0^8*f1^7*f′0^9*f′1^4*f′′′′1 + 2300*f0^8*f1^7*f′0^9*f′1^3*f′′1*f′′′1 + 3450*f0^8*f1^7*f′0^9*f′1^2*f′′1^2 + 145*f0^8*f1^7*f′0^4*f′1^9*f′′′′0 + 1450*f0^8*f1^7*f′0^3*f′1^9*f′′0*f′′′0 + 2175*f0^8*f1^7*f′0^2*f′1^9*f′′0^2 - 2180*f0^8*f1^6*f′0^9*f′1^5*f′′′1 + 6540*f0^8*f1^6*f′0^9*f′1^4*f′′1^2 + 560*f0^8*f1^6*f′0^5*f′1^9*f′′′0 - 1680*f0^8*f1^6*f′0^4*f′1^9*f′′0^2 + 9900*f0^8*f1^5*f′0^9*f′1^6*f′′1 + 19800*f0^8*f1^4*f′0^9*f′1^8 - 19800*f0^8*f1^3*f′0^9*f′1^9*x1 + 5*f0^7*f1^9*f′0^9*f′1^3*f′′′′′1 - 75*f0^7*f1^9*f′0^9*f′1^2*f′′1*f′′′′1 - 50*f0^7*f1^9*f′0^9*f′1^2*f′′′1^2 + 525*f0^7*f1^9*f′0^9*f′1*f′′1^2*f′′′1 - 525*f0^7*f1^9*f′0^9*f′′1^4 - 10*f0^7*f1^9*f′0^3*f′1^9*f′′′′′0 + 150*f0^7*f1^9*f′0^2*f′1^9*f′′0*f′′′′0 + 100*f0^7*f1^9*f′0^2*f′1^9*f′′′0^2 - 1050*f0^7*f1^9*f′0*f′1^9*f′′0^2*f′′′0 + 1050*f0^7*f1^9*f′1^9*f′′0^4 - 145*f0^7*f1^8*f′0^9*f′1^4*f′′′′1 - 1450*f0^7*f1^8*f′0^9*f′1^3*f′′1*f′′′1 - 2175*f0^7*f1^8*f′0^9*f′1^2*f′′1^2 - 230*f0^7*f1^8*f′0^4*f′1^9*f′′′′0 - 2300*f0^7*f1^8*f′0^3*f′1^9*f′′0*f′′′0 - 3450*f0^7*f1^8*f′0^2*f′1^9*f′′0^2 + 1840*f0^7*f1^7*f′0^9*f′1^5*f′′′1 - 5520*f0^7*f1^7*f′0^9*f′1^4*f′′1^2 - 1840*f0^7*f1^7*f′0^5*f′1^9*f′′′0 + 5520*f0^7*f1^7*f′0^4*f′1^9*f′′0^2 - 12240*f0^7*f1^6*f′0^9*f′1^6*f′′1 - 5040*f0^7*f1^6*f′0^6*f′1^9*f′′0 - 39600*f0^7*f1^5*f′0^9*f′1^8 + 39600*f0^7*f1^4*f′0^9*f′1^9*x1 - f0^6*f1^10*f′0^9*f′1^3*f′′′′′1 + 15*f0^6*f1^10*f′0^9*f′1^2*f′′1*f′′′′1 + 10*f0^6*f1^10*f′0^9*f′1^2*f′′′1^2 - 105*f0^6*f1^10*f′0^9*f′1*f′′1^2*f′′′1 + 105*f0^6*f1^10*f′0^9*f′′1^4 + 5*f0^6*f1^10*f′0^3*f′1^9*f′′′′′0 - 75*f0^6*f1^10*f′0^2*f′1^9*f′′0*f′′′′0 - 50*f0^6*f1^10*f′0^2*f′1^9*f′′′0^2 + 525*f0^6*f1^10*f′0*f′1^9*f′′0^2*f′′′0 - 525*f0^6*f1^10*f′1^9*f′′0^4 + 35*f0^6*f1^9*f′0^9*f′1^4*f′′′′1 + 350*f0^6*f1^9*f′0^9*f′1^3*f′′1*f′′′1 + 525*f0^6*f1^9*f′0^9*f′1^2*f′′1^2 + 170*f0^6*f1^9*f′0^4*f′1^9*f′′′′0 + 1700*f0^6*f1^9*f′0^3*f′1^9*f′′0*f′′′0 + 2550*f0^6*f1^9*f′0^2*f′1^9*f′′0^2 - 560*f0^6*f1^8*f′0^9*f′1^5*f′′′1 + 1680*f0^6*f1^8*f′0^9*f′1^4*f′′1^2 + 2180*f0^6*f1^8*f′0^5*f′1^9*f′′′0 - 6540*f0^6*f1^8*f′0^4*f′1^9*f′′0^2 + 5040*f0^6*f1^7*f′0^9*f′1^6*f′′1 + 12240*f0^6*f1^7*f′0^6*f′1^9*f′′0 + 25200*f0^6*f1^6*f′0^9*f′1^8 - 25200*f0^6*f1^6*f′0^8*f′1^9 - 55440*f0^6*f1^5*f′0^9*f′1^9*x1 - f0^5*f1^11*f′0^3*f′1^9*f′′′′′0 + 15*f0^5*f1^11*f′0^2*f′1^9*f′′0*f′′′′0 + 10*f0^5*f1^11*f′0^2*f′1^9*f′′′0^2 - 105*f0^5*f1^11*f′0*f′1^9*f′′0^2*f′′′0 + 105*f0^5*f1^11*f′1^9*f′′0^4 - 55*f0^5*f1^10*f′0^4*f′1^9*f′′′′0 - 550*f0^5*f1^10*f′0^3*f′1^9*f′′0*f′′′0 - 825*f0^5*f1^10*f′0^2*f′1^9*f′′0^2 - 1100*f0^5*f1^9*f′0^5*f′1^9*f′′′0 + 3300*f0^5*f1^9*f′0^4*f′1^9*f′′0^2 - 9900*f0^5*f1^8*f′0^6*f′1^9*f′′0 + 39600*f0^5*f1^7*f′0^8*f′1^9 + 55440*f0^5*f1^6*f′0^9*f′1^9*x0 + 5*f0^4*f1^11*f′0^4*f′1^9*f′′′′0 + 50*f0^4*f1^11*f′0^3*f′1^9*f′′0*f′′′0 + 75*f0^4*f1^11*f′0^2*f′1^9*f′′0^2 + 220*f0^4*f1^10*f′0^5*f′1^9*f′′′0 - 660*f0^4*f1^10*f′0^4*f′1^9*f′′0^2 + 3300*f0^4*f1^9*f′0^6*f′1^9*f′′0 - 19800*f0^4*f1^8*f′0^8*f′1^9 - 39600*f0^4*f1^7*f′0^9*f′1^9*x0 - 20*f0^3*f1^11*f′0^5*f′1^9*f′′′0 + 60*f0^3*f1^11*f′0^4*f′1^9*f′′0^2 - 660*f0^3*f1^10*f′0^6*f′1^9*f′′0 + 6600*f0^3*f1^9*f′0^8*f′1^9 + 19800*f0^3*f1^8*f′0^9*f′1^9*x0 + 60*f0^2*f1^11*f′0^6*f′1^9*f′′0 - 1320*f0^2*f1^10*f′0^8*f′1^9 - 6600*f0^2*f1^9*f′0^9*f′1^9*x0 + 120*f0*f1^11*f′0^8*f′1^9 + 1320*f0*f1^10*f′0^9*f′1^9*x0 - 120*f1^11*f′0^9*f′1^9*x0)/(120*f′0^9*f′1^9*(f0^11 - 11*f0^10*f1 + 55*f0^9*f1^2 - 165*f0^8*f1^3 + 330*f0^7*f1^4 - 462*f0^6*f1^5 + 462*f0^5*f1^6 - 330*f0^4*f1^7 + 165*f0^3*f1^8 - 55*f0^2*f1^9 + 11*f0*f1^10 - f1^11))

end

function lmm(::LithBoonkkampIJzerman{3,5}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s)
    x0,x1,x2 = xs
    f0,f1,f2 = fs
    f′0,f′1,f′2 = f′s
    f′′0,f′′1,f′′2 = f′′s
    f′′′0,f′′′1,f′′′2 = f′′′s
    f′′′′0,f′′′′1,f′′′′2 = f′′′′s
    f′′′′′0,f′′′′′1,f′′′′′2 = f′′′′′s

    error("not computed")

end

## n = 6
function lmm(::LithBoonkkampIJzerman{1,6}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s, f′′′′′′s)

    x0 = xs[1]
    f0= fs[1]
    f′0 = f′s[1]
    f′′0 = f′′s[1]
    f′′′0 = f′′′s[1]
    f′′′′0 = f′′′′s[1]
    f′′′′′0 = f′′′′′s[1]
    f′′′′′′0 = f′′′′′′s[1]

    (f0^6*(-f′0^4*f′′′′′′0 + 21*f′0^3*f′′0*f′′′′′0 + 35*f′0^3*f′′′0*f′′′′0 - 210*f′0^2*f′′0^2*f′′′′0 - 280*f′0^2*f′′0*f′′′0^2 + 1260*f′0*f′′0^3*f′′′0 - 945*f′′0^5) + 6*f0^5*f′0^2*(f′0^3*f′′′′′0 - 15*f′0^2*f′′0*f′′′′0 - 10*f′0^2*f′′′0^2 + 105*f′0*f′′0^2*f′′′0 - 105*f′′0^4) - 30*f0^4*f′0^4*(f′0^2*f′′′′0 + 10*f′0*f′′0*f′′′0 + 15*f′′0^2) + 120*f0^3*f′0^6*(f′0*f′′′0 - 3*f′′0^2) - 360*f0^2*f′0^8*f′′0 - 720*f0*f′0^10 + 720*f′0^11*x0)/(720*f′0^11)

end

function lmm(::LithBoonkkampIJzerman{2,6}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s, f′′′′′′s)

    x0,x1 = xs
    f0,f1 = fs
    f′0,f′1 = f′s
    f′′0,f′′1 = f′′s
    f′′′0,f′′′1 = f′′′s
    f′′′′0,f′′′′1 = f′′′′s
    f′′′′′0,f′′′′′1 = f′′′′′s
    f′′′′′′0,f′′′′′′1 = f′′′′′′s

    (-f0^13*f1^6*f′0^11*f′1^4*f′′′′′′1 + 21*f0^13*f1^6*f′0^11*f′1^3*f′′1*f′′′′′1 + 35*f0^13*f1^6*f′0^11*f′1^3*f′′′1*f′′′′1 - 210*f0^13*f1^6*f′0^11*f′1^2*f′′1^2*f′′′′1 - 280*f0^13*f1^6*f′0^11*f′1^2*f′′1*f′′′1^2 + 1260*f0^13*f1^6*f′0^11*f′1*f′′1^3*f′′′1 - 945*f0^13*f1^6*f′0^11*f′′1^5 + 6*f0^13*f1^5*f′0^11*f′1^5*f′′′′′1 - 90*f0^13*f1^5*f′0^11*f′1^4*f′′1*f′′′′1 - 60*f0^13*f1^5*f′0^11*f′1^4*f′′′1^2 + 630*f0^13*f1^5*f′0^11*f′1^3*f′′1^2*f′′′1 - 630*f0^13*f1^5*f′0^11*f′1^2*f′′1^4 - 30*f0^13*f1^4*f′0^11*f′1^6*f′′′′1 - 300*f0^13*f1^4*f′0^11*f′1^5*f′′1*f′′′1 - 450*f0^13*f1^4*f′0^11*f′1^4*f′′1^2 + 120*f0^13*f1^3*f′0^11*f′1^7*f′′′1 - 360*f0^13*f1^3*f′0^11*f′1^6*f′′1^2 - 360*f0^13*f1^2*f′0^11*f′1^8*f′′1 - 720*f0^13*f1*f′0^11*f′1^10 + 720*f0^13*f′0^11*f′1^11*x1 + 6*f0^12*f1^7*f′0^11*f′1^4*f′′′′′′1 - 126*f0^12*f1^7*f′0^11*f′1^3*f′′1*f′′′′′1 - 210*f0^12*f1^7*f′0^11*f′1^3*f′′′1*f′′′′1 + 1260*f0^12*f1^7*f′0^11*f′1^2*f′′1^2*f′′′′1 + 1680*f0^12*f1^7*f′0^11*f′1^2*f′′1*f′′′1^2 - 7560*f0^12*f1^7*f′0^11*f′1*f′′1^3*f′′′1 + 5670*f0^12*f1^7*f′0^11*f′′1^5 + f0^12*f1^7*f′0^4*f′1^11*f′′′′′′0 - 21*f0^12*f1^7*f′0^3*f′1^11*f′′0*f′′′′′0 - 35*f0^12*f1^7*f′0^3*f′1^11*f′′′0*f′′′′0 + 210*f0^12*f1^7*f′0^2*f′1^11*f′′0^2*f′′′′0 + 280*f0^12*f1^7*f′0^2*f′1^11*f′′0*f′′′0^2 - 1260*f0^12*f1^7*f′0*f′1^11*f′′0^3*f′′′0 + 945*f0^12*f1^7*f′1^11*f′′0^5 - 78*f0^12*f1^6*f′0^11*f′1^5*f′′′′′1 + 1170*f0^12*f1^6*f′0^11*f′1^4*f′′1*f′′′′1 + 780*f0^12*f1^6*f′0^11*f′1^4*f′′′1^2 - 8190*f0^12*f1^6*f′0^11*f′1^3*f′′1^2*f′′′1 + 8190*f0^12*f1^6*f′0^11*f′1^2*f′′1^4 + 390*f0^12*f1^5*f′0^11*f′1^6*f′′′′1 + 3900*f0^12*f1^5*f′0^11*f′1^5*f′′1*f′′′1 + 5850*f0^12*f1^5*f′0^11*f′1^4*f′′1^2 - 1560*f0^12*f1^4*f′0^11*f′1^7*f′′′1 + 4680*f0^12*f1^4*f′0^11*f′1^6*f′′1^2 + 4680*f0^12*f1^3*f′0^11*f′1^8*f′′1 + 9360*f0^12*f1^2*f′0^11*f′1^10 - 9360*f0^12*f1*f′0^11*f′1^11*x1 - 15*f0^11*f1^8*f′0^11*f′1^4*f′′′′′′1 + 315*f0^11*f1^8*f′0^11*f′1^3*f′′1*f′′′′′1 + 525*f0^11*f1^8*f′0^11*f′1^3*f′′′1*f′′′′1 - 3150*f0^11*f1^8*f′0^11*f′1^2*f′′1^2*f′′′′1 - 4200*f0^11*f1^8*f′0^11*f′1^2*f′′1*f′′′1^2 + 18900*f0^11*f1^8*f′0^11*f′1*f′′1^3*f′′′1 - 14175*f0^11*f1^8*f′0^11*f′′1^5 - 6*f0^11*f1^8*f′0^4*f′1^11*f′′′′′′0 + 126*f0^11*f1^8*f′0^3*f′1^11*f′′0*f′′′′′0 + 210*f0^11*f1^8*f′0^3*f′1^11*f′′′0*f′′′′0 - 1260*f0^11*f1^8*f′0^2*f′1^11*f′′0^2*f′′′′0 - 1680*f0^11*f1^8*f′0^2*f′1^11*f′′0*f′′′0^2 + 7560*f0^11*f1^8*f′0*f′1^11*f′′0^3*f′′′0 - 5670*f0^11*f1^8*f′1^11*f′′0^5 + 300*f0^11*f1^7*f′0^11*f′1^5*f′′′′′1 - 4500*f0^11*f1^7*f′0^11*f′1^4*f′′1*f′′′′1 - 3000*f0^11*f1^7*f′0^11*f′1^4*f′′′1^2 + 31500*f0^11*f1^7*f′0^11*f′1^3*f′′1^2*f′′′1 - 31500*f0^11*f1^7*f′0^11*f′1^2*f′′1^4 - 48*f0^11*f1^7*f′0^5*f′1^11*f′′′′′0 + 720*f0^11*f1^7*f′0^4*f′1^11*f′′0*f′′′′0 + 480*f0^11*f1^7*f′0^4*f′1^11*f′′′0^2 - 5040*f0^11*f1^7*f′0^3*f′1^11*f′′0^2*f′′′0 + 5040*f0^11*f1^7*f′0^2*f′1^11*f′′0^4 - 2340*f0^11*f1^6*f′0^11*f′1^6*f′′′′1 - 23400*f0^11*f1^6*f′0^11*f′1^5*f′′1*f′′′1 - 35100*f0^11*f1^6*f′0^11*f′1^4*f′′1^2 + 9360*f0^11*f1^5*f′0^11*f′1^7*f′′′1 - 28080*f0^11*f1^5*f′0^11*f′1^6*f′′1^2 - 28080*f0^11*f1^4*f′0^11*f′1^8*f′′1 - 56160*f0^11*f1^3*f′0^11*f′1^10 + 56160*f0^11*f1^2*f′0^11*f′1^11*x1 + 20*f0^10*f1^9*f′0^11*f′1^4*f′′′′′′1 - 420*f0^10*f1^9*f′0^11*f′1^3*f′′1*f′′′′′1 - 700*f0^10*f1^9*f′0^11*f′1^3*f′′′1*f′′′′1 + 4200*f0^10*f1^9*f′0^11*f′1^2*f′′1^2*f′′′′1 + 5600*f0^10*f1^9*f′0^11*f′1^2*f′′1*f′′′1^2 - 25200*f0^10*f1^9*f′0^11*f′1*f′′1^3*f′′′1 + 18900*f0^10*f1^9*f′0^11*f′′1^5 + 15*f0^10*f1^9*f′0^4*f′1^11*f′′′′′′0 - 315*f0^10*f1^9*f′0^3*f′1^11*f′′0*f′′′′′0 - 525*f0^10*f1^9*f′0^3*f′1^11*f′′′0*f′′′′0 + 3150*f0^10*f1^9*f′0^2*f′1^11*f′′0^2*f′′′′0 + 4200*f0^10*f1^9*f′0^2*f′1^11*f′′0*f′′′0^2 - 18900*f0^10*f1^9*f′0*f′1^11*f′′0^3*f′′′0 + 14175*f0^10*f1^9*f′1^11*f′′0^5 - 540*f0^10*f1^8*f′0^11*f′1^5*f′′′′′1 + 8100*f0^10*f1^8*f′0^11*f′1^4*f′′1*f′′′′1 + 5400*f0^10*f1^8*f′0^11*f′1^4*f′′′1^2 - 56700*f0^10*f1^8*f′0^11*f′1^3*f′′1^2*f′′′1 + 56700*f0^10*f1^8*f′0^11*f′1^2*f′′1^4 + 246*f0^10*f1^8*f′0^5*f′1^11*f′′′′′0 - 3690*f0^10*f1^8*f′0^4*f′1^11*f′′0*f′′′′0 - 2460*f0^10*f1^8*f′0^4*f′1^11*f′′′0^2 + 25830*f0^10*f1^8*f′0^3*f′1^11*f′′0^2*f′′′0 - 25830*f0^10*f1^8*f′0^2*f′1^11*f′′0^4 + 6060*f0^10*f1^7*f′0^11*f′1^6*f′′′′1 + 60600*f0^10*f1^7*f′0^11*f′1^5*f′′1*f′′′1 + 90900*f0^10*f1^7*f′0^11*f′1^4*f′′1^2 + 1080*f0^10*f1^7*f′0^6*f′1^11*f′′′′0 + 10800*f0^10*f1^7*f′0^5*f′1^11*f′′0*f′′′0 + 16200*f0^10*f1^7*f′0^4*f′1^11*f′′0^2 - 34320*f0^10*f1^6*f′0^11*f′1^7*f′′′1 + 102960*f0^10*f1^6*f′0^11*f′1^6*f′′1^2 + 102960*f0^10*f1^5*f′0^11*f′1^8*f′′1 + 205920*f0^10*f1^4*f′0^11*f′1^10 - 205920*f0^10*f1^3*f′0^11*f′1^11*x1 - 15*f0^9*f1^10*f′0^11*f′1^4*f′′′′′′1 + 315*f0^9*f1^10*f′0^11*f′1^3*f′′1*f′′′′′1 + 525*f0^9*f1^10*f′0^11*f′1^3*f′′′1*f′′′′1 - 3150*f0^9*f1^10*f′0^11*f′1^2*f′′1^2*f′′′′1 - 4200*f0^9*f1^10*f′0^11*f′1^2*f′′1*f′′′1^2 + 18900*f0^9*f1^10*f′0^11*f′1*f′′1^3*f′′′1 - 14175*f0^9*f1^10*f′0^11*f′′1^5 - 20*f0^9*f1^10*f′0^4*f′1^11*f′′′′′′0 + 420*f0^9*f1^10*f′0^3*f′1^11*f′′0*f′′′′′0 + 700*f0^9*f1^10*f′0^3*f′1^11*f′′′0*f′′′′0 - 4200*f0^9*f1^10*f′0^2*f′1^11*f′′0^2*f′′′′0 - 5600*f0^9*f1^10*f′0^2*f′1^11*f′′0*f′′′0^2 + 25200*f0^9*f1^10*f′0*f′1^11*f′′0^3*f′′′0 - 18900*f0^9*f1^10*f′1^11*f′′0^5 + 510*f0^9*f1^9*f′0^11*f′1^5*f′′′′′1 - 7650*f0^9*f1^9*f′0^11*f′1^4*f′′1*f′′′′1 - 5100*f0^9*f1^9*f′0^11*f′1^4*f′′′1^2 + 53550*f0^9*f1^9*f′0^11*f′1^3*f′′1^2*f′′′1 - 53550*f0^9*f1^9*f′0^11*f′1^2*f′′1^4 - 510*f0^9*f1^9*f′0^5*f′1^11*f′′′′′0 + 7650*f0^9*f1^9*f′0^4*f′1^11*f′′0*f′′′′0 + 5100*f0^9*f1^9*f′0^4*f′1^11*f′′′0^2 - 53550*f0^9*f1^9*f′0^3*f′1^11*f′′0^2*f′′′0 + 53550*f0^9*f1^9*f′0^2*f′1^11*f′′0^4 - 7590*f0^9*f1^8*f′0^11*f′1^6*f′′′′1 - 75900*f0^9*f1^8*f′0^11*f′1^5*f′′1*f′′′1 - 113850*f0^9*f1^8*f′0^11*f′1^4*f′′1^2 - 4590*f0^9*f1^8*f′0^6*f′1^11*f′′′′0 - 45900*f0^9*f1^8*f′0^5*f′1^11*f′′0*f′′′0 - 68850*f0^9*f1^8*f′0^4*f′1^11*f′′0^2 + 60600*f0^9*f1^7*f′0^11*f′1^7*f′′′1 - 181800*f0^9*f1^7*f′0^11*f′1^6*f′′1^2 - 14400*f0^9*f1^7*f′0^7*f′1^11*f′′′0 + 43200*f0^9*f1^7*f′0^6*f′1^11*f′′0^2 - 257400*f0^9*f1^6*f′0^11*f′1^8*f′′1 - 514800*f0^9*f1^5*f′0^11*f′1^10 + 514800*f0^9*f1^4*f′0^11*f′1^11*x1 + 6*f0^8*f1^11*f′0^11*f′1^4*f′′′′′′1 - 126*f0^8*f1^11*f′0^11*f′1^3*f′′1*f′′′′′1 - 210*f0^8*f1^11*f′0^11*f′1^3*f′′′1*f′′′′1 + 1260*f0^8*f1^11*f′0^11*f′1^2*f′′1^2*f′′′′1 + 1680*f0^8*f1^11*f′0^11*f′1^2*f′′1*f′′′1^2 - 7560*f0^8*f1^11*f′0^11*f′1*f′′1^3*f′′′1 + 5670*f0^8*f1^11*f′0^11*f′′1^5 + 15*f0^8*f1^11*f′0^4*f′1^11*f′′′′′′0 - 315*f0^8*f1^11*f′0^3*f′1^11*f′′0*f′′′′′0 - 525*f0^8*f1^11*f′0^3*f′1^11*f′′′0*f′′′′0 + 3150*f0^8*f1^11*f′0^2*f′1^11*f′′0^2*f′′′′0 + 4200*f0^8*f1^11*f′0^2*f′1^11*f′′0*f′′′0^2 - 18900*f0^8*f1^11*f′0*f′1^11*f′′0^3*f′′′0 + 14175*f0^8*f1^11*f′1^11*f′′0^5 - 246*f0^8*f1^10*f′0^11*f′1^5*f′′′′′1 + 3690*f0^8*f1^10*f′0^11*f′1^4*f′′1*f′′′′1 + 2460*f0^8*f1^10*f′0^11*f′1^4*f′′′1^2 - 25830*f0^8*f1^10*f′0^11*f′1^3*f′′1^2*f′′′1 + 25830*f0^8*f1^10*f′0^11*f′1^2*f′′1^4 + 540*f0^8*f1^10*f′0^5*f′1^11*f′′′′′0 - 8100*f0^8*f1^10*f′0^4*f′1^11*f′′0*f′′′′0 - 5400*f0^8*f1^10*f′0^4*f′1^11*f′′′0^2 + 56700*f0^8*f1^10*f′0^3*f′1^11*f′′0^2*f′′′0 - 56700*f0^8*f1^10*f′0^2*f′1^11*f′′0^4 + 4590*f0^8*f1^9*f′0^11*f′1^6*f′′′′1 + 45900*f0^8*f1^9*f′0^11*f′1^5*f′′1*f′′′1 + 68850*f0^8*f1^9*f′0^11*f′1^4*f′′1^2 + 7590*f0^8*f1^9*f′0^6*f′1^11*f′′′′0 + 75900*f0^8*f1^9*f′0^5*f′1^11*f′′0*f′′′0 + 113850*f0^8*f1^9*f′0^4*f′1^11*f′′0^2 - 48600*f0^8*f1^8*f′0^11*f′1^7*f′′′1 + 145800*f0^8*f1^8*f′0^11*f′1^6*f′′1^2 + 48600*f0^8*f1^8*f′0^7*f′1^11*f′′′0 - 145800*f0^8*f1^8*f′0^6*f′1^11*f′′0^2 + 297000*f0^8*f1^7*f′0^11*f′1^8*f′′1 + 118800*f0^8*f1^7*f′0^8*f′1^11*f′′0 + 926640*f0^8*f1^6*f′0^11*f′1^10 - 926640*f0^8*f1^5*f′0^11*f′1^11*x1 - f0^7*f1^12*f′0^11*f′1^4*f′′′′′′1 + 21*f0^7*f1^12*f′0^11*f′1^3*f′′1*f′′′′′1 + 35*f0^7*f1^12*f′0^11*f′1^3*f′′′1*f′′′′1 - 210*f0^7*f1^12*f′0^11*f′1^2*f′′1^2*f′′′′1 - 280*f0^7*f1^12*f′0^11*f′1^2*f′′1*f′′′1^2 + 1260*f0^7*f1^12*f′0^11*f′1*f′′1^3*f′′′1 - 945*f0^7*f1^12*f′0^11*f′′1^5 - 6*f0^7*f1^12*f′0^4*f′1^11*f′′′′′′0 + 126*f0^7*f1^12*f′0^3*f′1^11*f′′0*f′′′′′0 + 210*f0^7*f1^12*f′0^3*f′1^11*f′′′0*f′′′′0 - 1260*f0^7*f1^12*f′0^2*f′1^11*f′′0^2*f′′′′0 - 1680*f0^7*f1^12*f′0^2*f′1^11*f′′0*f′′′0^2 + 7560*f0^7*f1^12*f′0*f′1^11*f′′0^3*f′′′0 - 5670*f0^7*f1^12*f′1^11*f′′0^5 + 48*f0^7*f1^11*f′0^11*f′1^5*f′′′′′1 - 720*f0^7*f1^11*f′0^11*f′1^4*f′′1*f′′′′1 - 480*f0^7*f1^11*f′0^11*f′1^4*f′′′1^2 + 5040*f0^7*f1^11*f′0^11*f′1^3*f′′1^2*f′′′1 - 5040*f0^7*f1^11*f′0^11*f′1^2*f′′1^4 - 300*f0^7*f1^11*f′0^5*f′1^11*f′′′′′0 + 4500*f0^7*f1^11*f′0^4*f′1^11*f′′0*f′′′′0 + 3000*f0^7*f1^11*f′0^4*f′1^11*f′′′0^2 - 31500*f0^7*f1^11*f′0^3*f′1^11*f′′0^2*f′′′0 + 31500*f0^7*f1^11*f′0^2*f′1^11*f′′0^4 - 1080*f0^7*f1^10*f′0^11*f′1^6*f′′′′1 - 10800*f0^7*f1^10*f′0^11*f′1^5*f′′1*f′′′1 - 16200*f0^7*f1^10*f′0^11*f′1^4*f′′1^2 - 6060*f0^7*f1^10*f′0^6*f′1^11*f′′′′0 - 60600*f0^7*f1^10*f′0^5*f′1^11*f′′0*f′′′0 - 90900*f0^7*f1^10*f′0^4*f′1^11*f′′0^2 + 14400*f0^7*f1^9*f′0^11*f′1^7*f′′′1 - 43200*f0^7*f1^9*f′0^11*f′1^6*f′′1^2 - 60600*f0^7*f1^9*f′0^7*f′1^11*f′′′0 + 181800*f0^7*f1^9*f′0^6*f′1^11*f′′0^2 - 118800*f0^7*f1^8*f′0^11*f′1^8*f′′1 - 297000*f0^7*f1^8*f′0^8*f′1^11*f′′0 - 570240*f0^7*f1^7*f′0^11*f′1^10 + 570240*f0^7*f1^7*f′0^10*f′1^11 + 1235520*f0^7*f1^6*f′0^11*f′1^11*x1 + f0^6*f1^13*f′0^4*f′1^11*f′′′′′′0 - 21*f0^6*f1^13*f′0^3*f′1^11*f′′0*f′′′′′0 - 35*f0^6*f1^13*f′0^3*f′1^11*f′′′0*f′′′′0 + 210*f0^6*f1^13*f′0^2*f′1^11*f′′0^2*f′′′′0 + 280*f0^6*f1^13*f′0^2*f′1^11*f′′0*f′′′0^2 - 1260*f0^6*f1^13*f′0*f′1^11*f′′0^3*f′′′0 + 945*f0^6*f1^13*f′1^11*f′′0^5 + 78*f0^6*f1^12*f′0^5*f′1^11*f′′′′′0 - 1170*f0^6*f1^12*f′0^4*f′1^11*f′′0*f′′′′0 - 780*f0^6*f1^12*f′0^4*f′1^11*f′′′0^2 + 8190*f0^6*f1^12*f′0^3*f′1^11*f′′0^2*f′′′0 - 8190*f0^6*f1^12*f′0^2*f′1^11*f′′0^4 + 2340*f0^6*f1^11*f′0^6*f′1^11*f′′′′0 + 23400*f0^6*f1^11*f′0^5*f′1^11*f′′0*f′′′0 + 35100*f0^6*f1^11*f′0^4*f′1^11*f′′0^2 + 34320*f0^6*f1^10*f′0^7*f′1^11*f′′′0 - 102960*f0^6*f1^10*f′0^6*f′1^11*f′′0^2 + 257400*f0^6*f1^9*f′0^8*f′1^11*f′′0 - 926640*f0^6*f1^8*f′0^10*f′1^11 - 1235520*f0^6*f1^7*f′0^11*f′1^11*x0 - 6*f0^5*f1^13*f′0^5*f′1^11*f′′′′′0 + 90*f0^5*f1^13*f′0^4*f′1^11*f′′0*f′′′′0 + 60*f0^5*f1^13*f′0^4*f′1^11*f′′′0^2 - 630*f0^5*f1^13*f′0^3*f′1^11*f′′0^2*f′′′0 + 630*f0^5*f1^13*f′0^2*f′1^11*f′′0^4 - 390*f0^5*f1^12*f′0^6*f′1^11*f′′′′0 - 3900*f0^5*f1^12*f′0^5*f′1^11*f′′0*f′′′0 - 5850*f0^5*f1^12*f′0^4*f′1^11*f′′0^2 - 9360*f0^5*f1^11*f′0^7*f′1^11*f′′′0 + 28080*f0^5*f1^11*f′0^6*f′1^11*f′′0^2 - 102960*f0^5*f1^10*f′0^8*f′1^11*f′′0 + 514800*f0^5*f1^9*f′0^10*f′1^11 + 926640*f0^5*f1^8*f′0^11*f′1^11*x0 + 30*f0^4*f1^13*f′0^6*f′1^11*f′′′′0 + 300*f0^4*f1^13*f′0^5*f′1^11*f′′0*f′′′0 + 450*f0^4*f1^13*f′0^4*f′1^11*f′′0^2 + 1560*f0^4*f1^12*f′0^7*f′1^11*f′′′0 - 4680*f0^4*f1^12*f′0^6*f′1^11*f′′0^2 + 28080*f0^4*f1^11*f′0^8*f′1^11*f′′0 - 205920*f0^4*f1^10*f′0^10*f′1^11 - 514800*f0^4*f1^9*f′0^11*f′1^11*x0 - 120*f0^3*f1^13*f′0^7*f′1^11*f′′′0 + 360*f0^3*f1^13*f′0^6*f′1^11*f′′0^2 - 4680*f0^3*f1^12*f′0^8*f′1^11*f′′0 + 56160*f0^3*f1^11*f′0^10*f′1^11 + 205920*f0^3*f1^10*f′0^11*f′1^11*x0 + 360*f0^2*f1^13*f′0^8*f′1^11*f′′0 - 9360*f0^2*f1^12*f′0^10*f′1^11 - 56160*f0^2*f1^11*f′0^11*f′1^11*x0 + 720*f0*f1^13*f′0^10*f′1^11 + 9360*f0*f1^12*f′0^11*f′1^11*x0 - 720*f1^13*f′0^11*f′1^11*x0)/(720*f′0^11*f′1^11*(f0^13 - 13*f0^12*f1 + 78*f0^11*f1^2 - 286*f0^10*f1^3 + 715*f0^9*f1^4 - 1287*f0^8*f1^5 + 1716*f0^7*f1^6 - 1716*f0^6*f1^7 + 1287*f0^5*f1^8 - 715*f0^4*f1^9 + 286*f0^3*f1^10 - 78*f0^2*f1^11 + 13*f0*f1^12 - f1^13))
end

function lmm(::LithBoonkkampIJzerman{3,6}, xs, fs, f′s, f′′s, f′′′s, f′′′′s, f′′′′′s, f′′′′′′s)
    error("not computed")
end


text:
#!/usr/bin/env julia

using Sandbox, Pkg.Artifacts, Scratch, SHA, ghr_jll

if isempty(ARGS)
    println("Usage: julia --project build_docker_image.jl <dir>")
    println("  Where <dir> is something like `debian-julia-python3`")
end

image_name = ARGS[1]
run(`docker build -t $(image_name) $(image_name)`)

artifact_hash = create_artifact() do dir
    @info("Building $(image_name)")
    Sandbox.export_docker_image(image_name, dir; verbose=true)
    @info("Hashing")
end

# Write out to a file
tarball_path = joinpath(@get_scratch!("archived"), "$(image_name).tar.gz")
@info("Archiving out to $(tarball_path)")
archive_artifact(artifact_hash, tarball_path)

# Hash the tarball
@info("Hashing tarball")
tarball_hash = open(io -> bytes2hex(sha256(io)), tarball_path)

# Upload to `staticfloat/Sandbox.jl`, create a tag based on this docker image
tag_name = "$(image_name)-$(bytes2hex(artifact_hash.bytes[end-3:end]))"
@info("Uploading to staticfloat/Sandbox.jl@$(tag_name)")
run(`$(ghr_jll.ghr()) -replace $(tag_name) $(tarball_path)`)

# Bind it into `Artifacts.toml`
tarball_url = "https://github.com/staticfloat/Sandbox.jl/releases/download/$(tag_name)/$(basename(tarball_path))"
bind_artifact!(
    joinpath(dirname(@__DIR__), "Artifacts.toml"),
    "$(image_name)-rootfs",
    artifact_hash;
    download_info=[(tarball_url, tarball_hash)],
    lazy=true,
    force=true,
)


text:
# - LOAD SYSTEM PACKAGES ------------------------------------------- #
using MAT
using RowEchelon
using LinearAlgebra
using SparseArrays
using Distributions
using JSON
using GLPK
using Gurobi
using DelimitedFiles
using Optim
using Logging
using ProgressMeter
# ------------------------------------------------------------------ #

# constants -
const path_to_package = dirname(pathof(@__MODULE__))
#const path_to_package = pwd()
const growth_rate_reaction_index = 13

# local code -
include("./flux/Flux.jl")
include("./flux/Data.jl")
include("./flux/Sample.jl")
include("./flux/Utility.jl")
include("./flux/Types.jl")
include("./flux/Export.jl")
include("./flux/Rules.jl")
include("./flux/Solve.jl")
include("./flux/MOMA.jl")
include("./flux/optimizeRegModel.jl")
include("./flux/solveBooleanReg.jl")
include("./flux/cobrafunctions.jl")


text:
@testset "Doctests" begin
    random_seed = 443591
    DocMeta.setdocmeta!(ExperimentalDesign,
                        :DocTestSetup,
                        :(using ExperimentalDesign, Distributions,
                          Random, StatsModels, DataFrames;
                          Random.seed!($random_seed););
                        recursive = true)
    doctest(ExperimentalDesign)
end


text:
using FastSpike
using ProgressMeter
using Plots
using Statistics

const time = 3600 #seconds
const N = 1000
const rec_every = 60 #seconds
const group_size = 50

if isfile("network.jld2") # Load Network
    net = load("network.jld2")
    excitatory = (sum(net.weight; dims=2).>0)[:, 1]
    inhibitory = (sum(net.weight; dims=2).<0)[:, 1]
else # Or construct a new network
    connection = randomConnection(N, 0.1; amplitude=5)
    excitatory, inhibitory, min_weight, max_weight = EI(
        N, 0.8;
        E_E=(0, 10), E_I=(0, 10), I_E=(-5, -5), I_I=(-5, -5),
        adjacency=connection.adjacency, shuffle=true, seed=0
    )
    #Izhikevich neurons parameters
    a = reshape(0.02 .* excitatory + 0.1 .* inhibitory, (1, N))
    b = 0.2
    c = -65
    d = reshape(8 .* excitatory + 2 .* inhibitory, (1, N))
    # Define Network
    net = Network(
        Izhikevich(1, a, b, c, d, 30.0),
        1,
        STDP(0.1, 0.12, 20, 20; min_weight, max_weight)
    )
    # Add euron groups
    US = add_group!(net, group_size; name="US")
    CS = add_group!(net, group_size; name="CS")
    UR = add_group!(net, group_size; name="UR")
    S = add_group!(net, N - group_size * 3; name="S")
    neurons = US + CS + UR + S
    # connections in neural group
    connect!(net, neurons, neurons, connection)
end
# Define Monitor to record network activities
monitor = Monitor()

net = net |> gpu
# Define training function
function train(time)
    dig_pad = Int(ceil(log10(time)))
    mkdir("plots")
    @showprogress 1 "training " for t = 1:time
        # Make thalamic radom input
        random_input = zeros(1000, N)
        for t_ = 1:1000
            random_input[t_, rand(1:N)] = 20
        end
        # Simulate for 1000 ms
        if t == 1 || t % rec_every == 0 # Record
            for ms = 1:1000 # each second
                run!(net; input_voltage=random_input[ms:ms, :] |> gpu)
                record!(monitor, net)
            end
            # Save recordings
            plot(
                scatter(
                    raster(monitor)[1], markersize=1, c=:black,
                    title="time: $(t)s, #Spikes= $(sum(sum(monitor.spikes))), Mean voltage= $(round(mean(mean(monitor.voltage));digits = 3))",),
                histogram(
                    net.weight[excitatory, :][net.adjacency[excitatory, :]],
                    title="Excitatory weight histogram", xlims=(0, 10));
                layout=grid(2, 1, heights=[0.75, 0.25]),
                legend=false
            )
            savefig("plots/$(lpad(t,dig_pad,"0"))s.png")
            # save (overwrite) the model
            save(net, "network.jld2")
            reset!(monitor)
        else
            for ms = 1:1000 # each second
                run!(net; input_voltage=random_input[ms:ms, :] |> gpu)
            end
        end
    end
    save(net, "network.jld2")
end

# Run
train(time)


text:
#######################################################################
# History: global In/Out and other history variables exported to Main
const In = Dict{Integer,UTF8String}()
const Out = Dict{Integer,Any}()
_ = __ = __ = ans = nothing
export In, Out, _, __, ___, ans

# execution counter
_n = 0

#######################################################################
# methods to clear history or any subset thereof

export clear_history # user-visible in IJulia

function clear_history(indices)
    for n in indices
        delete!(In, n, nothing)
        if haskey(Out, n)
            delete!(Out, n)
            # no way to undefine _{n} variable, but we can set it to nothing
            eval(Main, :($(symbol(string("_",n))) = nothing))
        end
    end
end

# since a range could be huge, intersect it with 1:_n first
clear_history{T<:Integer}(r::Ranges{T}) =
    invoke(clear_history, (Any,), intersect(r, 1:_n))

function clear_history()
    clear_history(1:_n)
    global _, __, ___, ans
    ans = _ = __ = ___ = nothing
end

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


text:
module SyntheticDatasets

using PyCall
using DataFrames
using Random

const datasets = PyNULL()

function __init__()
    copy!(datasets, pyimport("sklearn.datasets"))
end

include("sklearn.jl")
include("matlab.jl")

function convert(features::Array{T, 2}, labels::Array{D, 1})::DataFrame where {T <: Number, D <: Number}
    df = DataFrame()

    for i = 1:size(features)[2]
        df[!, Symbol("feature_$(i)")] = eltype(features)[]
    end

    df[!, :label] = eltype(labels)[]

    for label in unique(labels)
        for i in findall(r->r == label, labels)
            push!(df, (features[i, :]... , label))
        end
    end

    return df
end

function convert(features::Array{T, 2}, labels::Array{D, 2})::DataFrame where {T <: Number, D <: Number}
    df = DataFrame()

    for i = 1:size(features)[2]
        df[!, Symbol("feature_$(i)")] = eltype(features)[]
    end

    for i = 1:size(labels)[2]
        df[!, Symbol("label_$(i)")] = eltype(labels)[]
    end

    for row in 1:size(features)[1]
        push!(df, (features[row, :]... , labels[row, :]...))
    end

    return df
end

end # module


text:
struct SubgradientMethod <: UnconstrainedNonSmoothSolver
  step::StepSizeRule
  n_iter::Int
end

function solve(p::UnconstrainedNonSmoothProblem, params::SubgradientMethod, x0::Vector{Float64}; kwargs...)
  return _solve_sd(p, params, x0; kwargs...)
end

struct ProjectedSubgradientMethod <: ProjectedConstrainedNonSmoothSolver
  step::StepSizeRule
  n_iter::Int
end

function solve(p::ProjectedConstrainedNonSmoothProblem, params::ProjectedSubgradientMethod, x0::Vector{Float64}; kwargs...)
  return _solve_sd(p, params, x0; kwargs...)
end

function _solve_sd(p::NonSmoothProblem, params::NonSmoothSolver, x0::Vector{Float64}; info_callback::Union{Nothing, Function}=nothing)
  # Assumption for p' and params' types:
  # - either (UnconstrainedNonSmoothProblem, SubgradientMethod)
  # - or (ProjectedConstrainedNonSmoothProblem, ProjectedSubgradientMethod).
  @assert length(x0) == p.dimension

  x_best = copy(x0)
  if p isa ProjectedConstrainedNonSmoothProblem
    x0 = p.project(x0)
  end
  f_best = p.f(x0)

  dir = ((p.sense == Minimise) ? -1.0 : 1.0)::Float64

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

    # Perform a gradient step.
    g = p.g(x)
    x += dir * step(params.step, k, norm(g)) * g

    # Projection, if need be.
    if p isa ProjectedConstrainedNonSmoothProblem
      x = p.project(x)
    end

    # Only take the solution if it improves the objective function.
    f = p.f(x)
    if p.sense == Minimise
      if f < f_best
        f_best = f
        x_best = x
      end
    else
      if f > f_best
        f_best = f
        x_best = x
      end
    end

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

  return x_best
end


text:
function llhfit(fit::FitFunction{T, 1, NP}, h::Histogram)::Optim.MultivariateOptimizationResults where {T <: AbstractFloat, NP}
    first_bin::Int = StatsBase.binindex(h, first(fit.fitranges[1]))
    last_bin::Int  = StatsBase.binindex(h, last(fit.fitranges[1]))
    if first_bin < 1 first_bin = 1 end
    if (last_bin > length(h.weights)) last_bin = length(h.weights) end
    bin_centers::Vector{T} = StatsBase.midpoints(h.edges[1])[first_bin:last_bin]
    bin_widths::Vector{T} = [StatsBase.binvolume(h, StatsBase.binindex(h, mp)) for mp in bin_centers]
    counts::Vector{T} = h.weights[first_bin:last_bin]

    function log_likelihood(params)
        s = 0
        if in(0, in.(params, fit.parameter_bounds))
            return Inf
        end
        @inbounds for i in eachindex(counts)
            expected_counts = fit.model(bin_centers[i], params) * bin_widths[i]
            if isnan(expected_counts) || expected_counts < 0
                expected_counts = T(Inf)
            end
            s += -logpdf(Poisson(expected_counts), counts[i])
        end
        return s
    end
    optim_result = Optim.optimize( log_likelihood, fit.initial_parameters )

    uncertainties = try
        sqrt.(diag(inv(ForwardDiff.hessian(log_likelihood, Optim.minimizer(optim_result)))))
    catch e
        zeros(T, length(fit.initial_parameters))
    end

    _set_fit_backend_result!(fit, (optim_result, uncertainties))
    _set_fitted_parameters!(fit, optim_result.minimizer)
    _set_bin_widths!(fit, bin_widths)
    _set_bin_centers!(fit, bin_centers)
    _set_residuals!(fit, bin_centers, counts)
    _set_Χ²!(fit)
    optim_result
end

"""
    llhfit!(fit::FitFunction{T, 1, NP}, h::Histogram)::Nothing where {T <: AbstractFloat, NP}

Performs a log-likelihood fit of the model function `fit.model` and the initial parameters `fit.initial_parameters`
on the histogram `h` in the range `fit.fitranges[1]`. The determined parameters are stored in `fit.fitted_parameters`.

The likelihood for each individual bin is the Poission distribution.
"""
function llhfit!(fit::FitFunction{T, 1, NP}, h::Histogram) where {T <: AbstractFloat, NP}
    llhfit(fit, h)
    fit
end

get_fit_backend_result(fr::Tuple{<:Optim.MultivariateOptimizationResults,<:Any}) = fr[1]

_get_standard_deviations(f::FitFunction, fr::Tuple{<:Optim.MultivariateOptimizationResults,<:Any}) = fr[2]


text:
"""
    collate(samples)

Collates a vector of samples into a single batch. See [collating](../documents/docs/collate.md).
"""
function collate(samples::AbstractVector{<:Tuple})
    !isempty(samples) || return samples
    return tuple([collate([s[key] for s in samples]) for key in keys(samples[1])]...)
end

function collate(samples::AbstractVector{<:Dict})
    !isempty(samples) || return samples
    return Dict(key => collate([s[key] for s in samples]) for key in keys(samples[1]))
end

function collate(samples::AbstractVector{<:NamedTuple})
    !isempty(samples) || return samples
    return (;(key => collate([s[key] for s in samples]) for key in keys(samples[1]))...)
end

collate(samples::AbstractVector{<:AbstractArray{T, N}}) where {T, N} = cat(samples...; dims = N + 1)
collate(samples::AbstractVector) = samples

collate(batch::Tuple) = batch


text:
include("update_beta.jl")

function assumefields!(to::T, from::T) where T
  foreach(fn -> setfield!(to, fn, getfield(from, fn)), fieldnames(T))
end

function update_theta!(state::State, data::Data, prior::Prior,
                       tuners::Tuners; rep_aux::Integer=1,
                       fix::Vector{Symbol}=Symbol[],
                       flags::Vector{Flag}=Flag[])
  update_state!(state, data, prior, tuners,
                fix=[fix; [:p, :beta]], flags=flags)
  for _ in 1:rep_aux
    update_state!(state, data, prior, tuners, flags=flags,
                  fix=[fix; [:p, :beta, :lambda, :gamma, :eta]])
  end
end

function update_state_via_pseudo_prior!(state::State,
                                        state0::State, state1::State,
                                        data::Data, prior::Prior,
                                        tuners::Tuners,
                                        tuners0::Tuners, tuners1::Tuners;
                                        rep_aux::Integer=1,
                                        fix::Vector{Symbol}=Symbol[],
                                        flags::Vector{Flag}=Flag[])
  # Update (β, θᵦ). Though it technically updates β only.
  update_beta_via_pseudo_prior!(state, state0, state1, data, prior)

  # Update θ₀ and θ₁.
  update_theta!(state0, data, prior, tuners0, rep_aux=rep_aux, fix=fix,
                flags=flags)
  update_theta!(state1, data, prior, tuners1, rep_aux=rep_aux, fix=fix,
                flags=flags)

  # Save θᵦ.
  if state.beta
    assumefields!(state, deepcopy(state1))
  else
    assumefields!(state, deepcopy(state0))
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation