text:
using InteractiveUtils

struct RegisteredPlugin
    type::Type
    deps::Vector{<:Type}
end

const registry = IdDict{Type,RegisteredPlugin}()

function register(plugin::Type, dependencies = deps(plugin))
    registry[plugin] = RegisteredPlugin(plugin, dependencies)
    return nothing
end

"""
    Plugins.deps(::Type{T}) = Type[] # where T is your plugin type

Add a method to declare your dependencies.

The plugin type must have a constructor accepting an instance of every of their dependencies.

# Examples

abstract type InterfaceLeft end
struct ImplLeft <: InterfaceLeft end

abstract type InterfaceRight end
struct ImplRight <: InterfaceRight end

Plugins.deps(::Type{ImplLeft}) = [ImplRight]
"""
deps(t) = Type[]

"""
    function autoregister(base=Plugin)

Find and register every concrete subtype of 'base' as a Plugin
"""
function autoregister(base=Plugin)
    for t in subtypes(base) # TODO subtypes is extremely slow
        if isconcretetype(t)
            if !haskey(registry, t)
                @debug "Auto-registering plugin type $t"
                register(t)
            end
        else
            autoregister(t)
        end
    end
end

# t1 implements t2 or
# t1 is an interface that is more specific than t2, or
# t1 and t2 are implementations and t1's direct interface
# is a real subinterface of t2's
function ismorespecific(t1, t2)
    t1 <: t2 && return true
    if isconcretetype(t2) && isconcretetype(t2) &&
        supertype(t1) != supertype(t2) &&
        supertype(t2) != Plugin &&
        supertype(t1) <: supertype(t2)
        return true
    end
    return false
end

function getplugin(t::Type, throw_on_missing=true)::RegisteredPlugin
    found = get(registry, t, nothing)
    !isnothing(found) && return deepcopy(registry[t])
    for p in values(registry) # find the most specific implementation
        if ismorespecific(p.type, t) && (isnothing(found) || ismorespecific(p.type, found.type))
            found = p
        end
    end
    isnothing(found) && throw_on_missing && error("No implementing plugin found for $t")
    return deepcopy(found)
end
getplugin(p::RegisteredPlugin) = p
getplugin(p) = error("The $(typeof(p)) is instantiated. Please provide a type instead!")

allplugins() = values(registry)
plugintypes(plugins) = map(p->p.type, plugins)
plugindeps(plugins) = unique(Iterators.flatten(map(p->p.deps, plugins)))
findimplementation(req, impls) = findfirst(t -> t <: req, impls)

function missingdeps(plugins, dependencies)
    deptypes = plugintypes(dependencies)
    return filter(plugindeps(plugins)) do p
        isnothing(findimplementation(p, deptypes))
    end
end

function find_deps(plugins)
    m = missingdeps(plugins, allplugins())
    length(m) != 0 && error("Missing dependencies: $m")
    lastlength = 0
    retval = plugins
    while lastlength != length(retval)
        lastlength = length(retval)
        retval = vcat(retval, getplugin.(missingdeps(retval, retval)))
    end
    return retval
end

_selectone(plugins) = plugins[findlast(p -> isempty(p.deps), plugins)]

function _removedep(plugins, dep)
    foreach(p -> filter!(d -> !(dep.type <: d), p.deps), plugins)
    return filter!(p -> !(dep.type <: p.type), plugins)
end

function instantiation_order(plugintypes)
    plugins = getplugin.(plugintypes)
    deps = getplugin.(find_deps(plugins))
    remaining = unique(vcat(deps, plugins))
    sorted = RegisteredPlugin[]
    while !isempty(remaining)
        p = _selectone(remaining)
        push!(sorted, registry[p.type])
        _removedep(remaining, p)
    end
    return map(p -> p.type, sorted)
end

# Return a reordered copy of instances to reflect the order of plugintypes.
# Instances not represented in plugintypes will go to the end of the
# returned vector in their original order.
# Instances that are represented several times in plugintypes
# will appear at the place determined by the first representation
function order_instances(instances, plugintypes)
    cache = Dict([typeof(instance) => instance for instance in instances])
    result = []
    for t in plugintypes
        instance = pop!(cache, getplugin(t, false).type, nothing)
        !isnothing(instance) && push!(result, instance)
    end
    for instance in values(cache)
        push!(result, instance)
    end
    return result
end

function instantiate(plugintypes; options...)
    order = instantiation_order(plugintypes)
    instances = []
    cache = IdDict{Type, Any}()
    for t in order
        p = getplugin(t)
        injected_deps = (get(cache, getplugin(dep).type, nothing) for dep in p.deps)
        instance = t(injected_deps...; options...)
        cache[p.type] = instance
        push!(instances, instance)
    end
    return order_instances(instances, plugintypes)
end


text:
"""
xxx
# Examples
```@example

```
"""
struct LTIListener
   sourceList::Vector{Source}
   position::Vector{Float64}
end

function (Ψ::LTIListener)(t₀::Float64)
sourceList = Ψ.sourceList
𝐩ᵣ = Ψ.position
   val = 0.0
   for i = 1:length(sourceList)
      val+=sourceList[i](𝐩ᵣ,t₀)
   end
   return val
end

# DISPLAY
Base.show(io::IO, x::LTIListener) = print(io, "LTI omnidirectional listener")

ext: jl


text:
function testmesh()
  inductor = coilcraft1010vsgroup("1010VS-111ME")
  meshtest1 = mesh(inductor)
  meshtest2 = mesh(Node(1,1,1))
  (g1,junk) = groupfortests()
  meshtest3 = mesh(g1)
end
testmesh()


text:
using TableIO
using Test
using DataFrames
using Dates

using Parquet
using SQLite
using LibPQ

# the following packages are imported automatically when a corresponding file type is used
# using JDF
# using XLSX
# using StatFiles
# using ZipFile
# using JSONTables
# using Arrow

testpath = mktempdir()
println("Temporary directory for test files: ", testpath)

# defining Tables.jl compatible test data
df = DataFrame(a=1:10, b=rand(10), c="hello".* string.(1:10), d=Bool.((1:10) .% 2), e=Date("2020-08-15") .+ Day.(1:10), f="world!" .* string.(1:10))
nt = [(a=1, b=0.5, c="hello"), (a=2, b=0.9, c="world"), (a=3, b=5.5, c="!")]

@testset "TableIO.jl" begin

    include("file_io.jl")
    include("database_io.jl")
    include("plutoui_file_picker.jl")

end


text:
@testset "481.magical-string.jl" begin
    let
        res = [1, 1, 1, 2, 3, 3, 4, 4, 4, 5]
        for i in 1:10
            @test magical_string(i) == res[i]
        end
    end
end


text:
using SafeTestsets

@time begin
    @time @safetestset "Quantum Harmonic oscillator tests" begin
        include("harmonic_oscillator.jl")
    end
    @time @safetestset "Schwarzschild quasinormal modes tests" begin
        include("schwarzschild_qnm.jl")
    end
end


text:
@testset "81.search-in-rotated-sorted-array-ii.jl" begin
    @test search([2, 5, 6, 0, 0, 1, 2], 0) == true
    @test search([2, 5, 6, 0, 0, 1, 2], 3) == false
end


text:
include("..\\temperature_render.jl")

# choosing the range to animate over
# vars = 75.0:1.0:85.0 # changing angle small
vars = [-179.0:1.0:-1.0;1.0:1.0:179.0] # changing angle full
# vars = 5:2.5:85                 # smaller range of angles
# vars = 1.0:5.0:200 # changing mass
# vars = 0:0.1:0.998 # changing spin

# generating frames
n_frames = length(vars)
anim = @animate for (i, var) in enumerate(vars)
    print("Frame: $(i-1) / $n_frames complete")
    fvar = Float64(var)
    temperature_render(
                        mass=10,
                        spin=0.998,
                        obs_angle=fvar,
                        tolerance=1e-12,
                        resolution=1080,
                        size_multiplier=4,
                        dtmax=5
                        )
    # combined_plot(
    #                 mass=10,
    #                 spin=0.998,
    #                 obs_angle=fvar,
    #                 tolerance=1e-12,
    #                 size_multiplier=6,
    #                 fov=10,
    #                 dtmax=20
    # )
end

# saving animation
gif(anim, "full_gif.gif", fps=10)


text:
abstract type AbstractGenomicData end
abstract type AbstractMetabolism end

@units @description mutable struct PMetabolismCN{T<:AbstractGenomicData, D1<:AbstractArray, D2<:AbstractArray} <: AbstractMetabolism
   Isolate::T              | "Isolate defined by BaseGenome"
   y_EM::D1      | mol/mol | "yield of maintenance on reserve"
   y_EX::D1      | mol/mol | "yield of enzyme on reserve"
   y_EV::D1      | mol/mol | "yield of structure on reserve"
   k_E::D2       | d^-1    | "specific reserve turnover rate"
   k_M::D2       | d^-1    | "specific maintenance rate"
   α::D2                   | "fraction of growth flux used for enzyme production"
   α_n::D2                 | "normalized fraction of growth flux used for enzyme production"

   function PMetabolismCN{T, D1, D2}(Isolate::T) where {T<:AbstractGenomicData, D1<:AbstractArray, D2<:AbstractArray}
         y_EV = reserve_yield_constant(genome_size_to_cell_volume(Isolate.L_DNA))
         y_EM = similar(y_EV)
         y_EX = similar(y_EV)
         k_E  = rRNA_copy_number_to_translation_power(Isolate.L_DNA, Isolate.rRNA, Isolate.min_gentime)
         k_M  = cell_volume_to_specific_maintenance_rate(genome_size_to_cell_volume(Isolate.L_DNA))
         α    = similar(k_M)
         α_n  = enzyme_production(α, Isolate.enzyme_distr)
         new(Isolate, y_EM, y_EX, y_EV, k_E, k_M, α, α_n)
     end
end

@units @description mutable struct PMetabolismC{T<:AbstractGenomicData, D1<:AbstractArray} <: AbstractMetabolism
   Isolate::T              | "Isolate defined by BaseGenome"
   y_EM::D1      | mol/mol | "yield of maintenance on reserve"
   y_EX::D1      | mol/mol | "yield of enzyme on reserve"
   y_EV::D1      | mol/mol | "yield of structure on reserve"
   k_E::D1       | d^-1    | "specific reserve turnover rate"
   k_M::D1       | d^-1    | "specific maintenance rate"
   α::D1                   | "fraction of growth flux used for enzyme production"
   α_n::D1                 | "normalized fraction of growth flux used for enzyme production"

   function PMetabolismC{T, D1}(Isolate::T) where {T<:AbstractGenomicData, D1<:AbstractArray}
         y_EV = metabolic_partitioning(genome_size_to_cell_volume(Isolate.L_DNA))
         #y_EV = reserve_yield_constant(genome_size_to_cell_volume(Isolate.L_DNA))[1,:]
         y_EM = similar(y_EV)
         y_EX = similar(y_EV)
         k_E  = rRNA_copy_number_to_translation_power(Isolate.L_DNA, Isolate.rRNA, Isolate.min_gentime)
         k_M  = cell_volume_to_specific_maintenance_rate(genome_size_to_cell_volume(Isolate.L_DNA))
         α    = similar(k_M)
         α_n  = enzyme_production(α, Isolate.enzyme_distr)
         new(Isolate, y_EM, y_EX, y_EV, k_E, k_M, α, α_n)
     end
end

function growth!(r0::Array{Float64,1}, E::Array{Float64,1}, V::Array{Float64,1}, p::PMetabolismC)
   y_VM = p.y_EM./p.y_EV
   #
   function f(r, j)
     m_E     = max(1e-8, E[j]/V[j])
     j_EC    = m_E*(p.k_E[j] - r)
     j_EM    = p.k_M[j]*p.y_EM[j]
     jEM     = max(1e-8, min(j_EC, j_EM))
     jVM     = (j_EM - jEM)*y_VM[j]/p.y_EM[j]
     j_EG    = (j_EC - jEM)/(1.0 + p.α_n[j])
     j_G     = j_EG/p.y_EV[j]
     res     = r - (j_G - jVM)
   end
   #
   f1 = [r->f(r,j) for j in 1:size(r0,1)]
   r = Roots.find_zero.(f1, r0)
 end

 function growth!(r0::Array{Float64,1}, E::Array{Float64,2}, V::Array{Float64,1}, p::PMetabolismCN)
    y_VM = p.y_EM./p.y_EV
    #
    function f(r, j)
      m_E     = @. E[:,j]/V[j]
      j_EC    = @. m_E*(p.k_E - r)
      j_EM    = @. p.k_M*p.y_EM[:,j]
      jEM     = @. min(j_EC, j_EM)
      jVM     = @. (j_EM - jEM)*y_VM[:,j]/p.y_EM[:,j]
      j_EG    = @. (j_EC - jEM)/(1.0 + p.α_n[j])
      j_EG_V  = @. max(1e-6, j_EG./p.y_EV[:,j])
      j_G     = 1.0./(sum(1.0./j_EG_V) - 1.0./sum(j_EG_V))
      res     = r - (j_G - sum(jVM))
    end
    #
    f1 = [r->f(r,j) for j in 1:size(r0,1)]
    r = Roots.find_zero.(f1, r0)
  end

