text:
using Plots
using Plots.PlotMeasures
using LambertW
import PlotlyBase

plotly()
Plots.PlotlyBackend()

Plots.scalefontsizes()
Plots.scalefontsizes(1.3)

R = 0:0.01:1;

## No-repair reliability
reliability_single = R;
reliability_dr = @. 1 - (1 - R)^2;
reliability_cold = @.  R * (1 - log(R));
reliability_3r = @. 1 - (1 - R)^3;
reliability_tmr = @. 3 * R^2 - 2 * R^3

plot(
    R,
    [reliability_single, reliability_cold, reliability_dr, reliability_tmr, reliability_3r],
    label = ["Single component" "Dual redundancy (cold)" "Dual redundancy (hot)" "Triple modular redundancy (naive)" "Triple modular redundancy (resilient)"],
    # title = "Resilience to permanent failures",
    legend = :outerright,
    xlabel = "Component reliability",
    ylabel = "System reliability",
    lw = 2,
)
savefig("reliability_norepair_en.pdf")

## No-repair reliability with vulnerable voter
Rv = 0:0.01:1;
reliability_tmrv = reliability_tmr * permutedims(Rv)

heatmap(
    R,
    Rv,
    reliability_tmrv,
    xlabel = "Component reliability R",
    ylabel = "Voter reliability Rv",
    colorbar_title = "System reliability",
    aspect_ratio = 1,
    c = :inferno,
)
savefig("reliability_norepair_voter_en.pdf")

## Ρrepair reliability
λ = 10 .^ (LinRange(-11, -2, 1000))

λ_single = λ;
λ_dr = @. -lambertw(-(2λ * exp(-λ) - λ^2 * exp(-2λ)));
λ_tmr = @. -lambertw(-(3(λ * exp(-λ))^2 - 2(λ * exp(-λ))^3));

plot(
    λ,
    [λ_single, λ_single, λ_dr, λ_tmr],
    label = permutedims(
        [
            "Single component"
            "Dual redundancy, cold"
            "Dual redundancy, hot"
            "Triple modular redundancy"
        ],
    ),
    #title = "Resilience to temporary failures",
    legend = (0.6, 0.35),
    xlabel = "Component failure rate",
    ylabel = "System failure rate",
    yaxis = :log,
    xaxis = :log,
    left_margin = -150px,
    lw = 2,
    linecolor = [palette(:tab10)[1] palette(:tab10)[1] palette(:tab10)[2] palette(
        :tab10,
    )[3]],
)
savefig("reliability_repair_en.pdf")


text:
module FDModBiotCyl

    # the dependency of this module
    using FFTW, #rff,irfft (planewave.jl)
    Dierckx, #Spline1D (planewave.jl)
    SpecialFunctions #besselj,hankelh1 (planewave.jl)

    include("./FD/fdtd_export.jl")
    include("./planewave/planewave_export.jl")

end

#using Revise
#]activate .

ext: jl


text:
using Plots, VortexDistributions

@load "./examples/one_frame.jld2"