function growth_production(r::Array{Float64,1}, E::Array{Float64,1}, V::Array{Float64,1}, p::PMetabolismC)
   y_VM    = p.y_EM./p.y_EV
   m_E     = @. E/V
   j_EC    = @. m_E*(p.k_E - r)
   j_EM    = @. p.k_M*p.y_EM
   jEM     = @. min(j_EC, j_EM)
   jVM     = @. (j_EM - jEM)*y_VM/p.y_EM
   j_EG    = @. (j_EC - jEM)/(1.0 + p.α_n)
   x       = @. p.α_n*j_EG/p.y_EX
   #
   rG_CO2  = @. (1-1/p.y_EV)*j_EG*V
   rM_CO2  = @. (jEM + jVM)*V
   rX_CO2  = @. (1 - 1/p.y_EX)*p.α_n*j_EG*V
   #
   return x, rG_CO2, rM_CO2, rX_CO2
end

function growth_production(r::Array{Float64,1}, E::Array{Float64,2}, V::Array{Float64,1}, p::PMetabolismCN)
  y_VM    = p.y_EM./p.y_EV
  #
  x = zeros(size(p.α))
  j_Ex   = zeros(size(p.y_EM))
  rG_CO2 = zeros(size(p.y_EM))
  rM_CO2 = zeros(size(p.y_EM))
  rX_CO2 = zeros(size(p.y_EM))
  rEG_rej = zeros(size(p.y_EM))
  rEX_rej = zeros(size(p.y_EM))
  #
  for j in 1:size(r,1)
    m_E     = @. E[:,j]/V[j]
    j_EC    = @. m_E*(p.k_E - r)
    j_EM    = @. p.k_M*p.y_EM[:,j]
    jEM     = @. min(j_EC, j_EM)
    jVM     = @. (j_EM - jEM)*y_VM[:,j]/p.y_EM[:,j]
    j_EG    = @. (j_EC - jEM)/(1.0 + p.α_n[j])
    #
    j_Ex    = @. max(1e-8, p.α_n[j]*j_EG/p.y_EX[:,j])
    x[j]    = 1.0./(sum(1.0./j_Ex) - 1.0./sum(j_Ex))
    #x[isnan.(x)] .= 0.0
    #rG_CO2[:,j]  = @. (1-1/p.y_EV[:,j])*j_EG*V[j]
    rM_CO2[:,j]  = @. (jEM + jVM)*V[j]
    #rX_CO2[:,j]  = @. (1 - 1/p.y_EX[:,j])*p.α_n[j]*j_EG*V[j]
    rEG_rej[:,j]  = max.(0.0, (p.k_E - r).* m_E - jEM - p.y_EV[:,j]*(r[j] + sum(jVM)))
    rEX_rej[:,j] = max.(0.0, p.α_n[j].*((p.k_E - r).* m_E - jEM) - p.y_EX[:,j]*x[j])
  end
  return x, rM_CO2, rEG_rej, rEX_rej
end

function cell_volume_to_specific_maintenance_rate(V_c::Array{Float64,1})
    # Lynch & Marinov (2015), Eq. 1a
    n_0 = cell_volume_to_cell_number_density(V_c)
    V_cell = V_c./1e-18         #  V_cell in μm^3
    E_M = 0.39*V_cell.^0.88     # [1e9 ATP/(cell h)]
    E_G = 1e9*E_M/26            # 26 ATP per glucose molecule: [glucose/(cell h)]
    E_C = E_G*6*24              # [C atoms/(cell/d)]
    mol_C = E_C/(6.022e23)      # [mol C/(cell d)]
    k_M = mol_C./n_0            # 1/d
end

function reserve_yield_constant(V_c::Array{Float64,1})
    m_cell = cell_volume_to_dry_mass(V_c)
    m_baseline = cell_volume_to_dry_mass_baseline(V_c)
    z = m_baseline./m_cell
    X_z = [0.43, 0.143]
    m_growth = cell_volume_to_dry_mass_growth(V_c)
    g = m_growth./m_cell
    X_g = [0.37, 0.153]
    y_EV = zeros(size(X_z,1), size(z,1))
    for i in 1:size(z,1)
            y_EV[:,i] = @. g[i]*X_g/(z[i]*X_z)
    end
    return y_EV
end

function rRNA_copy_number_to_translation_power(rRNA::Array{Float64,1})
    log2_rna = log2.(rRNA)
    y = @. 0.8 + 0.61*log2_rna
    fg = y*24 #per day
end

function rRNA_copy_number_to_translation_power(L_DNA::Array{Float64,1}, rRNA::Array{Float64,1}, min_gentime::Array{Float64,1})
    log2_rna = log2.(rRNA)
    y = @. 0.8 + 0.61*log2_rna
    fg = y #per day
    V_c = genome_size_to_cell_volume(L_DNA)
    V_r = cell_volume_to_ribosome_volume(V_c)
    V_p = cell_volume_to_protein_volume(V_c)
    α = @. V_r./V_p
    k_E = fg.*α
    μ_max = @. log(2)/(min_gentime/24)
    k_E[k_E.>μ_max] = 0.95*μ_max[k_E.>μ_max]
    return k_E
end

function metabolic_partitioning(V_c::Array{Float64,1})
    μ_0 = 4e7
    β_B = 1.64
    μ = μ_0*V_c.^(β_B-1)
    b = 2.42e-6
    γ = @. 1/(1 + b/μ)
    return @. 1.0/γ
end

function enzyme_production(α::Array{Float64,1}, genome_distr::Array{Float64,1})
  genome_distr_norm = genome_distr./sum(genome_distr)
  closure = α.*genome_distr_norm
end


text:
using Documenter, ITensorMPOCompression, ITensors

makedocs(sitename="ITensorMPOCompression",
        format = Documenter.HTML(prettyurls = false),
        modules = [ITensorMPOCompression]
        )


text:
```
Adapt Zygote and Tracker modes

```
track(m) = fmap(x -> x isa AbstractArray ? Tracker.param(x) : x, m)

function device(isgpu, args...)
    if isgpu
        return gpu(args)
    else
        return cpu(args)
    end
end


text:
using Taxonomy
using Test

db = Taxonomy.DB("db/nodes.dmp", "db/names.dmp")

@testset "taxon.jl" begin
    @test Taxon <: AbstractTaxon
    @test UnclassifiedTaxon <: AbstractTaxon

    human = Taxon(9606,db)
    @test typeof(human) == Taxon
    @test human.name == "Homo sapiens"
    @test human.taxid == 9606
    @test human.rank == :species
    @test sprint(io -> show(io, human)) == "9606 [species] Homo sapiens"

    @test get(db, 9606, nothing) == human

    @test taxid(human) == human.taxid
    @test rank(human) == human.rank
    @test parent(human) == Taxon(9605,db)
    @test children(human) == [Taxon(741158,db), Taxon(63221,db)]
    denisova = Taxon(741158, db)
    @test children(denisova) == Taxon[]
    @test isempty(children(denisova))

    unclassified_human_subspecies = UnclassifiedTaxon(:subspecies, human)
    @test typeof(unclassified_human_subspecies) == UnclassifiedTaxon
    @test unclassified_human_subspecies.name == "unclassified Homo sapiens subspecies"
    @test unclassified_human_subspecies.rank == :subspecies
    @test unclassified_human_subspecies.source == human
    @test sprint(io -> show(io, unclassified_human_subspecies)) == "Unclassified [subspecies] unclassified Homo sapiens subspecies"

    @test rank(unclassified_human_subspecies) == unclassified_human_subspecies.rank

    @test_throws KeyError Taxon(99999999, db)
    @test get(db, 9999999, nothing) === nothing
end

@testset "lineage.jl" begin
    human = Taxon(9606,db)
    lineage = Lineage(human)
    @test lineage[1] == Taxon(1,db)
    @test size(lineage) == (32,)
    @test lineage[32] == lineage[end] == lineage[:species] == human
    @test lineage[1:9] == lineage[Between(1, 9)]
    @test lineage[All()] == lineage
    @test lineage[Between(3, 29)] == lineage[Between(:superkingdom, 29)] == lineage[Between(3, :family)] == lineage[Between(:superkingdom, :family)]
    @test lineage[From(9)] == lineage[From(:phylum)] == lineage[9:32]
    @test lineage[Until(24)] == lineage[Until(:order)] == lineage[1:24]

    @test get(lineage, 2, nothing) == Taxon(131567,db)
    @test get(lineage, :class, nothing) == Taxon(40674, db)
    @test get(lineage, 99, nothing) === nothing
    @test get(lineage, :strain, nothing) === nothing
    @test sprint(io -> print_lineage(io, lineage)) == "root;cellular organisms;Eukaryota;Opisthokonta;Metazoa;Eumetazoa;Bilateria;Deuterostomia;Chordata;Craniata;Vertebrata;Gnathostomata;Teleostomi;Euteleostomi;Sarcopterygii;Dipnotetrapodomorpha;Tetrapoda;Amniota;Mammalia;Theria;Eutheria;Boreoeutheria;Euarchontoglires;Primates;Haplorrhini;Simiiformes;Catarrhini;Hominoidea;Hominidae;Homininae;Homo;Homo sapiens"
    @test sprint(io -> print_lineage(io, lineage; delim="+")) == "root+cellular organisms+Eukaryota+Opisthokonta+Metazoa+Eumetazoa+Bilateria+Deuterostomia+Chordata+Craniata+Vertebrata+Gnathostomata+Teleostomi+Euteleostomi+Sarcopterygii+Dipnotetrapodomorpha+Tetrapoda+Amniota+Mammalia+Theria+Eutheria+Boreoeutheria+Euarchontoglires+Primates+Haplorrhini+Simiiformes+Catarrhini+Hominoidea+Hominidae+Homininae+Homo+Homo sapiens"

    reformated_human_lineage = reformat(lineage,[:superkingdom,:phylum,:class,:order,:family,:genus,:species])
    @test reformated_human_lineage[1] == Taxon(2759, db)
    @test reformated_human_lineage[3] == Taxon(40674, db)
    @test reformated_human_lineage[7] == Taxon(9606, db)
    @test sprint(io -> print_lineage(io, reformated_human_lineage) ) == "Eukaryota;Chordata;Mammalia;Primates;Hominidae;Homo;Homo sapiens"

    primate = Taxon(9443,db)
    reformated_primate_lineage = reformat(Lineage(primate), [:superkingdom,:phylum,:class,:order,:family,:genus,:species])
    @test reformated_primate_lineage[7] == UnclassifiedTaxon(:species, primate)
    @test sprint(io -> print_lineage(io, reformated_primate_lineage)) == "Eukaryota;Chordata;Mammalia;Primates;;;"
    @test sprint(io -> print_lineage(io, reformated_primate_lineage; skip=true)) == "Eukaryota;Chordata;Mammalia;Primates"
    @test sprint(io -> print_lineage(io, reformated_primate_lineage; fill=true)) == "Eukaryota;Chordata;Mammalia;Primates;unclassified Primates family;unclassified Primates genus;unclassified Primates species"

    co = Taxon(131567,db)
    reformated_co_lineage = reformat(Lineage(co), [:superkingdom,:phylum,:class,:order,:family,:genus,:species])
    @test reformated_co_lineage[7] == UnclassifiedTaxon(:species, co)
    @test sprint(io -> print_lineage(io, reformated_co_lineage)) == ";;;;;;"
    @test sprint(io -> print_lineage(io, reformated_co_lineage; skip=true)) == ""
    @test sprint(io -> print_lineage(io, reformated_co_lineage; fill=true)) == "unclassified cellular organisms superkingdom;unclassified cellular organisms phylum;unclassified cellular organisms class;unclassified cellular organisms order;unclassified cellular organisms family;unclassified cellular organisms genus;unclassified cellular organisms species"

    human = Taxon(9606,db)
    primate = Taxon(9443,db)
    @test isdescendant(human, primate)
    @test isancestor(primate, human)
end

@testset "lca.jl" begin
    human = Taxon(9606, db)
    gorilla = Taxon(9593, db)
    pan = Taxon(9598, db)
    @test lca([human,gorilla]) == lca(human,gorilla) == Taxon(207598, db)
    @test lca([human,gorilla,pan]) == lca(human,gorilla,pan) == Taxon(207598, db)
end

@testset "tree.jl" begin
    human = Taxon(9606, db)
    gorilla = Taxon(9593, db)
    tree = topolgoy([human, gorilla])
    @test tree.node == Taxon(207598, db)
end


text:
using ChangePointMean
using Documenter

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

makedocs(;
    modules=[ChangePointMean],
    authors="Daymond Ling",
    repo="https://github.com/DaymondLing/ChangePointMean.jl/blob/{commit}{path}#{line}",
    sitename="ChangePointMean.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DaymondLing.github.io/ChangePointMean.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DaymondLing/ChangePointMean.jl",
)


text:
abstract type UncertainElement{C, F} end

struct δ{C, F<:Real} <: UncertainElement{C, F}
    val::C
    radius::F
    name::Symbol
    # function δr(val = 0, interval = Interval(val-1, val+1))
    #     lower ≤ val ≤ upper || throw(ArgumentError("val must be between lower and upper"))
    #     val, lower, upper = promote(val, lower, upper)
    #     new{typeof(val)}(val, lower, upper)
    # end
    function δ(val, radius, name = gensym("δ"))
        radius ≥ 0 || throw(ArgumentError("radius must be positive"))
        new{typeof(val), typeof(radius)}(val, radius, name)
    end
end

Base.zero(::Type{δ{C, F}}) where {C,F} = δ(zero(C), zero(F))

"""
    δr(val::Real = 0.0, radius::Real = 1.0, name)

Create a real, uncertain parameter. If no name is given, a boring name will be generated automatically.
"""
function δr(val::Real = 0.0, radius::Real = 1.0, args...)
    radius ≥ 0 || throw(ArgumentError("radius must be positive"))
    val, radius = promote(val, radius)
    δ(val, radius, args...)
end

"""
    δc(val::Complex = complex(0.0), radius::Real = 1.0, name)

Create a complex, uncertain parameter. If no name is given, a boring name will be generated automatically.
"""
function δc(val = complex(0.0), radius=1.0, args...)
    radius ≥ 0 || throw(ArgumentError("radius must be positive"))
    val = complex(val)
    δ(val, radius, args...)
end

Base.size(::δ) = (1,)
Base.size(::δ, n) = 1
Base.length(::δ) = 1
Base.eltype(::UncertainElement{C}) where C = C
Base.eltype(::Type{<:UncertainElement{C}}) where C = C

function Base.promote_rule(::Type{N}, ::Type{δ{C, F}}) where {N <: Number, F, C}
    δ{promote_type(N, C), promote_type(N, F)}
end

Base.convert(::Type{δ{C, F}}, n::Number) where {C,F} = δ(C(n), zero(F))
Base.convert(::Type{δ{C, F}}, d::δ) where {C,F} = δ(C(d.val), F(d.radius), d.name)

#  + and - we define here, all other ops generate ExtendedStateSpace. THis is to avoid creating a complicated tagging mechanism
function Base.:(+)(d1::δ, d2::δ)
    if d1.name === d2.name
        δ(+(d1.val, d2.val), +(d1.radius, d2.radius), d1.name)
    else
        uss(d1) + uss(d2)
    end
end

function Base.:(+)(n::Number, d2::δ)
    δ(+(n, d2.val), d2.radius, d2.name)
end

Base.:(+)(d2::δ, n::Number) = +(n, d2)

function Base.:(*)(n::Number, d2::δ)
    δ(*(n, d2.val), abs(n)*d2.radius, d2.name)
end

function Base.:(-)(d2::δ)
    δ(-d2.val, d2.radius, d2.name)
end

function Base.:(-)(d1::δ, d2::δ)
    # δr(d1.val+d2.val, d1.lower+d2.lower, d1.upper+d2.upper)
    if d1.name === d2.name
        δ(-(d1.val, d2.val), -(d1.radius, d2.radius), d1.name)
    else
        # δ(-(d1.val, d2.val), +(d1.radius, d2.radius))
        uss(d1) + (-uss(d2))
    end
end

function Base.:(*)(d1::δ, d2::δ)
    uss(d1)*uss(d2)
end

function Base.:(/)(d1::δ, d2::δ)
    uss(d1)*inv(uss(d2))
end

function Base.:(/)(d1::Number, d2::δ)
    d1*inv(uss(d2))
end

function Base.inv(d2::δ)
    inv(uss(d2))
end

"""
    UncertainSS{TE} <: AbstractStateSpace{TE}

Represents LFT_u(M, Diagonal(Δ))
"""
struct UncertainSS{TE} <: AbstractStateSpace{TE}
    sys::ExtendedStateSpace{TE}
    Δ
end

function Base.getproperty(sys::UncertainSS, s::Symbol)
    s ∈ fieldnames(typeof(sys)) && return getfield(sys, s)
    if s === :nz
        # return foldl((l, d)->l + size(d, 1), sys.Δ, init=0)
        return size(sys.C1, 1)
    elseif s === :nw
        # return foldl((l, d)->l + size(d, 2), sys.Δ, init=0)
        return size(sys.B1, 2)
    elseif s === :zinds
        return 1:sys.nz
    elseif s === :yinds
        return sys.nz .+ (1:size(sys.C2, 1))
    elseif s === :winds
        return 1:sys.nw
    elseif s === :uinds
        return sys.nw .+ (1:size(sys.B2, 2))
    elseif s ===:M
        return sminreal(performance_mapping(sys.sys))
        # return sys.sys
    elseif s ===:delta
        return sys.Δ
    end
    getproperty(getfield(sys, :sys), s)
end

for f in [:system_mapping, :performance_mapping, :ss]
    @eval $f(s::UncertainSS, args...; kwargs...) = $f(s.sys, args...; kwargs...)
end

function Base.promote_rule(::Type{UncertainSS{TE}}, ::Type{δ{C, F}}) where {TE, F, C}
    UncertainSS{TE}
end

function Base.promote_rule(::Type{UncertainSS{TE}}, ::Type{L}) where {TE, L <: LTISystem}
    UncertainSS{TE}
end

function Base.promote_rule(::Type{UncertainSS{TE}}, ::Type{L}) where {TE, L <: AbstractArray}
    UncertainSS{TE}
end

function Base.convert(::Type{U}, d::δ) where {U <: UncertainSS}
    uss(d)
end

function Base.convert(::Type{U}, d::AbstractArray) where {U <: UncertainSS}
    uss(d)
end

function Base.convert(::Type{U}, s::LTISystem) where {U <: UncertainSS}
    s isa UncertainSS && return s
    sys = ss(s)
    esys = partition(sys, 0, 0)
    UncertainSS(esys, [])
end

Base.:+(sys::UncertainSS{TE}, n::δ) where TE <: ControlSystems.TimeEvolution = +(sys, uss(n))
Base.:+(n::δ, sys::UncertainSS{TE}) where TE <: ControlSystems.TimeEvolution = +(uss(n), sys)

Base.:+(sys::UncertainSS{TE}, n::Number) where TE <: ControlSystems.TimeEvolution = +(sys, uss(n))
Base.:+(n::Number, sys::UncertainSS{TE}) where TE <: ControlSystems.TimeEvolution = +(uss(n), sys)

Base.:*(G::LTISystem, d::UncertainElement) = uss(ss(G)) * uss(d)
Base.:*(d::UncertainElement, G::LTISystem) = uss(d) * uss(ss(G))

"""
    uss(D11, D12, D21, D22, Δ, Ts = nothing)

Create an uncertain statespace object with only gin matrices.
"""
function uss(D11, D12, D21, D22, Δ, Ts=nothing)
    D11, D12, D21, D22 = vcat.((D11, D12, D21, D22))
    sys = ExtendedStateSpace(zeros(0,0), zeros(0,size(D21,2)), zeros(0,size(D12,2)), zeros(size(D12,1),0), zeros(size(D21,1),0), D11, D12, D21, D22, Ts)
    #                                       B1                   B2                    C1                        C2
    UncertainSS(sys, Δ)
end

"""
    uss(d::δ{C, F}, Ts = nothing)

Convert a δ object to an UncertainSS
"""
function uss(d::δ{C,F}, Ts = nothing) where {C,F}
    # sys = partition(ss([d.val d.radius; 1 0]), u=2, y=2, w=1, z=1)
    uss(0, 1, d.radius, d.val, [normalize(d)], Ts)
end

"""
    uss(d::AbstractVector{<:δ}, Ts = nothing)

Create a diagonal uncertain statespace object with the uncertain elements `d` on the diagonal.
"""
function uss(d::AbstractVector{<:δ}, Ts = nothing)
    vals = getfield.(d, :val)
    rads = getfield.(d, :radius)
    nonzero = rads .> 0
    # vals = vals
    # rads = rads
    nd = count(nonzero)
    n = length(d)
    uss(0I(nd), I(n)[nonzero, :], diagm(rads)[:, nonzero], diagm(vals), normalize.(d[nonzero]), Ts)
end

uss(n::Number, Ts=nothing) = uss(zeros(0,0), zeros(0,1), zeros(1,0), 1, [], Ts)

"""
    uss(D::AbstractArray, Δ, Ts = nothing)

If only a single `D` matrix is provided, it's treated as `D11` if Δ is given, and as `D22` if no Δ is provided.
"""
function uss(D::AbstractArray, Δ=[], Ts=nothing)
    if length(Δ) == size(D,1)
        uss(D, zeros(size(D,1),0), zeros(0,size(D,2)), zeros(0,0), Δ, Ts)
    elseif isempty(Δ)
        uss(zeros(0,0), zeros(0,size(D,2)), zeros(size(D,1),0), D, Δ, Ts)
    else
        throw(DimensionMismatch("length(Δ) != size(D,1)"))
    end
end

uss(s::UncertainSS) = s
uss(s::AbstractStateSpace) = convert(UncertainSS, s)

# Lower LFT
# function Base.inv(s::UncertainSS)
#     s.nx == 0 || error("Haven't considered this case yet")
#     @unpack D11, D12, D21, D22 = s
#     iD11 = inv(D11)
#     sys = uss(iD11, -iD11*D12, D21*iD11, D22-D21*iD11*D12, s.Δ, s.timeevol)
# end

function Base.inv(s::UncertainSS) # Upper LFT
    s.nx == 0 || error("Haven't considered this case yet")
    @unpack D11, D12, D21, D22 = s
    iD22 = inv(D22)
    sys = uss(D11-D12*iD22*D21, D12*iD22, -iD22*D21, iD22, s.Δ, s.timeevol)
end

function Base.:*(s1::UncertainSS, s2::UncertainSS)
    sys1 = s1.sys
    sys2 = s2.sys
    sys = invert_mappings(invert_mappings(sys1)*invert_mappings(sys2))
    UncertainSS(sys, [s1.Δ; s2.Δ])
end

function Base.:*(s1::UncertainSS, n::Number)  # reverse case is handled by switching
    sys = invert_mappings(n*invert_mappings(s1.sys))
    UncertainSS(sys, s1.Δ)
end

function Base.:+(s1::UncertainSS, s2::UncertainSS)
    sys1 = s1.sys
    sys2 = s2.sys
    sys = invert_mappings(invert_mappings(sys1)+invert_mappings(sys2))
    UncertainSS(sys, [s1.Δ; s2.Δ])
end

function Base.:-(s::UncertainSS)
    UncertainSS(-s.sys, s.Δ)
end

function Base.hcat(systems::UncertainSS...)
    # Perform checks
    timeevol = common_timeevol(systems...)

    ny = systems[1].ny
    if !all(s.ny == ny for s in systems)
        error("All systems must have same ouput dimension")
    end

    A = blockdiag([s.A for s in systems]...)

    B1 = blockdiag([s.B1 for s in systems]...)
    B2 = blockdiag([s.B2 for s in systems]...)

    C2 = reduce(hcat, [s.C2 for s in systems])
    C1 = blockdiag([s.C1 for s in systems]...)

    D21 = reduce(hcat, [s.D21 for s in systems])
    D22 = reduce(hcat, [s.D22 for s in systems])
    D11 = blockdiag([s.D11 for s in systems]...)
    D12 = blockdiag([s.D12 for s in systems]...)

    sysnew = ss(A, B1, B2, C1, C2, D11, D12, D21, D22, timeevol)
    return UncertainSS(sysnew, reduce(vcat, [s.Δ for s in systems]))
end

function Base.vcat(systems::UncertainSS...)
    # Perform checks
    timeevol = common_timeevol(systems...)

    nu = systems[1].nu
    if !all(s.nu == nu for s in systems)
        error("All systems must have same input dimension")
    end

    A = blockdiag([s.A for s in systems]...)

    B2 = reduce(vcat, [s.B2 for s in systems])
    B1 = blockdiag([s.B1 for s in systems]...)

    C1 = blockdiag([s.C1 for s in systems]...)
    C2 = blockdiag([s.C2 for s in systems]...)

    D12 = reduce(vcat, [s.D12 for s in systems])
    D22 = reduce(vcat, [s.D22 for s in systems])
    D11 = blockdiag([s.D11 for s in systems]...)
    D21 = blockdiag([s.D21 for s in systems]...)

    sysnew = ss(A, B1, B2, C1, C2, D11, D12, D21, D22, timeevol)
    return UncertainSS(sysnew, reduce(vcat, [s.Δ for s in systems]))
end

function ControlSystems.lft(G::UncertainSS, Δ::AbstractArray=G.Δ, type=:u)
    if ndims(Δ) == 1
        Δ = append(ss.(Δ)...)
    end

    if type !== :u
        error("Invalid type of lft ($type), specify type=:u")
    else
        lft(ss(G.sys), ss(Δ), :u)
    end
end

function ControlSystems.lft(G::UncertainSS, K::LTISystem, type=:l)
    if type === :u
        sys = lft(ss(G.sys), ss(K), :u)
    else
        nu1, ny1 = G.nw, G.nz
        sys = lft(ss(G.sys), ss(K), :l)
        sys = partition(sys, nu1, ny1)
        UncertainSS(sys, G.Δ)
    end