heatmap(x,y,abs2.(ψ1'))

psi = Torus(ψ1,x,y)
vortices = findvortices(psi)

# zoom window
x0,y0 = 70,200
Lx,Ly = 100,100

# simple cartesian mask
xm = x[@. abs(x-x0)<Lx/2]
ym = y[@. abs(y-y0)<Ly/2]
psim = @. ψ1[abs(x-x0)<Lx/2,abs(y-y0)<Ly/2]

heatmap(xm,ym,angle.(psim),transpose=true)
psi2 = Torus(psim,xm,ym)
vort = findvortices(psi2)

scatter!([vort[1].xv],[vort[1].yv],alpha=0.4)
scatter!([vort[2].xv],[vort[2].yv],alpha=0.4)


text:
struct Colors
    RED::ST0
    GREEN::ST1
    BLUE::ST2
end

RED = auto()
GREEN = auto()
BLUE = auto()
struct Permissions
    R::ST0
    W::ST1
    X::ST2
end

R = 1
W = 2
X = 16


text:
##############################################################################
##
## read formula
##
##############################################################################

function decompose!(rf::Formula)
	(has_absorb, absorb_formula, absorb_terms) = decompose_absorb!(rf)
	(has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose_iv!(rf)
	return (has_absorb, absorb_formula, absorb_terms, has_iv, iv_formula, iv_terms, endo_formula, endo_terms)
end

# decompose formula into normal + iv vs absorbpart
function decompose_absorb!(rf::Formula)
	has_absorb = false
	absorb_formula = nothing
	absorb_terms = nothing
	if typeof(rf.rhs) == Expr && rf.rhs.args[1] == :(|>)
		has_absorb = true
		absorb_formula = Formula(nothing, rf.rhs.args[3])
		rf.rhs = rf.rhs.args[2]
		absorb_terms = Terms(absorb_formula)
	end
	return(has_absorb, absorb_formula, absorb_terms)
end

# decompose formula into normal vs iv part
function decompose_iv!(rf::Formula)
	has_iv = false
	iv_formula = nothing
	iv_terms = nothing
	endo_formula = nothing
	endo_terms = nothing

	if typeof(rf.rhs) == Expr
		if rf.rhs.head == :(=)
			has_iv = true
			if isa(rf.rhs.args[2], Expr) && rf.rhs.args[2].head == :block
				# happens when several endogeneous variable
				iv_formula = Formula(nothing,  rf.rhs.args[2].args[2])
			else
				iv_formula = Formula(nothing,  rf.rhs.args[2])
			end
			endo_formula = Formula(nothing, rf.rhs.args[1])
			rf.rhs = :1
		elseif rf.rhs.head == :call
			i = 1
			while !has_iv && i <= length(rf.rhs.args)
				if isa(rf.rhs.args[i], Expr) && rf.rhs.args[i].head == :(=)
					has_iv = true
					iv_vars = rf.rhs.args[2]
					if isa(rf.rhs.args[i].args[2], Expr) && rf.rhs.args[i].args[2].head == :block
						# happens when several endogeneous variable
						iv_formula = Formula(nothing,  rf.rhs.args[i].args[2].args[2])
					else
						iv_formula = Formula(nothing,  rf.rhs.args[i].args[2])
					end
					endo_formula = Formula(nothing, rf.rhs.args[i].args[1])
					splice!(rf.rhs.args, i)
				else
					i += 1
				end
			end
		else
			error("formula not correct")
		end
	end

	if has_iv
		iv_terms = Terms(iv_formula)
		iv_terms.intercept = false
		endo_terms = Terms(endo_formula)
		endo_terms.intercept = false
	end
	return(has_iv, iv_formula, iv_terms, endo_formula, endo_terms)
end

function secondstage!(rf::Formula)
	(has_absorb, absorb_formula, absorb_terms, has_iv, iv_formula, iv_terms, endo_formula, endo_terms) = decompose!(rf)
	if has_iv
	    if typeof(rf.rhs) == Symbol
	        rf.rhs = endo_formula.rhs
	    else
	        push!(rf.rhs.args, endo_formula.rhs)
	    end
	end
end

##############################################################################
##
## build model
##
##############################################################################

function ModelFrame2(trms::Terms, d::AbstractDataFrame, esample; contrasts::Dict = Dict())
	mf = ModelFrame(trms, d; contrasts = contrasts)
	mf.msng = esample
	return mf
end

#  remove observations with negative weights
function isnaorneg{T <: Real}(a::Vector{T})
	BitArray(a .> zero(eltype(a)))
end
function isnaorneg{T <: Real}(a::DataVector{T})
	out = !a.na
	@simd for i in 1:length(a)
		if out[i]
			@inbounds out[i] = a[i] > zero(Float64)
		end
	end
	BitArray(out)
end

# Directly from DataFrames.jl

function dropresponse(trms::Terms)
    if trms.response
        ckeep = 2:size(trms.factors, 2)
        rkeep = vec(any(trms.factors[:, ckeep], 2))
        Terms(trms.terms, trms.eterms[rkeep], trms.factors[rkeep, ckeep],
              trms.is_non_redundant[rkeep, ckeep], trms.order[ckeep], false, trms.intercept)
    else
        trms
    end
end

function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    vcat([allvars(a) for a in ex.args[2:end]]...)
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(::Any) = Array(Symbol, 0)

# used when removing certain rows in a dataset
# NA always removed
function dropUnusedLevels!(f::PooledDataVector)
	uu = unique(f.refs)
	length(uu) == length(f.pool) && return f
	sort!(uu)
	T = reftype(length(uu))
	dict = Dict{eltype(uu), T}(zip(uu, 1:convert(T, length(uu))))
	@inbounds @simd  for i in 1:length(f.refs)
		 f.refs[i] = dict[f.refs[i]]
	end
	f.pool = f.pool[uu]
	return f
end

dropUnusedLevels!(f::DataVector) = f


text:
IndexedTables.table(s::SelectedData) = isempty(orderby(s)) ? s.table : reindex(s.table, orderby(s))
IndexedTables.table(a::Analysis) = table(a.data)

process(::Type{<:IndexedTable}, a::Analysis) = table(a)


text:
# TODO make this stuff work again

function graphviz_identifier(x)
    s = string(x)
    result = replace(s, '"' => "\\\"")
    if endswith(result, '\\')
        result = result * " "
    end
    return result
end

function show(io::IO, m::MIME"text/vnd.graphviz", g::AbstractValueGraph)
    graph_symbol = is_directed(g) ? "digraph" : "graph"
    edge_symbol = is_directed(g) ? "->" : "--"
    println(io, "$graph_symbol {")
    for e in edges(g)
        s, d, v = src(e), dst(e), val(e)
        println(io, "$s $edge_symbol $d [value = \"$(graphviz_identifier(v))\"];")
    end
    println(io, "}")
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
include("compile.jl")

# env var to force compilation from source, for testing purposes
const forcecompile = get(ENV, "FORCE_COMPILE_BLOSC", "no") == "yes"

# Parse some basic command-line arguments
const verbose = ("--verbose" in ARGS) || forcecompile
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libblosc"], :libblosc),
]
verbose && forcecompile && Compat.@info("Forcing compilation from source.")

# Download binaries from hosted location
bin_prefix = "https://github.com/stevengj/BloscBuilder/releases/download/v1.14.3+5"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/Blosc.v1.14.3.aarch64-linux-gnu.tar.gz", "2c2af6cc54f9101420d2652a12d5da64dbd7526af8be7d748d4f1fd7686d5ecc"),
    Linux(:aarch64, :musl) => ("$bin_prefix/Blosc.v1.14.3.aarch64-linux-musl.tar.gz", "44da7bf787ae3992f619781e82d0546b156c1f4926c387201779184b710933f2"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/Blosc.v1.14.3.arm-linux-gnueabihf.tar.gz", "d151548cd0d4f2c74f22bb9810baa6957bc91f7770c2952a3f98b9f61761d831"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/Blosc.v1.14.3.arm-linux-musleabihf.tar.gz", "01216681178d3ae7d01dba6670774e321364affb7c3077b8e31fcf367693b924"),
    Linux(:i686, :glibc) => ("$bin_prefix/Blosc.v1.14.3.i686-linux-gnu.tar.gz", "8b833b6544485c906729081b2f862798feb860b4744b467e8ef520a940da0aa6"),
    Linux(:i686, :musl) => ("$bin_prefix/Blosc.v1.14.3.i686-linux-musl.tar.gz", "9830099baae78d859c9e88d8ab2a93af5a953f56d46a122504eb9a27f909b8a6"),
    Windows(:i686) => ("$bin_prefix/Blosc.v1.14.3.i686-w64-mingw32.tar.gz", "97fd5573b324bdfc5a1a3f6b83095be8e2660e5e5ac0e33d865a9c7804f2bcd0"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/Blosc.v1.14.3.powerpc64le-linux-gnu.tar.gz", "b25b60e907cba82d5f35387d6d2eb438d9898a4f6789bd702b114f3c0667b330"),
    MacOS(:x86_64) => ("$bin_prefix/Blosc.v1.14.3.x86_64-apple-darwin14.tar.gz", "d038eec02df45e3d83450b87182667d619c64b06a189d0d9be3be29877356c41"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/Blosc.v1.14.3.x86_64-linux-gnu.tar.gz", "eb98760d5c54592b82f71aa11f28a1d6a1d588eb48c9d23977a74a513544cda0"),
    Linux(:x86_64, :musl) => ("$bin_prefix/Blosc.v1.14.3.x86_64-linux-musl.tar.gz", "ec25976201fba46cd8ce9674c9b25f87ad5d8e1399cf002cab94fa242a17aa3d"),
    FreeBSD(:x86_64) => ("$bin_prefix/Blosc.v1.14.3.x86_64-unknown-freebsd11.1.tar.gz", "d58104e836b8b39e665b957b27ccba41829cbbd8c803a49af9faffe797ae5f49"),
    Windows(:x86_64) => ("$bin_prefix/Blosc.v1.14.3.x86_64-w64-mingw32.tar.gz", "0eb9f0efd47bdaba7a506770c876575e2c49e9d1648aa8a41d083b7c996fe2b7"),
)

# source code tarball and hash for fallback compilation
source_url = "https://github.com/Blosc/c-blosc/archive/v1.14.3.tar.gz"
source_hash = "7217659d8ef383999d90207a98c9a2555f7b46e10fa7d21ab5a1f92c861d18f7"

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key()) && !forcecompile
    url, tarball_hash = download_info[platform_key()]
    if !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)

        # check again whether the dependency is satisfied, which
        # may not be true if dlopen fails due to a libc++ incompatibility (#50)
        unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
    end
end

if unsatisfied || forcecompile
    # Fall back to building from source, giving the library a different name
    # so that it is not overwritten by BinaryBuilder downloads or vice-versa.
    libname = "libblosc_from_source"
    products = [ LibraryProduct(prefix, [libname], :libblosc) ]
    source_path = joinpath(prefix, "downloads", "src.tar.gz")
    if !isfile(source_path) || !verify(source_path, source_hash; verbose=verbose) || !satisfied(products[1]; verbose=verbose)
        compile(libname, source_url, source_hash, prefix=prefix, verbose=verbose)
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
using Random
include("bike.jl")

function initialize_random_bike() :: Bike
    make_num = rand((1,100))
    model_num = rand((1,10))
    weight = rand(Float64) * 5 + 15
    drag_coefficient = 1 - rand(Float64)
    bike_attributes = BikeAttributes("Make$(make_num)", "Model$(model_num)", weight, drag_coefficient)
    x = rand((1,100))
    y = rand((1,100))
    bike = Bike(bike_attributes, x, y)
    return bike

end


text:
struct MountainCarOnlineModelSearchAgent
    mountaincar::MountainCar
    model::MountainCar
    num_eval_samples::Int64
end

function run(
    agent::MountainCarOnlineModelSearchAgent,
    rng::MersenneTwister;
    max_steps = 3e3,
    debug = false,
    max_likelihood = false,
    optimistic = false,
    epsilon = 0.0,
    model_evaluation = false,
)
    # Start with initial policy
    start_params = MountainCarParameters(-0.0025, 3)
    params = vec(start_params)
    policy, _, _ = value_iteration(agent.model, params)
    num_steps = 0.0
    state = init_random(agent.mountaincar, rng)
    actions = getActions(agent.mountaincar)
    data::Array{MountainCarContTransition} = []
    while !checkGoal(agent.mountaincar, state) && num_steps < max_steps
        num_steps += 1
        toss = rand(rng, Uniform(0, 1))
        action = actions[policy[cont_state_to_idx(agent.mountaincar, state)]]
        if toss < epsilon
            # random action
            action = actions[rand(rng, [1, 2])]
        end
        next_state, cost = step(agent.mountaincar, state, action, true_params)
        push!(data, MountainCarContTransition(state, action, cost, next_state))
        if num_steps % agent.num_eval_samples == 0
            # Update policy
            if max_likelihood
                params = get_least_squares_fit(agent.model, params, data)
            else
                optimization_params = MountainCarOptimizationParameters(
                    [0.0012, 0.5], vec(start_params), 75, false,
                )
                horizon = 500
                params = return_based_model_search(
                    agent.model,
                    data,
                    optimization_params,
                    horizon;
                    debug = debug,
                    optimistic = optimistic,
                    model_evaluation = model_evaluation,
                )
            end
            if debug
                println("Moving to params ", params)
            end
            policy, _, _ = value_iteration(agent.model, params)
        end
        state = next_state
    end
    num_steps
end

function run_cmax(
    agent::MountainCarOnlineModelSearchAgent,
    rng::MersenneTwister;
    max_steps = 3e3,
    debug = false,
    cmax_threshold = 0.01,
)
    start_params = MountainCarParameters(-0.0025, 3)
    cost_matrix = generate_cost_matrix(agent.model)
    transition_matrix = generate_transition_matrix(agent.model, vec(start_params))
    policy, _, _ = value_iteration(agent.model, transition_matrix, cost_matrix)
    num_steps = 0.0
    state = init_random(agent.mountaincar, rng)
    actions = getActions(agent.mountaincar)
    while !checkGoal(agent.mountaincar, state) && num_steps < max_steps
        num_steps += 1
        action = actions[policy[cont_state_to_idx(agent.mountaincar, state)]]
        next_state, cost = step(agent.mountaincar, state, action, true_params)
        next_predicted_state, _ = step(agent.model, state, action, true_params)
        if (abs(next_state.position - next_predicted_state.position) +
            abs(next_state.speed - next_predicted_state.speed)) > cmax_threshold
            # Discrepancy; inflate cost
            cost_matrix[cont_state_to_idx(agent.mountaincar, state),
                        policy[cont_state_to_idx(agent.mountaincar, state)]] =
                            1e5
            # Update policy
            policy, _, _ = value_iteration(agent.model, transition_matrix, cost_matrix)
        end
        state = next_state
    end
    num_steps
end

function run_true(
    agent::MountainCarOnlineModelSearchAgent,
    rng::MersenneTwister;
    max_steps = 3e3,
    debug = false,
)
    start_params = true_params
    policy, _, converged = value_iteration(agent.mountaincar, vec(start_params))
    if !converged && debug
        println("Not converged with params ", vec(start_params))
    end
    num_steps = 0.0
    state = init_random(agent.mountaincar, rng)
    actions = getActions(agent.mountaincar)
    while !checkGoal(agent.mountaincar, state) && num_steps < max_steps
        num_steps += 1
        action = actions[policy[cont_state_to_idx(agent.mountaincar, state)]]
        state, cost = step(agent.mountaincar, state, action, true_params)
    end
    num_steps
end


text:
using Test
using NQCDynamics
using FiniteDiff
using LinearAlgebra: norm

function test_motion!(sim::Simulation{<:eCMM}, u)
    f(x) = DynamicsUtils.classical_hamiltonian(sim, x)

    grad = FiniteDiff.finite_difference_gradient(f, u)

    du = zero(u)
    DynamicsMethods.motion!(du, u, sim, 0.0)

    @test DynamicsUtils.get_positions(du) ≈ DynamicsUtils.get_velocities(grad) ./ sim.atoms.masses' rtol=1e-3
    @test DynamicsUtils.get_velocities(du) ≈ -DynamicsUtils.get_positions(grad) ./ sim.atoms.masses' rtol=1e-3
    @test DynamicsMethods.MappingVariableMethods.get_mapping_positions(du) ≈ DynamicsMethods.MappingVariableMethods.get_mapping_momenta(grad) rtol=1e-3
    @test DynamicsMethods.MappingVariableMethods.get_mapping_momenta(du) ≈ -DynamicsMethods.MappingVariableMethods.get_mapping_positions(grad) rtol=1e-3
end

sim = Simulation{eCMM}(Atoms(1), DoubleWell(); γ=0.0)
sim1 = Simulation{eCMM}(Atoms(1), DoubleWell(); γ=0.5)

v = randn(1,1)
r = randn(1,1)
u = DynamicsVariables(sim, v, r, PureState(1))
u1 = DynamicsVariables(sim1, v, r, PureState(1))

test_motion!(sim, u)
test_motion!(sim1, u1)

sol = run_trajectory(u, (0, 100.0), sim; output=(:hamiltonian, :position, :u), reltol=1e-10, abstol=1e-10)
@test sol.hamiltonian[1] ≈ sol.hamiltonian[end] rtol=1e-3
qmap = [u.qmap for u in sol.u]
pmap = [u.pmap for u in sol.u]
total_population = sum.(DynamicsMethods.MappingVariableMethods.mapping_kernel.(qmap, pmap, sim.method.γ))
@test all(isapprox.(total_population, 1, rtol=1e-3))

sol = run_trajectory(u1, (0, 100.0), sim1; output=(:hamiltonian, :position, :u), reltol=1e-10, abstol=1e-10)
@test sol.hamiltonian[1] ≈ sol.hamiltonian[end] rtol=1e-3
total_population = sum.(DynamicsMethods.MappingVariableMethods.mapping_kernel.(qmap, pmap, sim.method.γ))
@test all(isapprox.(total_population, 1, rtol=1e-3))

@testset "generate_random_points_on_nsphere" begin
    points = DynamicsMethods.MappingVariableMethods.generate_random_points_on_nsphere(10, 1)
    @test norm(points) ≈ 1
    points = DynamicsMethods.MappingVariableMethods.generate_random_points_on_nsphere(10, 10)
    @test norm(points) ≈ 10
end

@testset "Population correlation" begin
    # Tests that the initial state correlated with itself is 1 and correlated with other state is 0.
    gams = -0.5:0.5:1.5
    for γ in gams[2:end]
        @testset "Gamma = $γ" begin
            sim = Simulation{eCMM}(Atoms(1), DoubleWell(); γ=γ)
            out = zeros(2, 2)
            n = 1e4
            correlation = TimeCorrelationFunctions.PopulationCorrelationFunction(sim, Diabatic())
            normalisation = TimeCorrelationFunctions.evaluate_normalisation(sim, correlation)

            for i=1:n
                u = DynamicsVariables(sim, 0, 0, PureState(1))
                K = Estimators.initial_diabatic_population(sim, u)
                Kinv = Estimators.diabatic_population(sim, u)
                out .+= normalisation * K * Kinv'
            end
            @test out ./ n ≈ [1 0; 0 1] atol=0.1
        end
    end
end


text:
include("test_Vector2D.jl")
include("test_BilliardModels.jl")

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


text:
type GoldbergFitnessScheme{T} <: ParetoFitnessScheme{T}
  maximise::Vector{Bool}
end

"""
Implements a pareto-based multiple-objective fitness scheme through the use
of Goldberg fitness ranking, wherein individuals are assigned a rank based
upon the number of individuals by whom they are dominated, as given by
the formula:

  Rank(i) = 1 + DominatedBy(i)

**Parameters:**

* `of:Type`, the base fitness type (default is Float).
* `maximise:Bool`, a flag, indicating whether fitness values are to be maximised
or minimised.
"""
function goldberg(s::Dict{Any, Any})
   s["of"] = eval(Base.parse(Base.get(s, "of", "Float")))
   GoldbergFitnessScheme{s["of"]}(s["maximise"])
end

function scale!{F}(s::GoldbergFitnessScheme, inds::Vector{Tuple{Int, F}})
  n = length(inds)
  j = k = rank = 1

  # Keep calculating each of the pareto fronts until all individuals have
  # been handled.
  while j <= n

    # Check each remaining member for inclusion in the current pareto front.
    # If the individual belongs to the pareto front then swap it with the
    # first unsorted individual.
    for i in j:n
      (p1_id, p1_f) = inds[i]
      if all((p2_id, p2_f) -> p1_id == p2_id || !dominates(p2_f.scores, p1_f.scores, s.maximise), inds[j:end])
        p1_f.rank = rank
        inds[k], inds[i] = inds[i], inds[k]
        k += 1
      end
    end

    j = k
    rank += 1
  end
end


text:
"""
    AbstractRegion

`AbstractRegion` is an abstract type representing region of space which will be
used to perform a spatial query.

An implementation `Region <: AbstractRegion` must provide implementations for;
- `in(::SVector{3, T}, ::Region{T})`
- `boundingbox{T}(::Region{T})`

An implementation `Region <: AbstractRegion` may provided optimized
implementations for;
- `might_intersect(::BoundingBox{T}, ::Region{T})`

See also `BoundingBox`, `Sphere`, `Cylinder`, `Circle`, `TriangularPrism`, and `Triangle`.
"""
abstract type AbstractRegion{T <: Real} end

"""
    in(p, region)

Determine whether a point `p` is within the given `region`.
"""
@inline function in(p::AbstractVector{T}, region::AbstractRegion) where T
    return in(convert(SVector{3,T}, p), region)
end
function in(p::StaticVector{3, T}, region::Region) where {T <: Real, Region <: AbstractRegion}
    error("Base.in not implemented for region type $Region")
end

"""
    eltype(region)

Determine the type of Real used to describe the `region` geometry
"""
eltype(::AbstractRegion{T}) where {T} = SVector{3, T}
eltype(::Type{<:AbstractRegion{T}}) where {T} = SVector{3, T}

"""
    boundingbox(region)

Construct a `BoundingBox` which encapsulates all the `region` which is of a
sub-type of AbstractRegion.
"""
function boundingbox(::Region) where Region <: AbstractRegion
    error("bounding box not implemented for region type $(Region)")
end


text:
using RvSpectMLBase
using Test

@testset "ChunkList" begin
    num_pixels = 100
    num_orders = 10
    resolution = 1e5
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    inst = TheoreticalInstrument.TheoreticalInstrument2D(lambda2d)
    spectrum = TheoreticalInstrument.demo_generate_spectrum_line(inst)
    @test_nowarn make_orders_into_chunks(spectrum, inst)
    @test_nowarn make_orders_into_chunks(spectrum, orders_to_use=1:5, pixels_to_use=fill(20:80,5) )
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
## Note: this file is supposed to be run from the TensorBoardLogger/gen folder as
## $ cd .../TensorBoardLogger/gen
## $ julia --project=. compile_proto.jl

# The proto folder is generated by the following script
# and then manually cleaning the output from useless stuff
#
# mkdir -p tmp/proto && cd tmp
# git clone --recurse-submodules git@github.com:tensorflow/tensorboard.git
# git clone --recurse-submodules git@github.com:tensorflow/tensorflow.git
# rsync -zarvm --include="*/" --include="*.proto" --exclude="*" tensorboard proto
# rsync -zarvm --include="*/" --include="*.proto" --exclude="*" tensorflow proto

using ProtoBuf
using Glob
using FilePaths; using FilePathsBase: /

pbpath =dirname(dirname(PosixPath(pathof(ProtoBuf))))/p"gen"

## Setup input and output paths
cur_path = cwd()
TBL_root = dirname(cur_path)

src_dir = cur_path/"proto"
out_dir   = cur_path/"protojl"

## Clean the output directory
rm(out_dir, force=true, recursive=true)

## First module
function process_module(cur_module::AbstractString; base_module::AbstractString=cur_module, input_path=cur_module)
    # Include search paths
    includes   = [src_dir, src_dir/base_module]

    # Output folder
    module_out_dir    = out_dir/cur_module

    # Input files
    infiles = glob("*.proto", src_dir/input_path)

    mkpath(module_out_dir)
    includes_str=["--proto_path=$path" for path=includes]
    run(ProtoBuf.protoc(`$includes_str --julia_out=$module_out_dir $infiles`))

    nothing
end

#process_module("tensorflow", input_path="tensorflow/core/protobuf")

process_module("tensorboard", input_path="tensorboard/compat/proto")

#plugins = ["audio", "custom_scalar", "hparams", "histogram", "image", "scalar", "text"]
plugins = ["custom_scalar", "hparams", "text"]
for plugin in plugins
    process_module("tensorboard/plugins/$plugin", base_module="tensorboard")
end

## this fails but would be better
#cur_module  = "tensorboard"
#base_module = cur_module
#
## Include search paths
#includes   = [src_dir, src_dir/base_module]
#
## Output folder
#module_out_dir    = out_dir/("$cur_module"*"2")
#
## Input files
#infiles = glob("*.proto", src_dir/cur_module/"compat/proto")
#
#for plugin in plugins
#    plugin_proto_files = glob("*.proto", src_dir/cur_module/"plugins/$plugin")
#    append!(infiles, plugin_proto_files)
#end
#
#mkpath(module_out_dir)
#includes_str=["--proto_path=$path" for path=includes]
#run(ProtoBuf.protoc(`$includes_str --julia_out=$module_out_dir $infiles`))

# Finally move the output directory to the src folder
mv(out_dir, TBL_root/"src"/"protojl")


text:
using Base.Test, QuantumInfo

import Base.isapprox

rtoldefault = Base.rtoldefault

#function isapprox(m1::Matrix,m2::Matrix; rtol::Real=rtoldefault(abs(m1[1,1]),abs(m2[1,1])), atol::Real=0)
#  all(x->isapprox(x,0.0,rtol=rtol,atol=atol),abs(m1-m2))
#end

function isapprox{T1,T2}(v1::Vector{Matrix{T1}},v2::Vector{Matrix{T2}}; rtol::Real=rtoldefault(abs(v1[1][1,1]),abs(v2[1][1,1])), atol::Real=0)
  local v
  if length(v1) != length(v2)
    false
  else
    v = Matrix{T1}[ abs(v1[i]-v2[i]) for i in 1:length(v1) ]
    all(x->isapprox(x,zeros(size(x,1),size(x,2)),rtol=rtol,atol=atol),v)
  end
end

#function ishermitian(m::Matrix; rtol::Real=rtoldefault(abs((m[1,1]+conj(m[1,1]))/2),abs(m[1,1])), atol::Real=0)
#  isapprox((m+m')/2,m,rtol=rtol,atol=atol)
#end

#function ispositive(m::Matrix;rtol::Real=-1,atol::Real=0)
#  if rtol<0
#      ishermitian(m) && all(x->isapprox(x,0,atol=atol),filter(x->x<0,eigvals(Hermitian(m))))
#  else
#      ishermitian(m) && all(x->isapprox(x,0,atol=atol,rtol=rtol),filter(x->x<0,eigvals(Hermitian(m))))
#  end
#end

let

  for i=1:100
    da = round(Int,round(rand()*10))+1
    db = round(Int,round(rand()*10))+1
    ra = partialtrace(projector(randn(da*3)+1im*randn(da*3)),[da,3],2)
    rb = partialtrace(projector(randn(db*3)+1im*randn(db*3)),[db,3],2)
    @test isapprox(partialtrace(kron(ra,rb),[da,db],1),rb,atol=1e-15)
    @test isapprox(partialtrace(kron(ra,rb),[da,db],2),ra,atol=1e-15)
  end

  id_choi = 1/2.*[1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1];

  @test liou(eye(2),eye(2)) == eye(4)
  @test liou(eye(2)) == liou(eye(2),eye(2))
  @test mat(vec(eye(2))) == eye(2)
  @test mat([1,2,3,4]) == [1 3;2 4]
  @test vec([1 2; 3 4]) == [1,3,2,4]

  @test choi_liou_involution(choi_liou_involution(eye(4)))==eye(4)

  @test isapprox(choi2liou(id_choi),eye(4))
  @test isapprox(liou2choi(eye(4)),id_choi)
  @test isapprox(choi2kraus(id_choi),Matrix{Complex128}[zeros(2,2), zeros(2,2), zeros(2,2), eye(2)],atol=1e-7)
  @test isapprox(kraus2liou(Matrix{Complex128}[eye(2)]),eye(4))
  @test isapprox(liou2kraus(eye(4)),Matrix{Complex128}[zeros(2,2), zeros(2,2), zeros(2,2), eye(2)],atol=1e-7)
  @test isapprox(kraus2choi(Matrix{Complex128}[eye(2)]),id_choi)

  @test isapprox(depol(2),projector(vec(eye(2))),atol=1e-15)
  @test isapprox(depol(2),depol(2,1.))
  @test isapprox(depol(2,0.),eye(4),atol=1e-15)

  for i=1:100
    rrho = partialtrace(projector(randn(3^3)+1im*randn(3^3)),[3,9],2)
    @test isapprox(trace(rrho),1.)
    @test ishermitian(rrho)
    @test ispossemidef(rrho)

    ru = svd(randn(3,3)+1im*randn(3,3))[1]
    @test isapprox(ru*ru',eye(3),atol=1e-13)
    @test isapprox(ru'*ru,eye(3),atol=1e-13)

    rv = svd(randn(3,3)+1im*randn(3,3))[1]
    re = (liou(ru)+liou(rv))/2
    @test ispossemidef(liou2choi(re),tol=1e-13)
    @test isapprox(partialtrace(liou2choi(re),[3,3],2),eye(3)/3,atol=1e-15)
  end

  @test isapprox(liou2pauliliou(eye(4)), eye(4))
  @test isapprox(pauliliou2liou(eye(4)), eye(4))

end


text:
using CLEARSWI
using Documenter

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

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

deploydocs(;
    repo="github.com/korbinian90/CLEARSWI.jl",
    devbranch="master",
)


text:
"""
    AbstractSolver

Abstract supertype for all direct integration methods.
"""
abstract type AbstractSolver end

"""
    step_size(alg::AbstractSolver)

Return the step size of the given algorithm.

### Input

- `alg` -- algorithm

### Output

The step size of the algorithm, or `nothing` if the step-size is not fixed.
"""
step_size(alg::AbstractSolver) = alg.Δt

"""
    AbstractSolution

Abstract supertype that holds the solution of a numerical integration.
"""
abstract type AbstractSolution end

"""
    Solution{T<:AbstractSolver, UT, VT, AT} <: AbstractSolution

### Fields

- `alg` -- algorithm used in the integration
- `U`   -- displacements
- `U′`  -- velocities
- `U′′` -- accelerations
- `t`   -- vector of time values
"""
struct Solution{T<:AbstractSolver, UT, VT, AT, ST} <: AbstractSolution
    alg::T
    U::UT
    U′::VT
    U′′::AT
    t::ST
end

# constructor with missing fields
Solution(alg, U, t) = Solution(alg, U, nothing, nothing, t)

"""
    dim(sol::Solution)

Return the ambient dimension of the state space of the solution.
"""
dim(sol::Solution) = length(first(sol.U))

"""
    displacements(sol::Solution)

Return the vector of displacements of the given solution.
"""
displacements(sol::Solution) = sol.U

"""
    displacements(sol::Solution, i::Int)

Return the vector of displacements of the given solution along coordinate `i`.
"""
function displacements(sol::Solution, i::Int)
    1 ≤ i ≤ dim(sol) || throw(ArgumentError("expected the coordinate to be between 1 and $(dim(sol)), got $i"))
    U = displacements(sol)
    return [u[i] for u in U]
end

"""
    velocities(sol::Solution)

Return the vector of velocities of the given solution.
"""
velocities(sol::Solution) = sol.U′

"""
    velocities(sol::Solution, i::Int)

Return the vector of velocities of the given solution along coordinate `i`.
"""
function velocities(sol::Solution, i::Int)
    1 ≤ i ≤ dim(sol) || throw(ArgumentError("expected the coordinate to be between 1 and $(dim(sol)), got $i"))
    U′ = velocities(sol)
    return [u′[i] for u′ in U′]
end

"""
    accelerations(sol::Solution)

Return the vector of accelerations of the given solution.
"""
accelerations(sol::Solution) = sol.U′′

"""
    accelerations(sol::Solution, i::Int)

Return the vector of accelerations of the given solution along coordinate `i`.
"""
function accelerations(sol::Solution, i::Int)
    1 ≤ i ≤ dim(sol) || throw(ArgumentError("expected the coordinate to be between 1 and $(dim(sol)), got $i"))
    U′′ = accelerations(sol)
    return [u′′[i] for u′′ in U′′]
end

"""
    times(sol::Solution)

Return the vector of times of the given solution.
"""
times(sol::Solution) = sol.t

struct SolutionExtrema{T<:AbstractSolver, VT, ST, DT, RT} <: AbstractSolution
    alg::T
    U_min::VT
    U_max::VT
    U′_min::VT
    U′_max::VT
    idx::Int
    hist::DT
    times::RT
    t::ST
end

# ===============
# Problem types
# ===============

abstract type AbstractProblem end

struct StructuralDynamicsProblem{T, PT} <: AbstractProblem
    alg::T
    ivp::PT
    NSTEPS::Int
end

"""
    solve(ivp::InitialValueProblem, alg::AbstractSolver, args..; kwargs...)

Solve an initial-value problem.

### Input

- `ivp` -- initial-value problem
- `alg` -- algorithm

### Output

A solution structure (`Solution`) that holds the result and the algorithm used
to obtain it.
"""
function solve(ivp::IVP{<:AbstractContinuousSystem}, alg::AbstractSolver, args...; kwargs...)
    sdprob = init(ivp, alg, args...; kwargs...)
    return _solve(sdprob.alg, sdprob.ivp, sdprob.NSTEPS; kwargs...)
end

const SOACS = SecondOrderConstrainedLinearControlContinuousSystem
const SOCLCCS  = SecondOrderConstrainedLinearControlContinuousSystem

function init(ivp::InitialValueProblem{ST, XT},
              alg::AbstractSolver;
              kwargs...) where {N, VT,
                             ST, # FIXME restrict to SOACS and SOCLCCS
                             XT<:Tuple{VT, VT}}

    if haskey(kwargs, :NSTEPS)
        NSTEPS = kwargs[:NSTEPS]
    elseif haskey(kwargs, :T) || haskey(kwargs, :finalTime)
        Δt = step_size(alg)
        if haskey(kwargs, :T)
            NSTEPS = ceil(Int, kwargs[:T] / Δt)
        else haskey(kwargs, :finalTime)
            NSTEPS = ceil(Int, kwargs[:finalTime] / Δt)
        end
    elseif haskey(kwargs, :tspan)
        Δt = step_size(alg)
        tsp = kwargs[:tspan]
        @assert iszero(tsp[1]) "expected that the initial time is zero, got $(tsp[1])"
        T = tsp[2]
        NSTEPS = ceil(Int, T / Δt)
    else
        throw(ArgumentError("please define `NSTEPS`, `T`, `finalTime` or `tspan`"))
    end
    return StructuralDynamicsProblem(alg, ivp, NSTEPS)
end

# lazily extend the vector to the required number of steps
function _init_input(R::AbstractVector{N}, IMAX) where {N<:Number}
    return Fill(R, IMAX)
end

# no-op, checking that the number of forcing terms has the correct length
function _init_input(R::AbstractVector{VT}, IMAX) where {N, VT<:AbstractVector{N}}
    @assert length(R) == IMAX "expected the forcing term to be an array of length $IMAX, got $(length(R))"
    return R
end

# TODO dispatch on B (eg. IdentityMultiple)
function _init_input(R::AbstractVector{VT}, B::AbstractMatrix, IMAX) where {N, VT<:AbstractVector{N}}
    @assert length(R) == IMAX "expected the forcing term to be an array of length $IMAX, got $(length(R))"
    return [B*Ri for Ri in R]
end

# unwrap a second order system into its each component
function _unwrap(sys::SecondOrderAffineContinuousSystem, IMAX)
    M = mass_matrix(sys)
    C = viscosity_matrix(sys)
    K = stiffness_matrix(sys)
    R = affine_term(sys)
    R = _init_input(R, IMAX)
    return M, C, K, R
end

function _unwrap(sys::SecondOrderConstrainedLinearControlContinuousSystem, IMAX)
    M = mass_matrix(sys)
    C = viscosity_matrix(sys)
    K = stiffness_matrix(sys)
    R = inputset(sys)
    B = input_matrix(sys)
    R = _init_input(R, B, IMAX)
    return M, C, K, R
end

# ================
# Plot recipes
# ================

function _check_vars(vars)
    if vars == nothing
        throw(ArgumentError("default ploting variables not implemented yet; you need " *
              "to pass the `vars=(...)` option, e.g. `vars=(0, 1)` to plot variable with " *
              "index 1 vs. time, or `vars=(1, 2)` to plot variable with index 2 vs. variable with index 1`"))
    end
    D = length(vars)
    @assert (D == 1) || (D == 2) "can only plot in one or two dimensions, " *
                                 "but received $D variable indices where `vars = ` $vars"
end

# plot displacements of the solution for the given vars tuple, eg. vars=(0, 1) for x1(t) vs t
@recipe function plot_solution(sol::Solution; vars=nothing, func=displacements)

   seriestype -->  :path # :scatter
   markershape --> :circle

   _check_vars(vars)

   if vars[1] == 0 && vars[2] != 0
       x = times(sol)
       y = func(sol, vars[2])
       x, y
    else
       x = func(sol, vars[1])
       y = func(sol, vars[2])
    end
    return x, y
end


text:
# invert all the ITMIX 2

include("itmix-setup.jl") # run this first to make sure all precompilation this through
                          # before the parallel run starts

if nprocs()<2
    if Sys.CPU_CORES>=17
        addprocs(17) # /2 to get to physical cores
    else
        addprocs(Sys.CPU_CORES÷2 + 1) # /2 to get to physical cores
    end
end

@everywhere begin
    const debugging = true
    include("itmix-setup.jl")
end
# whether doing a testing, semi-testing or production run
runtyp = ["test", "testmid", "prodlow", "prod"][1]
runtyp == :test && println("\n\nTEST RUN !!!!!!!!\n\n")

fit_target = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.length][2]  # making this a global to avoid below odd bug.
pl_kws = Dict{Symbol,Any}()

dir = "results_Werder/tmp"

for priority=1:4
    println("\n\n++++++++++++++++ PRIORITY $priority ++++++++++++++++++++++++++\n\n")
    for gid in I2.priorities[priority] # Note: parallel execution will block here until all processes have finished
                                       # the pmap.
        println("\n\n============================  Running Glacier $(BM.getname(gid))")
        #    try
        gl,gb,pp,pm,pn,pl = BM.init_forward(gid; pl_kws...);
        # I get an odd bug in pmap with ERROR: LoadError: On worker 2: UndefVarError: h_iv not defined
        if nprocs()>1
            pmap(exp->run_one_exp(exp, gl, gb, pp, pm, pn, pl, runtyp, fit_target; dir=dir), 0:I2.nexps)
        else
            map(exp->run_one_exp(exp, gl, gb, pp, pm, pn, pl, runtyp, fit_target; dir=dir), 0:I2.nexps)
        end
    end
end


text:
module GBGraph

using Graphs
using SuiteSparseGraphBLAS
using LinearAlgebra

include("GBEdge.jl")
include("GBGraph.jl")
end


text:
function WilkinsonTicks(k_ideal::Int; k_min = 2, k_max = 10,
        Q = [(1.0,1.0), (5.0, 0.9), (2.0, 0.7), (2.5, 0.5), (3.0, 0.2)],
        granularity_weight = 1/4,
        simplicity_weight = 1/6,
        coverage_weight = 1/3,
        niceness_weight = 1/4,
        min_px_dist = 50.0)

    if !(0 < k_min <= k_ideal <= k_max)
        error("Invalid tick number specifications k_ideal $k_ideal, k_min $k_min, k_max $k_max")
    end

    WilkinsonTicks(k_ideal, k_min, k_max, Q, granularity_weight, simplicity_weight,
        coverage_weight, niceness_weight, min_px_dist)
end

get_tick_labels(ticks::WilkinsonTicks, tickvalues) = linearly_spaced_tick_labels(tickvalues)

function compute_tick_values(ticks::WilkinsonTicks, vmin, vmax, pxwidth)

    min_px_dist = ticks.min_px_dist
    n_max_allowed_ticks = max(ticks.k_min, round(Int, pxwidth / min_px_dist))

    tickvalues, _ = PlotUtils.optimize_ticks(vmin, vmax;
        extend_ticks = false, strict_span=true, span_buffer = nothing,
        k_min = ticks.k_min,
        k_max = min(ticks.k_max, n_max_allowed_ticks),
        k_ideal = min(ticks.k_ideal, n_max_allowed_ticks),
        Q = ticks.Q,
        granularity_weight = ticks.granularity_weight,
        simplicity_weight = ticks.simplicity_weight,
        coverage_weight = ticks.coverage_weight,
        niceness_weight = ticks.niceness_weight)

    tickvalues
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 9f18ccd0-8e25-11eb-20e6-ab941f663550
begin
	import Pkg
    Pkg.activate(".")
    Pkg.add("PlutoUI")
	Pkg.add("Plots")
	Pkg.add("Distributions")

	using PlutoUI
	using Plots
	using Distributions
end

# ╔═╡ 10a9374c-8e1a-11eb-10b4-fb8272dd46d1
md"""
## Random SIRD Model with two Populations

#### Intro

The **SIRD Model** has been developed to simulate an epidemic over time. The model consists of a system of 4 differential equations that express the rates of change of 4 variables over time. The 4 variables are:

${ \begin{align*}
	S & - \text{the susceptibles of getting the infection} \\
	I & - \text{the infected} \\
	R & - \text{the recovered from the infection} \\
	D & - \text{the number of dead}
\end{align*} }$

#### The Model
Define the parameters of the dynamics as follows:

${ \begin{align*}
	\beta & \to \text{ rate of infection} \\
	\gamma & \to \text{ rate of recovery} \\
	\delta & \to \text{ rate of immunity loss} \\
	\rho & \to \text{ mortality}
\end{align*} }$

Then the following 4 equations govern the dynamis of the model:

${ \begin{align*}
\frac{\mathrm{d}S}{\mathrm{d}t} & = -\beta \cdot \frac{S \cdot I}{N} + \delta \cdot R \\
\frac{\mathrm{d}I}{\mathrm{d}t} & = \beta \cdot \frac{S \cdot I}{N} - (\gamma+\rho) \cdot I \\
\frac{\mathrm{d}R}{\mathrm{d}t} & = \gamma \cdot I - \delta \cdot R \\
\frac{\mathrm{d}D}{\mathrm{d}t} & = \rho \cdot I
\end{align*} }$

To model two mixing populations with independent parameters double the system and add resp. substract the mixing term

${ \pm \beta_{ij} \frac{S_j\cdot I_i}{N_i + N_j} }$

to the number of infected resp. susceptible individuals. This term models the infection from Population \\( i \\) to Population \\( j \\), for \\( i,j = 1,2\\).

Finally to randomize the dynamics of the system in every step, all quantities that move between within the system get drawn from a Poisson distribution with the respective rate as mean.

#### Model simulation
The following is a simulation of the model described above. First choose the model parameters
"""

# ╔═╡ 43970fa0-8e1b-11eb-3052-1701966a4478
begin
	T = 300 	#period of 300 daysplot
	Δt = 1/4 	#time interval of 6 hours (1/4 of a day)
end;

# ╔═╡ 742d1fd0-907d-11eb-27d1-59ca2af59f97
parameter = Dict(
	"rate of infection" => 0.14,
	"rate of recovery" => 0.07,
	"rate of immuity loss" => 0.1,
	"rate of death of infected" => 0.001,
	"total population" => 1*10^4,
	"initial number of infected" => 10
);

# ╔═╡ d59335ac-913b-11eb-228f-273be8d7a86a
md"""
---
**Additional parameter and functions for simulations with two subgroups**
"""

# ╔═╡ 56b9df4c-913b-11eb-180a-eb784c2b5181
begin
	parameter_subpopulation = Dict(
		"rate of infection" => 0.14,
		"rate of recovery" => 0.035,
		"rate of immuity loss" => 0.2,
		"rate of death of infected" => 0.01,
		"total population" => 1*10^2,
		"initial number of infected" => 0.0
		)

	mixing_parameter = Dict(
		"rate of infection 1 -> 2" => 0.14,
		"rate of infection 2 -> 1" => 0.14
	)
end;

# ╔═╡ f8e23e02-9142-11eb-1700-6974da1a09c0
md"""
---
#### Worksheet
Here all the calculations start. Set the initial state, choose the set of parameters and if needed the parameter changing function.
"""

# ╔═╡ da008842-9142-11eb-388c-d7e1f662d9db
md"""
---
#### Plots
"""

# ╔═╡ 93b1fcc4-9145-11eb-3d5e-4b9bd577a545
md"""
---
#### Model Functions

**Euler Algorithm to solve system. With constant or changing parameters**
"""

# ╔═╡ 2680b820-8e32-11eb-28de-6f039360524b
begin
	function euler(f::Function,change_par!::Function,t0,tn,Δt,x₀,par)
		T = t0:Δt:tn
		F = Matrix{eltype(x₀)}(undef,(length(T),length(x₀)))
		F[1,:] = x₀
		for (n,t) ∈ enumerate(T[2:end])
			change_par!(par,t,F,n)
			F[n+1,:] = max.(F[n,:] + f(t,F,n,par) .* Δt,0)
		end
		return F
	end

	function no_change!(par,t,F,n) end

	euler(f::Function,t0,tn,Δt,x₀,par) = euler(f,no_change!,t0,tn,Δt,x₀,par)
end;

# ╔═╡ 58f42866-9146-11eb-29ea-87830bc8db66
md"""
**Defines the infinitesimal change of system**
"""

# ╔═╡ 3ab68bec-8e2f-11eb-203e-f9530bc39073
begin
	function SIR_step(t,x,β,δ,γ,ρ,N)
		S, I, R, D = x[1], x[2], x[3], x[4]
		nI = rand(Poisson(β*I*S/N))
		nIL = rand(Poisson(δ*R))
		nR = rand(Poisson(γ*I))
		nD = rand(Poisson(ρ*I))
		return [
			-nI + nIL,
			nI - nR - nD,
			nR - nIL,
			nD
			]
	end

	SIR_step(t,x,par::Dict) = SIR_step(t,x,
				par["rate of infection"],par["rate of immuity loss"],
				par["rate of recovery"],par["rate of death of infected"],
				par["total population"]
			)

	SIR_step(t,H,n,par::Dict) = SIR_step(t,H[n,:],par)
end;

# ╔═╡ a7ac774a-913b-11eb-267b-8f9841c44cb4
begin
	function SIR_step_2D(t,H,n,par₁::Dict,par₂::Dict,par_mix::Dict)
		x₁,x₂ = view(H,n,1:4), view(H,n,5:8)
		x = vcat(SIR_step(t,x₁,par₁),SIR_step(t,x₂,par₂))
		return x + SIR_interaction(t,x₁,x₂,
					par_mix["rate of infection 1 -> 2"],
					par_mix["rate of infection 2 -> 1"],
					par₁["total population"]+par₂["total population"])
	end

	SIR_step_2D(t,H,n,par) = SIR_step_2D(t,H,n,par[1],par[2],par[3])

	function SIR_interaction(t,x₁,x₂,β₃,β₄,N)
		S₁, I₁, S₂, I₂ = x₁[1], x₁[2], x₂[1], x₂[2]
		nI₁ = rand(Poisson(β₄*S₁*I₂/N))
		nI₂ = rand(Poisson(β₃*S₂*I₁/N))
		return [
			-nI₁, nI₁, 0, 0,
			-nI₂, nI₂, 0, 0
		]
	end
end;

# ╔═╡ c6bbb678-9080-11eb-2047-a3dcd34e8c77
begin
	N₁,N₂ = parameter["total population"], parameter_subpopulation["total population"]
	I₁,I₂ = parameter["initial number of infected"], parameter_subpopulation["initial number of infected"]

	x₀ = [
		N₁-I₁,
		I₁,
		0.0,
		0.0,
		N₂-I₂,
		I₂,
		0.0,
		0.0
	]

	result = euler(
		SIR_step_2D,
		0.0,T,Δt,
		x₀,
		(parameter,parameter_subpopulation,mixing_parameter)
		)
end;

# ╔═╡ 3d9185a6-908c-11eb-301c-358ed4a3eb37
let
	p1 = plot(framestlye=:zerolines,title="Main Population")

	plot!(p1, 0.0:Δt:T,result[:,1:4],
		label = ["S" "I" "R" "D"],
		color = [:blue :red :green :yellow],
		xlabel = "Days",
		ylabel = "Number of Individuals"
	)

	p2 = plot(framestlye=:zerolines,title="Subpopulation")

	plot!(p2, 0.0:Δt:T,result[:,5:8],
		label = ["S" "I" "R" "D"],
		color = [:blue :red :green :yellow],
		xlabel = "Days",
		ylabel = "Number of Individuals"
	)

	plot(p1,p2,layout=(2,1),size=(680,680))
end

# ╔═╡ de101aba-9145-11eb-2457-7b1831145eeb
md"""

**Supporting function for calculating the Incidence**
"""

# ╔═╡ 16a50518-906d-11eb-1e2b-3faa1765f30c
begin
	#Calculates at one time within the simulation
	function incidence(nI,t,pop_size;ndays=7,scale=100_000)
		#set starting time in range
		t₀ = max(t-round(Int,ndays/Δt),1)
		#return Incidence
		return (nI[t] - nI[t₀]) .* (scale/pop_size)
	end

	#Calculates for all times after simulation
	function incidence(nI,pop_size;ndays=7,scale=100_000)
		#rescale days
		ndays = round(Int,ndays/Δt)
		#sum new cases in first days
		I = nI[1:ndays]
		#iterate and add incidence
		for t₀ ∈ 1:length(nI)-ndays
			append!(I,nI[t₀+ndays] - nI[t₀])
		end
		#scale and return
		return I .* scale/pop_size
	end
end;

# ╔═╡ ceedca8e-9148-11eb-189c-a3a088520c2c
md"""
---
Package Management
"""

# ╔═╡ Cell order:
# ╟─10a9374c-8e1a-11eb-10b4-fb8272dd46d1
# ╠═43970fa0-8e1b-11eb-3052-1701966a4478
# ╠═742d1fd0-907d-11eb-27d1-59ca2af59f97
# ╟─d59335ac-913b-11eb-228f-273be8d7a86a
# ╠═56b9df4c-913b-11eb-180a-eb784c2b5181
# ╟─f8e23e02-9142-11eb-1700-6974da1a09c0
# ╠═c6bbb678-9080-11eb-2047-a3dcd34e8c77
# ╟─da008842-9142-11eb-388c-d7e1f662d9db
# ╟─3d9185a6-908c-11eb-301c-358ed4a3eb37
# ╟─93b1fcc4-9145-11eb-3d5e-4b9bd577a545
# ╠═2680b820-8e32-11eb-28de-6f039360524b
# ╟─58f42866-9146-11eb-29ea-87830bc8db66
# ╠═3ab68bec-8e2f-11eb-203e-f9530bc39073
# ╠═a7ac774a-913b-11eb-267b-8f9841c44cb4
# ╟─de101aba-9145-11eb-2457-7b1831145eeb
# ╠═16a50518-906d-11eb-1e2b-3faa1765f30c
# ╟─ceedca8e-9148-11eb-189c-a3a088520c2c
# ╠═9f18ccd0-8e25-11eb-20e6-ab941f663550

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


text:
abstract type Kernel end

struct PolynomialKernel <: Kernel
    # return (x^T @ y)^M
    _degree::Int64 # degree of polynomial feature
    _intercept::Float64 # intercept
    function PolynomialKernel(degree::Int64, intercept = 0.0)
        new(degree, intercept)
    end
end

"""
    kernel(kernel::PolynomialKernel, x::AbstractArray{Float64, 2}, y::AbstractArray{Float64, 2})

For x = [x_1,... x_N] of size (n_feature, n_samples1) and y = [y_1,... y_N] of size (n_feature, n_samples2),
 compute the gram matrix K of size(n_samples1, n_samples2), where K_ij =(x_i * y_j)^M.
"""
function kernel(
    kernel::PolynomialKernel,
    x::AbstractArray{Float64,2},
    y::AbstractArray{Float64,2},
)
    n_samples1 = size(x)[2]
    n_samples2 = size(y)[2]
    n_feature = size(x)[1]

    intercept = kernel._intercept
    degree = kernel._degree

    n_feature == size(y)[1] || error("feature dimension does not match")
    ret = zeros(n_samples1, n_samples2)
    for i = 1:n_samples1
        for j = 1:n_samples2
            x_i = view(x, :, i)
            y_j = view(y, :, j)
            ret[i, j] = (sum(x_i .* y_j) + intercept)^degree
        end
    end

    return ret
end

"Kernel to compute c0 * exp(-0.5 * c1^T * (x1 - x2)^2) where c0 is scalar and c1 is Array{Float64, 1}"
mutable struct RBFKernel <: Kernel
    _params::Array{Float64,1}
    _n_params::Int64
    function RBFKernel(params::Array{Float64,1})
        new(copy(params), size(params)[1])
    end
end

function rbf(
    c0::Float64,
    c1::AbstractArray{Float64,1},
    x1::AbstractArray{Float64,1},
    x2::AbstractArray{Float64,1},
)
    d = (x1 - x2) .^ 2
    return c0 * exp(-0.5 * transpose(c1) * d)
end

"""
    kernel(kernel::PolynomialKernel, x::AbstractArray{Float64, 2}, y::AbstractArray{Float64, 2})

For x = [x_1,... x_N] of size (n_feature, n_samples1) and y = [y_1,... y_N] of size (n_feature, n_samples2),
 compute the gram matrix K of size(n_samples1, n_samples2), where K_ij = rbf(x_i, x_j).
"""
function kernel(kernel::RBFKernel, x::AbstractArray{Float64,2}, y::AbstractArray{Float64,2})
    n_feature = size(x)[1]
    n_samples1 = size(x)[2]
    n_samples2 = size(y)[2]
    n_feature == size(y)[1] || error("size does not match")
    n_feature == kernel._n_params - 1 || error("kernel size does not match")

    ret = zeros(n_samples1, n_samples2)
    for i = 1:n_samples1
        x_i = x[:, i]
        for j = 1:n_samples2
            y_j = y[:, j]
            ret[i, j] = rbf(kernel._params[1], kernel._params[2:end], x_i, y_j)
        end
    end

    return ret
end

"""
    derivatives(kernel::RBFKernel, x::AbstractArray{Float64, 2}, y::AbstractArray{Float64, 2})

Return the derivative of K(gram matrix) with respect to c0 and c1 as Array{Float64, 3}
"""
function derivatives(
    kernel::RBFKernel,
    x::AbstractArray{Float64,2},
    y::AbstractArray{Float64,2},
)
    n_feature = size(x)[1]
    n_samples1 = size(x)[2]
    n_samples2 = size(y)[2]
    n_feature == size(y)[1] || error("size does not match")
    n_feature == kernel._n_params - 1 || error("kernel size does not match")
    n_params = kernel._n_params

    ret = zeros(n_params, n_samples1, n_samples2)
    c0 = kernel._params[1]
    c1 = kernel._params[2:end]
    for i = 1:n_samples1
        for j = 1:n_samples2
            x_i = x[:, i]
            y_j = y[:, j]
            d = (x_i - y_j) .^ 2
            delta_c0 = exp(-0.5 * sum(c0 * d))
            delta_c1 = -0.5 * d * c0 * delta_c0
            ret[:, i, j] = vcat([delta_c0], delta_c1)
        end
    end

    return ret
end

function update_parameters(kernel::RBFKernel, updates_to_add::AbstractArray{Float64,1})
    kernel._params += updates_to_add
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _mag2geo(lat::T, long::T, pole_lat::T, pole_long::T) where {T<:AbstractFloat}
    r    = 1 # Distance from planet center.  Value unimportant -- just need
             # a length for conversion to rectangular coordinates
    sin_lat, cos_lat = sincos(deg2rad(lat))
    sin_long, cos_long = sincos(deg2rad(long))

    # convert to rectangular coordinates
    #   x-axis: defined by the vector going from Earth's center towards
    #        the intersection of the equator and Greenwich's meridian.
    #   z-axis: axis of the geographic poles
    #   y-axis: defined by y=z^x
    x = r * cos_lat * cos_long
    y = r * cos_lat * sin_long
    z = r * sin_lat

    # First rotation: in the plane of the current meridian from magnetic pole to
    # geographic pole.
    s, c = sincos(T(pi) / 2 - pole_lat)
    togeolat = SMatrix{3,3}(c,       zero(T),      -s,
                            zero(T),  one(T), zero(T),
                            s,       zero(T),       c)

    # Second rotation matrix: rotation around plane of the equator, from the
    # meridian containing the magnetic poles to the Greenwich meridian.
    sin_pole_long, cos_pole_long = sincos(pole_long)
    maglong2geolong = SMatrix{3,3}( cos_pole_long, sin_pole_long, zero(T),
                                   -sin_pole_long, cos_pole_long, zero(T),
                                    zero(T),             zero(T),  one(T))
    out = maglong2geolong * togeolat * SVector(x, y, z)

    geolat  = rad2deg(atan(out[3], hypot(out[1], out[2])))
    geolong = rad2deg(atan(out[2], out[1]))
    # I don't care about that one...just put it there for completeness' sake
    # magalt  = norm(out) - r
    return geolat, geolong
end

"""
    mag2geo(latitude, longitude[, year]) -> geographic_latitude, geographic_longitude

### Purpose ###

Convert from geomagnetic to geographic coordinates.

### Explanation ###

Converts from geomagnetic (latitude, longitude) to geographic (latitude,
longitude).  Altitude is not involved in this function.

### Arguments ###

* `latitude`: geomagnetic latitude (North), in degrees.
* `longitude`: geomagnetic longitude (East), in degrees.
* `year` (optional numerical argument): the year in which to perform conversion.
  If omitted, defaults to current year.

The coordinates can be passed as arrays of the same length.

### Output ###

The 2-tuple of geographic (latitude, longitude) coordinates, in degrees.

If geomagnetic coordinates are given as arrays, a 2-tuple of arrays of the same
length is returned.

### Example ###

Find position of North Magnetic Pole in 2016

```jldoctest
julia> using AstroLib

julia> mag2geo(90, 0, 2016)
(86.395, -166.29000000000002)
```

### Notes ###

This function uses list of North Magnetic Pole positions provided by World
Magnetic Model (https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy).

`geo2mag` converts geographic coordinates to geomagnetic coordinates.

Code of this function is based on IDL Astronomy User's Library.
"""
mag2geo(lat::Real, long::Real, year::Real=Dates.year(Dates.now())) =
    _mag2geo(promote(float(lat), float(long),
                     deg2rad(POLELATLONG[year][1]::Float64),
                     deg2rad(POLELATLONG[year][2]::Float64))...)

function mag2geo(lat::AbstractArray{LA}, long::AbstractArray{LO},
                 year::Real=Dates.year(Dates.now())) where {LA<:Real, LO<:Real}
    @assert length(lat) == length(long)
    typela   = float(LA)
    geolat   = similar(lat, typela)
    geolong  = similar(lat, typela)
    polelat  = deg2rad(POLELATLONG[year][1])
    polelong = deg2rad(POLELATLONG[year][2])
    for i in eachindex(lat)
        geolat[i], geolong[i] =
            _mag2geo(promote(float(lat[i]), float(long[i]),
                             polelat, polelong)...)
    end
    return geolat, geolong
end


text:
# MIRT/z-all.jl

include("../data/z-list.jl")

include("algorithm/z-list.jl")
include("fbp/z-list.jl")
include("io/z-list.jl")
include("plot/z-list.jl")
include("regularize/z-list.jl")
include("system/z-list.jl")
include("utility/z-list.jl")

include("../test/test_all_mirt.jl")
export test_all_mirt

tmp = homedir() * "/l/src/julia/mirt/um/"
if isdir(tmp) # UM-only tools
	include(tmp * "z-list.jl")
end


text:
using ImageMagick, Images, ColorTypes, FixedPointNumbers, FileIO
using FactCheck
using Compat

ontravis = haskey(ENV, "TRAVIS")

facts("IO") do
    workdir = joinpath(tempdir(), "Images")
    isdir(workdir) && rm(workdir, recursive=true)
    mkdir(workdir)

    context("Binary png") do
        a = rand(Bool,5,5)
        fn = joinpath(workdir, "5by5.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        a8 = convert(Array{Gray{UFixed8}}, a) # IM won't read back as Bool
        @fact convert(Array, b) --> a8
        aim = grayim(a)
        ImageMagick.save(fn, aim)
        b = ImageMagick.load(fn)
        @fact b --> copyproperties(aim, a8)
        a = bitrand(5,5)
        fn = joinpath(workdir, "5by5.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        a8 = convert(Array{Gray{UFixed8}}, a)
        @fact convert(Array, b) --> a8
        aim = grayim(a)
        ImageMagick.save(fn, aim)
        b = ImageMagick.load(fn)
        @fact b --> copyproperties(aim, a8)
    end

    context("Gray png") do
        a = [0 1/2^16 1/2^8; 1-1/2^8 1-1/2^16 1]
        aa = convert(Array{UFixed8}, a)
        fn = joinpath(workdir, "2by3.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        @fact convert(Array, b) --> aa
        ImageMagick.save(fn, aa)
        b = ImageMagick.load(fn)
        @fact convert(Array, b) --> aa
        open(fn, "w") do io
            @compat show(io, MIME("image/png"), b; minpixels=0)
        end
        bb = ImageMagick.load(fn)
        @fact bb.data --> b.data
        aaimg = Images.grayim(aa)
        ImageMagick.save(fn, aaimg)
        b = ImageMagick.load(fn)
        @fact b --> aaimg
        aa = convert(Array{UFixed16}, a)
        ImageMagick.save(fn, aa)
        b = ImageMagick.load(fn)
        @fact eltype(eltype(b)) --> UFixed16
        @fact convert(Array, b) --> aa
    end

    context("Color") do
        fn = joinpath(workdir, "2by2.png")
        A = rand(3,2,2)
        A[1] = 1
        img = Images.colorim(A)
        img24 = convert(Images.Image{RGB24}, img)
        ImageMagick.save(fn, img24)
        b = ImageMagick.load(fn)
        imgrgb8 = convert(Images.Image{RGB{UFixed8}}, img)
        @fact Images.data(imgrgb8) --> Images.data(b)

        open(fn, "w") do io
            @compat show(io, MIME("image/png"), imgrgb8; minpixels=0)
        end
        bb = ImageMagick.load(fn)
        @fact data(bb) --> data(imgrgb8)
    end

    context("Colormap usage") do
        datafloat = reshape(linspace(0.5, 1.5, 6), 2, 3)
        dataint = round(UInt8, 254*(datafloat .- 0.5) .+ 1)  # ranges from 1 to 255
        # build our colormap
        b = RGB(0,0,1)
        w = RGB(1,1,1)
        r = RGB(1,0,0)
        cmaprgb = Array(RGB{Float64}, 255)
        f = linspace(0,1,128)
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = Images.ImageCmap(dataint, cmaprgb)
        ImageMagick.save(joinpath(workdir,"cmap.jpg"), img)
        cmaprgb = Array(RGB, 255) # poorly-typed cmap, issue #336
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = Images.ImageCmap(dataint, cmaprgb)
        ImageMagick.save(joinpath(workdir, "cmap.png"), img)
    end

    context("Alpha") do
        c = reinterpret(BGRA{UFixed8}, [0xf0884422]'')
        fn = joinpath(workdir, "alpha.png")
    	ImageMagick.save(fn, c)
        C = ImageMagick.load(fn)
        if !ontravis
            # disabled on Travis because it has a weird, old copy of
            # ImageMagick for which this fails (see Images#261)
            @fact C[1] --> c[1]
        end
        ImageMagick.save(fn, reinterpret(ARGB32, [0xf0884422]''))
        D = ImageMagick.load(fn)
        if !ontravis
            @fact D[1] --> c[1]
        end

        # Images#396
        c = colorim(reshape(collect(0x00:0x11:0xff), 2, 2, 4), "RGBA")
        ImageMagick.save(fn, c)
        D = ImageMagick.load(fn)
        D8 = convert(Image{base_colorant_type(eltype(D)){UFixed8}}, D)
        C = permutedims(convert(Image{eltype(D8)}, c), spatialorder(D8))
        if !ontravis
            for i = 1:length(D8)
                @fact D8[i] --> C[i]
            end
        end
    end

    context("3D TIFF (issue #307)") do
        Ar = rand(0x00:0xff, 2, 2, 4)
        Ar[1] = 0xff
        A = Image(map(x->Gray(UFixed8(x,0)), Ar); colorspace="Gray", spatialorder=["x", "y"], timedim=3) # grayim does not use timedim, but load does...
        fn = joinpath(workdir, "3d.tif")
        ImageMagick.save(fn, A)
        B = ImageMagick.load(fn)

        @fact A --> B
    end

    context("Clamping (issue #256)") do
        Ar = rand(2,2)
        Ar[1] = 1
        A = grayim(Ar)
        A[1,1] = -0.4
        fn = joinpath(workdir, "2by2.png")
        println("The following warning is a sign of normal operation:")
        @fact_throws InexactError ImageMagick.save(fn, A)
        ImageMagick.save(fn, A, mapi=mapinfo(Images.Clamp, A))
        B = ImageMagick.load(fn)
        A[1,1] = 0
        @fact B --> map(Gray{UFixed8}, A)
    end

    is_unix() && context("Reading from a stream (issue #312)") do
        fn = joinpath(workdir, "2by2.png")
        img = open(fn) do io
            ImageMagick.load(io)
        end
        @fact isa(img, Images.Image) --> true
    end

    is_unix() && context("Writing to a stream (PR #22)") do
        orig_img = ImageMagick.load(joinpath(workdir, "2by2.png"))
        fn = joinpath(workdir, "2by2_fromstream.png")
        open(fn, "w") do f
            ImageMagick.save(Stream(format"PNG", f), orig_img)
        end
        img = ImageMagick.load(fn)
        @fact img --> orig_img
    end

    is_unix() && context("Reading from a byte array (issue #279)") do
        fn = joinpath(workdir, "2by2.png")
        io = open(fn)
        arr = read(io)
        close(io)
        img = readblob(arr)
        @fact isa(img, Images.Image) --> true
    end

    context("show(MIME)") do
        Ar = rand(UInt8, 3, 2, 2)
        Ar[1] = typemax(eltype(Ar))
        a = colorim(Ar)
        fn = joinpath(workdir, "2by2.png")
        open(fn, "w") do file
            @compat show(file, MIME("image/png"), a, minpixels=0)
        end
        b = ImageMagick.load(fn)
        @fact data(b) --> data(a)

        Ar = rand(UInt8, 3, 1021, 1026)
        Ar[1] = typemax(eltype(Ar))
        abig = colorim(Ar)
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            @compat show(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = ImageMagick.load(fn)
        @fact data(b) --> convert(Array{RGB{UFixed8},2}, data(restrict(abig, (1,2))))

        # Issue #269
        Ar = rand(UInt16, 3, 1024, 1023)
        Ar[1] = typemax(eltype(Ar))
        abig = colorim(Ar)
        open(fn, "w") do file
            @compat show(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = ImageMagick.load(fn)
        @fact data(b) --> convert(Array{RGB{UFixed8},2}, data(restrict(abig, (1,2))))
    end

end


text:
using .Sobol: SobolSeq
import .Sobol
using ..TupleVectors

export SobolHypercube

struct SobolHypercube{k} <: Hypercube{k}
    seq :: SobolSeq{k}
    value :: Vector{Float64}
    index :: Ref{Int}  # start at zero

    function SobolHypercube(k::Int)
        seq = SobolSeq(k)
        value = Sobol.next!(seq)
        return new{k}(seq, value, 0)
    end
end

export next!

function next!(s::SobolHypercube)
    s.index[] = 0
    Sobol.next!(s.seq, s.value)
end

function Base.rand(ω::SobolHypercube{k}) where {k}
    ω.value[ω.index[] += 1]
end

# function makeplot()
#     ω = SobolHypercube(100)
#     x = range(-2,2,length=100)
#     d = For(x) do xj Normal(2*xj, 1/(1 + xj^2)) end
#     y = rand(ω, d)

#     plt = scatter(x,y,markersize=2, color=colorant"rgba(0,0,0,0.1)")
#     for j in 1:100
#         next!(ω)
#         y = rand(ω, d)
#         scatter!(x, y, markersize=2, color=colorant"rgba(0,0,0,0.1)")
#     end
#     plt
# end


text:
###############################################################################
#
# Canopy radiation
#
###############################################################################
"""
    mutable struct CanopyRads{FT}

A struct for canopy radiation information

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct CanopyRads{FT}
    # local storage of dimension information
    "Number of azimuth angles"
    nAzi  ::Int = 36
    "Number of inclination agles"
    nIncl ::Int = 9
    "Number of canopy layers"
    nLayer::Int = 5
    "Number of canopy levels"
    nLevel::Int = nLayer+1
    "Number of wave lengths"
    nWL   ::Int = 10
    "Number of wave lengths for SIF"
    nWLF  ::Int = 10

    # Scalars
    "Integrated TOC outgoing flux `[W m⁻²]`"
    intEout            ::FT = 0
    "Incident spectrally integrated total PAR `[mol m⁻² s⁻¹]`"
    incomingPAR        ::FT = 0
    "Incident spectrally integrated direct PAR `[mol m⁻² s⁻¹]`"
    incomingPAR_direct ::FT = 0
    "Incident spectrally integrated diffuse PAR `[mol m⁻² s⁻¹]`"
    incomingPAR_diffuse::FT = 0
    "Net radiation of shaded soil `[W m⁻²]`"
    RnSoil_diffuse     ::FT = 0
    "Net Short-wave radiation of sunlit soil `[W m⁻²]`"
    RnSoil_direct      ::FT = 0
    "Net Short-wave radiation of soil (shaded + sunlit) `[W m⁻²]`"
    RnSoil             ::FT = 0
    "Net long-wave radiation of soil (shaded + sunlit) `[W m⁻²]`"
    RnSoilLW           ::FT = 0

    # Dim of nLayer
    "Net PAR of shaded leaves `[mol m⁻² s⁻¹]`"
    absPAR_shade   ::Array{FT,1} = zeros(FT, nLayer)
    "Net PAR by Cab+Car of shaded leaves `[moles m⁻² s⁻¹]`"
    absPAR_shadeCab::Array{FT,1} = zeros(FT, nLayer)
    "Spectrally integrated net absorbed direct radiation in each layer per leaf area `[W m⁻²]`"
    intNetSW_sunlit::Array{FT,1} = zeros(FT, nLayer)
    "Spectrally integrated net absorbed diffuse radiation in each layer per leaf area `[W m⁻²]`"
    intNetSW_shade ::Array{FT,1} = zeros(FT, nLayer)
    "Spectrally integrated net absorbed direct radiation in each layer per leaf area `[W m⁻²]`"
    intNetLW_sunlit::Array{FT,1} = zeros(FT, nLayer)
    "Spectrally integrated net absorbed diffuse radiation in each layer per leaf area `[W m⁻²]`"
    intNetLW_shade ::Array{FT,1} = zeros(FT, nLayer)
    "Leaf temperature (sunlit) `[K]`"
    T_sun          ::Array{FT,1} = zeros(FT, nLayer) .+ FT(298.15)
    "Leaf temperature (shaded) `[K]`"
    T_shade        ::Array{FT,1} = zeros(FT, nLayer) .+ FT(298.15)
    "Fluorescence yield for shaded leaves"
    ϕ_shade        ::Array{FT,1} =  ones(FT, nLayer) .* FT(0.01)
    "Sensible Heat flux H of shaded leaves `[W m⁻²]`"
    H_shade        ::Array{FT,1} = zeros(FT, nLayer)
    "Latent Heat flux LE of shaded leaves `[W m⁻²]`"
    LE_shade       ::Array{FT,1} = zeros(FT, nLayer)
    "NPQ of shaded leaves"
    NPQ_shade      ::Array{FT,1} = zeros(FT, nLayer)
    # TODO remove these?
    "GPP of shaded leaves `[μmol m⁻² s⁻¹]`"
    GPP_shade      ::Array{FT,1} = zeros(FT, nLayer)
    "gs of shaded leaves `[mol m⁻² s⁻¹]`"
    gs_shade       ::Array{FT,1} = zeros(FT, nLayer)
    "Leaf water potential of shaded leaves `[MPa]`"
    ψl_shade       ::Array{FT,1} = zeros(FT, nLayer)
    "Cc of shaded leaves `[µmol/mol]`"
    Cc_shade       ::Array{FT,1} = zeros(FT, nLayer)
    "internal CO₂ concentration of shaded leaves `[µmol/mol]`"
    Pi_shade       ::Array{FT,1} = zeros(FT, nLayer)

    # Dimension of wavelength
    "Short-wave TOC outgoing radiance in observation direction `[mW m⁻² nm⁻¹ sr⁻¹]`"
    Lo         ::Array{FT,1} = zeros(FT, nWL)
    "Short-wave TOC outgoing radiation `[mW m⁻² nm⁻¹]`"
    Eout       ::Array{FT,1} = zeros(FT, nWL)
    "Short-wave Albedo in viewing direction"
    alb_obs    ::Array{FT,1} = zeros(FT, nWL)
    "Short-wave Albedo for direct incoming radiation"
    alb_direct ::Array{FT,1} = zeros(FT, nWL)
    "Short-wave Albedo for diffuse incoming radiation"
    alb_diffuse::Array{FT,1} = zeros(FT, nWL)

    # Dimension of nLevel * nWavelengths
    "Upwelling diffuse short-wave radiation within canopy `[mW m⁻² nm⁻¹]`"
    E_up  ::Array{FT,2} = zeros(FT, (nWL,nLevel))
    "Downwelling diffuse short-wave radiation within canopy `[mW m⁻² nm⁻¹]`"
    E_down::Array{FT,2} = zeros(FT, (nWL,nLevel))

    # Dimension of nLayer * nWavelengths
    "Net absorbed direct radiation in each layer `[mW m⁻² nm⁻¹]`"
    netSW_sunlit   ::Array{FT,2} = zeros(FT, (nWL,nLayer))
    "net absorbed diffuse radiation in each layer `[mW m⁻² nm⁻¹]`"
    netSW_shade    ::Array{FT,2} = zeros(FT, (nWL,nLayer))

    # Dimension of nLeafInclination * nLeafAzimuth * nLayer
    "net PAR of sunlit leaves `[mol m⁻² s⁻¹]`"
    absPAR_sun   ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "net PAR by Cab+Car of sunlit leaves `[mol m⁻² s⁻¹]`"
    absPAR_sunCab::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "Leaf temperature (sunlit) `[K]`"
    T_sun3D      ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer)) .+ FT(298.15)
    "Fluorescence yield for sunlit leaves"
    ϕ_sun        ::Array{FT,3} =  ones(FT, (nIncl,nAzi,nLayer)) .* FT(0.01)
    "Sensible Heat flux H of sunlit leaves `[W m⁻²]`"
    H_sun        ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "Latent Heat flux LE of sunlit leaves `[W m⁻²]`"
    LE_sun       ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "NPQ of sunlit leaves"
    NPQ_sun      ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    # TODO remove these?
    "GPP of sunlit leaves `[μmol m⁻² s⁻¹]`"
    GPP_sun      ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "gs of sunlit leaves `[mol m⁻² s⁻¹]`"
    gs_sun       ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "Leaf water potential of sunlit leaves `[MPa]`"
    ψl_sun       ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "Cc of sunlit leaves `[µmol/mol]`"
    Cc_sun       ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))
    "Internal CO₂ concentration of sunlit leaves `[µmol/mol]`"
    Pi_sun       ::Array{FT,3} = zeros(FT, (nIncl,nAzi,nLayer))

    # Fluorescence Output:
    "Hemispheric total outgoing SIF flux `[mW m⁻² nm⁻¹]`)"
    SIF_hemi         ::Array{FT,1} = zeros(FT, nWLF)
    "Observer-direction outgoing SIF radiance (mW m⁻² nm⁻¹ sr⁻¹))"
    SIF_obs          ::Array{FT,1} = zeros(FT, nWLF)
    "Observer-direction outgoing SIF radiance, sunlit leaves (mW m⁻² nm⁻¹ sr⁻¹)"
    SIF_obs_sunlit   ::Array{FT,1} = zeros(FT, nWLF)
    "Observer-direction outgoing SIF radiance, shaded leaves (mW m⁻² nm⁻¹ sr⁻¹)"
    SIF_obs_shaded   ::Array{FT,1} = zeros(FT, nWLF)
    "Observer-direction outgoing SIF radiance, scattered (mW m⁻² nm⁻¹ sr⁻¹)"
    SIF_obs_scattered::Array{FT,1} = zeros(FT, nWLF)
    "Observer-direction outgoing SIF radiance, soil-reflected (mW m⁻² nm⁻¹ sr⁻¹)"
    SIF_obs_soil     ::Array{FT,1} = zeros(FT, nWLF)
    "Total SIF sum of layer sources  `[mW m⁻² nm⁻¹]`)"
    SIF_sum          ::Array{FT,1} = zeros(FT, nWLF)
end

max_stars_repo_name: climate-machine/Land


text:
include("multi_thread_env.jl")

"""
Many policy gradient based algorithms require that the `env` is a
`MultiThreadEnv` to increase the diversity during training. So the training
pipeline is different from the default one in `RLCore`.
"""
function RLCore._run(
    policy::AbstractPolicy,
    env::MultiThreadEnv,
    stop_condition,
    hook::AbstractHook = EmptyHook(),
)

    while true
        reset!(env)  # this is a soft reset!, only environments reached the end will be reset.
        action = policy(env)
        policy(PRE_ACT_STAGE, env, action)
        hook(PRE_ACT_STAGE, policy, env, action)

        env(action)
        policy(POST_ACT_STAGE, env)
        hook(POST_ACT_STAGE, policy, env)

        if stop_condition(policy, env)
            break
        end
    end
    action = policy(env)
    policy(PRE_ACT_STAGE, env, action)  # let the policy see the last observation
    hook(PRE_ACT_STAGE, policy, env, action)
    nothing
end


text:
module AccretionDiscProfiles

using Parameters

import StaticArrays: @SVector

# this should be an optional dependency
# as should CarterBoyerLindquist
using ComputedGeodesicEquations

import GeodesicBase: AbstractMetricParams
import GeodesicTracer: tracegeodesics
import AccretionGeometry: AbstractAccretionGeometry, tracegeodesics

import GeometricalPredicates: Point2D, getx, gety, geta, getb, getc

include("sky-geometry.jl")
include("corona-models.jl")

function tracegeodesics(
    m::AbstractMetricParams{T},
    model::AbstractCoronaModel{T},
    time_domain::Tuple{T,T};
    n_samples = 1024,
    sampler = WeierstrassSampler(res = 100.0),
    kwargs...,
) where {T}
    us = sample_position(m, model, n_samples)
    vs = sample_velocity(m, model, sampler, us, n_samples)
    tracegeodesics(m, us, vs, time_domain; kwargs...)
end

function tracegeodesics(
    m::AbstractMetricParams{T},
    model::AbstractCoronaModel{T},
    time_domain::Tuple{T,T},
    d::AbstractAccretionGeometry{T},
    ;
    n_samples = 1024,
    sampler = WeierstrassSampler(res = 100.0),
    kwargs...,
) where {T}
    us = sample_position(m, model, n_samples)
    vs = sample_velocity(m, model, sampler, us, n_samples)
    tracegeodesics(m, us, vs, time_domain, d; kwargs...)
end

function renderprofile(
    m::AbstractMetricParams{T},
    model::AbstractCoronaModel{T},
    max_time::T,
    d::AbstractAccretionGeometry{T};
    n_samples = 2048,
    sampler = WeierstrassSampler(res = 100.0),
    kwargs...,
) where {T}
    __renderprofile(m, model, d, n_samples, (0.0, max_time); kwargs...)
end

export AbstractCoronaModel, tracegeodesics, renderprofile

end # module


text:
module PopMemberModule

import ..CoreModule: Options, Dataset, Node, copy_node
import ..UtilsModule: get_birth_order
import ..LossFunctionsModule: score_func

# Define a member of population by equation, score, and age
mutable struct PopMember{T<:Real}
    tree::Node
    score::T  # Inludes complexity penalty, normalization
    loss::T  # Raw loss
    birth::Int

    # For recording history:
    ref::Int
    parent::Int
end

generate_reference() = abs(rand(Int))

"""
    PopMember(t::Node, score::T, loss::T)

Create a population member with a birth date at the current time.

# Arguments

- `t::Node`: The tree for the population member.
- `score::T`: The score (normalized to a baseline, and offset by a complexity penalty)
- `loss::T`: The raw loss to assign.
"""
function PopMember(
    t::Node, score::T, loss::T; ref::Int=-1, parent::Int=-1, deterministic=false
) where {T<:Real}
    if ref == -1
        ref = generate_reference()
    end
    return PopMember{T}(
        t, score, loss, get_birth_order(; deterministic=deterministic), ref, parent
    )
end

"""
    PopMember(dataset::Dataset{T}, baseline::T,
              t::Node, options::Options)

Create a population member with a birth date at the current time.
Automatically compute the score for this tree.

# Arguments

- `dataset::Dataset{T}`: The dataset to evaluate the tree on.
- `baseline::T`: The baseline loss.
- `t::Node`: The tree for the population member.
- `options::Options`: What options to use.
"""
function PopMember(
    dataset::Dataset{T},
    baseline::T,
    t::Node,
    options::Options;
    ref::Int=-1,
    parent::Int=-1,
    deterministic=nothing,
) where {T<:Real}
    score, loss = score_func(dataset, baseline, t, options)
    return PopMember(t, score, loss; ref=ref, parent=parent, deterministic=deterministic)
end

function copy_pop_member(p::PopMember{T}) where {T<:Real}
    tree = copy_node(p.tree)
    score = copy(p.score)
    loss = copy(p.loss)
    birth = copy(p.birth)
    ref = copy(p.ref)
    parent = copy(p.parent)
    return PopMember{T}(tree, score, loss, birth, ref, parent)
end

end


text:
abstract type AbstractChargeDriftModel{T <: SSDFloat} end
abstract type AbstractTemperatureModel{T <: SSDFloat} end

"""
    get_electron_drift_field(ef::Array{SVector{3, T},3}, chargedriftmodel::AbstractChargeDriftModel)::Array{SVector{3,T},3} where {T <: SSDFloat}

Applies the charge drift model onto the electric field vectors. The field vectors have to be in cartesian coordinates.
"""
function get_electron_drift_field(ef::Array{SVector{3, T},3}, chargedriftmodel::AbstractChargeDriftModel;
            use_nthreads::Int = Base.Threads.nthreads())::Array{SVector{3,T},3} where {T <: SSDFloat}
    df = Array{SVector{3,T}, 3}(undef, size(ef))
    @onthreads 1:use_nthreads for i3 in workpart(axes(df, 3), 1:use_nthreads, Base.Threads.threadid())
        for i2 in axes(df, 2)
            for i1 in axes(df, 1)
                @inbounds df[i1, i2, i3] = getVe(ef[i1, i2, i3], chargedriftmodel)
            end
        end
    end
    return df
end

"""
    get_hole_drift_field(ef::Array{SVector{3, T},3}, chargedriftmodel::AbstractChargeDriftModel)::Array{SVector{3,T},3} where {T <: SSDFloat}

Applies the charge drift model onto the hole field vectors. The field vectors have to be in cartesian coordinates.
"""
function get_hole_drift_field(ef::Array{SVector{3,T},3}, chargedriftmodel::AbstractChargeDriftModel;
            use_nthreads::Int = Base.Threads.nthreads())::Array{SVector{3,T},3} where {T <: SSDFloat}
    df = Array{SVector{3,T}, 3}(undef, size(ef))
    @onthreads 1:use_nthreads for i3 in workpart(axes(df, 3), 1:use_nthreads, Base.Threads.threadid())
        for i2 in axes(df, 2)
            for i1 in axes(df, 1)
                @inbounds df[i1, i2, i3] = getVh(ef[i1, i2, i3], chargedriftmodel)
            end
        end
    end
    return df
end

include("Vacuum/Vacuum.jl")
include("ADL/ADL.jl")


text:
# Quadrotor in Maze

function QuadrotorMaze(::Type{Rot}; use_rot=true, costfun=:Quadratic,
        normcon=false) where Rot<:Rotation

# model
model = Dynamics.Quadrotor2{Rot}(use_rot=use_rot)
n,m = size(model)

N = 101 # number of knot points
tf = 5.0
dt = tf/(N-1) # total time

x0 = Dynamics.build_state(model, [0,-15,1], I(UnitQuaternion), zeros(3), zeros(3))
xf = Dynamics.build_state(model, [0, 15,1], I(UnitQuaternion), zeros(3), zeros(3))
utrim = Dynamics.trim_controls(model)

# cost
if n == 13
    rm_quat = @SVector [1,2,3,4,5,6,8,9,10,11,12,13]
else
    rm_quat = @SVector [1,2,3,4,5,6,7,8,9,10,11,12]
end
costfun == :QuatLQR ? sq = 0 : sq = 1
Q_diag = Dynamics.fill_state(model, 1e-3, 1e-2*sq, 1e-3, 1e-3)
R_diag = @SVector fill(1e-4,m)
Q = Diagonal(Q_diag)
R = Diagonal(R_diag)
Qf = Diagonal(@SVector fill(1e3,n))

if costfun == :Quadratic
    cost = LQRCost(Q, R, xf)
    obj = Objective(cost, N)
    # obj = LQRObjective(Q, R, Qf, xf, N) # objective with same stagewise costs
elseif costfun == :QuatLQR
    cost = QuatLQRCost(Q, R, xf, w=1e-1)
    cost_term = QuatLQRCost(Qf, R, xf, w=1.0)
    obj = Objective(cost, cost_term, N)
elseif costfun == :ErrorQuad
    cost = ErrorQuadratic(model, Diagonal(Q_diag[rm_quat]), R, xf)
    cost_term = ErrorQuadratic(model, Diagonal(diag(Qf)[rm_quat]), R, xf)
    obj = Objective(cost, cost_term, N)
end

# constraints
r_quad_maze = 1.0
r_cylinder_maze = 1.2
maze_cylinders = []
zh = 3
l1 = 5
l2 = 4
l3 = 7
l4 = 10

d = 3  # 0.5 door width
w = 10  # y location of wall
mid = 3  # 0.5 middle width

x_enter=-10
x_mid=0
x_exit=10

for i = range(-w,stop=-d,length=l1) # enter wall
    push!(maze_cylinders,(i, x_enter,r_cylinder_maze))
end

for i = range(d,stop=w,length=l1) # enter wall
    push!(maze_cylinders,(i, x_enter, r_cylinder_maze))
end

for i = range(-mid,stop=mid,length=l3) # middle obstacle
    push!(maze_cylinders,(i, x_mid, r_cylinder_maze))
end

for i = range(-w,stop=-d,length=l1) # exit wall
    push!(maze_cylinders,(i, x_exit, r_cylinder_maze))
end

for i = range(d,stop=w,length=l1) # exit wall
    push!(maze_cylinders,(i, x_exit, r_cylinder_maze))
end

# Top and bottom walls
for i = range(x_enter+2*r_cylinder_maze,stop=x_exit-2*r_cylinder_maze,length=l4)
    push!(maze_cylinders,(-w, i, r_cylinder_maze))
end

for i = range(x_enter+2*r_cylinder_maze,stop=x_exit-2*r_cylinder_maze,length=l4)
    push!(maze_cylinders,(w, i, r_cylinder_maze))
end

n_maze_cylinders = length(maze_cylinders)
maze_xyr = collect(zip(maze_cylinders...))

cx = SVector{n_maze_cylinders}(maze_xyr[1])
cy = SVector{n_maze_cylinders}(maze_xyr[2])
cr = SVector{n_maze_cylinders}(maze_xyr[3])

obs = CircleConstraint(n, cx, cy, cr .+ r_quad_maze)

u_min = 0.
u_max = 50.
x_max = Inf*ones(n)
x_min = -Inf*ones(n)

x_max[1:3] = [25.0; Inf; 20]
x_min[1:3] = [-25.0; -Inf; 0.]

if n == 13
    noquat = @SVector [1,2,3,8,9,10,11,12,13]
    x_max[8:10] .= 30.
else
    noquat = @SVector [1,2,3,7,8,9,10,11,12]
    x_max[7:9] .= 30.
end

bnd1 = BoundConstraint(n,m,u_min=u_min,u_max=u_max)
bnd2 = BoundConstraint(n,m,u_min=u_min,u_max=u_max,x_min=x_min,x_max=x_max)
goal = GoalConstraint(xf, noquat)

u0 = @SVector fill(0.5*9.81/4.0, m)
U_hover = [copy(u0) for k = 1:N-1] # initial hovering control trajectory

conSet = ConstraintSet(n,m,N)
add_constraint!(conSet, obs, 1:N-1)
add_constraint!(conSet, bnd2, 2:N-1)
add_constraint!(conSet, goal, N:N)

# Quaternion norm constraint
if normcon
    if use_rot == :slack
        add_constraint!(conSet, QuatSlackConstraint(), 1:N-1)
    else
        add_constraint!(conSet, QuatNormConstraint(), 1:N-1)
        u0 = [u0; (@SVector [1.])]
    end
end

quadrotor_maze = Problem(model, obj, xf, tf, x0=x0, constraints=conSet)
initial_controls!(quadrotor_maze,U_hover); # initialize problem with controls

X_guess = zeros(n,7)
X_guess[:,1] = x0
X_guess[:,7] = xf
wpts = [2 -7 1;
        5 -5 1;
        7  0 1;
        5  5 1;
        2  7 1;]
X_guess[1:3,2:6] .= wpts'

if n == 13
    X_guess[4:7,:] .= q0
end
X0 = interp_rows(N,tf,X_guess);
initial_states!(quadrotor_maze, X0)

quadrotor_maze_objects = (cx,cy,cr)

return quadrotor_maze
end


text:
############################################################################################
### models w.o. meanstructure
############################################################################################

# observed ---------------------------------------------------------------------------------
observed = SemObservedData(specification = spec, data = dat)

# imply
imply_ram = RAM(specification = spec)

imply_ram_sym = RAMSymbolic(specification = spec)

# loss functions ---------------------------------------------------------------------------
ml = SemML(observed = observed)

wls = SemWLS(observed = observed)

ridge = SemRidge(α_ridge = .001, which_ridge = 16:20, n_par = 31)

constant = SemConstant(constant_loss = 3.465)

# loss -------------------------------------------------------------------------------------
loss_ml = SemLoss(ml)

loss_wls = SemLoss(wls)

# optimizer -------------------------------------------------------------------------------------
optimizer_obj = semoptimizer()

# models -----------------------------------------------------------------------------------

model_ml = Sem(observed, imply_ram, loss_ml, optimizer_obj)

model_ls_sym = Sem(observed, RAMSymbolic(specification = spec, vech = true), loss_wls, optimizer_obj)

model_ml_sym = Sem(observed, imply_ram_sym, loss_ml, optimizer_obj)

model_ridge = Sem(observed, imply_ram, SemLoss(ml, ridge), optimizer_obj)

model_constant = Sem(observed, imply_ram, SemLoss(ml, constant), optimizer_obj)

model_ml_weighted = Sem(observed, imply_ram, SemLoss(ml; loss_weights = [n_obs(model_ml)]), optimizer_obj)

############################################################################################
### test gradients
############################################################################################

models = [model_ml, model_ls_sym, model_ridge, model_constant, model_ml_sym, model_ml_weighted]
model_names = ["ml", "ls_sym", "ridge", "constant", "ml_sym", "ml_weighted"]

for (model, name) in zip(models, model_names)
    try
        @testset "$(name)_gradient" begin
            @test test_gradient(model, start_test; rtol = 1e-9)
        end
    catch
    end
end

############################################################################################
### test solution
############################################################################################

models = [model_ml, model_ls_sym, model_ml_sym, model_constant]
model_names = ["ml", "ls_sym", "ml_sym", "constant"]
solution_names = Symbol.("parameter_estimates_".*["ml", "ls", "ml", "ml"])

for (model, name, solution_name) in zip(models, model_names, solution_names)
    try
        @testset "$(name)_solution" begin
            solution = sem_fit(model)
            update_estimate!(partable, solution)
            @test compare_estimates(partable, solution_lav[solution_name]; atol = 1e-2)
        end
    catch
    end
end

@testset "ridge_solution" begin
    solution_ridge = sem_fit(model_ridge)
    solution_ml = sem_fit(model_ml)
    # solution_ridge_id = sem_fit(model_ridge_id)
    @test abs(solution_ridge.minimum - solution_ml.minimum) < 1
end

# test constant objective value
@testset "constant_objective_and_gradient" begin
    @test (objective!(model_constant, start_test) - 3.465) ≈ objective!(model_ml, start_test)
    grad = similar(start_test); grad2 = similar(start_test)
    gradient!(grad, model_constant, start_test)
    gradient!(grad2, model_ml, start_test)
    @test grad ≈ grad2
end

@testset "ml_solution_weighted" begin
    solution_ml = sem_fit(model_ml)
    solution_ml_weighted = sem_fit(model_ml_weighted)
    @test isapprox(solution(solution_ml), solution(solution_ml_weighted), rtol = 1e-3)
    @test isapprox(n_obs(model_ml)*StructuralEquationModels.minimum(solution_ml),
        StructuralEquationModels.minimum(solution_ml_weighted), rtol = 1e-6)
end

############################################################################################
### test fit assessment
############################################################################################

@testset "fitmeasures/se_ml" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml];
        atol = 1e-3))

    update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable, solution_lav[:parameter_estimates_ml];
        atol = 1e-3, col = :se, lav_col = :se)
end

@testset "fitmeasures/se_ls" begin
    solution_ls = sem_fit(model_ls_sym)
    fm = fit_measures(solution_ls)
    @test all(test_fitmeasures(fm, solution_lav[:fitmeasures_ls]; atol = 1e-3,
        fitmeasure_names = fitmeasure_names_ls))
    @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing)

    update_partable!(partable, identifier(model_ls_sym), se_hessian(solution_ls), :se)
    @test compare_estimates(partable, solution_lav[:parameter_estimates_ls]; atol = 1e-2,
        col = :se, lav_col = :se)
end

############################################################################################
### test hessians
############################################################################################

if semoptimizer == SemOptimizerOptim
    using Optim, LineSearches

    optimizer_obj = SemOptimizerOptim(
        algorithm = Newton(
            ;linesearch = BackTracking(order=3),
            alphaguess = InitialHagerZhang()
        )
    )

    imply_sym_hessian_vech = RAMSymbolic(specification = spec, vech = true, hessian = true)

    imply_sym_hessian = RAMSymbolic(specification = spec, hessian = true)

    model_ls = Sem(observed, imply_sym_hessian_vech, loss_wls, optimizer_obj)

    model_ml = Sem(observed, imply_sym_hessian, loss_ml, SemOptimizerOptim(algorithm = Newton()))

    @testset "ml_hessians" begin
        @test test_hessian(model_ml, start_test; atol = 1e-4)
    end

    @testset "ls_hessians" begin
        @test test_hessian(model_ls, start_test; atol = 1e-4)
    end

    @testset "ml_solution_hessian" begin
        solution = sem_fit(model_ml)
        update_estimate!(partable, solution)
        @test compare_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3)
    end

    @testset "ls_solution_hessian" begin
        solution = sem_fit(model_ls)
        update_estimate!(partable, solution)
        @test compare_estimates(partable, solution_lav[:parameter_estimates_ls]; atol = 1e-3)
    end

end

############################################################################################
### meanstructure
############################################################################################

# observed ---------------------------------------------------------------------------------
observed = SemObservedData(specification = spec_mean, data = dat, meanstructure = true)

# imply
imply_ram = RAM(specification = spec_mean, meanstructure = true)

imply_ram_sym = RAMSymbolic(specification = spec_mean, meanstructure = true)

# loss functions ---------------------------------------------------------------------------
ml = SemML(observed = observed, meanstructure = true)

wls = SemWLS(observed = observed, meanstructure = true)

# loss -------------------------------------------------------------------------------------
loss_ml = SemLoss(ml)

loss_wls = SemLoss(wls)

# optimizer -------------------------------------------------------------------------------------
optimizer_obj = semoptimizer()

# models -----------------------------------------------------------------------------------
model_ml = Sem(observed, imply_ram, loss_ml, optimizer_obj)

model_ls = Sem(
    observed,
    RAMSymbolic(specification = spec_mean, meanstructure = true, vech = true),
    loss_wls,
    optimizer_obj)

model_ml_sym = Sem(observed, imply_ram_sym, loss_ml, optimizer_obj)

############################################################################################
### test gradients
############################################################################################

models = [model_ml, model_ls, model_ml_sym]
model_names = ["ml", "ls_sym", "ml_sym"]

for (model, name) in zip(models, model_names)
    try
        @testset "$(name)_gradient_mean" begin
            @test test_gradient(model, start_test_mean; rtol = 1e-9)
        end
    catch
    end
end

############################################################################################
### test solution
############################################################################################

solution_names = Symbol.("parameter_estimates_".*["ml", "ls", "ml"].*"_mean")

for (model, name, solution_name) in zip(models, model_names, solution_names)
    try
        @testset "$(name)_solution_mean" begin
            solution = sem_fit(model)
            update_estimate!(partable_mean, solution)
            @test compare_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2)
        end
    catch
    end
end

############################################################################################
### test fit assessment
############################################################################################

@testset "fitmeasures/se_ml_mean" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml_mean];
        atol = 1e-3))

    update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_ml_mean];
        atol = 1e-3, col = :se, lav_col = :se)
end

@testset "fitmeasures/se_ls_mean" begin
    solution_ls = sem_fit(model_ls)
    fm = fit_measures(solution_ls)
    @test all(test_fitmeasures(fm,
        solution_lav[:fitmeasures_ls_mean];
        atol = 1e-3,
        fitmeasure_names = fitmeasure_names_ls))
    @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing)

    update_partable!(partable_mean, identifier(model_ls), se_hessian(solution_ls), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_ls_mean]; atol = 1e-2, col = :se, lav_col = :se)
end

############################################################################################
### fiml
############################################################################################

observed = SemObservedMissing(specification = spec_mean, data = dat_missing)

fiml = SemFIML(observed = observed, specification = spec_mean)

loss_fiml = SemLoss(fiml)

model_ml = Sem(observed, imply_ram, loss_fiml, optimizer_obj)

model_ml_sym = Sem(observed, imply_ram_sym, loss_fiml, optimizer_obj)

############################################################################################
### test gradients
############################################################################################

@testset "fiml_gradient" begin
    @test test_gradient(model_ml, start_test_mean; atol = 1e-6)
end

@testset "fiml_gradient_symbolic" begin
    @test test_gradient(model_ml_sym, start_test_mean; atol = 1e-6)
end

############################################################################################
### test solution
############################################################################################

@testset "fiml_solution" begin
    solution = sem_fit(model_ml)
    update_estimate!(partable_mean, solution)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
end

@testset "fiml_solution_symbolic" begin
    solution = sem_fit(model_ml_sym)
    update_estimate!(partable_mean, solution)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
end

############################################################################################
### test fit measures
############################################################################################

@testset "fitmeasures/se_fiml" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_fiml];
        atol = 1e-3))

    update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml];
        atol = 1e-3, col = :se, lav_col = :se)
end


text:
using Nextion
using Test

include("test_constants.jl")
include("test_object.jl")
include("test_events.jl")
include("test_range_number.jl")
# include("test_examples.jl")  # can't be run without display