end

## Uncertain dynamics

struct δDyn{F, TE} <: UncertainElement{Complex{F}, F}
    ny::Int
    nu::Int
    nx::Int
    bound::F
    timeevol::TE
    name::Symbol
end

function Base.getproperty(d::δDyn, s::Symbol)
    s ∈ fieldnames(typeof(d)) && return getfield(d, s)
    if s === :Ts
        return d.timeevol.Ts
    else
        throw(ArgumentError("$(typeof(d)) has no property named $s"))
    end
end

function δDyn(ny::Int, nu::Int; bound=1, nx::Int = 2, Ts=0, name=gensym("δDyn"))
    timeevol = Ts > 0 ? Discrete(Ts) : Continuous()
    δDyn(ny,nu,nx,bound,timeevol,name)
end

Base.size(d::δDyn, n=:) = (d.ny, d.nu)[n]
# Base.length(d::δDyn) = d.ny + d.nu # NOTE: feels dangerous not to have prod(size)

function δss(ny, nu; bound=1, kwargs...)
    Δ = δDyn(ny, nu; bound, kwargs...)
    b = sqrt(bound)
    # D11 = Matrix([zeros(nu, ny) b*I(nu); b*I(ny) zeros(ny, nu)])
    # D22 = zeros(ny,nu)
    # D12 = zeros(ny+nu, nu)
    # D21 = zeros(ny, nu+ny)
    D11 = zeros(nu, ny)
    D22 = zeros(ny, nu)
    D12 = b*I(nu)
    D21 = b*I(ny)
    uss(D11, D12, D21, D22, [normalize(Δ)], Δ.timeevol)
end

## Sampling =======================================================
MonteCarloMeasurements.nominal(d::UncertainElement) = d.val

function Base.rand(d::δ{V}, N::Int) where V <: Real
    T = float(V)
    d.radius == 0 && return Particles{T, N}(fill(T(d.val), N))
    Particles(N, Uniform(d.val - d.radius, d.val + d.radius))
end

function Base.rand(D::AbstractArray{<:δ}, N::Int)
    M = rand.(D, N)
    T = eltype(M[1])
    T.(M)
end

function Base.rand(d::δ{Complex{F}}, N::Int) where F
    T = float(F)
    d.radius == 0 && return Particles(fill(T(d.val), N))
    r = Particles(N, Uniform(0, d.radius))
    θ = Particles(N, Uniform(0, 2pi))
    c = r*cis(θ)
    d.val == 0 ? c : c + d.val
end

function Base.rand(d::δDyn, N::Int)
    G = map(1:N) do i
        while true
            G = ssrand(d.ny, d.nu, d.nx; Ts = d.timeevol isa Continuous ? nothing : d.Ts)
            try
                n = norm(G, Inf)
                return (1/n) * G
            catch
            end
        end
    end
    ss2particles(G)
end

normalize(d::δ) = δ(0*d.val, one(d.radius), d.name)
normalize(d::δDyn) = δDyn(d.ny, d.nu, d.nx, 1, d.timeevol, d.name)

function Base.rand(s::UncertainSS, N::Int)
    Δ = rand.(normalize.(s.Δ), N)
    lft(s, Δ, :u)
end

## MCM

# function LinearAlgebra.det(D::AbstractMatrix{<:Complex{<:AbstractParticles}})
#     D0 = similar(D, ComplexF64)
#     parts = map(1:nparticles(D[1].re)) do i
#         for j in eachindex(D0)
#             D0[j] = Complex(D[j].re.particles[i], D[j].im.particles[i])
#         end
#         det(D0)
#     end
#     Complex(StaticParticles(getfield.(parts, :re)), StaticParticles(getfield.(parts, :im)))
# end

function ControlSystems.tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where T <: AbstractParticles
    bymap(tzeros, A, B, C, D)
end

using MonteCarloMeasurements: vecindex
function sys_from_particles(P, i)
    A,B,C,D = ssdata(P)
    ss(vecindex(A, i), vecindex(B, i), vecindex(C, i), vecindex(D, i))
end

# function any0det(D::Matrix{<:Complex{<:AbstractParticles}})
#     D0 = similar(D, ComplexF64)
#     maxre = maxim = -1
#     minre = minim = 1
#     for i = 1:nparticles(D[1].re)
#         for j in eachindex(D0)
#             D0[j] = Complex(D[j].re.particles[i], D[j].im.particles[i])
#         end
#         d = det(D0)
#         maxre = max(maxre, d.re)
#         minre = min(minre, d.re)
#         maxim = max(maxim, d.im)
#         minim = min(minim, d.im)
#         if maxre > 0 && minre < 0 && maxim > 0 && minim < 0
#             return true
#         end
#     end
#     false
# end

function ss2particles(G::Vector{<:AbstractStateSpace})
    pdp(x) = Particles(permutedims(x))
    A = reduce(hcat, vec.(getproperty.(G, :A))) |> pdp
    B = reduce(hcat, vec.(getproperty.(G, :B))) |> pdp
    C = reduce(hcat, vec.(getproperty.(G, :C))) |> pdp
    D = reduce(hcat, vec.(getproperty.(G, :D))) |> pdp
    @unpack nx,ny,nu = G[1]
    A = reshape(A, nx, nx)
    B = reshape(B, nx, nu)
    C = reshape(C, ny, nx)
    D = reshape(D, ny, nu)
    ss(A,B,C,D, G[1].timeevol)
end

function δr(N::Int)
    Particles(N, Uniform(-1, 1))
end

function δc(N::Int)
    r = Particles(N, Uniform(0, 1))
    θ = Particles(N, Uniform(0, 2pi))
    c = r*cis(θ)
end

Δ(n, δ) = Diagonal([δ() for _ in 1:n])

## Intervals
IntervalArithmetic.Interval(d::δ{R}) where R <: Real =
    Interval(d.val-d.radius, d.val+d.radius)

function IntervalArithmetic.Interval(d::δ{C}) where C <: Complex
    re = Interval(d.val.re-d.radius, d.val.re+d.radius)
    im = Interval(d.val.im-d.radius, d.val.im+d.radius)
    Complex(re, im)
end

"""
    block_structure(Δ)

Take a vector of uncertain elements and return a vector of vectors with the block
structure of perturbation blocks as described by μ-tools, i.e.
- `[-N, 0]` denotes a repeated real block of size `N`
- `[N, 0]` denotes a repeated complex block of size `N`
- `[ny, nu]` denotes a full complex block of size `ny × nu`
"""
function block_structure(D)
    perm = sortperm(D, by=d->d.name)
    D = D[perm]
    blocks = Vector{Int}[]
    blockstart = 1
    for i in 2:length(D)
        if D[i].name == D[i-1].name
            if blockstart == 0
                blockstart = i-1
            end
        else
            push!(blocks, makeblock(D[blockstart:i-1]))
            blockstart = i
        end
    end
    push!(blocks, makeblock(D[blockstart:length(D)]))
    blocks, perm
end

function makeblock(d) # TODO: full uncertain block that is not a dynamical system is not supported
    N = length(d)
    if N > 1 || (N==1 && d[] isa δ) # Δ = δI -> μ(M) = ρ(M) 8.80
        if eltype(d[1]) <: Complex
            b = [N, 0] # complex*I
        else
            b = [-N, 0] # real*I
        end
    else
        b = [size(d[], 1), size(d[], 2)]
    end
end

"""
    blocks, M = blocksort(P::UncertainSS)

Returns the block structure of `P.Δ` as well as `P.M` permuted according to the sorted block structure.
`blocks` is a vector of vectors with the block
structure of perturbation blocks as described by μ-tools, i.e.
- `[-N, 0]` denotes a repeated real block of size `N`
- `[N, 0]` denotes a repeated complex block of size `N`
- `[ny, nu]` denotes a full complex block of size `ny × nu`
"""
function blocksort(P::UncertainSS)
    blocks, perm = block_structure(P.Δ)
    Poutinds = []
    Pininds = []
    # Build indices into M that corresponds to the Δ blocks and permute those indices with the same permutation as was applied to Δ
    for (i, d) in enumerate(P.Δ)
        if i == 1
            push!(Poutinds, 1:size(d, 2))
            push!(Pininds, 1:size(d, 1))
        else
            push!(Poutinds, (1:size(d, 2)) .+ Poutinds[i-1][end])
            push!(Pininds, (1:size(d, 1)) .+ Pininds[i-1][end])
        end
    end
    Poutinds = reduce(vcat, Poutinds[perm])
    Pininds = reduce(vcat, Pininds[perm])
    M = P.M[Poutinds, Pininds]
    blocks, M
end

function getinds(blocks::Vector{Vector{Int}}, d::Int)
    lengths = [b[2] <= 1 ? b[1] : b[d] for b in blocks]
    endinds = cumsum(lengths)
    startinds = [1; endinds[1:end-1] .+ 1]
    inds = UnitRange.(0 .+ startinds, lengths .+ startinds .- 1)
end

# function getinds(D::AbstractVector, d::Int)
#     lengths = size.(D, d)
#     endinds = cumsum(lengths)
#     startinds = [1; endinds[1:end-1] .+ 1]
#     inds = UnitRange.(0 .+ startinds, lengths .+ startinds .- 1)
# end

# struct BlockDiag{T}
#     D::Vector{T}
#     rinds::Vector{UnitRange{Int}}
#     cinds::Vector{UnitRange{Int}}
# end
# BlockDiag(D::AbstractVector) = BlockDiag(D, getinds(D, 1), getinds(D, 2))

# Base.size(B::BlockDiag) = (B.rinds[end][end], B.cinds[end][end])

# function Base.:*(A, B::BlockDiag)
#     rinds = B.rinds
#     submats = map(enumerate(rinds)) do (i,inds)
#         A[:, inds] * B.D[i]
#     end
#     reduce(hcat, submats)
# end

# function quadform_minus(Q::BlockDiag, M, b, op)
#     # con = M'Q*M - b^2*Q ⪯ I(n)
#     rinds = Q.rinds
#     cinds = Q.cinds
#     map(eachindex(rinds)) do i
#         r = rinds[i]
#         c = cinds[i]
#         q = Q.D[i]
#         op(M[r,:]'q*M[c,:] .- b^2*q, I(size(M,2)))
#     end
# end

#=
D_from_Delta: I-full, full-I, diag-diag
Version that constructs from initial guess
=#

# See Skogestad 8.8.3
# for complex δ:
# Δ = δI -> μ(M) = ρ(M) 8.80
# Full block complex -> μ(M) = σ̄(M) = opnorm(M)
# In general ρ(M) <= μ(M) <= σ̄(M)

# function make_D(blocks)
#     D = []
#     for b in blocks
#         if b[2] == 0
#             n = abs(b[1])
#             T = b[1] > 0 ? ComplexF64 : Float64
#             push!(D, zeros(T, n, n))
#         else
#             push!(D, 0*I)
#         end
#     end
#     BlockDiagonal(D)
# end


text:
# Brusselator figure
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

function brusselator_f(x, p, t)
    u, v = x
    a, b = p
    du = 1 - (b+1)*u + a*u^2*v
    dv = b*u - a*u^2*v
    return SVector(du, dv)
end

a = 0.3; b = 1.5
p0 = [a, b]
u0 = [2.0, 2.0]

ds = ContinuousDynamicalSystem(brusselator_f, u0, p0)

Δt = 0.001
T = 100
tr = trajectory(ds, T; Δt, Ttr = 100)
t = 0:Δt:T
u, v = columns(tr)

period = estimate_period(u, :zerocrossing, t)
ip = round(Int, period/Δt)

fig = figure(figsize = (figx/2, figy))
ax = gca()

sc = ax.scatter(u[1:ip], v[1:ip]; c = range(0, 2π, length = ip), cmap = :twilight, s = 20)
cb = colorbar(sc)
cb.set_ticks([0, 2π])
cb.set_ticklabels(["\$0\$", "\$2\\pi\$"])
cb.set_label("\$\\phi\$", labelpad = -30)
ax.set_xticks([0.5, 2.0])
ax.set_yticks([3, 5])
ax.set_xlabel("\$u\$", labelpad = -20)
ax.set_ylabel("\$w\$", labelpad = -10)
ax.plot(1, b/a; marker = "o", mec = "C2", mew = 2,
markersize = 10, mfc = :white)

r = ip ÷ 10
ax.scatter(u[1:r:ip], v[1:r:ip]; color = "C2", marker = "D", s = 150, edgecolors="white")

fig.tight_layout(pad=0.3)

# wsave(plotsdir("2", "brusselator"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Adjacency{P,T}

The adjacency relation of rank `P` for topology of type `T`.
"""
struct Adjacency{P,T<:Topology} <: TopologicalRelation
  topology::T
end

Adjacency{P}(topology::T) where {P,T} = Adjacency{P,T}(topology)

# -------------------
# HALF-EDGE TOPOLOGY
# -------------------

function (𝒜::Adjacency{0,T})(vert::Integer) where {T<:HalfEdgeTopology}
  e = half4vert(vert, 𝒜.topology)

  # initialize result
  vertices = [e.half.head]

  # search in CCW orientation
  p = e.prev
  h = p.half
  while !isnothing(h.elem) && h != e
    push!(vertices, p.head)
    p = h.prev
    h = p.half
  end

  # if border edge is hit
  if isnothing(h.elem)
    # add last arm manually
    push!(vertices, p.head)

    # search in CW orientation
    h = e.half
    while !isnothing(h.elem)
      n = h.next
      h = n.half
      pushfirst!(vertices, h.head)
    end
  end

  vertices
end


text:
using HetaSimulator, Plots

p = load_platform("./alco");
m = p.models[:nameless]

res = sim(m; tspan = (0.0, 100.0), parameters = [:dose=>10.])
plot(res)

# load data

conditions_df = read_conditions("./alco/conditions.csv")
add_conditions!(p, conditions_df)

measurements_df = read_measurements("./alco/saturation_C1.csv")
add_measurements!(p, measurements_df)

p |> sim |> plot

# fitting

fit_params = [
    :Vd => 5.3,
    :kabs => 0.1,
    :Vmax => 0.2,
    :Km => 0.1,
    :sigma => 0.05
]

fit1 = fit(p, fit_params)

sim(p) |> plot # default
sim(p, parameters_upd = fit1.optim) |> plot # best fit

max_stars_repo_name: hetalang/qsp-benchmark

max_issues_repo_name: hetalang/qsp-benchmark

max_forks_repo_name: hetalang/qsp-benchmark


text:
module TestIR

using Test

if Base.JLOptions().check_bounds == 1
    # Run IR tests in a different process as `--check-bounds=yes` is
    # specified.
    path = joinpath(@__DIR__, "__test_ir.jl")
    code = """
    $(Base.load_path_setup_code())
    include($(repr(path)))
    """
    cmd = Base.julia_cmd()
    cmd = `$cmd --check-bounds=no`
    if Base.JLOptions().color == 1
        cmd = `$cmd --color=yes`
    end
    @test success(run(pipeline(`$cmd -e $code`; stdout=stdout, stderr=stderr)))
else
    include("__test_ir.jl")
end

end  # module


text:
type randomForest
    min_split::Int64
    min_gain::Float64
    max_depth::Integer
    n_estimators::Int64
    max_features::Union{Int64, String}
    feature_index::Dict{Int64, Vector}
    trees::Vector{ClassificationTree}
end

function randomForest(;
                      min_split::Int64 = 2,
                      min_gain::Float64 = 1e-7,
                      max_depth::Integer = 10000000,
                      n_estimators::Int64 = 200,
                      max_features::Union{Int64, String} = "nothing",
                      feature_index::Dict{Int64,Vector} = Dict{Int64,Vector}(),
                      trees::Vector{ClassificationTree} = Vector{ClassificationTree}()
                      )
    for i = 1:n_estimators
        push!(trees, ClassificationTree(min_gain = min_gain, min_samples_split = min_split, max_depth = max_depth))
    end
    return randomForest(min_split, min_gain, max_depth, n_estimators, max_features, feature_index, trees)
end

function train!(model::randomForest, X::Matrix, y::Vector)
    n_sample, n_feature = size(X)
    if model.max_features == "nothing"
        model.max_features = trunc(Int64, sqrt(n_feature))
    end
    sets = get_random_subsets(X, y, model.n_estimators)
    #println("Start training Random Forests")
    for i = 1:model.n_estimators
        #temp = i/model.n_estimators*100
        #println("Progress: $(temp) %")
        temp_data = sets[1,:,:]
        X = temp_data[:, 1:(end-1)]
        y = temp_data[:, end]
        idx = sample(1:n_feature, model.max_features, replace = false)
        model.feature_index[i] = idx
        X = X[:, idx]
        train!(model.trees[i], X, y)
    end
end

function predict(model::randomForest,
                 x::Matrix)
    n_sample = size(x,1)
    res = zeros(n_sample, model.n_estimators)
    for i = 1:model.n_estimators
        idx = model.feature_index[i]
        res[:,i] = predict(model.trees[i], x[:, idx])
    end
    ans = zeros(n_sample)
    for k = 1:n_sample
        most = 0
        max_time = 0
        for i in unique(res[k,:])
            n_times = 0
            for j = 1:model.n_estimators
                if i == res[j]
                    n_times += 1
                end
            end
            if n_times > max_time
                most = i
            end
        end
        ans[k] = most
    end
    return ans
end

function test_randomForest()
    X_train, X_test, y_train, y_test = make_iris()
    model = randomForest()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    println("classification accuracy ", accuracy(y_test, predictions))
    #pca
    pcamodel = PCA()
    train!(pcamodel, X_test)
    plot_in_2d(pcamodel, X_test, predictions, "RadomForest")
end


text:
## 2D Panels
#==========================================================================================#

abstract type AbstractPanel2D <: AbstractPanel end

struct Panel2D{T <: Real} <: AbstractPanel2D
    p1 :: SVector{2,T}
    p2 :: SVector{2,T}
end

struct WakePanel2D{T <: Real} <: AbstractPanel2D
    p1 :: SVector{2,T}
    p2 :: SVector{2,T}
end

Base.length(:: Panel2D) = 1

Panel2D(p1 :: FieldVector{2,T}, p2 :: FieldVector{2,T}) where T <: Real = Panel2D{T}(p1, p2)
WakePanel2D(p1 :: FieldVector{2,T}, p2 :: FieldVector{2,T}) where T <: Real = WakePanel2D{T}(p1, p2)

a :: AbstractPanel2D + b :: AbstractPanel2D = Panel2D(p1(a) + p1(b), p2(a) + p2(b))
a :: AbstractPanel2D - b :: AbstractPanel2D = Panel2D(p1(a) - p1(b), p2(a) - p2(b))

collocation_point(panel :: AbstractPanel2D, a = 0.5) = a * (p1(panel) + p2(panel))
panel_vector(panel :: AbstractPanel2D) = p2(panel) - p1(panel)
panel_length(panel :: AbstractPanel2D) = (norm ∘ panel_vector)(panel)

function transform_panel_points(panel_1 :: AbstractPanel2D, panel_2 :: AbstractPanel2D)
    x1, y1, x2, y2 = p1(panel_2), p2(panel_2)
    xs, ys = p1(panel_1)

    xp1, yp1 = affine_2D(x1, y1, xs, ys, panel_angle(panel_1))
    xp2, yp2 = affine_2D(x2, y2, xs, ys, panel_angle(panel_1))

    xp1, yp1, xp2, yp2
end

function transform_panel(panel :: AbstractPanel2D, point :: SVector{2,<: Real})
    xs, ys = p1(panel)
    affine_2D(first(point), last(point), xs, ys, panel_angle(panel))
end

panel_angle(panel :: AbstractPanel2D) = let (x, y) = panel_vector(panel); atan(y, x) end
panel_tangent(panel :: AbstractPanel2D) = rotation(1., 0., -panel_angle(panel))
panel_normal(panel :: AbstractPanel2D) = inverse_rotation(0., 1., panel_angle(panel))
panel_location(panel :: AbstractPanel2D) = let angle = panel_angle(panel); ifelse((π/2 <= angle <= π) || (-π <= angle <= -π/2), "lower", "upper") end

panel_points(panels) = [ p1.(panels); [(p2 ∘ last)(panels)] ]

reverse_panel(panel :: AbstractPanel2D) = Panel2D(panel.p2, panel.p1)

trailing_edge_panel(panels) = Panel2D((p2 ∘ last)(panels), (p1 ∘ first)(panels))

function wake_panel(panels, bound, α)
    firstx, firsty   = (p1 ∘ first)(panels)
    lastx, lasty     = (p2 ∘ last)(panels)
    y_mid            = (firsty + lasty) / 2
    y_bound, x_bound = bound .* sincos(α)
    WakePanel2D(SVector(lastx, y_mid), SVector(x_bound * lastx, y_bound * y_mid))
end

function wake_panels(panels, chord, length, num)
    _, firsty  = (p1 ∘ first)(panels)
    _, lasty   = (p2 ∘ last)(panels)
    y_mid      = (firsty + lasty) / 2
    bounds     = cosine_spacing(chord + length / 2, length, num)
    @. WakePanel2D(SVector(bounds[1:end-1], y_mid), SVector(bounds[2:end], y_mid))
end

function panel_scalar(scalar_func, strength, panel :: AbstractPanel2D, x, y)
    # Transform point to local panel coordinates
    xp, yp = transform_panel(panel, SVector(x, y))
    scalar_func(strength, xp, yp, 0., panel_length(panel))
end

function panel_scalar(scalar_func, strength, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D, point = 0.5)
    x, y = collocation_point(panel_i, point)
    panel_scalar(scalar_func, strength, panel_j, x, y)
end

function panel_velocity(velocity_func, strength, panel :: AbstractPanel2D, x, y)
    # Transform point to local panel coordinates
    xp, yp = transform_panel(panel, SVector(x, y))

    # Compute velocity in local panel coordinates
    u, w = velocity_func(strength, xp, yp, 0., panel_length(panel))

    # Transform velocity to original coordinates
    inverse_rotation(u, w, panel_angle(panel))
end

function panel_velocity(velocity_func, strength, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D, point = 0.5)
    x, y = collocation_point(panel_i, point)
    u, w = panel_velocity(velocity_func, strength, panel_j, x, y)
end

panel_velocity(f1, f2, str1, str2, panel :: AbstractPanel2D, x, y) = panel_velocity(f1, str1, panel, x, y) .+ panel_velocity(f2, str2, panel, x, y)

panel_velocity(f1, f2, str_j1, str_j2, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(f1, str_j1, panel_j, panel_i) .+ panel_velocity(f2, str_j2, panel_j, panel_i)

get_surface_values(panels, vals, surf = "upper", points = false) = partition(x -> (panel_location ∘ first)(x) == surf, (collect ∘ zip)(panels[2:end], vals), x -> ((first ∘ ifelse(points, p1, collocation_point) ∘ first)(x), last(x)))


text:
direct_tr = (NWI(), Dogleg())
direct_globals = (HZAW(), Backtracking(), Backtracking(;interp=FFQuadInterp()), NWI(), Dogleg())
inverse_ls = (HZAW(), Backtracking(), Backtracking(; interp=FFQuadInterp()))
inverse_globals = (HZAW(), Backtracking(), Backtracking(;interp=FFQuadInterp()))
run_all(td_obj, x0, BFGS(Inverse()), inverse_globals)
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals)
#run_all(td_obj, x0, SR1(Inverse()), inverse_ls)
run_all(td_obj, x0, DFP(Inverse()), inverse_globals)
run_all(td_obj, x0, BFGS(Direct()), direct_globals)
run_all(td_obj, x0, DBFGS(Direct()), direct_globals)
run_all(td_obj, x0, SR1(Direct()), direct_globals)
run_all(td_obj, x0, DFP(Direct()), direct_globals)

run_all(td_obj, x0, Newton(), direct_globals)
run_all(td_obj, x0, Newton(), direct_globals)

run_all(od_obj, x0, BFGS(Inverse()), inverse_globals)
run_all(od_obj, x0, DBFGS(Inverse()), inverse_globals)
run_all(od_obj, x0, SR1(Inverse()), inverse_ls)
run_all(od_obj, x0, DFP(Inverse()), inverse_globals)
run_all(od_obj, x0, BFGS(Direct()), direct_globals)
run_all(od_obj, x0, DBFGS(Direct()), direct_globals)
run_all(od_obj, x0, SR1(Direct()), direct_globals)
run_all(od_obj, x0, DFP(Direct()), direct_globals)

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

examples["Polynomial"] = OptimizationProblem("Polynomial",
                                             polynomial,
                                             polynomial_gradient!,
                                             nothing,
                                             polynomial_hessian!,
                                             nothing, # Constraints
                                             [0.0, 0.0, 0.0],
                                             [10.0, 7.0, 108.0],
                                             polynomial([10.0, 7.0, 108.0]),
                                             true,
                                             true)

##########################################################################
###
### Powell (d=4)
###
### Problem 35 in [1]
### Difficult since the hessian is singular at the optimum
##########################################################################

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);

##########################################################################
###
### Rosenbrock (2D)
###
### Problem 38 in [1]
###
### Saddle point makes optimization difficult
##########################################################################

td_obj = TwiceDiffed(rosenbrock!)

run_all(td_obj, x0, BFGS(Inverse()), inverse_globals[2:end]);
run_all(td_obj, x0, DBFGS(Inverse()), inverse_globals[2:end]);
run_all(td_obj, x0, SR1(Inverse()), inverse_ls[2:end]);
run_all(td_obj, x0, DFP(Inverse()), inverse_globals[2:end]);
run_all(td_obj, x0, BFGS(Direct()), direct_globals);
run_all(td_obj, x0, DBFGS(Direct()), direct_globals);
run_all(td_obj, x0, SR1(Direct()), direct_globals);
run_all(td_obj, x0, DFP(Direct()), direct_globals);
run_all(td_obj, x0, Newton(), direct_globals);


text:
module WAVEFORM

using Printf, LinearAlgebra, Random, JOLI, SparseArrays

# package code goes here
include("Operators/Operators.jl")
include("Krylov/Krylov.jl")
include("Utilities/Utilities.jl")
include("Modeling/Modeling.jl")
include("Optimization/Optimization.jl")
include("Multigrid/Multigrid.jl")
include("Interface/Interface.jl")

end # module


text:
# Copyright (c) 2019 Mathieu Besançon, Oscar Dowson, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

using Test

@testset "$(file)" for file in readdir(@__DIR__)
    if file == "runtests.jl" || !endswith(file, ".jl")
        continue
    end
    include(file)
end


text:
import ScikitLearnBase: BaseClassifier, BaseRegressor, predict, predict_proba,
                        fit!, get_classes, @declare_hyperparameters

################################################################################
# Classifier

"""
    DecisionTreeClassifier(; pruning_purity_threshold=nothing,
                           nsubfeatures::Int=0, maxdepth::Int=-1,
                           rng=Base.GLOBAL_RNG)
Decision tree classifier. See [DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:

- `pruning_purity_threshold`: merge leaves having `>=thresh` combined purity (default: no pruning)
- `nsubfeatures`: number of features to select at random (default: keep all)
- `maxdepth`: maximum depth of the decision tree (default: no maximum)
- `rng`: the random number generator to use. Can be an `Int`, which will be used
  to seed and create a new random number generator.

Implements `fit!`, `predict`, `predict_proba`, `get_classes`
"""
type DecisionTreeClassifier <: BaseClassifier
    pruning_purity_threshold::Nullable{Float64} # no pruning if null
    # Does nsubfeatures make sense for a stand-alone decision tree?
    nsubfeatures::Int
    maxdepth::Int
    rng::AbstractRNG
    root::LeafOrNode
    # classes (in scikit-learn) === labels (in DecisionTree.jl)
    classes::Vector   # an arbitrary ordering of the labels
    # No pruning by default
    DecisionTreeClassifier(;pruning_purity_threshold=nothing, nsubfeatures=0,
                           maxdepth=-1, rng=Base.GLOBAL_RNG) =
        new(convert(Nullable{Float64}, pruning_purity_threshold), nsubfeatures,
            maxdepth, mk_rng(rng))
end

get_classes(dt::DecisionTreeClassifier) = dt.classes
@declare_hyperparameters(DecisionTreeClassifier,
                         [:pruning_purity_threshold, :nsubfeatures, :maxdepth,
                          :rng])

function fit!(dt::DecisionTreeClassifier, X, y)
    dt.root = build_tree(y, X, dt.nsubfeatures, dt.maxdepth; rng=dt.rng)
    if !isnull(dt.pruning_purity_threshold)
        dt.root = prune_tree(dt.root, get(dt.pruning_purity_threshold))
    end
    dt.classes = sort(unique(y))
    dt
end

predict(dt::DecisionTreeClassifier, X) = apply_tree(dt.root, X)

predict_proba(dt::DecisionTreeClassifier, X) =
    apply_tree_proba(dt.root, X, dt.classes)

predict_log_proba(dt::DecisionTreeClassifier, X) =
    log(predict_proba(dt, X)) # this will yield -Inf when p=0. Hmmm...

################################################################################
# Regression

"""
    DecisionTreeRegressor(; pruning_purity_threshold=nothing,
                          maxlabels::Int=5,
                          nsubfeatures::Int=0,
                          rng=Base.GLOBAL_RNG)
Decision tree regression. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:

- `pruning_purity_threshold`: merge leaves having `>=thresh` combined purity (default: no pruning)
- `maxlabels`: maximum number of samples per leaf, split leaf if exceeded
- `nsubfeatures`: number of features to select at random (default: keep all)
- `maxdepth`: maximum depth of the decision tree (default: no maximum)
- `rng`: the random number generator to use. Can be an `Int`, which will be used
  to seed and create a new random number generator.

Implements `fit!`, `predict`, `get_classes`
"""
type DecisionTreeRegressor <: BaseRegressor
    pruning_purity_threshold::Nullable{Float64}
    maxlabels::Int
    nsubfeatures::Int
    maxdepth::Int
    rng::AbstractRNG
    root::LeafOrNode
    # No pruning by default (I think purity_threshold=1.0 is a no-op, maybe
    # we could use that)
    DecisionTreeRegressor(;pruning_purity_threshold=nothing, maxlabels=5,
                          nsubfeatures=0, maxdepth=-1, rng=Base.GLOBAL_RNG) =
        new(convert(Nullable{Float64}, pruning_purity_threshold), maxlabels,
            nsubfeatures, maxdepth, mk_rng(rng))
end

@declare_hyperparameters(DecisionTreeRegressor,
                         [:pruning_purity_threshold, :maxlabels, :nsubfeatures,
                          :maxdepth, :rng])

function fit!{T<:Real}(dt::DecisionTreeRegressor, X::Matrix, y::Vector{T})
    # build_tree knows that its a regression problem by its argument types. I'm
    # not sure why X has to be Float64, but the method signature requires it
    # (as of April 2016).
    dt.root = build_tree(y, convert(Matrix{Float64}, X), dt.maxlabels,
                         dt.nsubfeatures, dt.maxdepth; rng=dt.rng)
    if !isnull(dt.pruning_purity_threshold)
        dt.root = prune_tree(dt.root, get(dt.pruning_purity_threshold))
    end
    dt
end

predict(dt::DecisionTreeRegressor, X) = apply_tree(dt.root, X)

################################################################################
# Random Forest Classification

"""
    RandomForestClassifier(; nsubfeatures::Int=0,
                           ntrees::Int=10,
                           partialsampling=0.7,
                           maxdepth=-1,
                           rng=Base.GLOBAL_RNG)
Random forest classification. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:

- `nsubfeatures`: number of features to select in each tree at random (default:
  keep all)
- `ntrees`: number of trees to train
- `partialsampling`: fraction of samples to train each tree on
- `maxdepth`: maximum depth of the decision trees (default: no maximum)
- `rng`: the random number generator to use. Can be an `Int`, which will be used
  to seed and create a new random number generator.

Implements `fit!`, `predict`, `predict_proba`, `get_classes`
"""
type RandomForestClassifier <: BaseClassifier
    nsubfeatures::Int
    ntrees::Int
    partialsampling::Float64
    maxdepth::Int
    rng::AbstractRNG
    ensemble::Ensemble
    classes::Vector
    RandomForestClassifier(; nsubfeatures=0, ntrees=10, partialsampling=0.7,
                           maxdepth=-1, rng=Base.GLOBAL_RNG) =
        new(nsubfeatures, ntrees, partialsampling, maxdepth, mk_rng(rng))
end

get_classes(rf::RandomForestClassifier) = rf.classes
@declare_hyperparameters(RandomForestClassifier,
                         [:nsubfeatures, :ntrees, :partialsampling, :maxdepth,
                          :rng])

function fit!(rf::RandomForestClassifier, X::Matrix, y::Vector)
    rf.ensemble = build_forest(y, X, rf.nsubfeatures, rf.ntrees,
                               rf.partialsampling, rf.maxdepth; rng=rf.rng)
    rf.classes = sort(unique(y))
    rf
end

predict_proba(rf::RandomForestClassifier, X) =
    apply_forest_proba(rf.ensemble, X, rf.classes)

predict(rf::RandomForestClassifier, X) = apply_forest(rf.ensemble, X)

################################################################################
# Random Forest Regression

"""
    RandomForestRegressor(; nsubfeatures::Int=0,
                          maxlabels::Int=5,
                          ntrees::Int=10,
                          partialsampling=0.7,
                          maxdepth=-1,
                          rng=Base.GLOBAL_RNG)
Random forest regression. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:

- `nsubfeatures`: number of features to select in each tree at random (default:
  keep all)
- `maxlabels`: maximum number of samples per leaf, split leaf if exceeded
- `ntrees`: number of trees to train
- `partialsampling`: fraction of samples to train each tree on
- `maxdepth`: maximum depth of the decision trees (default: no maximum)
- `rng`: the random number generator to use. Can be an `Int`, which will be used
  to seed and create a new random number generator.

Implements `fit!`, `predict`, `get_classes`
"""
type RandomForestRegressor <: BaseRegressor
    nsubfeatures::Int
    maxlabels::Int
    ntrees::Int
    partialsampling::Float64
    maxdepth::Int
    rng::AbstractRNG
    ensemble::Ensemble
    RandomForestRegressor(; nsubfeatures=0, ntrees=10, maxlabels=5, partialsampling=0.7, maxdepth=-1, rng=Base.GLOBAL_RNG) =
        new(nsubfeatures, maxlabels, ntrees, partialsampling, maxdepth,
            mk_rng(rng))
end

@declare_hyperparameters(RandomForestRegressor,
                         [:nsubfeatures, :ntrees, :maxlabels, :partialsampling,
                          # I'm not crazy about :rng being a hyperparameter,
                          # since it'll change throughout fitting, but it works
                          :maxdepth, :rng])

function fit!{T<:Real}(rf::RandomForestRegressor, X::Matrix, y::Vector{T})
    rf.ensemble = build_forest(y, convert(Matrix{Float64}, X), rf.nsubfeatures,
                               rf.ntrees, rf.maxlabels, rf.partialsampling,
                               rf.maxdepth; rng=rf.rng)
    rf
end

predict(rf::RandomForestRegressor, X) = apply_forest(rf.ensemble, X)

################################################################################
# AdaBoost Stump Classifier

"""
    AdaBoostStumpClassifier(; niterations::Int=0)

Adaboosted decision tree stumps. See
[DecisionTree.jl's documentation](https://github.com/bensadeghi/DecisionTree.jl)

Hyperparameters:

- `niterations`: number of iterations of AdaBoost
- `rng`: the random number generator to use. Can be an `Int`, which will be used
  to seed and create a new random number generator.

Implements `fit!`, `predict`, `predict_proba`, `get_classes`
"""
type AdaBoostStumpClassifier <: BaseClassifier
    niterations::Int
    rng::AbstractRNG
    ensemble::Ensemble
    coeffs::Vector{Float64}
    classes::Vector
    AdaBoostStumpClassifier(; niterations=10, rng=Base.GLOBAL_RNG) =
        new(niterations, mk_rng(rng))
end
@declare_hyperparameters(AdaBoostStumpClassifier, [:niterations, :rng])
get_classes(ada::AdaBoostStumpClassifier) = ada.classes

function fit!(ada::AdaBoostStumpClassifier, X, y)
    ada.ensemble, ada.coeffs =
        build_adaboost_stumps(y, X, ada.niterations; rng=ada.rng)
    ada.classes = sort(unique(y))
    ada
end

predict(ada::AdaBoostStumpClassifier, X) =
    apply_adaboost_stumps(ada.ensemble, ada.coeffs, X)

predict_proba(ada::AdaBoostStumpClassifier, X) =
    apply_adaboost_stumps_proba(ada.ensemble, ada.coeffs, X, ada.classes)

################################################################################
# Common functions

# Cannot use Union{} in Julia 0.3
depth(dt::DecisionTreeClassifier) = depth(dt.root)
depth(dt::DecisionTreeRegressor) = depth(dt.root)


text:
using newtonroot
using Test

@testset "newtonroot.jl" begin
    # Write your own tests here.
end


text:
module SobolSearches

export SobolSearch, setup_sobol_search, sobol_search_next!, sobol_search!, sobol_search

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES
using Parameters: @unpack
import Sobol

struct SobolSearch{TM,TF,TV,TS} <: AbstractVector{TM}
    minima::Vector{TM}
    minimand::TF
    lower_bounds::TV
    upper_bounds::TV
    sobol_seq::TS
end

function _transform_position(lower_bounds, upper_bounds, s)
    map((l, u, α) -> l*(1-α) + u*α, lower_bounds, upper_bounds, s)
end

function _evaluated_position(minimand, lower_bounds, upper_bounds, s)
    position = _transform_position(lower_bounds, upper_bounds, s)
    (value = minimand(position), position = position)
end

function setup_sobol_search(minimand, lower_bounds, upper_bounds, keep;
                            skip = 0)
    @argcheck length(lower_bounds) == length(upper_bounds)
    @argcheck all(lower_bounds .< upper_bounds)
    @argcheck keep ≥ 1
    sobol_seq = Sobol.SobolSeq(length(lower_bounds))
    skip > 0 && Sobol.skip(sobol_seq, skip)
    minima = [_evaluated_position(minimand, lower_bounds, upper_bounds, Sobol.next!(sobol_seq))
              for _ in 1:keep]
    sort!(minima, by = x -> x.value)
    SobolSearch(minima, minimand, lower_bounds, upper_bounds, sobol_seq)
end

function sobol_search_next!(search::SobolSearch)
    @unpack lower_bounds, upper_bounds, minimand, sobol_seq, minima = search
    candidate = _evaluated_position(minimand, lower_bounds, upper_bounds,
                                    Sobol.next!(sobol_seq))
    len = length(minima)
    if candidate.value < minima[len].value
        n = searchsortedfirst(minima, candidate, by = x -> x.value)
        for i in (len - 1):-1:n
            minima[i + 1] = minima[i]
        end
        minima[n] = candidate
        true
    else
        false
    end
end

function sobol_search!(search::SobolSearch, steps)
    for _ in 1:steps
        sobol_search_next!(search)
    end
    search
end

"""
$(SIGNATURES)

Find the `keep` lowers minima of `miniman` along a Sobol sequence of length `steps` in the
hypercube between `lower_bounds` and `upper_bounds`.

Return a [`SobolSearch`](@ref) object.

```julia
```
"""
function sobol_search(minimand, lower_bounds, upper_bounds, keep, steps; skip = 0)
    search = setup_sobol_search(minimand, lower_bounds, upper_bounds, keep;
                                skip = skip)
    sobol_search!(search, steps)
end

Base.IndexStyle(::Type{<:SobolSearch{TM}}) where {TM} = Base.IndexStyle(Vector{TM})

Base.size(search::SobolSearch) = size(search.minima)

Base.getindex(search::SobolSearch, I) = getindex(search.minima, I)

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using LinearAlgebra, BenchmarkTools

A = randn(10, 1000)
B = randn(1000, 1000)
C = randn(1000, 10000)
P = zeros(10, 1000)
Q = zeros(10, 10000)
R = zeros(1000, 10000)
S = zeros(10, 10000)

"""(A * B) * C"""
function f!(A, B, C, P, Q)
    mul!(P, A, B)
    mul!(Q, P, C)
end

"""A * (B * C)"""
function g!(A, B, C, R, S)
    mul!(R, B, C)
    mul!(S, A, R)
end

f!(A, B, C, P, Q) ≈ g!(A, B, C, R, S)

# %%
@benchmark f!($A, $B, $C, $P, $Q)

# %%
@benchmark g!($A, $B, $C, $R, $S)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module TestExtras

using Maybe.Extras
using Test

@testset "getnested" begin
    d = Dict(:a => Dict(:b => 1, :c => nothing))
    @test getnested(d, :a, :b) === Some(1)
    @test getnested(d, :a, :c) === Some(nothing)
    @test getnested(d, :a, :d) === nothing
end

@testset "ifnothing" begin
    @test nothing |> ifnothing() do
        0
    end == 0
    @test Some(1) |> ifnothing() do
        0
    end == 1
end

@testset "defaultto" begin
    @test nothing |> defaultto(0) == 0
    @test Some(1) |> defaultto(0) == 1
end

end  # module


text:
print("testing")
a = Base.prompt("Press Any key to coninue")

print("test1")
b = Base.prompt("press again")

print("aaaaaaaaaa")

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
const Symbolic = Union{Expr, Symbol}
const Numeric = Union{Number, Array}

⊕(ex::Symbolic, v::Numeric) = :($ex .+ $v)
⊕(v::Numeric, ex::Symbolic) = :($v .+ $ex)
⊕(ex1::Symbolic, ex2::Symbolic) = :($ex1 .+ $ex2)
⊕(x, y) = x .+ y

⊗(ex::Symbolic, v::Numeric) = :($ex .* $v)
⊗(v::Numeric, ex::Symbolic) = :($v .* $ex)
⊗(ex1::Symbolic, ex2::Symbolic) = :($ex1 .* $ex2)
⊗(x, y) = x .* y


text:
function pow(base::Number, exp::Integer)
  r = one(base)
  for i = 1:exp
    r *= base
  end
  return r
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SimpleML

using Optim
using OptionsMod

export predictLogistic, predictLogisticAll, learnLogistic, learnLogisticAll, predict2nn, learn2nn

sigmoid(x) = (1 + exp (-1 * x)) .^ -1

sigmoidGradient(z) = sigmoid(z) .* (1-sigmoid(z))

#Cost function for a logistic regression, with regularisation
function logisticCost(theta, X, y, grad, lambda)
    m = length(y) #size of training set
    n = length(theta) #number of features

    h = sigmoid(X * theta)

    sum =0
    for i = 1:m
        sum = sum + ( -1 * y[i] * log(h[i]) - (1-y[i]) * log (1-h[i]))
    end
    for j=2:n
        sum = sum + (lambda/2) * (theta[j] ^ 2)
    end

    if !(grad === nothing)
        gptr = pointer(grad)
        for j=1:n
            grad[j]=0
            for i = 1:m
                grad[j] =  grad[j]  + (h[i] - y[i])*X[i,j]
            end

           if j==1
                grad[j] = grad[j]/m
           else
                grad[j] = ( grad[j] + lambda* theta[j]) / m
           end
        end
        @assert gptr == pointer(grad)
    end

    return sum / m
end

#Learn via a logistic regression of two states
#  X : features
#  y : lables (binary)
#  lambda : regularisation parameter
function learnLogistic(X, y, lambda)
    initial_theta = zeros(size(X,2))
    ops = @options itermax=50 tol=1e-5

    results, fval, fcount, converged = cgdescent((g,t)->logisticCost(t, X, y, g, lambda),  initial_theta, ops)
    #println("Converged: $converged @ $(fval[end]) with coefficients $results in $fcount iterations. ")

    return results
end

#Learn via a logistic regression
#  X : features
#  y : lables
#  nlables : number of lables
#  lambda : regularisation parameter
function learnLogisticAll(X, y, lambda, nlables)

    m = size(X, 1); #number of data points
    n = size(X, 2); #number of features

    all_theta = zeros(nlables, n)

    for c=1:nlables
        initial_theta = zeros(n , 1)
        all_theta[c,:]=learnLogistic(X,int((y.==c)),lambda)[:]
    end

     return all_theta

end

function predictLogistic(theta, X )
    m = size(X, 1)
    p = zeros(m, 1)
    h = sigmoid(X * theta)

    for i=1:m
        if h[i] < 0.5
            p[i] = 0
        else
            p[i] = 1
        end
    end
    return p
end

function predictLogisticAll(all_theta, X)
    m = size(X, 1) #Number of training examples
    n = size(X, 2) #Number of features
    nlables = size(all_theta, 1)
    all_h = zeros(nlables, m)
    for c = 1: nlables
        theta = all_theta[c,:]
        h = sigmoid(X * theta')
        all_h[c,:] = h[:]
    end

    p = zeros(m)

    for i = 1:m
        _ , p[i] = findmax(all_h[:,i])
    end

    return p
end

#Learn a 2 layer neural network
#  X : features
#  y : lables
#  input_sz: size of input layer
#  hidden_sz: size of hidden layer
#  nlables : number of lables
#  lambda : regularisation parameter
function learn2nn(X, y, input_sz, hidden_sz, nlables, lambda)

    initial_theta1 = SimpleML.randomInitialWeights(input_sz, hidden_sz)
    initial_theta2 = SimpleML.randomInitialWeights(hidden_sz, nlables)
    initial_params = [initial_theta1[:] ; initial_theta2[:]]

    ops = @options itermax=50 tol=1e-5

    results, fval, fcount, converged = cgdescent((g,t)->nn2Cost(t, X, y,input_sz, hidden_sz, nlables, lambda, g),  initial_params, ops)

    #println("Converged: $converged @ $(fval[end]) with coefficients $results in $fcount iterations. ")

    theta1 = reshape(results[1:hidden_sz * (input_sz + 1)], hidden_sz, (input_sz + 1))
    theta2 = reshape(results[(1 + (hidden_sz * (input_sz + 1))):end], nlables, (hidden_sz + 1))
    return theta1 ,theta2

end

#Predict a 2 layer neural network
function predict2nn(theta1, theta2, X)
    m = size(X, 1)
    nlables = size(theta2, 1)
    p = zeros(size(X, 1), 1)
    X = [ones(m, 1) X]
    a2 = sigmoid(X * theta1');
    a2 = [ones(m, 1) a2]
    a3 = sigmoid (a2 * theta2')

    for i=1:m
        _, p[i] = findmax(a3[i,:])
    end

    return p
end

#Cost function for a 2 layer neural network
function nn2Cost(params, X, y, input_sz,  hidden_sz, nlables,  lambda, grad)

    theta1 = reshape(params[1:hidden_sz * (input_sz + 1)], hidden_sz, (input_sz + 1))
    theta2 = reshape(params[(1 + (hidden_sz * (input_sz + 1))):end], nlables, (hidden_sz + 1))

    m = size(X, 1)

    J = 0
    theta1_grad = zeros(size(theta1))
    theta2_grad = zeros(size(theta2))

    X = [ones(m, 1) X]

    z2 = X * theta1'
    a2 = sigmoid(z2)

    a2 = [ones(m, 1) a2]
    a3 = sigmoid (a2 * theta2')
    s=0
    for i=1:m
        yvector = int([1:nlables] .== y[i])
        s=s+sum(-yvector' .* log(a3[i,:]) - (1-yvector') .* log (1-a3[i,:]));
    end

    temp1 = theta1
    temp2 = theta2
    temp1[:, 1] = zeros(size(temp1, 1), 1)
    temp2[:, 1] = zeros(size(temp2, 1), 1)

    if !(grad === nothing)
        gptr = pointer(grad)
        for t=1:m
            yvector = int([1:nlables] .== y[t])

            delta3 = zeros(1,nlables)
            delta3 = a3[t,:] - yvector'

            delta2 = zeros (1, hidden_sz)
            delta2 = (theta2' * delta3' )[2:end] .* sigmoidGradient(z2[t,:])'

            theta2_grad = theta2_grad + delta3' * a2[t,:]
            theta1_grad = theta1_grad + delta2 * X[t,:]
        end

        theta1_grad = theta1_grad / m
        theta2_grad = theta2_grad / m

        theta1_grad = theta1_grad + (lambda/m) * temp1
        theta2_grad = theta2_grad + (lambda/m) * temp2

        grad[:] = [theta1_grad[:] ; theta2_grad[:]]
        @assert gptr == pointer(grad)

    end
    return s/m  + (lambda/(2*m))*(sum(sum(temp1 .* temp1)) + sum(sum(temp2 .* temp2)))
end

function randomInitialWeights(l_in, l_out)
    ep = 0.12
    rand(l_out, 1 + l_in) * 2 * ep - ep
end

end


text:
struct TimestampArrayTableRowsIterator{D, T, TP}
    m::TimestampArray{D, T}
    names::TP
end

Base.length(m::TimestampArrayTableRowsIterator) = length(m.m.tsps)
function Base.iterate(m::TimestampArrayTableRowsIterator, i = 1)
    i > length(m) && return nothing
    tsp = @inbounds m.m.tsps[i]
    return (NamedTuple{m.names}((tsp.ts, tsp.val...)), i+1)
end

struct TimestampArrayTableColsIterator{D, T}
    mapping::Dict{Symbol, Int}
    m::TimestampArray{D, T}
end

function Base.iterate(m::TimestampArrayTableColsIterator, i = 1)
    i > length(m.m.names) + 1 && return nothing
    return (Tables.getcolumn(m, i), i+1)
end

Tables.istable(::Type{<:TimestampArray}) = true
Tables.columnnames(m::TimestampArray) = (:timestamp, m.names...)
Tables.rowaccess(::Type{<:TimestampArray}) = true
Tables.rows(m::TimestampArray) = TimestampArrayTableRowsIterator(m, Tables.columnnames(m))
Tables.columnaccess(::Type{<:TimestampArray}) = true
Tables.columns(m::TimestampArray) = TimestampArrayTableColsIterator(Dict(:timestamp => 1, (m.names .=> 2:length(m.names)+1)...), m)
Tables.schema(m::TimestampArray) = Tables.Schema((:timestamp, m.names...), nothing)
Tables.schema(m::TimestampArray{D, T}) where {D, T <: Tuple} = Tables.Schema((:timestamp, m.names...), Tuple{D, T.types...})
Tables.schema(m::TimestampArrayTableRowsIterator) = Tables.schema(m.m)

function Tables.getcolumn(m::TimestampArray{D, T}, i::Int) where {D, T <: Tuple}
    i == 1 && return map(x -> x.ts, m.tsps)
    return map(x -> x.val[i - 1], m.tsps)
end
function Tables.getcolumn(m::TimestampArray, nm::Symbol) where {D, T <: Tuple}
    nm == :timestamp && return map(x -> x.ts, m.tsps)
    i = findfirst(==(nm), m.names)
    return map(x -> x.val[i], m.tsps)
end

Tables.columnnames(mci::TimestampArrayTableColsIterator) = Tables.columnnames(mci.m)
Tables.getcolumn(mci::TimestampArrayTableColsIterator, i::Int) = Tables.getcolumn(mci.m, i)
Tables.getcolumn(mci::TimestampArrayTableColsIterator, nm::Symbol) = Tables.getcolumn(mci.m, mci.mapping[nm])

########################################
# Constructor
########################################

# Ok, this one is experimental and probably should be moved to a different package
istimetable(::T) where T = istimetable(T)
istimetable(::Type{T}) where T = false
timeaxis(::T) where T = timeaxis(T)
timeaxis(::Type{T}) where T = :timestamp

istimetable(::Type{<:TimestampArray}) = true
timeaxis(x::TimestampArray) = :timestamp

function TimestampArray(x; timestamp = :timestamp)
    names = Tables.columnnames(x)
    tscol = istimetable(x) ? timeaxis(x) : timestamp # fallback
    indx = Vector{Int}(undef, length(names) - 1)
    k = 0
    for (i, nm) in pairs(names)
        nm == tscol && continue
        k += 1
        indx[k] = i
    end
    rows = Tables.rows(x)
    tsps = map(rows) do row
        rt = ntuple(i -> Tables.getcolumn(row, indx[i]), length(row) - 1)
        Timestamp(Tables.getcolumn(row, tscol), rt)
    end

    return TimestampArray(tsps, collect(names[indx]))
end


text:
#module A
#export Astruct

struct Astruct
end

#end


text:
"Tableau for the explicit 1-stage stochastic Euler method"
function TableauStochasticEuler()
  a = zeros(Float64, 1, 1)
  b = [1.0]
  c = [0.0]

  TableauSERK(:Stochastic_Euler_explicit_method, 1, a, b, c, 1, a, b, c)
end

"Tableau for the explicit 2-stage stochastic Heun method"
function TableauStochasticHeun()

    a = [[0.0 0.0]
         [1.0 0.0]]
    b = [0.5, 0.5]
    c = [0.,  1.]

    TableauSERK(:Stochastic_Heun_explicit_method, 2, a, b, c, 2, a, b, c)
end

"""
Tableau for the explicit Platen method
  Platen's method cited in Eq. (52) in
   K. Burrage, P. Burrage (1996) "High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations".
  According to the paper, the method has strong order 1.0 for one-dimensional Brownian motion.
  Appears to have a rather poor long-time performance.
"""
function TableauPlaten()

    a_drift = [[0.0 0.0]
               [1.0 0.0]]
    b_drift = [1.0, 0.0]
    c_drift = [0.0, 1.0]

    a_diff = [[0.0 0.0]
              [1.0 0.0]]
    b_diff = [0.5, 0.5]
    c_diff = [0.0, 1.0]

    TableauSERK(:Platen_explicit_method, 1, a_drift, b_drift, c_drift, 2, a_diff, b_diff, c_diff)
end

"""
Tableau for the explicit 2-stage R2 method due to K. Burrage and P. Burrage
  Method cited in Eq. (51) in
   K. Burrage, P. Burrage (1996) "High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations".
  According to the paper, the method has strong order 1.0 for one-dimensional Brownian motion
"""
function TableauBurrageR2()

    a = [[0.0 0.0]
         [2. / 3. 0.0]]
    b = [0.25, 0.75]
    c = [0., 2. / 3.]

    TableauSERK(:BurrageR2_explicit_method, 2, a, b, c, 2, a, b, c)
end

"""
Tableau for the explicit 4-stage CL method due to K. Burrage and P. Burrage
  Method cited in Eq. (56) in
   K. Burrage, P. Burrage (1996) "High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations".
  According to the paper, the method has strong order 1.5 for one-dimensional Brownian motion. Reduces to the classical R-K method
  of order 4 when noise is zero.
"""
function TableauBurrageCL()

    a_drift = [[0.0 0.0 0.0 0.0]
               [0.5 0.0 0.0 0.0]
               [0.0 0.5 0.0 0.0]
               [0.0 0.0 1.0 0.0]]
    b_drift =  [1. / 6., 1. / 3., 1. / 3., 1. / 6.]
    c_drift =  [0.0, 0.5, 0.5, 1.0]

    a_diff =  [[ 0.0           0.0         0.0         0.0]
               [-0.7242916356  0.0         0.0         0.0]
               [ 0.4237353406 -0.199443705 0.0         0.0]
               [-1.578475506   0.840100343 1.738375163 0.0]]
    b_diff =   [-0.7800788474, 0.0736376824, 1.486520013, 0.2199211524]
    c_diff =   [0.0, -0.7242916356, 0.2242916356, 1.0]

    a_diff2=  [[ 0.0           0.0         0.0         0.0]
               [ 2.70200041    0.0         0.0         0.0]
               [ 1.757261649   0.0         0.0         0.0]
               [-2.918524118   0.0         0.0         0.0]]
    b_diff2=   [ 1.693950844, 1.636107882, -3.024009558, -0.3060491602]
    c_diff2=   [ 0.0, 2.70200041, 1.757261649, -2.918524118]

    TableauSERK(:BurrageCL_explicit_method, 4, a_drift, b_drift, c_drift, 4, a_diff, b_diff, c_diff, 0, a_diff2, b_diff2, c_diff2)
end

"""
Tableau for the explicit 4-stage E1 method due to K. Burrage and P. Burrage
  Method cited in Eq. (4.2)-(4.3) in
   K. Burrage, P. Burrage (2000) "Order conditions for stochastic Runge-Kutta methods by B-series".
  According to the paper, the method has strong order 1.0 for one-dimensional Brownian motion.
"""
function TableauBurrageE1()

    a_drift = [[0.0      0.0     0.0 0.0]
               [2. / 3.  0.0     0.0 0.0]
               [3. / 2. -1. / 3. 0.0 0.0]
               [7. / 6.  0.0     0.0 0.0]]
    b_drift =  [0.25, 0.75, -0.75, 0.75]
    c_drift =  [0.0, 2. / 3., 7. / 6., 7. / 6.]

    a_diff =  [[ 0.0      0.0      0.0    0.0]
               [ 2. / 3.  0.0      0.0    0.0]
               [ 0.5      1. / 6.  0.0    0.0]
               [-0.5      0.0      0.5    0.0]]
    b_diff =   [-0.5, 1.5, -0.75, 0.75]
    c_diff =   [0.0, 2. / 3., 2. / 3., 0.0]

    a_diff2=  [[ 0.0    0.0    0.0    0.0]
               [ 0.0    0.0    0.0    0.0]
               [ -2. / 3. 0.0    0.0    0.0]
               [ 1. / 6.  0.5    0.0    0.0]]
    b_diff2=   [ 1.5, -1.5, 0.0, 0.0]
    c_diff2=   [ 0.0, 0.0, -2. / 3., 2. / 3.]

    TableauSERK(:BurrageE1_explicit_method, 4, a_drift, b_drift, c_drift, 4, a_diff, b_diff, c_diff, 0, a_diff2, b_diff2, c_diff2)
end

"""
Tableau for the explicit 5-stage G5 method due to K. Burrage and P. Burrage
  Method cited in Section 4 of
   K. Burrage, P. Burrage (2000) "Order conditions for stochastic Runge-Kutta methods by B-series".
  According to the paper, the method has strong order 1.5 for one-dimensional Brownian motion.
"""
function TableauBurrageG5()

    a_drift = [[ 0.0               0.0               0.0               0.0              0.0]
               [ 0.52494822322232  0.0               0.0               0.0              0.0]
               [ 0.07167584568902  0.27192330512685  0.0               0.0              0.0]
               [ 0.13408162649312  0.24489042208103 -0.02150276857782  0.0              0.0]
               [-0.07483338680171 -0.07276896351874  0.55202897082453 -0.50752343840006 0.0]]

    b_drift =  [-5.60958180689351, -0.67641638321828, -5.44025143434789, 8.76396506407891, 3.96228456038077]
    c_drift =  [ 0.0, 0.52494822322232, 0.34359915081587, 0.35746927999633, -0.10309681789598002]

    a_diff =  [[ 0.0               0.0               0.0              0.0              0.0]
               [ 0.52494822322232  0.0               0.0              0.0              0.0]
               [ 0.49977623528582 -0.14576793502675  0.0              0.0              0.0]
               [ 0.60871134749146  0.58291821365556 -0.94596532788804 0.0              0.0]
               [-0.04005606091567 -0.22719654397712 -0.1292628422212  0.42881625288868 0.0]]

    b_diff =   [ 6.68050246229861, 0.0, 4.28273528343281, -3.25408735237225, -6.7091503933593]
    c_diff =   [ 0.0, 0.52494822322232, 0.35400830025907, 0.24566423325898, 0.03230080577469002]

    a_diff2=  [[ 0.0               0.0              0.0               0.0              0.0]
               [ 0.0               0.0              0.0               0.0              0.0]
               [-0.23101439602069  0.59278042710702 0.0               0.0              0.0]
               [-0.54946055077234  0.86811263829203 0.06772607159055  0.0              0.0]
               [ 0.03847082280344 -0.16953882944054 0.88387761274601 -0.85833118389518 0.0]]

    b_diff2=   [ 1.90494977554482, -1.90494977554482, 0.0, 0.0, 0.0]
    c_diff2=   [ 0.0, 0.0, 0.36176603108633, 0.38637815911024, -0.10552157778627003]

    TableauSERK(:BurrageG5_explicit_method, 4, a_drift, b_drift, c_drift, 4, a_diff, b_diff, c_diff, 0, a_diff2, b_diff2, c_diff2)
end


text:
#=
Setup for estimation of the probability for the Schlögl chemical reaction network
with X(0)= x₀ to reach X(T) ∈ A, the target set.

This reaction network corresponds to:
A + 2S → 3S
3S → A + 2S
B → S
S → B

Parameters were chosen so that the system is bistable, as in:
https://doi.org/10.1063/1.5017955
=#

using LinearAlgebra
using Random
using DiffEqBiological

V = 25
a = 1
b = 2
c1 = 3 * 2
c2 = 0.6 * 6
c3 = 0.25
c4 = 2.95

params = (c1, c2, c3, c4, a, b, V);
Tmax = 1.0;

schlogl = @reaction_network begin
    c1 * a/ V, 2*S --> 3*S
    c2 / V^2, 3*S --> 2*S
    c3 * b * V, ∅ --> S
    c4, S --> ∅
end c1 c2 c3 c4 a b V

# preallocate data structure
dprob = DiscreteProblem(schlogl, [1], (0.0, 1.0), params);

# define rxn integrator
function integrate_schlogl(x, t)
    jprob = JumpProblem(remake(dprob, u0=x, tspan=(0,t)), Direct(),
        schlogl, save_positions = (false,false));
    jsol = solve(jprob, SSAStepper());
    return jsol.u[end];
end

x₀ = [5];       # initial state
A = (80,100);   # target set, (80, 100)
f = x-> Int(A[1]<x[1]<A[2]); # indicator function for A


text:
#some schemes
#Emb 4/3 AK p
embedded_scheme_43 = EmbeddedScheme(
          ( 0.125962888700250514,  0.333588446797901933,
            0.751193431379145450, -0.338296598434303506,
            0.127551831557005609,  0.127551831557005609,
           -0.338296598434303506,  0.751193431379145450,
            0.333588446797901933,  0.125962888700250514 ),

          ( 0.125962888700250514,  0.333588446797901933,
            0.751193431379145450, -0.338296598434303506,
            0.0                 ,  0.261153550449697153,
           -0.242703571757396124,  0.596114052266110425,
            0.365547251678000160,  0.147440548920593995 ),

            4 )
#PP 3/4 A
palindromic_scheme_34 = PalindromicScheme(
          ( 0.268330095781759925,  0.919661523017399857,
           -0.187991618799159782, -0.187991618799159782,
            0.919661523017399857,  0.268330095781759925 ),
            3 )
#PP 5/6 A
palindromic_scheme_56 = PalindromicScheme(
          ( 0.201651044312324230,   0.578800656272664932,
            0.562615975356569200,   0.273128836056524479,
            0.253874038247554845,  -0.102733803148432142,
           -0.835351693190370636,   0.068014946093165092,
            0.068014946093165092,  -0.835351693190370636,
           -0.102733803148432142,   0.253874038247554845,
            0.273128836056524479,   0.562615975356569200,
            0.578800656272664932,   0.201651044312324230 ),
            5 )
#defectbased
embedded_scheme_43_D = DefectBasedScheme(
            embedded_scheme_43.scheme2,
            embedded_scheme_43.order-1 )
palindromic_scheme_34_D = DefectBasedScheme(
            palindromic_scheme_34.scheme,
            palindromic_scheme_34.order )
palindromic_scheme_56_D = DefectBasedScheme(
            palindromic_scheme_56.scheme,
            palindromic_scheme_56.order )
schemes=[embedded_scheme_43  ,palindromic_scheme_34  ,palindromic_scheme_56,
         embedded_scheme_43_D,palindromic_scheme_34_D,palindromic_scheme_56_D];


text:
module GRASS

# import external modules
using CSV
using FITSIO
using Random
using StatsBase
using DataFrames
using Statistics
using SharedArrays
using Interpolations
import Glob.glob
import Dates.DateTime

# abbreviations for commonly used types
const AA = AbstractArray
const AF = AbstractFloat

# configure directories
include("config.jl")

# ancillary functions + constants
include("utils.jl")
include("constants.jl")

# user-defined composite types
include("structures.jl")

# star geometry + thermal/RT physics
include("stargeometry.jl")
include("starphysics.jl")

# data read-in + calculations
include("lineIO.jl")
include("bisectors.jl")

# star simulation
include("trim.jl")
include("synthesize.jl")
include("disk_sim.jl")

# processing spectra
include("velocities.jl")

# simulating observations
include("observing/convolutions.jl")
include("observing/signaltonoise.jl")
include("observing/ObservationPlan.jl")

# export some stuff
export SpecParams, DiskParams, synthesize_spectra, calc_ccf, calc_rvs_from_ccf, calc_rms

end # module

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS