text:
function deserialize_tensor_summary(summary)
    metadata    = summary.metadata
    plugin_data = metadata.plugin_data

    plugin_name = plugin_data.plugin_name

    dim_szs, dim_nms = deserialize_tensor_shape(summary.tensor)

    if plugin_name == "text"
        val = deserialize_text_tensor(summary, dim_szs, dim_nms)
    end
    return val
end

function deserialize_text_tensor(summary, sizes, names)
    version_data = readproto(PipeBuffer(summary.metadata.plugin_data.content),
                             TextPluginData())
    tensor = summary.tensor
    strings = Array{String, length(sizes)}(undef, sizes...)

    # Ignore dimension names
    if version_data.version == 0
        for (i, in) = enumerate(tensor.string_val)
            #TODO should undo `markdown_repr`, because all strings
            # will be displayed with a starting and ending \"
            strings[i] = String(in)
            delimiter = '\"'
            if delimiter === first(strings[i]) === last(strings[i])
                strings[i] = strings[i][2:end-1]
            end
        end
    else
        @warn "Unknown string version $(version_data.version)"
    end

    return strings
end

function deserialize_tensor_shape(tensor)
    t_shape = tensor.tensor_shape

    dim_names = Union{String,Nothing}[]
    dim_sz    = Int[]
    for dim = t_shape.dim
        push!(dim_sz, dim.size)
        isempty(dim.name) ? push!(dim_names, nothing) : push!(dim_names, dim.name)
    end

    return tuple(dim_sz...), tuple(dim_names...)
end


text:
function wrap_api_call(spec::SpecFunc, args; with_func_ptr = false)
    ex = :($(spec.name)($((with_func_ptr ? [args; first(func_ptrs(spec))] : args)...)))
    ex = if with_func_ptr
        ex
    else
        maybe_dispatch(spec, ex)
    end
    wrap_return(
        ex,
        spec.return_type,
        idiomatic_julia_type(spec.return_type),
    )
end

function maybe_dispatch(spec::SpecFunc, ex)
    maybe_handle = !isempty(children(spec)) ? innermost_type(first(children(spec)).type) : nothing
    use_dispatch_macro = spec.name ∉ (:vkGetInstanceProcAddr, :vkGetDeviceProcAddr)
    use_dispatch_macro || return ex
    if maybe_handle in spec_handles.name
        handle = handle_by_name(maybe_handle)
        handle_id = wrap_identifier(handle)
        hierarchy = parent_hierarchy(handle)
        if handle.name == :VkDevice || handle.name == :VkInstance
            # to avoid name conflicts
            :(@dispatch $handle_id $ex)
        elseif :VkDevice in hierarchy
            :(@dispatch device($handle_id) $ex)
        elseif :VkInstance in hierarchy
            :(@dispatch instance($handle_id) $ex)
        end
    else
        :(@dispatch nothing $ex)
    end
end

function wrap_enumeration_api_call(spec::SpecFunc, exs::Expr...; free = [])
    if must_repeat_while_incomplete(spec)
        if !isempty(free)
            free_block = quote
                if _return_code == VK_INCOMPLETE
                    $(map(x -> :(Libc.free($(x.name))), free)...)
                end
            end
            [:(@repeat_while_incomplete $(Expr(:block, exs..., free_block)))]
        else
            [:(@repeat_while_incomplete $(Expr(:block, exs...)))]
        end
    else
        exs
    end
end

function APIFunction(spec::SpecFunc, with_func_ptr)
    p = Dict(
        :category => :function,
        :name => nc_convert(SnakeCaseLower, remove_vk_prefix(spec.name)),
        :relax_signature => is_promoted,
    )

    count_ptr_index = findfirst(x -> (is_length(x) || is_size(x)) && x.requirement == POINTER_REQUIRED, children(spec))
    queried_params = getindex(children(spec), findall(is_implicit_return, children(spec)))
    if !isnothing(count_ptr_index)
        count_ptr = children(spec)[count_ptr_index]
        queried_params =
            getindex(children(spec), findall(x -> x.len == count_ptr.name && !x.is_constant, children(spec)))

        first_call_args = map(@λ(begin
            &count_ptr => count_ptr.name
            GuardBy(in(queried_params)) => :C_NULL
            x => vk_call(x)
        end), children(spec))

        i = 0
        second_call_args = map(@λ(begin
            :C_NULL && Do(i += 1) => queried_params[i].name
            x => x
        end), first_call_args)

        p[:body] = concat_exs(
            initialize_ptr(count_ptr),
            wrap_enumeration_api_call(
                spec,
                wrap_api_call(spec, first_call_args; with_func_ptr),
                (is_length(count_ptr) ? initialize_array : initialize_ptr).(queried_params, count_ptr)...,
                wrap_api_call(spec, second_call_args; with_func_ptr),
                ;
                free = filter(is_data, queried_params),
            )...,
            wrap_implicit_return(spec, queried_params; with_func_ptr),
        )

        args = filter(!in(vcat(queried_params, count_ptr)), children(spec))

        ret_type = @match length(queried_params) begin
            if any(is_data_with_retrievable_size, queried_params)
            end => Expr(:curly, :Tuple, idiomatic_julia_type.([unique(len.(queried_params)); queried_params])...)
            1 => idiomatic_julia_type(first(queried_params))
            _ => Expr(
                :curly,
                :Tuple,
                (idiomatic_julia_type(param) for param ∈ queried_params)...,
            )
        end
    elseif !isempty(queried_params)
        call_args = map(@λ(begin
            x && GuardBy(in(queried_params)) => x.name
            x => vk_call(x)
        end), children(spec))

        p[:body] = concat_exs(
            map(initialize_ptr, queried_params)...,
            wrap_api_call(spec, call_args; with_func_ptr),
            wrap_implicit_return(spec, queried_params; with_func_ptr),
        )

        args = filter(!in(filter(x -> x.requirement ≠ POINTER_REQUIRED, queried_params)), children(spec))

        ret_type = @match length(queried_params) begin
            1 => idiomatic_julia_type(first(queried_params))
            _ => Expr(:curly, :Tuple, (idiomatic_julia_type(param) for param ∈ queried_params)...)
        end
    else
        p[:short] = true
        p[:body] = :($(wrap_api_call(spec, map(vk_call, children(spec)); with_func_ptr)))

        args = children(spec)

        ret_type = idiomatic_julia_type(spec.return_type)
    end

    add_func_args!(p, spec, args; with_func_ptr)
    p[:return_type] = wrap_return_type(spec, ret_type)
    APIFunction(spec, with_func_ptr, p)
end

"""
Extend functions that create (or allocate) one or several handles,
by exposing the parameters of the associated CreateInfo structures.
`spec` must have one or several CreateInfo arguments.
"""
function APIFunction(spec::CreateFunc, with_func_ptr)
    @assert !isnothing(spec.create_info_param) "Cannot extend handle constructor with no create info parameter."
    def = APIFunction(spec.func, false)
    p_func = def.p
    p_info = Constructor(StructDefinition{false}(spec.create_info_struct)).p

    args = [p_func[:args]; p_info[:args]]
    kwargs = [p_func[:kwargs]; p_info[:kwargs]]

    info_expr = reconstruct_call(p_info; is_decl = false)
    info_index = findfirst(==(spec.create_info_param), filter(!is_optional, children(spec.func)))
    deleteat!(args, info_index)

    func_call_args::Vector{ExprLike} = name.(p_func[:args])
    func_call_args[info_index] = info_expr

    if with_func_ptr
        append!(args, func_ptr_args(spec.func))
        append!(func_call_args, func_ptrs(spec.func))
    end

    body = reconstruct_call(Dict(:name => name(def), :args => func_call_args, :kwargs => name.(p_func[:kwargs])))

    p = Dict(
        :category => :function,
        :name => p_func[:name],
        :args => args,
        :kwargs => kwargs,
        :short => true,
        :body => body,
        :relax_signature => is_promoted,
    )
    APIFunction(spec, with_func_ptr, p)
end

function contains_api_structs(def::Union{APIFunction,Constructor})
    any(x -> x ≠ promote_hl(x), def.p[:args])
end

is_promoted(ex) = ex == promote_hl(ex)

function promote_hl(def::APIFunction)
    APIFunction(def, def.with_func_ptr, promote_hl(def.p))
end

function promote_hl(def::Constructor)
    Constructor(def, promote_hl(def.p))
end

function promote_hl(arg::ExprLike)
    id, type = @match arg begin
        :($id::$t) => (id, t)
        _ => return arg
    end
    type = postwalk(type) do ex
        if ex isa Symbol && startswith(string(ex), '_')
            Symbol(string(ex)[2:end]) # remove underscore prefix
        else
            ex
        end
    end
    :($id::$type)
end

function promote_hl(p::Dict)
    args = promote_hl.(p[:args])
    modified_args = [arg for (arg, new_arg) in zip(p[:args], args) if arg ≠ new_arg]
    !isempty(modified_args) || error("Cannot define high-level function for $(reconstruct_call(p)): there are no arguments to adjust, so the low-level function is enough.")
    call_args = map(p[:args]) do arg
        id, type = @match arg begin
            :($id::$t) => (id, t)
            id => (id, nothing)
        end
        if arg in modified_args
            T = @match type begin
                :(AbstractArray{<:$t}) => :(Vector{$t})
                t => t
            end
            :(convert($T, $id))
        else
            id
        end
    end
    p = Dict(
        :category => :function,
        :name => p[:name],
        :args => args,
        :kwargs => p[:kwargs],
        :body => reconstruct_call(Dict(:name => p[:name], :args => call_args, :kwargs => name.(p[:kwargs]))),
        :short => true,
        :relax_signature => true,
    )
end

ext: jl


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["DllMain"], :dllmain),
    LibraryProduct(prefix, String["EnvConst"], :envconst),
    LibraryProduct(prefix, String["AstroFunc"], :astrofunc),
    LibraryProduct(prefix, String["TimeFunc"], :timefunc),
    LibraryProduct(prefix, String["Tle"], :tle),
    LibraryProduct(prefix, String["Sgp4Prop"], :sgp4prop),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/benelsen/SAASGP4Builder.jl/releases/download/v7.8-0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    BinaryProvider.Linux(:x86_64, :glibc) => ("$bin_prefix/SAASGP4.x86_64-linux-gnu.tar.gz", "15369188833880f827e1d536ee30d6aabb8bd6a359699f544aa941ac83fdd1bb"),
)

# First, check to see if we're all satisfied
if any(!satisfied(p; verbose=verbose) for p in products)
    if haskey(download_info, platform_key())
        # Download and install binaries
        url, tarball_hash = download_info[platform_key()]
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    else
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something more even more ambitious here.
        error("Your platform $(Sys.MACHINE) is not supported by this package!")
    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:
a, b, c = parse.(Int, split(readline()))
if a < b < c
    println("Yes")
else
    println("No")
end


text:
# Firms' parameters
const A = 1
const N = 1
const α = 0.33
const β = 0.96
const δ = 0.05

"""
Compute wage rate given an interest rate, r
"""
function r_to_w(r::Real)
    return A * (1 - α) * (A * α / (r + δ)) ^ (α / (1 - α))
end

"""
Inverse demand curve for capital. The interest rate
associated with a given demand for capital K.
"""
function rd(K::Real)
    return A * α * (N / K) ^ (1 - α) - δ
end

"""
Map prices to the induced level of capital stock.

##### Arguments
- `am::Household` : Household instance for problem we want to solve
- `r::Real` : interest rate

##### Returns
- The implied level of aggregate capital
"""
function prices_to_capital_stock(am::Household, r::Real)

    # Set up problem
    w = r_to_w(r)
    setup_R!(am, r, w)
    aiyagari_ddp = DiscreteDP(am.R, am.Q, am.β)

    # Compute the optimal policy
    results = solve(aiyagari_ddp, PFI)

    # Compute the stationary distribution
    stationary_probs = stationary_distributions(results.mc)[:, 1][1]

    # Return K
    return dot(am.s_vals[:, 1], stationary_probs)
end

# Create an instance of Household
am = Household(β=β, a_max=20.0)

# Create a grid of r values at which to compute demand and supply of capital
num_points = 20
r_vals = linspace(0.005, 0.04, num_points)

# Compute supply of capital
k_vals = prices_to_capital_stock.(am, r_vals)

# Plot against demand for capital by firms
demand = rd.(k_vals)
labels =  ["demand for capital" "supply of capital"]
plot(k_vals, [demand r_vals], label=labels, lw=2, alpha=0.6)
plot!(xlabel="capital", ylabel="interest rate", xlim=(2, 14), ylim=(0.0, 0.1))


text:
using Test, GpABC, OrdinaryDiffEq, Distances, Distributions

@testset "SMC ABC Test" begin
    #
    # ABC settings
    #
    n_particles = 1000
    threshold_schedule = [3.0, 2.0, 1.0]
    priors = [Uniform(0., 5.), Uniform(0., 5.)]
    progress_every = 1000

    #
    # Emulation settings
    #
    n_design_points = 100

    #
    # True parameters
    #
    true_params =  [2.0, 1.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]

    #
    # ODE solver settings
    #
    Tspan = (0.0, 10.0)
    x0 = [3.0, 2.0, 1.0]
    solver = RK4()
    saveat = 0.1

    #
    # Returns the solution to the toy model as solved by OrdinaryDiffEq
    #
    GeneReg = function(params::AbstractArray{Float64,1},
        Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1},
        solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

      if size(params,1) != 10
        throw(ArgumentError("GeneReg needs 10 parameters, $(size(params,1)) were provided"))
      end

      function ODE_3GeneReg(dx, x, pars, t)
        dx[1] = pars[1]/(1+pars[7]*x[3]) - pars[4]*x[1]
        dx[2] = pars[2]*pars[8]*x[1]./(1+pars[8]*x[1]) - pars[5]*x[2]
        dx[3] = pars[3]*pars[9]*x[1]*pars[10]*x[2]./(1+pars[9]*x[1])./(1+pars[10]*x[2]) - pars[6]*x[3]
      end

      prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params)
      Obs = solve(prob, solver, saveat=saveat)

      return Array{Float64, 2}(Obs)
    end

    reference_data = GeneReg(true_params, Tspan, x0, solver, saveat)
    simulator_function(var_params) = GeneReg(vcat(var_params, true_params[length(var_params)+1:end]), Tspan, x0, solver, saveat)

    #
    # Test using keep all as summary statistic
    #
    sim_out = SimulatedABCSMC(reference_data, simulator_function, priors, threshold_schedule, n_particles)
    @test size(sim_out.population, 1) > 0
    emu_out = EmulatedABCSMC(reference_data, simulator_function, priors, threshold_schedule, n_particles, n_design_points)
    @test size(emu_out.population, 1) > 0

    #
    # Test using built-in summary statistics
    #
    sum_stats = ["mean", "variance", "max", "min", "range", "median", "q1", "q3", "iqr"]
    sim_out = SimulatedABCSMC(reference_data, simulator_function, priors, threshold_schedule, n_particles;
        summary_statistic=sum_stats)
    @test size(sim_out.population, 1) > 0
    emu_out = EmulatedABCSMC(reference_data, simulator_function, priors, threshold_schedule, n_particles, n_design_points;
        summary_statistic=sum_stats)
    @test size(emu_out.population, 1) > 0

    #
    # Test using custom summary statistic
    #
    function sum_stat(data::AbstractArray{Float64,2})
        return std(data, dims=2)[:]
    end
    sim_out = SimulatedABCSMC(reference_data, simulator_function, priors, 3.0 * threshold_schedule, n_particles;
        summary_statistic=sum_stat, write_progress=false)
    @test size(sim_out.population, 1) > 0
    emu_out = EmulatedABCSMC(reference_data, simulator_function, priors, 3.0 * threshold_schedule, n_particles, n_design_points;
        summary_statistic=sum_stat, write_progress=false)
    @test size(emu_out.population, 1) > 0

    #
    # Custom emulator (re-)training and particle selection
    #
    emu_out = EmulatedABCSMC(reference_data, simulator_function, priors, 3.0 * threshold_schedule, n_particles, n_design_points;
        summary_statistic=sum_stat,
        emulator_training=DefaultEmulatorTraining(SquaredExponentialIsoKernel()),
        emulator_retraining=IncrementalRetraining(10, 100),
        emulated_particle_selection=MeanVarEmulatedParticleSelection(2.0))
    @test size(emu_out.population, 1) > 0

    for selection in [
        MeanEmulatedParticleSelection(),
        MeanVarEmulatedParticleSelection(),
        MeanVarEmulatedParticleSelection(2.0),
        PosteriorSampledEmulatedParticleSelection(),
        PosteriorSampledEmulatedParticleSelection(true)]

        emu_out = EmulatedABCSMC(reference_data, simulator_function, priors, 3.0 * threshold_schedule, n_particles, n_design_points;
            summary_statistic=sum_stat,
            emulator_training=DefaultEmulatorTraining(SquaredExponentialIsoKernel()),
            emulator_retraining=IncrementalRetraining(10, 100),
            emulated_particle_selection=selection)
        @test size(emu_out.population, 1) > 0

    end

end


text:
module ConstrainedOptim

using Optim
import Optim: optimize, Optimizer, @add_generic_fields()

export EqualityConstraint, BoxConstraint, BallConstraint, ProjectedGradientDescent, optimize
include("types.jl")
include("objective_helpers.jl")
include("augmented_lagrangian.jl")
include("optimize.jl")
end # module

max_stars_repo_name: cortner/AugmentedLagrangian

max_issues_repo_name: cortner/AugmentedLagrangian

max_forks_repo_name: cortner/AugmentedLagrangian


text:
using BenchmarkTools
using EnhancedGJK
import GeometryTypes
using CoordinateTransformations: IdentityTransformation, Translation, AffineMap
using EnhancedGJK: projection_weights, transform_simplex, penetration_distance
using StaticArrays: SVector, SMatrix
using Random

Random.seed!(1)

const gt = GeometryTypes

function comparison_benchmark(c1, c2)
    group = BenchmarkGroup()
    cache = CollisionCache(c1, c2)
    group["geometrytypes"] = @benchmarkable gt.gjk($c1, $c2)
    group["enhanced"] = @benchmarkable gjk!($cache, IdentityTransformation(), IdentityTransformation())
    group["enhanced no cache"] = @benchmarkable gjk($c1, $c2)
    return group
end

let
    suite = BenchmarkGroup()

    suite["3d_simplex_1"] = @benchmarkable projection_weights(simplex) setup=(
        simplex = SVector{4}(SVector{3, Float64}[[1., -1, 0], [1., 1, 0], [0.4, 0, -0.1], [0.5, 0, 1]])
    )

    suite["3d_simplex_2"] = @benchmarkable projection_weights(simplex) setup=(
        simplex = SVector{4}(SVector{3, Float64}[[1., -1, 0], [1., 1, 0], [1.4, 0, -0.1], [1.5, 0, 1]])
    )

    suite["2d_simplex_1"] = @benchmarkable projection_weights(simplex) setup=(
        simplex = SVector{3}(SVector{2, Float64}[[1., -1], [1., 1], [0.5, 0]])
    )

    suite["2d_simplex_2"] = @benchmarkable projection_weights(simplex) setup=(
        simplex = SVector{3}(SVector{2, Float64}[[1., -1], [1., 1], [1.5, 0]])
    )

    comparisons = suite["GeometryTypes comparisons"] = BenchmarkGroup()
    comparisons["2d_simplex_simplex"] = let
        c1 = gt.Simplex(gt.Vec(0.,0), gt.Vec(0.,1), gt.Vec(1.,0),gt.Vec(1.,1))
        c2 = gt.Simplex(gt.Vec(4.,0.5))
        comparison_benchmark(c1, c2)
    end
    comparisons["3d_point_box"] = let
        c1 = gt.HyperRectangle(rand(gt.Vec{3}), rand(gt.Vec{3}))
        c2 = rand(gt.Vec{3})
        comparison_benchmark(c1, c2)
    end

    suite["transform_simplex"] = @benchmarkable transform_simplex(cache, poseA, poseB) setup = begin
        c1 = gt.HyperRectangle(rand(gt.Vec{3}), rand(gt.Vec{3}))
        c2 = gt.Simplex(rand(gt.Vec{3}))
        cache = CollisionCache(c1, c2)
        poseA = AffineMap(rand(SMatrix{3, 3}), rand(SVector{3}))
        poseB = AffineMap(rand(SMatrix{3, 3}), rand(SVector{3}))
        gjk!(cache, poseA, poseB)
    end

    suite["penetration_distance"] = @benchmarkable penetration_distance(simplex) setup = begin
        c1 = gt.HyperRectangle(rand(gt.Vec{3}), rand(gt.Vec{3}))
        c2 = rand(gt.Vec{3})
        result = gjk(c1, c2)
        simplex = result.simplex
    end

    tune!(suite)
    results = run(suite)
    show(IOContext(stdout, :compact => false), results)
end


text:
import Pkg
using PyCall, Conda, JSON, MacroTools
import MacroTools: splitdef, combinedef

## PYTHON MODULES
const Vosk   = PyNULL()
const Pynput = PyNULL()

function __init__()
    ENV["PYTHON"] = ""                                              # FORCE PyCall to use Conda.jl
    if !any(startswith.(PyCall.python, DEPOT_PATH))                 # Rebuild of PyCall if it has not been built with Conda.jl
        @info "Rebuilding PyCall for using Julia Conda.jl..."
        Pkg.build("PyCall")
        @info "...rebuild completed. Restart Julia and JustSayIt."
        exit()
    end
    copy!(Vosk,   pyimport_pip("vosk"))
    copy!(Pynput, pyimport_pip("pynput"))
end

## CONSTANTS

const SAMPLERATE = 44100                #[Hz]
const AUDIO_READ_MAX = 512              #[bytes]
const AUDIO_ALLOC_GRANULARITY = 1024^2  #[bytes]
const AUDIO_HISTORY_MIN = 1024^2        #[bytes]
const AUDIO_ELTYPE = Int16
const AUDIO_IN_CHANNELS = 1
Base.ENDIAN_BOM == 0x04030201 || error("only little-endian is supported")
const RECORDER_BACKEND = "arecord"
const RECORDER_ARGS = ["--rate=$SAMPLERATE", "--channels=$AUDIO_IN_CHANNELS", "--format=S16_LE", "--quiet"] # (S16_LE: signed 16 bit little endian)
const COMMAND_NAME_EXIT  = "terminus"
const COMMAND_NAME_SLEEP = "sleep"
const COMMAND_NAME_AWAKE = "awake"
const COMMAND_ABORT = "abortus"
const VARARG_END = "terminus"
const VALID_VOICEARGS_KWARGS = Dict(:model=>String, :valid_input=>AbstractArray{String}, :valid_input_auto=>Bool, :interpret_function=>Function, :use_max_accuracy=>Bool, :vararg_end=>String, :vararg_max=>Integer, :vararg_timeout=>AbstractFloat)
const DEFAULT_MODEL_NAME = "default"
const DEFAULT_RECORDER_ID = "default"
const TYPE_MODEL_NAME = "type"
const COMMAND_RECOGNIZER_ID = "" # NOTE: This is a safe ID as it cannot be taken by any model (raises error).
const NOISES_ENGLISH = ["huh"]
const DIGITS_ENGLISH = Dict("zero"=>"0", "one"=>"1", "two"=>"2", "three"=>"3", "four"=>"4", "five"=>"5", "six"=>"6", "seven"=>"7", "eight"=>"8", "nine"=>"9", "dot"=>".")
const UNKNOWN_TOKEN = "[unk]"

## FUNCTIONS

function pyimport_pip(modulename::AbstractString)
    try
        pyimport(modulename)
    catch e
        if isa(e, PyCall.PyError)
            Conda.pip_interop(true)
            Conda.pip("install", modulename)
            pyimport(modulename)
        else
            rethrow(e)
        end
    end
end

let
    global tic, toc
    t0::Union{Float64,Nothing} = nothing

    tic()::Float64            = ( t0 = time() )
    toc()::Float64            = ( time() - t0 )
    toc(t1::Float64)::Float64 = ( time() - t1 )
end

function pretty_dict_string(dict::Dict{String,<:Any})
    key_length_max = maximum(length.(keys(dict)))
    return join([map(sort([keys(dict)...])) do x
                    join((x, dict[x]), " "^(key_length_max+1-length(x)) * "=> ")
                end...
                ], "\n")
end


text:
using TrajectoryOptimization
using BenchmarkTools
using TrajOptPlots
using Plots
using LinearAlgebra
using MeshCat

vis = Visualizer()
open(vis)

solver = DIRCOLSolver(Problems.Quadrotor()...,integration=HermiteSimpson)
solve!(solver)
cost(solver)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Double Integrator
solver = ALTROSolver(Problems.DoubleIntegrator()...)
benchmark_solve!(solver)
iterations(solver) # 8
plot(solver)

# Pendulum
solver = ALTROSolver(Problems.Pendulum()...)
benchmark_solve!(solver)
iterations(solver) # 19
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Cartpole
solver = ALTROSolver(Problems.Cartpole()...)
benchmark_solve!(solver)
iterations(solver) # 40
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Acrobot
solver = ALTROSolver(Problems.Acrobot()...)
benchmark_solve!(solver)
iterations(solver) # 50
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Parallel Park
solver = ALTROSolver(Problems.DubinsCar(:parallel_park)...)
benchmark_solve!(solver)
iterations(solver)  # 13
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Three Obstacles
solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)...)
benchmark_solve!(solver)
iterations(solver) # 20
delete!(vis)
add_cylinders!(vis, solver, robot_radius=get_model(solver).radius)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Escape
solver = ALTROSolver(Problems.DubinsCar(:escape)..., infeasible=true, R_inf=0.1)
benchmark_solve!(solver, samples=1, evals=1)
iterations(solver) # 13
delete!(vis)
set_mesh!(vis, get_model(solver))
add_cylinders!(vis, solver, robot_radius=get_model(solver).model.radius, height=0.2)
visualize!(vis, solver)
A = rand(10,10)
typeof(view(A,1:0,1:0))

# Zig-zag
solver = ALTROSolver(Problems.Quadrotor(:zigzag)...)
benchmark_solve!(solver)
iterations(solver) # 15
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)

# Barrell Roll
solver = ALTROSolver(Problems.YakProblems()...)
benchmark_solve!(solver)
iterations(solver) # 17
max_violation(solver)
delete!(vis)
set_mesh!(vis, get_model(solver))
visualize!(vis, solver)


text:
module SkyCoords
using Compat

export AbstractSkyCoords,
       ICRSCoords,
       GalCoords,
       FK5Coords

import Base: convert, *, transpose

# -----------------------------------------------------------------------------
# Types

abstract AbstractSkyCoords

immutable ICRSCoords <: AbstractSkyCoords
    ra::Float64
    dec::Float64
    ICRSCoords(ra::Real, dec::Real) = @compat new(mod(Float64(ra), 2pi),
                                                  Float64(dec))
end

immutable GalCoords <: AbstractSkyCoords
    l::Float64
    b::Float64
    GalCoords(l::Real, b::Real) = @compat new(mod(Float64(l), 2pi),
                                              Float64(b))
end

# FK5 is parameterized by equinox (e)
immutable FK5Coords{e} <: AbstractSkyCoords
    ra::Float64
    dec::Float64
    FK5Coords(ra::Real, dec::Real) = @compat new(mod(Float64(ra), 2pi),
                                                 Float64(dec))
end

# -----------------------------------------------------------------------------
# Helper functions: Immutable array operations

# We use immutable array operations to avoid allocating memory for
# small arrays. Eventually, base Julia should support immutable arrays, at
# which point we should switch to using that functionality in base.

immutable Matrix33
    a11::Float64
    a12::Float64
    a13::Float64
    a21::Float64
    a22::Float64
    a23::Float64
    a31::Float64
    a32::Float64
    a33::Float64
end

immutable Vector3
    x::Float64
    y::Float64
    z::Float64
end

function *(x::Matrix33, y::Matrix33)
    Matrix33(x.a11 * y.a11 + x.a12 * y.a21 + x.a13 * y.a31,
             x.a11 * y.a12 + x.a12 * y.a22 + x.a13 * y.a32,
             x.a11 * y.a13 + x.a12 * y.a23 + x.a13 * y.a33,
             x.a21 * y.a11 + x.a22 * y.a21 + x.a23 * y.a31,
             x.a21 * y.a12 + x.a22 * y.a22 + x.a23 * y.a32,
             x.a21 * y.a13 + x.a22 * y.a23 + x.a23 * y.a33,
             x.a31 * y.a11 + x.a32 * y.a21 + x.a33 * y.a31,
             x.a31 * y.a12 + x.a32 * y.a22 + x.a33 * y.a32,
             x.a31 * y.a13 + x.a32 * y.a23 + x.a33 * y.a33)
end

transpose(m::Matrix33) = Matrix33(m.a11, m.a21, m.a31,
                                  m.a12, m.a22, m.a32,
                                  m.a13, m.a23, m.a33)

function *(m::Matrix33, v::Vector3)
    Vector3(m.a11 * v.x + m.a12 * v.y + m.a13 * v.z,
            m.a21 * v.x + m.a22 * v.y + m.a23 * v.z,
            m.a31 * v.x + m.a32 * v.y + m.a33 * v.z)
end

# -----------------------------------------------------------------------------
# Helper functions: Create rotation matrix about a given axis (x, y, z)

function xrotmat(angle)
    s = sin(angle)
    c = cos(angle)
    Matrix33(1., 0., 0.,
             0.,  c,  s,
             0., -s,  c)
end

function yrotmat(angle)
    s = sin(angle)
    c = cos(angle)
    Matrix33(c,  0., -s,
             0., 1., 0.,
             s,  0.,  c)
end

function zrotmat(angle)
    s = sin(angle)
    c = cos(angle)
    Matrix33(c,   s,  0.,
             -s,  c,  0.,
             0., 0.,  1.)
end

# (lon, lat) -> [x, y, z] unit vector
function coords2cart(lon, lat)
    coslat = cos(lat)
    Vector3(coslat*cos(lon), coslat*sin(lon), sin(lat))
end

# [x, y, z] unit vector -> (lon, lat)
cart2coords(v) = atan2(v.y, v.x), atan2(v.z, sqrt(v.x*v.x + v.y*v.y))

# -----------------------------------------------------------------------------
# Constant rotation matricies and precession matrix function

# ICRS --> FK5 at J2000 (See USNO Circular 179, section 3.5)
eta0 = deg2rad(-19.9 / 3600000.)
xi0 = deg2rad(9.1 / 3600000.)
da0 = deg2rad(-22.9 / 3600000.)
const ICRS_TO_FK5J2000 = xrotmat(-eta0) * yrotmat(xi0) * zrotmat(da0)
const FK5J2000_TO_ICRS = ICRS_TO_FK5J2000'

# FK5J2000 --> Gal
# Note that galactic pole and zeropoint of l are somewhat arbitrary
# and not officially defined (as far as I know). The values below are
# from astropy.coordinates, which includes the following comment:
# | "This gives better consistency with other codes than using the values
# |  from Reid & Brunthaler 2004 and the best self-consistency between FK5
# |  -> Galactic and FK5 -> FK4 -> Galactic. The lon0 angle was found by
# |  optimizing the self-consistency."
ngp_fk5j2000_ra = deg2rad(192.8594812065348)
ngp_fk5j2000_dec = deg2rad(27.12825118085622)
lon0_fk5j2000 = deg2rad(122.9319185680026)
const FK5J2000_TO_GAL = (zrotmat(pi - lon0_fk5j2000) *
                         yrotmat(pi/2. - ngp_fk5j2000_dec) *
                         zrotmat(ngp_fk5j2000_ra))
const GAL_TO_FK5J2000 = FK5J2000_TO_GAL'

# Gal --> ICRS: simply chain through FK5J2000
const GAL_TO_ICRS = FK5J2000_TO_ICRS * GAL_TO_FK5J2000
const ICRS_TO_GAL = GAL_TO_ICRS'

# FK5J2000 --> FK5{epoch}
# Computes the precession matrix from J2000 to the given Julian equinox.
# Expression from from Capitaine et al. 2003 as expressed in the USNO
# Circular 179.  This should match the IAU 2006 standard from SOFA.
const pzeta = [2.650545, 2306.083227, 0.2988499, 0.01801828, -0.000005971,
               -0.0000003173]
const pz = [-2.650545, 2306.077181, 1.0927348, 0.01826837, -0.000028596,
            -0.0000002904]
const ptheta = [0.0, 2004.191903, -0.4294934, -0.04182264, -0.000007089,
                -0.0000001274]
function precess_from_j2000(equinox)
    t = (equinox - 2000.0) / 100.0
    tn = 1.0
    zeta = pzeta[1]
    z = pz[1]
    theta = ptheta[1]
    for i = 2:6
        tn *= t
        zeta += pzeta[i] * tn
        z += pz[i] * tn
        theta += ptheta[i] * tn
    end
    zeta = deg2rad(zeta/3600.0)
    z = deg2rad(z/3600.0)
    theta = deg2rad(theta/3600.0)
    zrotmat(-z) * yrotmat(theta) * zrotmat(-zeta)
end

# -----------------------------------------------------------------------------
# Type-dependent methods

# Abstract away specific field names (ra, dec vs l, b)
coords2cart(c::ICRSCoords) = coords2cart(c.ra, c.dec)
coords2cart(c::GalCoords) = coords2cart(c.l, c.b)
coords2cart{e}(c::FK5Coords{e}) = coords2cart(c.ra, c.dec)

# Rotation matrix between coordinate systems: `rotmat(to, from)`
rotmat(::Type{GalCoords}, ::Type{ICRSCoords}) = ICRS_TO_GAL
rotmat(::Type{ICRSCoords}, ::Type{GalCoords}) = GAL_TO_ICRS
rotmat{e}(::Type{FK5Coords{e}}, ::Type{ICRSCoords}) =
    precess_from_j2000(e) * ICRS_TO_FK5J2000
rotmat{e}(::Type{FK5Coords{e}}, ::Type{GalCoords}) =
   precess_from_j2000(e) * GAL_TO_FK5J2000
rotmat{e}(::Type{ICRSCoords}, ::Type{FK5Coords{e}}) =
    FK5J2000_TO_ICRS * precess_from_j2000(e)'
rotmat{e}(::Type{GalCoords}, ::Type{FK5Coords{e}}) =
    FK5J2000_TO_GAL * precess_from_j2000(e)'
rotmat{e1, e2}(::Type{FK5Coords{e1}}, ::Type{FK5Coords{e2}}) =
    precess_from_j2000(e1) * precess_from_j2000(e2)'

# Scalar coordinate conversions
convert{T<:AbstractSkyCoords}(::Type{T}, c::T) = c
function convert{T<:AbstractSkyCoords, S<:AbstractSkyCoords}(::Type{T}, c::S)
    r = rotmat(T, S) * coords2cart(c)
    lon, lat = cart2coords(r)
    T(lon, lat)
end

# Vector coordinate conversions
# This is useful for FK5Coords so that we only compute the rotation
# matrix (a function of epoch) once and apply it to the entire
# vector. The compiler doesn't seem to figure out that rotmat(T, S) is
# a constant for FK5Coords{e} types.
convert{T<:AbstractSkyCoords,n}(::Type{Array{T,n}}, c::Array{T,n}) = c
function convert{T<:AbstractSkyCoords, n, S<:AbstractSkyCoords}(
    ::Type{Array{T,n}}, c::Array{S, n})
    m = rotmat(T, S)
    result = similar(c, T)
    for i in 1:length(c)
        r = m * coords2cart(c[i])
        lon, lat = cart2coords(r)
        result[i] = T(lon, lat)
    end
    result
end

end # module


text:
abstract AbstractMesh
abstract AbstractMesh1d <: AbstractMesh

typealias BCnames Vector{String}
"""
Stores an ordered 1D mesh
"""
type Mesh1d{T <: Real} <: AbstractMesh1d
    "Number of nodes"
    nb::Int
    "Number of elements"
    nel::Int
    "Coordinates of the nodes"
    x::Vector{T}
    "Index of neighbor elements"
    neigh::Array{Int,2}
    bcs::BCnames
end

function check_elem_order{T<:Real}(x::AbstractVector{T})
    n = length(x)

    for i = 2:n
        if x[i] <= x[i-1]
            return false
        end
    end
    return true
end

"""
Outer constructer for `Mesh1d` objects

It builds a mesh from the nodes

 * `x`: Vector with node positions
 * `bcs`: Vector with boundary condition tags

Alternatively, the mesh can be built from the end points
and number of elements:

 * `a`: Left endpoint of the 1d domain
 * `b`: Right endpoint of the 1d domain
 * `nel`: Number of elements that the domain will be divided into
 * `bcs`: Vector with boundary condition tags
"""
function Mesh1d{T<:Real}(x::AbstractVector{T}, bcs=[])
    nb::Int = length(x)
    if nb < 2
        throw(DimensionMismatch("There should be at least 2 nodes in the 1d mesh"))
    end

    nel = nb-1

    if !check_elem_order(x)
        throw(ArgumentError("Nodes should be sequential from smaller to larger values"))
    end

    y = [xx for xx in x]

    if length(bcs) > 2
        throw(DimensionMismatch("There should be at most one BC for each end of the domain"))
    end

    bcs2 = String[utf8(b) for b in bcs]

    neigh = zeros(Int, 2, nel)

    for e = 1:nel-1
        neigh[1,e] = e-1
        neigh[2,e] = e+1
    end
    neigh[1,nel] = nel-1
    neigh[2,nel] = 0

    return Mesh1d(nb, nel, y, neigh, bcs2)

end

function Mesh1d{T<:Real}(a::T, b::T, nel::Int, bcs=[])

    if a >= b
        throw(ArgumentError("Nodes should be sequential from smaller to larger values"))
    end

    nb = nel + 1

    x = linspace(a, b, nb)

    return Mesh1d(x, bcs)

end

"""
Stores information about an element.
"""
immutable mshElem1d{T<:Real}
    id::Int
    a::T
    b::T
    neigh::Tuple{Int,Int}
end

import Base.getindex
getindex(msh::Mesh1d, e::Int) = mshElem1d(e, msh.x[e], msh.x[e+1], (msh.neigh[1,e], msh.neigh[2,e]))


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a78a6056-49dc-4315-b12a-f290a658fc2b
begin
	using Pluto
	using PlutoUI

	using LinearAlgebra
	using Statistics
	using Plots
	plotly()

	# Some useable definitions for the presentation
	# Taken from https://andreaskroepelin.de/blog/plutoslides/
	struct Foldable{C}
    	title::String
    	content::C
	end

	function Base.show(io, mime::MIME"text/html", fld::Foldable)
    	write(io,"$(fld.title)
")
    	show(io, mime, fld.content)
    	write(io,"
")
	end

	struct TwoColumn{L, R}
    	left::L
    	right::R
	end

	function Base.show(io, mime::MIME"text/html", tc::TwoColumn)
    	write(io, """

""")
    	show(io, mime, tc.left)
    	write(io, """

""")
    	show(io, mime, tc.right)
    	write(io, """

""")
	end
	nothing
end

# ╔═╡ 11152d22-899f-4a01-bec9-8ffe887306ad
html"<button onclick=present()>Present</button>"

# ╔═╡ 7a5e88a0-cdc0-11eb-23d1-8b1a18ca5072
md""" # Julia

Outline
+ Julia - A Fresh Approach to Numerical Computation
+ Julia - A *Fast* Approach to Numerical Computation
+ Programming Paradigms and Multiple Dispatch
+ Hands on Matching Pursuit
"""

# ╔═╡ afdf3c0f-8f09-425a-9ac9-94d5944c6e91
md""" # Preface
![](https://imgs.xkcd.com/comics/real_programmers.png)

Everyone is entitled to their personal opinion, workflow and preferences. The focus should be on the scientific progress.
"""

# ╔═╡ b0425c95-3eee-4461-9a48-3727fc2f6832
md""" # Julia : Raison d'être and Goals [[5]](https://julialang.org/blog/2012/02/why-we-created-julia/)

*In short, because we are greedy.*

+ The speed of C
+ The dynamism of Ruby
+ The obvious syntax of Matlab
+ The generalizability of Python
+ The adhesivity of the Shell
+ ...

"""

# ╔═╡ 1488e403-7d51-4bc5-b914-19b98df9ed51
md""" # Julia - A fresh approach to numerical computation [[1]](https://arxiv.org/abs/1411.1607)

Draft in high level language $\mapsto$ Reimplement in low level Language

Additional effort:
+ Map datastructures correctly
+ Ensure composability of datastructures and functions or methods
+ Implement interfaces to Open Source Projects
+ Add functionality to these
+ ...

This is known as the **Two Language Problem** in computer science.
"""

# ╔═╡ b9fd4165-9e51-48b1-81f0-4c5812ad26b9
begin
	text_c = md""" _C_

	+ Staticly Typed
	+ Compiled
	+ Highly performant

	""";
	text_cpp = md""" _C++_

	+ Staticly Typed
	+ Object Oriented
	+ Compiled
	+ Highly performant

	""";

	text_python = md""" _Python_

	+ Object Oriented
	+ Dynamicly Interpreted
	+ Performant (?)

	"""

	text_julia = md""" _Julia_
	+ Dynamicly Typed
	+ Just-In-Time Compiled
	+ Highly performant
	"""

	# Create a dict
	comp_lang = Dict(
		"text_c" => text_c,
		"text_cpp" => text_cpp,
		"text_python" => text_python
		)
	nothing
end

# ╔═╡ d50f89d9-06fe-402f-8613-d71591bdf9d2
md""" Let's compare Julia with $(@bind comparision_language_1 html"<select><option value='text_c'>C</option><option 		value='text_cpp'>C++</option><option 	value='text_python'>Python</option</select>")
"""

# ╔═╡ 19815c4f-b3ec-4708-aab8-3219b4ac5a77
TwoColumn(comp_lang[comparision_language_1], text_julia)

# ╔═╡ b7d2a5f4-c6f9-43ef-9beb-38f833b5d13c
md""" # The need for speed
"""

# ╔═╡ 2790a41c-cb5b-4eed-bb6e-7687b7c01824
# Consider the function
f(x,y) = exp(-(x-y)^2)

# ╔═╡ 014a0e4c-abaf-445d-9de8-18cc1fcf7e1d
begin
	# For two arrays
	x = fill(1, 100000)
	x̂ = randn(100000)
end

# ╔═╡ e158522d-fe91-49fb-bf98-1fcf909f0223
md""" The error above is due to the missing vectorization of the function `f`.
We can add this quite easily.
"""

# ╔═╡ bd0fe52b-444e-442c-89c5-ba84d8bad709
begin

	md""" ## Benchmarks
	![Julia Microbenchmarks](https://julialang.org/assets/benchmarks/benchmarks.svg)
	Microbenchmarks of Julia vs. different Languages as currently available [here](https://julialang.org/benchmarks/)

	**BUT** This is an old plot. Be careful.

	"""

end

# ╔═╡ bdb0db4b-ac6d-4311-b506-55e82b7974e0
Foldable("Just In Time (JIT) Compilation",
	md"""
	Julia is just in time compiled, using [LLVM](https://llvm.org/) under the hood. Similar approaches exist for [Matlab](https://de.mathworks.com/products/matlab/matlab-execution-engine.html) and [Python](https://doc.pypy.org/en/latest/architecture.html#jit-compiler). However, these are optional oppose to Julia, where it is the standard.
	"""
	)

# ╔═╡ b9958f53-c487-4237-b84b-b4e35e86a656
md""" ### Some Performance Pitfalls """

# ╔═╡ 4ccb9374-f661-47b2-8f25-3818c78d5cf3
Foldable(
	"Benchmarking",
md"""
Most commonly, new users consider Julia to be slow due to a first function call. This is due to the compilation of each method happening at compile level.

```
f(x) = sum(x.^2)
x = randn(Float64, 10000);

# First time call, the compiler sorts of the method to use and compiles it
@time f(x)
  0.070206 seconds (291.00 k allocations: 17.744 MiB, 99.93% compilation time)
10105.128612847464 # <- Result

# Second time call
@time f(x)
  0.000035 seconds (3 allocations: 78.219 KiB) # <- See the drastic reduction here
10105.128612847464 # <- Result

```

Use [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl) and do as I tell you and not as I do.
"""
	)

# ╔═╡ 363eebd3-0c8d-49c3-9d8a-9e6aaf6fce8b
Foldable("Memory Layout",
	md"""
	Julia is **column major** oriented. Meaning internally, the entries of an array get stored

	```
	x[1,1] x[2,1] x[3,1] x[1,2] x[2,2] x[3,2] ...
	```

	for speed ups always try to keep this in mind.
	"""
	)

# ╔═╡ fd8dd3db-aeb6-4105-b4c2-c9eca716ebed
Foldable("Referecing",
	md"""
	By default, we reference rather than copy. Consider

	```
	x = randn(10);
	x[1] = 1; # Set the first element to 1
	v = x[1] # Set v <- x[1]
	x[1] = 2 # Set x[1] to 2
	v # Will return 2
	```
	"""
	)

# ╔═╡ f51666f4-2896-4ba6-8f92-63795fb43d3a
md""" # Batteries included

Within the [LinearAlgebra](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) package - included in the standard library - we have BLAS and LAPACK at our fingertips.
"""

# ╔═╡ c78dee8d-9ff2-441f-9f5b-fc60f413869a
# Simple dot
BLAS.dot(10, fill(1.0, 10), 1, fill(1.0, 20), 2)

# ╔═╡ 3d110ae2-c8d0-4355-bccc-7d7a89224453
# a*x + y
BLAS.axpy!(2, [1;2;3], [4;5;6])

# ╔═╡ 606e95e0-f8d2-47dd-964c-fc33b72a91d1
# Solve A X = B via LU(A) and overwrites B with the solution
begin
	A = randn(3,3)
	B = A * [1 0 0; 0 0 1; 0 1 0]
	LAPACK.gesv!(A, B)
end

# ╔═╡ ff06c5fe-500b-496e-8801-d399f1b2f394
md""" # Classical Dispatch
"""

# ╔═╡ 62d14388-0616-4ab3-af97-82e31d380157
TwoColumn(
	md""" **Functional**

	_Global_ Namespace

	$f \in \mathcal{F}$

	_Unique_ Functions

	$f(x_1, x_2, x_3, \dots)$

	Results in

	```
	z = 3+i4
	r = 0.1
	complex_real_add(z, r)
	```
	""",

	md""" **Object Oriented**

	_Local_ Namespace

	$f \in \mathcal{O}$

	_Unique_ Functions in Namespace

	$o.f(x_1, x_2, x_3, \dots)$

	Can result in
	```
	z = 3+i4
	r = 0.1
	z + r
	```
	"""
)

# ╔═╡ 148805bf-e3c1-497b-a776-fd997a6d7acd
md""" # Multiple Dispatch

Global namespace

$f \in \mathcal{F}$

With unique arguments

$f : \mathcal{X}_1 \times \mathcal{X}_2 \dots \mapsto \mathcal{Y}$

```
	z = 3+i4
	r = 0.1
	z + r
```
"""

# ╔═╡ 1d53fe22-000f-4ee4-9c65-c4f755dd3ce6
begin

	_xx = md"Consider the expressive power [[2]](https://www.youtube.com/watch?v=kc9HwsxE1OY)

| | Functional | OOP | Multiple Dispatch
:--- | --- | ----- | ---:
Dispatch arg. | 0 | 1   |   N
Expressive Power | const. | linear | exponential
"
	Foldable("Is this useful?", _xx)
end

# ╔═╡ d1f967e2-bea9-45f4-8760-eb22d1908bf8

	Foldable("But why exactly is this useful?", md" **EFFECTIVE REUSE AND EXTENABILITY OF CODE** [[3]](https://www.oxinabox.net/2020/02/09/whycompositionaljulia.html)

	```
using StatsBase  # exports predict
using Foo  # overloads `StatsBase.predict(::FooModel)
using Bar  # overloads `StatsBase.predict(::BarModel)
training_data, test_data = ...
mbar = BarModel(training_data)
mfoo = FooModel(training_data)
evaluate(predict(mbar), test_data)
evaluate(predict(mfoo), test_data)
```
	")

# ╔═╡ 66079375-24c7-4e27-ba41-08fd5634fcd7
md""" ## Dual Numbers

$z = a + b \epsilon$

with

$\epsilon^2 = 0$
"""

# ╔═╡ 7b751a69-72aa-4543-b89a-0075a48f3bf0
# Type definition
struct DualNumber{T} <: Number
	a::T
	b::T
end

# ╔═╡ 1d55314e-7cb4-4be7-8e44-94c4a95739dd
# Addition
Base.:+(x::DualNumber{T}, y::DualNumber{T}) where T = DualNumber(x.a + y.a, x.b + y.b)

# ╔═╡ 55dfc8e4-efc5-4eb4-8534-c482d43059ae
z = DualNumber(1.0, 0.2)

# ╔═╡ 481e1757-cbd9-4b52-b6d7-b828934e0a74
sum([z for i in 1:100])

# ╔═╡ 1fb93892-2ab8-4e49-9076-1517afb2299d
# Multiplication
Base.:*(x::DualNumber, y::DualNumber) = DualNumber(
	x.a*y.a, x.a*y.b+y.a*x.b
	)

# ╔═╡ 93bca010-ef25-4672-a601-b42fbae73bd3
z*z # Works

# ╔═╡ 987772c6-35e6-4340-b8ad-65eb6c5521dc
prod([(z^3 + z^2) for i in 1:10])

# ╔═╡ ebc190c1-aa57-4ab0-8f6f-958eaae4008b
begin
	# Add multiplication
	Base.:*(x::Number, y::DualNumber{T}) where T = DualNumber(convert(T, x)*y.a, convert(T,x)*y.b)
	Base.:*(y::DualNumber{T}, x::Number) where T = DualNumber(convert(T, x)*y.a, convert(T,x)*y.b)
	# Add subtraction
	Base.:-(x::Number, y::DualNumber{T}) where T = DualNumber(convert(T,x)-y.a, -y.b)
	Base.:-(y::DualNumber{T}, x::Number) where T = DualNumber(y.a-convert(T,x), y.b)
	Base.:-(x::DualNumber, y::DualNumber) = x+(-1*y)
end

# ╔═╡ 892e2207-b7e0-410c-a22e-93aa185ea858
f(x) = x^2 - 3*(x-2)^4

# ╔═╡ 7584eb8c-9e6a-43cc-859a-f2940adadad0
# Lets see what happens
f(1, 5)

# ╔═╡ d4ec78eb-10ee-4b0f-812f-0c79f5dfce91
# And under the hood
@code_lowered f(0.2, 3.0)

# ╔═╡ e3033c49-9891-47eb-9d78-e1bd629a944a
with_terminal() do
	@time f(x, x̂)
end

# ╔═╡ 70764379-8132-4ca3-aee6-87871c2132a2
# Can we get it down? Hide at first
function fastf(x::AbstractVector{X}, y::AbstractVector{Y}) where {X, Y}
	# We want a common type
	T = promote_type(X, Y)
	# Convert the input to a common type
	x = convert.(T, x)
	y = convert.(T, y)
	# Preallocation of the result
	res = Vector{T}(undef, length(x))
	# Vectorize
	# We assume length(x) == length(y) -> Errors are not catched!
	# Additionally we add the @simd macro to instruct the compiler
	# to parallelize
	@simd for i in 1:length(x)
		# If we want, we could also use the @fastmath macro here
		# which does not improve our performance ( in this case )
		res[i] = f(x[i], y[i])
	end
	# The return argument
	return res
end

# ╔═╡ 76a0d6b9-f6e6-4e68-9607-17be0d8b8c19
with_terminal() do
	@time fastf(x, x̂)
end

# ╔═╡ 0c909529-b296-4a61-a369-e0d129b37501
# Mutating function and even faster because we assume the same type
function fastf!(res::AbstractVector{T}, x::AbstractVector{T}, y::AbstractVector{T}) where {T}
	# Vectorize
	# We assume length(x) == length(y) -> Errors are not catched!
	# Additionally we add the @simd macro to instruct the compiler
	# to parallelize
	@inbounds @simd for i in 1:length(x)
		# If we want, we could also use the @fastmath macro here
		# which does not improve our performance ( in this case )
		res[i] = f(x[i], y[i])
	end
	# The return argument
	return
end

# ╔═╡ e70fd29d-6ee0-4011-9f45-df4d4e72d6d4
with_terminal() do
	# Benchmarking
	mybench(f::Function, args...) = @time f(args...)

	mybench(f, x, x̂)
	mybench(fastf, x, x̂)

	res = similar(x̂)
	mybench(fastf!, res, convert.(eltype(x̂), x), x̂)
end

# ╔═╡ 97306279-df1f-4bb2-9d51-8e6345808100
f(z)

# ╔═╡ 4c45e54f-7d0b-4f3d-acfc-b4cd90e7afb6
f(0.3)

# ╔═╡ 422e28e6-c51e-4144-b0a6-970cdee830ed
f(5)

# ╔═╡ 5b6189c3-b2be-4487-aac4-0d182b8bd5bb
f(0.3+im*3)

# ╔═╡ fb3490f1-0b3a-44b0-9bea-50cb5e530c15
md""" # Matchig Pursuit [[4]](https://en.wikipedia.org/wiki/Matching_pursuit)

$\min_\Xi \quad \lVert \Xi \rVert_0, \quad \text{s.t.} \quad Y = \Psi(X) \Xi$

![Example](https://upload.wikimedia.org/wikipedia/commons/3/38/Orthogonal_Matching_Pursuit.gif)

**Pseudo-Code**

+ Compare each element of the normalized dictionary $\Psi(X)$ to the signal $Y$ via the inner product
+ Use the largest resemblance as a coefficient $\xi = \Xi_{i,j}$
+ Subtract the recovered signal and repeat until converged

We will use a more sophisticated version [based on this paper](https://asp-eurasipjournals.springeropen.com/track/pdf/10.1186/s13634-020-00680-9.pdf) where our goal is to find the right support $\Lambda$ of the coefficient vector.

*Note*

`LinearAlgebra` is already imported
"""

# ╔═╡ 34d14afe-2194-43f0-a829-15276ecca883
# Define the algorithm
function gOMP(y0::AbstractVector, Ψ::AbstractMatrix, K::Int = 2, S::Int = 1; max_iter::Int = 100, ϵ::Real = eps())
	# Get the dimensions
	m = length(y0)
	m_psi, n = size(Ψ)
	# Assert the dimensionality
	@assert m == m_psi "Please provide consistent input sizes"
	# Assert the selector
	@assert S <=  min(K, m/K) "S <= min(K, m/K)"

	# Normalize
	ψ = deepcopy(Ψ)
	normalize!.(eachcol(ψ), 2)

	# Preconditioning
	P = ψ'pinv(ψ*ψ')
	# New matrix
	ψ = P*ψ
	y = P*y0
	# Iteration
	iters = 0

	# Support
	Λ = zeros(Bool, n)
	u = zeros(eltype(y), n)
	r = y
	amps = zeros(eltype(y), n)

	# Find the magnitudes
	for i in 1:max_iter
		# Compute the similarities via magnitude
		amps .= ψ'r
		# Get the largest entry
		idx = sortperm(abs.(amps), rev = true)[1:S]
		# Update the support
		Λ[idx] .= true
		# Update the coefficients
		u[Λ] .= (y' / ψ')[1, Λ]
		# Update r
		r .= y - ψ*u
		# Convergence
		if norm(r,2) < ϵ || sum(Λ) >= K
			# Just for debug
			#@debug "Early break after $i iterations with $(norm(r,2))"
			break
		end
	end
	# Last time to get the right coefficients
	u[Λ] .= (y0' / Ψ')[1, Λ]
	return u

end

# ╔═╡ 2ebc2b89-27f2-4abc-816e-7950bce2267f
# Generate some test data
begin
	# Independent variable
	t = 0.0:0.1:10.0;
	# Signal
	y = 3.0*sin.(t).*exp.(-t./50.0) + 4.0*cos.(t) - t;
	# Dictionary
	ψ = [
		t t.^2 t.^3 sin.(t) ones(eltype(y), length(y)) cos.(t) exp.(t) sin.(t).*exp.(-t./50.0)
	];
	nothing
end

# ╔═╡ f7451fdb-b5ac-4bf7-82f2-fed46ed98e99
md""" Lets explore the data by adding a slider for the sparsity and a noise level.

Sparsity : $(@bind k Slider(1:size(ψ, 2)))

Noise : $(@bind n Slider(0.0:0.005:0.1))

"""

# ╔═╡ dac7e4e0-2532-45ff-8b0f-1d8a92e90ccf

# ╔═╡ 1679b86e-13ef-48c4-9a05-613387c213f9
# What about parallism?

# Simple

function gOMP(Y::AbstractMatrix, Ψ::AbstractMatrix, args...; kwargs...)

	# we know that the coefficients can be derived independent
	A = zeros(size(Ψ, 2), size(Y, 1)) # Init
	#Threads.@threads for i in 1:size(Y, 1)
	for i in 1:size(Y, 1)
		A[:, i] .= gOMP(Y[i, :], Ψ, args...; kwargs...)
	end
	return A
end

# ╔═╡ fb971894-a833-4955-93b0-da68c4a71db0
begin
	ynoise = y .+ n*mean(y, dims = 2).*randn(size(y))
	a = gOMP(ynoise, ψ, k)
	plot(t, ynoise, label = "Original Signal with $(n*100) percent noise", title = "Fit")
	plot!(t, ψ*a, label = "Recovered with  $(sum(abs.(a) .> 1e-10)) / $k elements ")
	annotate!([(0.0, -14.0, Plots.text("Coefficients : $(round.(a, digits = 2))", 10, :black, :left))])

end

# ╔═╡ 6e8805fe-05e8-49b5-a231-1efb995da173
Y = vcat([(y+1e-2*randn(size(y)))' for i in 1:100]...)

# ╔═╡ 05cccc05-89a8-48c9-b1ac-a8ff6b5089a9
with_terminal() do
	@time gOMP(Y, ψ, 4)
end

# ╔═╡ 253e609d-2fd1-4ec0-af3f-cefa10abe168
sum(gOMP(Y, ψ, 4), dims = 2)/100

# ╔═╡ 1d3939e1-bcec-46fd-84e0-b84abfbb3dc9
md""" # Ecosystem & Package Development

[JuliaHub](https://juliahub.com/ui/Packages) provides us with a nice, searchable database for all registered packages.

As an example for Package Development, we can have a look at [DataDrivenDiffEq.jl](https://github.com/SciML/DataDrivenDiffEq.jl).

And : [JuliaCon 2021](https://juliacon.org/2021/) is around the corner!
"""

# ╔═╡ Cell order:
# ╟─a78a6056-49dc-4315-b12a-f290a658fc2b
# ╟─11152d22-899f-4a01-bec9-8ffe887306ad
# ╟─7a5e88a0-cdc0-11eb-23d1-8b1a18ca5072
# ╟─afdf3c0f-8f09-425a-9ac9-94d5944c6e91
# ╟─b0425c95-3eee-4461-9a48-3727fc2f6832
# ╟─1488e403-7d51-4bc5-b914-19b98df9ed51
# ╟─b9fd4165-9e51-48b1-81f0-4c5812ad26b9
# ╟─d50f89d9-06fe-402f-8613-d71591bdf9d2
# ╟─19815c4f-b3ec-4708-aab8-3219b4ac5a77
# ╟─b7d2a5f4-c6f9-43ef-9beb-38f833b5d13c
# ╠═2790a41c-cb5b-4eed-bb6e-7687b7c01824
# ╠═7584eb8c-9e6a-43cc-859a-f2940adadad0
# ╠═d4ec78eb-10ee-4b0f-812f-0c79f5dfce91
# ╠═014a0e4c-abaf-445d-9de8-18cc1fcf7e1d
# ╠═e3033c49-9891-47eb-9d78-e1bd629a944a
# ╠═76a0d6b9-f6e6-4e68-9607-17be0d8b8c19
# ╟─e158522d-fe91-49fb-bf98-1fcf909f0223
# ╠═70764379-8132-4ca3-aee6-87871c2132a2
# ╠═0c909529-b296-4a61-a369-e0d129b37501
# ╠═e70fd29d-6ee0-4011-9f45-df4d4e72d6d4
# ╟─bd0fe52b-444e-442c-89c5-ba84d8bad709
# ╟─bdb0db4b-ac6d-4311-b506-55e82b7974e0
# ╟─b9958f53-c487-4237-b84b-b4e35e86a656
# ╟─4ccb9374-f661-47b2-8f25-3818c78d5cf3
# ╟─363eebd3-0c8d-49c3-9d8a-9e6aaf6fce8b
# ╟─fd8dd3db-aeb6-4105-b4c2-c9eca716ebed
# ╟─f51666f4-2896-4ba6-8f92-63795fb43d3a
# ╠═c78dee8d-9ff2-441f-9f5b-fc60f413869a
# ╠═3d110ae2-c8d0-4355-bccc-7d7a89224453
# ╠═606e95e0-f8d2-47dd-964c-fc33b72a91d1
# ╟─ff06c5fe-500b-496e-8801-d399f1b2f394
# ╟─62d14388-0616-4ab3-af97-82e31d380157
# ╟─148805bf-e3c1-497b-a776-fd997a6d7acd
# ╟─1d53fe22-000f-4ee4-9c65-c4f755dd3ce6
# ╟─d1f967e2-bea9-45f4-8760-eb22d1908bf8
# ╟─66079375-24c7-4e27-ba41-08fd5634fcd7
# ╠═7b751a69-72aa-4543-b89a-0075a48f3bf0
# ╠═1d55314e-7cb4-4be7-8e44-94c4a95739dd
# ╠═55dfc8e4-efc5-4eb4-8534-c482d43059ae
# ╠═481e1757-cbd9-4b52-b6d7-b828934e0a74
# ╠═1fb93892-2ab8-4e49-9076-1517afb2299d
# ╠═93bca010-ef25-4672-a601-b42fbae73bd3
# ╠═987772c6-35e6-4340-b8ad-65eb6c5521dc
# ╠═ebc190c1-aa57-4ab0-8f6f-958eaae4008b
# ╠═892e2207-b7e0-410c-a22e-93aa185ea858
# ╠═97306279-df1f-4bb2-9d51-8e6345808100
# ╠═4c45e54f-7d0b-4f3d-acfc-b4cd90e7afb6
# ╠═422e28e6-c51e-4144-b0a6-970cdee830ed
# ╠═5b6189c3-b2be-4487-aac4-0d182b8bd5bb
# ╟─fb3490f1-0b3a-44b0-9bea-50cb5e530c15
# ╠═34d14afe-2194-43f0-a829-15276ecca883
# ╠═2ebc2b89-27f2-4abc-816e-7950bce2267f
# ╟─fb971894-a833-4955-93b0-da68c4a71db0
# ╟─f7451fdb-b5ac-4bf7-82f2-fed46ed98e99
# ╟─dac7e4e0-2532-45ff-8b0f-1d8a92e90ccf
# ╠═1679b86e-13ef-48c4-9a05-613387c213f9
# ╠═6e8805fe-05e8-49b5-a231-1efb995da173
# ╠═05cccc05-89a8-48c9-b1ac-a8ff6b5089a9
# ╠═253e609d-2fd1-4ec0-af3f-cefa10abe168
# ╟─1d3939e1-bcec-46fd-84e0-b84abfbb3dc9


text:
"""
    Exo3Real{T <: Real} = Parabolic{Exo2Real{T}}

A 3-exo-real number is a parabolic Cayley-Dickson construct with bi-exoplex numbers.
"""
const Exo3Real{T <: Real} = Parabolic{Exo2Real{T}}

function Exo3Real(z::Exo2Real{T}) where T <: Real
    Exo3Real{T}(z, zero(Exo2Real{T}))
end

function Exo3Real(z::Exo1Real{T}) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(z, zero(Exo1Real{T})),
        zero(Exo2Real{T})
    )
end

function Exo3Real(a::T, b::T, c::T, d::T, f::T, g::T, h::T, j::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, d)
        ),
        Exo2Real{T}(
            Exo1Real{T}(f, g),
            Exo1Real{T}(h, j)
        )
    )
end

function Exo3Real(a::T, b::T, c::T, d::T, f::T, g::T, h::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, d)
        ),
        Exo2Real{T}(
            Exo1Real{T}(f, g),
            Exo1Real{T}(h, zero(T))
        )
    )
end

function Exo3Real(a::T, b::T, c::T, d::T, f::T, g::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, d)
        ),
        Exo2Real{T}(
            Exo1Real{T}(f, g),
            zero(Exo1Real{T})
        )
    )
end

function Exo3Real(a::T, b::T, c::T, d::T, f::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, d)
        ),
        Exo2Real{T}(
            Exo1Real{T}(f, zero(T)),
            zero(Exo1Real{T})
        )
    )
end

function Exo3Real(a::T, b::T, c::T, d::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, d)
        ),
        zero(Exo2Real{T})
    )
end

function Exo3Real(a::T, b::T, c::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            Exo1Real{T}(c, zero(T))
        ),
        zero(Exo2Real{T})
    )
end

function Exo3Real(a::T, b::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, b),
            zero(Exo1Real{T})
        ),
        zero(Exo2Real{T})
    )
end

function Exo3Real(a::T) where T <: Real
    Exo3Real{T}(
        Exo2Real{T}(
            Exo1Real{T}(a, zero(T)),
            zero(Exo1Real{T})
        ),
        zero(Exo2Real{T})
    )
end

function Exo3Real(a::Real, b::Real)
    Exo3Real(promote(a, b)...)
end

function Exo3Real(a::Real, b::Real, c::Real)
    Exo3Real(promote(a, b, c)...)
end

function Exo3Real(a::Real, b::Real, c::Real, d::Real)
    Exo3Real(promote(a, b, c, d)...)
end

function Exo3Real(a::Real, b::Real, c::Real, d::Real, f::Real)
    Exo3Real(promote(a, b, c, d, f)...)
end

function Exo3Real(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real)
    Exo3Real(promote(a, b, c, d, f, g)...)
end

function Exo3Real(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real)
    Exo3Real(promote(a, b, c, d, f, g, h)...)
end

function Exo3Real(a::Real, b::Real, c::Real, d::Real, f::Real, g::Real, h::Real, j::Real)
    Exo3Real(promote(a, b, c, d, f, g, h, j)...)
end

function selfstar(z::Exo3Real)
    (z + star(z)) / 2
end

function anti_selfstar(z::Exo3Real)
    (z - star(z)) / 2
end

"""
    unreal(z::Exo3Real)

The unreal part of a 3-exo-real number is a 7-vector.
"""
function unreal(z::Exo3Real)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::Exo3Real)
    print(io, "[1: ")
    print(io, z.l.l.l)
    print(io, ", a: ")
    print(io, z.l.l.r)
    print(io, ", b: ")
    print(io, z.l.r.l)
    print(io, ", ab: ")
    print(io, z.l.r.r)
    print(io, ", c: ")
    print(io, z.r.l.l)
    print(io, ", ac: ")
    print(io, z.r.l.r)
    print(io, ", bc: ")
    print(io, z.r.r.l)
    print(io, ", (ab)c: ")
    print(io, z.r.r.r)
    print(io, "]")
end

function random(::Type{Exo3Real{T}}) where T <: Real
    Exo3Real{T}(random(Exo2Real{T}), random(Exo2Real{T}))
end


text:
function initSignals(window::Gtk.GtkWindow, canvas::Gtk.GtkCanvas)
    @guarded draw(canvas) do widget
        if loadedState.map !== nothing && isa(loadedState.map, Map)
            startTime = time()

            lastDrawDeltaTime = startTime - lastCanvasDraw
            global lastCanvasDraw = time()

            debug.log("Canvas was last drawn: $lastDrawDeltaTime seconds ago", "DRAWING_CANVAS_DRAW_DURATION")

            drawMap(canvas, camera, loadedState.map)

            stopTime = time()
            deltaTime = stopTime - startTime
            debug.log("Drawing canvas took $deltaTime seconds, approx $(1 / deltaTime) fps", "DRAWING_CANVAS_DRAW_DURATION")

            persistence["camera_position_x"] = camera.x
            persistence["camera_position_y"] = camera.y
            persistence["camera_scale"] = camera.scale
        end
    end

    add_events(canvas,
        GConstants.GdkEventMask.SCROLL |
        GConstants.GdkEventMask.BUTTON_PRESS |
        GConstants.GdkEventType.DOUBLE_BUTTON_PRESS |
        GConstants.GdkEventMask.BUTTON_RELEASE |
        GConstants.GdkEventMask.BUTTON1_MOTION |
        GConstants.GdkEventMask.BUTTON3_MOTION |
        GConstants.GdkEventMask.LEAVE_NOTIFY
    )

    signal_connect(ExitWindow.exitAhorn, window, "delete_event")
    signal_connect(
        function(window::GtkWindow, event::Gtk.GdkEventAny)
            persistence["start_maximized"] = get_gtk_property(window, :is_maximized, Bool)

            return false
        end,
        window,
        "window-state-event"
    )

    # signal_connect(resize_event, window, "resize")
    signal_connect(key_press_event, window, "key-press-event")
    signal_connect(key_release_event, window, "key-release-event")

    signal_connect(focus_out_event, window, "focus-out-event")

    signal_connect(scroll_event, canvas, "scroll-event")
    signal_connect(motion_notify_event, canvas, "motion-notify-event")
    signal_connect(button_press_event, canvas, "button-press-event")
    signal_connect(button_release_event, canvas, "button-release-event")
    signal_connect(leave_notify_event, canvas, "leave-notify-event")
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
const HILLCOEFFICIENT = Param(0.1; bounds=(0.0, 100))
const MEDIAN_LOGLOGISTIC = Param(0.1; bounds=(0.0, 100))
const THRESHOLD_EXPOSURE = Param(0.0; bounds=(0.0, 1e9))

"""
    Mortality <: CellRule

Abstract super type for rules of survival effect
For best performance these should be chained with other
CellRule or following an NeighborhoodRule.
"""
abstract type Mortality{R,W} <: CellRule{R,W} end

"""
    LoglogisticMortality <: Mortality

    LoglogisticMortality(; median, hillcoefficient, timestep, [nsteps_type])
    LoglogisticMortality{R}(; median, hillcoefficient, timestep, [nsteps_type])
    LoglogisticMortality{R,W}(; median, hillcoefficient, timestep, [nsteps_type])

Loglogistic mortality based on median, ``α`` and hill coefficient ``β``.

Cumulative function of loglogistic:
```math
F(x; α, β) =  x^β/(α^β + x^β)
```
where ``α>0`` is the scale and ``β>0`` is the shape.

# Keywords

- `median`: Median of the loglogistic function
- `hillcoefficient`: Hill's coefficient, a measure of ultrasensitivity (i.e. how steep is the response curve).
    May be a `Number`, an `Aux` array or another `Grid`.
- `timestep`: Time step for the mortality rate, in a type compatible with the simulation `tspan`.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
`R` is a 2 Grids `NamedTuple` like `Tuple{:population,:exposure}` and `W` return only the first grid `:population`.
"""
struct LoglogisticMortality{R,W,MD,HC,TS,S} <: Mortality{R,W}
    median::MD
    hillcoefficient::HC
    timestep::TS
    nsteps::S
end
function LoglogisticMortality{R,W}(;
    median=MEDIAN_LOGLOGISTIC,
    hillcoefficient=HILLCOEFFICIENT,
    timestep,
    nsteps_type=Float64,
) where {R,W}
    LoglogisticMortality{R,W}(median, hillcoefficient, timestep, zero(nsteps_type))
end

modifyrule(rule::LoglogisticMortality, data) = precalc_nsteps(rule, data)

@inline function applyrule(data, rule::LoglogisticMortality, (N, X), I)
    N > zero(N) || return zero(N)
    β = get(data, rule.hillcoefficient, I...)
    α = get(data, rule.median, I...)
    # NOTE: 1 - x^β / (α^β + x^β) = α^β / (α^β + x^β)
    return N * α^β / (α^β + X^β)
end

"""
    ExponentialMortality <: Mortality

    ExponentialMortality(; rate, threshold, timestep, [nsteps_type])
    ExponentialMortality{R}(; rate, threshold, timestep, [nsteps_type])
    ExponentialMortality{R,W}(; rate, threshold, timestep, [nsteps_type])

Exponential mortality based on exposure threshold and mortality rate parameter, using exact solution.

Exponential mortality based on exposure grid ``X``, an exposure threshold parameter ``z``
    and a mortality rate ``r`` using exact solution between time ``t`` and ```t+1``:

```math
N_{t+1} = N_{t}e^{-r t(X-z)}
```

# Keywords

- `rate`: Mortality rate.
- `threshold`: Exposure threshold under which there is no effect.
- `timestep`: Time step for the growth rate, in a type compatible with the simulation `tspan`.
- `nsteps_type`: Specify the floating point type to use when `nsteps` is generated from the
    timestep, if it is required for type-stability or performance. The default is `Float64`.

`rate` and `threshold` can be a `Number`, an `Aux` array or another Grid`.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
`R` is a 2 Grids `NamedTuple` like `Tuple{:population,:exposure}` and `W` return only the first grid `:population`.
"""
struct ExponentialMortality{R,W,MR,ZX,TS,S} <: GrowthRule{R,W}
    rate::MR
    threshold::ZX
    timestep::TS
    nsteps::S
end
function ExponentialMortality{R,W}(;
    rate=INTRINSICRATE_PARAM,
    threshold=THRESHOLD_EXPOSURE,
    timestep,
    nsteps_type=Float64,
) where {R,W}
    ExponentialMortality{R,W}(rate, threshold, timestep, zero(nsteps_type))
end

modifyrule(rule::ExponentialMortality, data) = precalc_nsteps(rule, data)

@inline function applyrule(data, rule::ExponentialMortality, (N,X), I)
    N > zero(N) || return zero(N)
    rt = get(data, rule.rate, I...) * rule.nsteps
    z = get(data, rule.threshold, I...)
    return N * exp(- rt * max(0, X-z))
end


text:
# Translating DFM's python version:
include("terms.jl")
using LinearAlgebra

mutable struct Celerite
    kernel::Term
    computed::Bool
    D::Vector{Real}
    W::Array{Real}
    up::Array{Real}
    phi::Array{Real}
    x::Vector{Real}
    logdet::Real
    n::Int
    J::Int

    Q::Array{Real, 2}

#    Celerite(kernel) = new(kernel, false, [], [], [], [], [])
#    Celerite(kernel) = new(kernel, false, zeros(Float64,0),zeros(Float64,0),zeros(Float64,0,0), zeros(Float64,0,0), zeros(Float64,0,0), zeros(Float64,0,0), zeros(Float64,0),0.0,0,0)
    Celerite(kernel, Q=ones(1, 1)) = new(kernel, false, zeros(Float64,0),zeros(Float64,0,0),zeros(Float64,0,0),zeros(Float64,0,0),zeros(Float64,0),0.0,0,0,Q)
end

function cholesky_ldlt!(a_real::Vector{<:Real}, c_real::Vector{<:Real},
                       a_comp::Vector{<:Real}, b_comp::Vector{<:Real},
                       c_comp::Vector{<:Real}, d_comp::Vector{<:Real},
                       t::Vector{<:Real}, var::Vector{<:Real}, X::Array{<:Real,2},
                       phi::Array{<:Real,2}, u::Array{<:Real,2}, D::Vector{<:Real})
#
# Fast LDLT Cholesky solver based on low-rank decomposition due to Sivaram, plus
# real implementation of celerite term.
#
# Compute the dimensions of the problem:
    N = length(t)
# Number of real components:
    J_real = length(a_real)
# Number of complex components:
    J_comp = length(a_comp)
# Rank of semi-separable components:
    J = J_real + 2*J_comp
# phi is used to stably compute exponentials between time steps:
    phi::Array{Real} = zeros(J, N-1)
# u, X & D are low-rank matrices and diagonal component:
    u::Array{Real} = zeros(J, N)
    X::Array{Real} = zeros(J, N)
    D::Array{Real} = zeros(N)
    diag::Array{Real} = zeros(N)

# Sum over the diagonal kernel amplitudes:
    a_sum = sum(a_real) + sum(a_comp)
# Compute the first element:
#    D[1] = sqrt(var[1] + a_sum)
    diag[1] = var[1]+a_sum
    D[1] = diag[1]
    value = 1.0 / D[1]
    for j in 1:J_real
        u[j, 1] = a_real[j]
        X[j, 1] = value
    end
# We are going to compute cosine & sine recursively - allocate arrays for each complex
# component:
    cd::Vector{Float64} = zeros(J_comp)
    sd::Vector{Float64} = zeros(J_comp)
# Initialize the computation of X:
    for j in 1:J_comp
        cd[j] = cos(d_comp[j]*t[1])
        sd[j] = sin(d_comp[j]*t[1])
        u[J_real+2*j-1, 1] = a_comp[j]*cd[j] + b_comp[j]*sd[j]
        u[J_real+2*j  , 1] = a_comp[j]*sd[j] - b_comp[j]*cd[j]
        X[J_real+2*j-1, 1] = cd[j]*value
        X[J_real+2*j, 1]   = sd[j]*value
    end
# Allocate array for recursive computation of low-rank matrices:
    S::Array{Float64, 2} = zeros(J, J)
#    for j in 1:J
#      for k in 1:j
#        S[k,j] = X[k,1]*X[j,1]
#      end
#    end
# Allocate temporary variables:
    phij = 0.0 ; dx = 0.0 ; dcd = 0.0 ; dsd = 0.0 ; cdtmp= 0.0 ; uj = 0.0 ;
    Xj = 0.0 ; Dn = 0.0 ; Sk = 0.0 ; uk = 0.0 ; tmp = 0.0 ; tn = 0.0 ;
# Loop over remaining indices:
    @inbounds for n in 2:N
        # Update phi
        tn = t[n]
# Using time differences stabilizes the exponential component and speeds
# up cosine/sine computation:
        dx = tn - t[n-1]
# Compute real components of the low-rank matrices:
        for j in 1:J_real
            phi[j, n-1] = exp(-c_real[j]*dx)
            u[j, n] = a_real[j]
            X[j, n] = 1.0
        end
# Compute complex components:
        for j in 1:J_comp
            value = exp(-c_comp[j]*dx)
            phi[J_real+2*j-1, n-1] = value
            phi[J_real+2*j, n-1]   = value
            cdtmp = cd[j]
            dcd = cos(d_comp[j]*dx)
            dsd = sin(d_comp[j]*dx)
            cd[j] = cdtmp*dcd-sd[j]*dsd
            sd[j] = sd[j]*dcd+cdtmp*dsd
        # Update u and initialize X
            u[J_real+2*j-1, n] = a_comp[j]*cd[j] + b_comp[j]*sd[j]
            u[J_real+2*j  , n] = a_comp[j]*sd[j] - b_comp[j]*cd[j]
            X[J_real+2*j-1, n  ] = cd[j]
            X[J_real+2*j  , n  ] = sd[j]
        end

        # Update S
        Dn = D[n-1]
        for j in 1:J
            phij = phi[j,n-1]
            Xj = X[j,n-1]
            for k in 1:j
                S[k, j] = phij*phi[k, n-1]*(S[k, j] + Dn*Xj*X[k,n-1])
            end
        end

        # Update D and X
        Dn = 0.0
        for j in 1:J
            uj = u[j,n]
            Xj = X[j,n]
            for k in 1:j-1
                Sk = S[k, j]
                tmp = uj * Sk
                uk = u[k,n]
                Dn += uk * tmp
                Xj -= uk*Sk
                X[k, n] -= tmp
            end
            tmp = uj * S[j, j]
            Dn += .5*uj * tmp
            X[j, n] = Xj - tmp
        end
# Finalize computation of D:
        diag[n] = var[n]+a_sum
        Dn = diag[n]-2.0*Dn
        D[n] = Dn
# Finalize computation of X:
        for j in 1:J
            X[j, n] /= Dn
        end
        # Update S
#        Xj = 0.0
#        for j in 1:J
#            Xj = X[j,n]
#            for k in 1:j
#                S[k, j] += Xj*X[k, n]
#            end
#        end
    end
# Finished looping over n.  Now return components to the calling routine
# so that these may be used in arithmetic:
    return D,X,u,phi
end

function cholesky!(a_real::Vector{<:Real}, c_real::Vector{<:Real},
                       a_comp::Vector{<:Real}, b_comp::Vector{<:Real},
                       c_comp::Vector{<:Real}, d_comp::Vector{<:Real},
                       t::Vector{<:Real}, diag::Vector{<:Real}, X::Array{<:Real,2},
                       phi::Array{<:Real,2}, u::Array{<:Real,2}, D::Vector{<:Real}, Q::Array{<:Real, 2})

    a_sum = sum(a_real) + sum(a_comp)
    M = length(Q[1,:])
    t = kron(t, ones(M))
    N = length(t)
    J_real = length(a_real) * M
    J_comp = length(a_comp) * M
    J = J_real + 2*J_comp

    # initialize the first row of variables (n=1):
    S::Array{Real} = zeros(J, J)
    V::Array{Real} = zeros(N, J)
    W::Array{Real} = zeros(N, J)
    U::Array{Real} = zeros(N, J)
    phi::Array{Real} = zeros(N, J)
    D::Vector{Real} = zeros(N)

    D[1] = sqrt(diag[1] + Q[1, 1]*a_sum)

    l = 0
    for j in 1:J_real
        p = 1
        q = ((j-1) % M) + 1
        l += (q == 1)
        U[1, (l-1)*M+q] = a_real[l]*Q[p, q]
        if p == q
            V[1, (l-1)*M+q] = 1
        end
        W[1, (l-1)*M+q] = V[1, (l-1)*M+q]/D[1]
    end

    l = 0
    for j in 1:J_comp
        p = 1
        q = ((j-1) % M) + 1
        l += (q == 1)
        minus = J_real + (l-1)*2*M + q
        plus = J_real + (2*l-1)*M + q
        sd = sin(d_comp[l]*t[1])
        cd = cos(d_comp[l]*t[1])
        U[1, minus] = (a_comp[l]*cd + b_comp[l]*sd)*Q[p, q]
        U[1, plus] = (a_comp[l]*sd - b_comp[l]*cd)*Q[p, q]
        if p == q
            V[1, minus] = cd
            V[1, plus] = sd
        end
        W[1, minus] = V[1, minus]/D[1]
        W[1, plus] = V[1, plus]/D[1]
    end

    # iterate over n
    r=1
    for i in 2:N
        p = ((i-1) % M) + 1
        r += (p == 1)

        # real terms
        l = 0
        for j in 1:J_real
            q = ((j-1) % M) + 1
            l += (q == 1)
            phi[i-1, (l-1)*M+q] = exp(-c_real[l]*(t[i] - t[i-1]))
            U[i, (l-1)*M+q] = a_real[l]*Q[p, q]
            if p == q
                V[i, (l-1)*M+q] = 1
            end
        end

        # complex terms
        l = 0
        for j in 1:J_comp
            q = ((j-1) % M) + 1
            l += (q == 1)
            minus = J_real + (l-1)*2*M + q
            plus = J_real + (2*l-1)*M + q
            phi[i-1, minus] = exp(-c_comp[l]*(t[i] - t[i-1]))
            phi[i-1, plus] = phi[i-1, minus]
            sd = sin(d_comp[l]*t[i])
            cd = cos(d_comp[l]*t[i])
            U[i, minus] = (a_comp[l]*cd + b_comp[l]*sd)*Q[p, q]
            U[i, plus] = (a_comp[l]*sd - b_comp[l]*cd)*Q[p, q]
            if p == q
                V[i, minus] = cd
                V[i, plus] = sd
            end
        end

        # compute the factorization
        sum_usu = 0
        sum_us::Array{Real} = zeros(J)

        for j in 1:J
            Uj = U[i, j]
            phij = phi[i-1, j]
            for k in 1:j-1
                Uk = U[i, k]
                S[j, k] = phij*phi[i-1, k] * (S[j, k] + W[i-1, j]*W[i-1, k])
                tmp = Uj*S[j, k]
                sum_usu += Uk*tmp
                sum_us[j] += Uk*S[j, k]
                sum_us[k] += tmp
            end
            S[j, j] = phi[i-1,j]*phi[i-1,j] *(S[j, j] + W[i-1,j]*W[i-1,j])
            sum_usu += 0.5*U[i, j]*U[i, j]*S[j, j]
            sum_us[j] += U[i, j]*S[j, j]
        end
        D[i] = sqrt(diag[i] + Q[p, p]*a_sum - 2*sum_usu)
        for j in 1:J
            W[i, j] = (V[i, j] - sum_us[j])/D[i]
        end
    end
    return D, W', U', phi'
end

function compute_ldlt!(gp::Celerite, x, yerr = 0.0)
# Call the choleksy function to decompose & update
# the components of gp with X,D,V,U,etc.
  coeffs = get_all_coefficients(gp.kernel)
  var = yerr.^2 + zeros(Float64, length(x))
  gp.n = length(x)
  gp.D,gp.W,gp.up,gp.phi = cholesky_ldlt!(coeffs..., x, var, gp.W, gp.phi, gp.up, gp.D)
  gp.J = size(gp.W)[1]
# Compute the log determinant (square the determinant of the Cholesky factor):
#  gp.logdet = sum(log(gp.D))
  logdet=0.0
  for i=1:gp.n
    logdet += log(gp.D[i])
  end
#  gp.logdet = sum(log(gp.D))
  gp.logdet = logdet
  gp.x = x
  gp.computed = true
  return gp.logdet
end

function compute!(gp::Celerite, x, yerr = 0.0)
# Call the choleksy function to decompose & update
# the components of gp with X,D,V,U,etc.
  coeffs = get_all_coefficients(gp.kernel)
  m = length(gp.Q[1,:])
  var::Array{Real} = yerr.^2 .+ zeros(m*length(x))
  gp.n = length(x)*length(gp.Q[1,:])
  gp.D,gp.W,gp.up,gp.phi = cholesky!(coeffs..., x, var, gp.W, gp.phi, gp.up, gp.D, gp.Q)
  gp.J = size(gp.W)[1]
# Compute the log determinant (square the determinant of the Cholesky factor):
  gp.logdet = 2 * sum(log.(gp.D))
#  gp.logdet = sum(log(gp.D))
  gp.x = x
  gp.computed = true
  return gp.logdet
end

function invert_lower_ldlt(gp::Celerite,y)
# Applies just the lower inverse:  L^{-1}.y:
  @assert(gp.computed)
  N = gp.n
  @assert(length(y)==N)
  z = zeros(Float64,N)
# The following lines solve L.z = y for z:
  z[1] = y[1]
  f = zeros(Float64,gp.J)
  for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* z[n-1])
#    f .= gp.phi[:,n] .* (f .+ gp.W[:,n-1] .* z[n-1])
    z[n] = y[n] - dot(gp.up[:,n], f)
  end
  return z
end

function invert_lower(gp::Celerite,y)
# Applies just the lower inverse:  L^{-1}.y:
  @assert(gp.computed)
  N = gp.n
  @assert(length(y)==N)
  z = zeros(Float64,N)
# The following lines solve L.z = y for z:
  z[1] = y[1]/gp.D[1]
  f = zeros(Float64,gp.J)
  for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* z[n-1])
    z[n] = (y[n] - sum(gp.up[:,n].*f))/gp.D[n]
  end
  return z
end

function apply_inverse_ldlt(gp::Celerite, y)
# Solves for K.b=y for b with LDLT decomposition.
  @assert(gp.computed)
  N = gp.n
  @assert(length(y)==N)
  z = zeros(Float64,N)
# The following lines solve L.z = y for z:
#  z[1] = y[1]/gp.D[1]
  z[1] = y[1]
  f = zeros(Float64,gp.J)
  for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* z[n-1])
    z[n] = (y[n] - dot(gp.up[:,n], f))
  end
# The following solves L^T.z = y for z:
  y = copy(z)
  fill!(z, zero(Float64))
  z[N] = y[N] / gp.D[N]
  fill!(f, zero(Float64))
  for n=N-1:-1:1
    f .= gp.phi[:,n] .* (f .+  gp.up[:,n+1] .* z[n+1])
    z[n] = y[n]/ gp.D[n] - dot(gp.W[:,n], f)
  end
# The result is the solution of L.L^T.z = y for z,
# or z = {L.L^T}^{-1}.y = L^{T,-1}.L^{-1}.y
  return z
end

function apply_inverse(gp::Celerite, y)
# Solves for K.b=y for b.
  @assert(gp.computed)
  N = gp.n
  @assert(length(y)==N)
  z::Array{Real} = zeros(N)
# The following lines solve L.z = y for z:
  z[1] = y[1]/gp.D[1]
  f::Array{Real} = zeros(gp.J)
  for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* z[n-1])
    z[n] = (y[n] - dot(gp.up[:,n], f))/gp.D[n]
  end
# The following solves L^T.z = y for z:
  y = copy(z)
  z = zeros(N)
  z[N] = y[N] / gp.D[N]
  f = zeros(gp.J)
  for n=N-1:-1:1
    f = gp.phi[:,n] .* (f +  gp.up[:,n+1].*z[n+1])
    z[n] = (y[n] - dot(gp.W[:,n], f)) / gp.D[n]
  end
# The result is the solution of L.L^T.z = y for z,
# or z = {L.L^T}^{-1}.y = L^{T,-1}.L^{-1}.y
  return z
end

function simulate_gp_ldlt(gp::Celerite, z)
# Multiplies lower Cholesky factor times random Gaussian vector (z is N(1,0) ) to simulate
# a Gaussian process.
# Check that Cholesky factor has been computed
@assert(gp.computed)
# Check for length mismatch:
N=gp.n
@assert(length(z)==N)
# If iid is zeros, then draw from random normal deviates:
if maximum(abs.(z)) == 0.0
  z .= randn.()
  #randn!(z) # In Julia 0.6 this can be replaced with "z .= randn.()"
end
y = zeros(Float64,N)
# Carry out multiplication
tmp = sqrt(gp.D[1])*z[1]
y[1] = tmp
f = zeros(Float64,gp.J)
for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* tmp)
    tmp = sqrt(gp.D[n])*z[n]
    y[n] = tmp + dot(gp.up[:,n], f)
end
# Return simulated correlated noise vector, y=L.z
return y
end

function multiply_ldlt(gp::Celerite, x, z, yerr=0.0)
# Multiplies the full matrix times a vector, z.
# Need to compute  onal, so get coefficients:
a1 , c1, a2, b2, c2, d2= get_all_coefficients(gp.kernel)
# Number of real components:
J1= length(a1)
# Number of complex components:
J2= length(a2)
# Rank of semi-separable components:
J = J1+ 2*J2
N = length(x)
@assert(length(z)==N)
# Carry out multiplication
asum = sum(a1)+sum(a2)
phi=zeros(J, N-1)
u=zeros(J, N-1)
v=zeros(J, N-1)
cd = cos.(d2 .* x[1])
sd = sin.(d2 .* x[1])
dx = 0.0
for n=1:N-1
  dx = x[n+1] - x[n]
  if J1 > 0
    v[1:J1,n] .= 1.0
    u[1:J1,n] .= a1
    phi[1:J1,n] .= exp.(-c1 .* dx)
  end
  if J2 > 0
    v[J1+1:J1+J2,n] .= cd
    v[J1+J2+1:J,n]  .= sd
    cd .= cos.(d2 .* x[n+1])
    sd .= sin.(d2 .* x[n+1])
    u[J1+1:J1+J2,n] .= a2 .* cd .+ b2.* sd
    u[J1+J2+1:J,n]  .= a2 .* sd .- b2 .* cd
    phi[J1+1:J1+J2,n] .= exp.(-c2 .* dx)
    phi[J1+J2+1:J,n] .= phi[J1+1:J1+J2,n]
  end
end
# This is A_{n,n} z_n from paper:
y = (yerr .^ 2 .+ asum) .* z
# sweep upwards in n:
f = zeros(Float64,gp.J)
for n =2:N
  f .= phi[:,n-1] .* (f .+ v[:,n-1] .* z[n-1])
  # This is \sum_{j=1}^J \tilde U_{n,j} f^-_{n,j}
  y[n] += dot(u[:,n-1],f)
end
# sweep downwards in n:
fill!(f, zero(Float64))
for n = N-1:-1:1
  f .= phi[:,n] .* (f .+  u[:,n] .* z[n+1])
  # This is \sum_{j=1}^J \tilde U_{n,j} f^-_{n,j}
  y[n] += dot(v[:,n],f)
end
# Return result of multiplication:
return y
end

function simulate_gp(gp::Celerite, y)
# Multiplies Cholesky factor times random Gaussian vector (y is N(1,0) ) to simulate
# a Gaussian process.
# If iid is zeros, then draw from random normal deviates:
# Check that Cholesky factor has been computed
# Carry out multiplication
# Return simulated correlated noise vector
N=gp.n
@assert(length(y)==N)
z = zeros(Float64,N)
z[1] = gp.D[1]*y[1]
f = zeros(Float64,gp.J)
for n =2:N
    f .= gp.phi[:,n-1] .* (f .+ gp.W[:,n-1] .* y[n-1])
    z[n] = gp.D[n]*y[n] + sum(gp.up[:,n] .* f)
end
# Returns z=L.y
n = length(gp.Q[1,:])
if n > 1
    ret = zeros(n, length(gp.x))
    for i in 1:n
        ret[i,:] = z[i:n:end]
    end
    return ret
else
    return z
end
end

function log_likelihood_ldlt(gp::Celerite, y)
# O(N) log likelihood computation once the low-rank Cholesky decomposition is completed
    @assert(gp.computed)
    if size(y, 2) != 1
        error("y must be 1-D")
    end
    alpha = apply_inverse_ldlt(gp, y)
    nll = gp.logdet + gp.n * log(2*pi)
    for i in 1:gp.n
        nll = nll + alpha[i] * y[i]
    end
    return -0.5 * nll
end

function log_likelihood(gp::Celerite, y)
# O(N) log likelihood computation once the low-rank Cholesky decomposition is completed
    @assert(gp.computed)
    y = vec(y)
    alpha = apply_inverse(gp, y)
    nll = gp.logdet + gp.n * log(2*pi)
    for i in 1:gp.n
        nll = nll + alpha[i] * y[i]
    end
    return -0.5 * nll
end

function full_solve(t::Vector,y0::Vector,aj::Vector,bj::Vector,cj::Vector,dj::Vector,yerr::Vector)
# This carries out the full GP solver using linear algebra on full covariance matrix.
# WARNING: do not use this with large datasets.
  N = length(t)
  J = length(aj)
  @assert(length(y0)==length(t))
  u = zeros(Float64,N,J*2)
  v = zeros(Float64,N,J*2)
# Compute the full U/V matrices:
  for i=1:N
    for j=1:J
      expcjt = exp(-cj[j]*t[i])
      cosdjt = cos(dj[j]*t[i])
      sindjt = sin(dj[j]*t[i])
      u[i,j*2-1]=aj[j]*expcjt*cosdjt+bj[j]*expcjt*sindjt
      u[i,j*2  ]=aj[j]*expcjt*sindjt-bj[j]*expcjt*cosdjt
      v[i,j*2-1]=cosdjt/expcjt
      v[i,j*2  ]=sindjt/expcjt
    end
  end
# Diagonal components:
#  diag = fill(yerr.^2 + sum(aj),N)
  diag = yerr.^2 .+ sum(aj)

# Compute the kernel:
  K = zeros(Float64,N,N)
  for i=1:N
    for j=1:N
      abs_ti_tj = abs(t[i] - t[j])
      K[i,j] = sum(aj .* exp.(-cj .* abs_ti_tj) .* cos.(dj .* abs_ti_tj) .+ bj .* exp.(-cj .* abs_ti_tj) .* sin.(dj .* abs_ti_tj))
    end
    K[i,i]=diag[i]
  end

# Check that equation (1) holds:
  K0 = tril(*(u, v'), -1) + triu(*(v, u'), 1)
  for i=1:N
    K0[i,i] = diag[i]
  end
  println("Semiseparable error: ", maximum(abs.(K .- K0)))
  return logdet(K),K
end

function predict_ldlt!(gp::Celerite, t, y, x, u=nothing, v=nothing, kern=nothing)
# Predict future times, x, based on a 'training set' of values y at times t.
# Runs in O((M+N)J^2) (variance is not computed, though)
    a_real, c_real, a_comp, b_comp, c_comp, d_comp = get_all_coefficients(gp.kernel)
    NQ = length(gp.Q[1,:])
    N = length(t)
    M = length(x)
    J_real = length(a_real)
    J_comp = length(a_comp)
    J = J_real + 2*J_comp

    b = apply_inverse_ldlt(gp, y)

    if kern != nothing
        NQ = length(u)
        MQ = length(v)
    else
        MQ = NQ
    end
    pred = zeros(M, MQ)
    z = zeros(N, MQ)
    Q = zeros(J, MQ)
    X = zeros(J)

    # Forward pass
    m = 1
    while m < M && x[m] <= t[1]
        m += 1
    end

    for n = 1:N
        if n < N
            tref = t[n+1]
        else
            tref = t[N]
        end
        if kern == nothing
            z[n,:] = gp.Q*b[(n-1)*NQ+1:n*NQ]
        else
            z[n,:] = predict1d!(kern, u, b[(n-1)*NQ+1:n*NQ], v)
        end
        for p = 1:MQ
            Q[1:J_real, p] .= (Q[1:J_real, p] .+ z[n,p]) .* exp.(-c_real .* (tref .- t[n]))
            Q[J_real+1:J_real+J_comp, p] .= (Q[J_real+1:J_real+J_comp, p] .+ z[n,p] .* cos.(d_comp .* t[n])) .*
                exp.(-c_comp .* (tref - t[n]))
            Q[J_real+J_comp+1:J, p] .= (Q[J_real+J_comp+1:J, p] .+ z[n,p] .* sin.(d_comp .* t[n])) .*
                exp.(-c_comp .* (tref - t[n]))
        end

        while m < M+1 && (n == N || x[m] <= t[n+1])
            X[1:J_real] .= a_real .* exp.(-c_real .* (x[m] - tref))
            X[J_real+1:J_real+J_comp] .= a_comp .* exp.(-c_comp .* (x[m] - tref)) .* cos.(d_comp .* x[m]) .+
                b_comp .* exp.(-c_comp .* (x[m] - tref)) .* sin.(d_comp .* x[m])
            X[J_real+J_comp+1:J] .= a_comp .* exp.(-c_comp .* (x[m] - tref)) .* sin.(d_comp .* x[m]) .-
                b_comp .* exp.(-c_comp .* (x[m] - tref)) .* cos.(d_comp .* x[m])

            for p = 1:MQ
                pred[m, p] = dot(X, Q[:, p])
            end
            m += 1
        end
    end

    # Backward pass
    m = M
    while m >= 1 && x[m] > t[N]
        m -= 1
    end
    fill!(Q, 0.0)
    for n=N:-1:1
        if n > 1
            tref = t[n-1]
        else
            tref = t[1]
        end
        for p=1:MQ
            Q[1:J_real, p] .= (Q[1:J_real, p] .+ z[n, p] .* a_real) .*
                exp.(-c_real .* (t[n] - tref))
            Q[J_real+1:J_real+J_comp, p] .= ((Q[J_real+1:J_real+J_comp, p] .+ z[n, p] .* a_comp .* cos.(d_comp .* t[n])) .+
                                          z[n, p] .* b_comp .* sin.(d_comp .* t[n])) .*
                                          exp.(-c_comp .* (t[n] - tref))
            Q[J_real+J_comp+1:J, p] .= (Q[J_real+J_comp+1:J, p] .+ z[n, p] .* a_comp .* sin.(d_comp .* t[n]) .-
                                     z[n, p] .* b_comp .* cos.(d_comp .* t[n])) .*
                                     exp.(-c_comp .* (t[n] - tref))
        end

        while m >= 1 && (n == 1 || x[m] > t[n-1])
            X[1:J_real] .= exp.(-c_real .* (tref - x[m]))
            X[J_real+1:J_real+J_comp] .= exp.(-c_comp .* (tref - x[m])) .* cos.(d_comp .* x[m])
            X[J_real+J_comp+1:J] .= exp.(-c_comp .* (tref - x[m])) .* sin.(d_comp .* x[m])

            for p = 1:MQ
                pred[m, p] += dot(X, Q[:, p])
            end
            m -= 1
        end
    end
    return pred
end

function predict!(gp::Celerite, y, t, x, u=nothing, v=nothing, kern=nothing)
# Predict future times, x, based on a 'training set' of values y at times t.
# Runs in O((M+N)J^2) (variance is not computed, though)
    a_real, c_real, a_comp, b_comp, c_comp, d_comp = get_all_coefficients(gp.kernel)
    NQ = length(gp.Q[1,:])
    N = length(t)
    M = length(x)
    J_real = length(a_real)
    J_comp = length(a_comp)
    J = J_real + 2*J_comp

    b = apply_inverse(gp, y)

    if kern != nothing
        NQ = length(u)
        MQ = length(v)
    else
        MQ = NQ
    end
    pred = zeros(M, MQ)
    z = zeros(N, MQ)
    Q = zeros(J, MQ)
    X = zeros(J)

    # Forward pass
    m = 1
    while m < M && x[m] <= t[1]
        m += 1
    end

    for n = 1:N
        if n < N
            tref = t[n+1]
        else
            tref = t[N]
        end
        if kern == nothing
            z[n,:] = gp.Q*b[(n-1)*NQ+1:n*NQ]
        else
            z[n,:] = predict1d!(kern, u, b[(n-1)*NQ+1:n*NQ], v)
        end
        for p = 1:MQ
            Q[1:J_real, p] .= (Q[1:J_real, p] .+ z[n,p]) .* exp.(-c_real .* (tref .- t[n]))
            Q[J_real+1:J_real+J_comp, p] .= (Q[J_real+1:J_real+J_comp, p] .+ z[n,p] .* cos.(d_comp .* t[n])) .*
                exp.(-c_comp .* (tref - t[n]))
            Q[J_real+J_comp+1:J, p] .= (Q[J_real+J_comp+1:J, p] .+ z[n,p] .* sin.(d_comp .* t[n])) .*
                exp.(-c_comp .* (tref - t[n]))
        end

        while m < M+1 && (n == N || x[m] <= t[n+1])
            X[1:J_real] .= a_real .* exp.(-c_real .* (x[m] - tref))
            X[J_real+1:J_real+J_comp] .= a_comp .* exp.(-c_comp .* (x[m] - tref)) .* cos.(d_comp .* x[m]) .+
                b_comp .* exp.(-c_comp .* (x[m] - tref)) .* sin.(d_comp .* x[m])
            X[J_real+J_comp+1:J] .= a_comp .* exp.(-c_comp .* (x[m] - tref)) .* sin.(d_comp .* x[m]) .-
                b_comp .* exp.(-c_comp .* (x[m] - tref)) .* cos.(d_comp .* x[m])

            for p = 1:MQ
                pred[m, p] = dot(X, Q[:, p])
            end
            m += 1
        end
    end

    # Backward pass
    m = M
    while m >= 1 && x[m] > t[N]
        m -= 1
    end
    fill!(Q, 0.0)
    for n=N:-1:1
        if n > 1
            tref = t[n-1]
        else
            tref = t[1]
        end
        for p=1:MQ
            Q[1:J_real, p] .= (Q[1:J_real, p] .+ z[n, p] .* a_real) .*
                exp.(-c_real .* (t[n] - tref))
            Q[J_real+1:J_real+J_comp, p] .= ((Q[J_real+1:J_real+J_comp, p] .+ z[n, p] .* a_comp .* cos.(d_comp .* t[n])) .+
                                          z[n, p] .* b_comp .* sin.(d_comp .* t[n])) .*
                                          exp.(-c_comp .* (t[n] - tref))
            Q[J_real+J_comp+1:J, p] .= (Q[J_real+J_comp+1:J, p] .+ z[n, p] .* a_comp .* sin.(d_comp .* t[n]) .-
                                     z[n, p] .* b_comp .* cos.(d_comp .* t[n])) .*
                                     exp.(-c_comp .* (t[n] - tref))
        end

        while m >= 1 && (n == 1 || x[m] > t[n-1])
            X[1:J_real] .= exp.(-c_real .* (tref - x[m]))
            X[J_real+1:J_real+J_comp] .= exp.(-c_comp .* (tref - x[m])) .* cos.(d_comp .* x[m])
            X[J_real+J_comp+1:J] .= exp.(-c_comp .* (tref - x[m])) .* sin.(d_comp .* x[m])

            for p = 1:MQ
                pred[m, p] += dot(X, Q[:, p])
            end
            m -= 1
        end
    end
    return pred'
end

function predict1d!(kern, u, b, v)
# straightforward prediction on second dimension covariance
# used by predict!
    if u == nothing || v == nothing
        error("requires u and v vectors")
        return
    else
        M = length(v)
        P = length(u)
        x = zeros(M)
        for p in 1:P
            for m in 1:M
                x[m] += kern(v[m], u[p])*b[p]
            end
        end
        return x
    end
end

function get_matrix(gp::Celerite, xs...)
# Gets the full covariance matrix.  Can provide autocorrelation or cross-correlation.
# WARNING: Do not use with large datasets.
    if length(xs) > 2
        error("At most 2 arguments can be provided")
    end
    local x1::Array
    local x2::Array
    if length(xs) >= 1
        x1 = xs[1]
    else
        if !gp.computed
            error("You must compute the GP first")
        end
        x1 = gp.x
    end
    if length(xs) == 2
        x2 = xs[2]
    else
        x2 = x1
    end

    if size(x1, 2) != 1 || size(x2, 2) != 1
        error("Inputs must be 1D")
    end

    tau = broadcast(-, reshape(x1, length(x1), 1), reshape(x2, 1, length(x2)))
    return kron(get_value(gp.kernel, tau), gp.Q)
end

function predict_full_ldlt(gp::Celerite, y, t; return_cov=true, return_var=false)
# Prediction with covariance using full covariance matrix
# WARNING: do not use this with large datasets!
    alpha = apply_inverse_ldlt(gp, y)
    Kxs = get_matrix(gp, t, gp.x)
    mu = Kxs * alpha
    if !return_cov && !return_var
        return mu
    end

    KxsT = Kxs'
    if return_var
        v = -sum(KxsT .* apply_inverse_ldlt(gp, KxsT), 1)
        v = v + get_value(gp.kernel, [0.0])[1]
        return mu, v[1, :]
    end

    cov = get_matrix(gp, t)
    cov = cov - Kxs * apply_inverse_ldlt(gp, KxsT)
    return mu, cov
end

function predict_full(gp::Celerite, y, t; return_cov=true, return_var=false)
# Prediction with covariance using full covariance matrix
# WARNING: do not use this with large datasets!
    alpha = apply_inverse(gp, y)
    Kxs = get_matrix(gp, t, gp.x)
    mu = Kxs * alpha
    if !return_cov && !return_var
        return mu
    end

    KxsT = transpose(Kxs)
    if return_var
        v = -sum(KxsT .* apply_inverse(gp, KxsT), 1)
        v = v + get_value(gp.kernel, [0.0])[1]
        return mu, v[1, :]
    end

    cov = get_matrix(gp, t)
    cov = cov - Kxs * apply_inverse(gp, KxsT)
    return mu, cov
end

function _reshape!(A::Array{<:Real}, dims...)
# Allocates arrays if size is not correct
    if size(A) != dims
        A = Array{Real}(undef, dims...)
    end
    return A
end


text:
@testset "cosine" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    k = CosineKernel()
    @test eltype(k) == Any
    @test kappa(k, 1.0) ≈ -1.0 atol=1e-5
    @test kappa(k, 2.0) ≈ 1.0 atol=1e-5
    @test kappa(k, 1.5) ≈ 0.0 atol=1e-5
    @test kappa(k, x) ≈ cospi(x) atol=1e-5
    @test k(v1, v2) ≈ cospi(sqrt(sum(abs2.(v1 - v2)))) atol=1e-5
    @test repr(k) == "Cosine Kernel"

    # Standardised tests.
    TestUtils.test_interface(k, Vector{Float64})
    test_ADs(CosineKernel)
end


text:
__precompile__()

module FNVHash

export fnv32, fnv16

# References:
#   http://www.isthe.com/chongo/tech/comp/fnv/
#   https://en.wikipedia.org/wiki/Fowler–Noll–Vo_hash_function
#   http://tools.ietf.org/html/draft-eastlake-fnv-03

function fnv32(data::Array{UInt8,1})
    hash = 0x811C9DC5
    for c in data
        hash = (hash * 0x01000193) $ c
    end
    return hash
end

fnv32(s::AbstractString) = fnv32(Array{UInt8,1}(s))

fnv16(hash::UInt32) = UInt16(hash >>> 16) $ UInt16(hash & 0xFFFF)
fnv16(data) = fnv16(fnv32(data))

end # module FNVHash


text:
"""
	subplot(fim=nothing; kwargs...)

Manage figure subplot configuration and selection.

Full option list at [`subplot`](http://docs.generic-mapping-tools.org/latest/subplot.html)

Parameters
----------

- **grid** :: [Type => Str | Tuple

    Specifies the number of rows and columns of subplots. Ex grid=(2,3)
- **F** | **dims** | **dimensions** | **size** | **sizes** :: [Type => Str | Tuple, NamedTuple]

    Specify the dimensions of the figure.
    [`-F`](http://docs.generic-mapping-tools.org/latest/subplot.html#f)
- **A** | **autolabel** | **fixedlabel** :: [Type => Str | number]

    Specify automatic tagging of each subplot. This sets the tag of the first, top-left subplot and others follow sequentially.
    [`-A`](http://docs.generic-mapping-tools.org/latest/subplot.html#a)
- $(GMT.opt_B)
- **C** | **clearance** :: [Type => Str | number]

    Reserve a space of dimension clearance between the margin and the subplot on the specified side. Settings specified under **begin** directive apply to all panels.
    [`-C`](http://docs.generic-mapping-tools.org/latest/subplot.html#c)
- $(GMT.opt_J)
- **M** | **margins** :: [Type => Str]

    The margin space that is added around each subplot beyond the automatic space allocated for tick marks, annotations, and labels.
    [`-M`](http://docs.generic-mapping-tools.org/latest/subplot.html#m)
- $(GMT.opt_R)
- **SC** | **SR** | **row_axes** | **col_axes** :: [Type => Str | NamedTuple]

    Set subplot layout for shared axes. Set separately for rows (SR) and columns (SC).
    [`-S`](http://docs.generic-mapping-tools.org/latest/subplot.html#s)
- **T** | **title** :: [Type => Str]

    While individual subplots can have titles, the entire figure may also have a overarching title.
    [`-T`](http://docs.generic-mapping-tools.org/latest/subplot.html#t)
"""
function subplot(fim=nothing; stop=false, kwargs...)

	global IamModern, IamSubplot

	d = KW(kwargs)
	# In case :title exists we must use and delete it to avoid double parsing
	cmd = ((val = find_in_dict(d, [:T :title], true)[1]) !== nothing) ? " -T\"" * val * "\"" : ""
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, "", false, " ")
	cmd = parse_common_opts(d, cmd, [:c :params], true)
	cmd = parse_these_opts(cmd, d, [[:M :margins]])
	cmd = add_opt(cmd, "A", d, [:A :autolabel :fixedlabel],
                  (Anchor=("+J", arg2str), anchor=("+j", arg2str), label="", clearance=("+c", arg2str), justify="+j", fill=("+g", add_opt_fill), pen=("+p", add_opt_pen), offset=("+o", arg2str), roman="_+r", Roman="_+R", vertical="_+v"))
	cmd = add_opt(cmd, "SC", d, [:SC :col_axes],
	              (top=("t", nothing, 1), bott=("b", nothing, 1), bottom=("b", nothing, 1), label="+l"))
	cmd = add_opt(cmd, "SR", d, [:SR :row_axes],
	              (left=("l", nothing, 1), right=("r", nothing, 1), label="+l", parallel="_+p", row_title="_+t", top_row_title="_+tc"))
	cmd = add_opt(cmd, "", d, [:C :clearance],
				  (left=(" -Cw", arg2str), right=(" -Ce", arg2str), bott=(" -Cs", arg2str), bottom=(" -Cs", arg2str), top=(" -Cn", arg2str)))

	if ((val = find_in_dict(d, [:F :dims :dimensions :size :sizes])[1]) !== nothing)
		if (isa(val, NamedTuple) && haskey(nt2dict(val), :width))	# Preferred way
			cmd *= " -F" * helper_sub_F(val)		# VAL = (width=x, height=x, fwidth=(...), fheight=(...))
		else
			cmd = add_opt(cmd, "F", d, [:F :dims :dimensions :size :sizes],
			              (panels=("s", nothing, 1), size=("", helper_sub_F), sizes=("", helper_sub_F), frac=("+f", helper_sub_F), fractions=("+f", helper_sub_F)))
		end
	end

	do_set = false;		do_show = false
	if (fim !== nothing)
		t = lowercase(string(fim))
		if     (t == "end" || t == "stop")  stop = true
		elseif (t == "show")  stop = true;  do_show = true
		elseif (t == "set")   do_set = true
		end
	end
	# ------------------------------ End parsing inputs --------------------------------

	if (!stop && !do_set)
		if (!haskey(d, :grid))  error("SUBPLOT: 'grid' keyword is mandatory")  end
		cmd = arg2str(d[:grid], 'x') * " " * cmd
		if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end		# Vd=2 cause this return

		if (!IamModern)			# If we are not in modern mode, issue a gmt("begin") first
			fname = ""			# Default name (GMTplot.ps) is set in gmt_main()
			if ((val = find_in_dict(d, [:name :savefig])[1]) !== nothing)
				fname = get_format(string(val), nothing, d)		# Get the fig name and format.
			end
			gmt("begin " * fname)
		end
		gmt("subplot begin " * cmd);
		IamSubplot = true
	elseif (do_set)
		if (!IamSubplot)  error("Cannot call subplot(set, ...) before setting dimensions")  end
		lix, pane = parse_c(cmd, d)
		cmd = pane * cmd				# Here we don't wanr the "-c" part
		if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end		# Vd=2 cause this return
		gmt("subplot set " * cmd)
	else
		(do_show || haskey(d, :show)) ? show = " show" : show = ""
		try
			gmt("subplot end");		gmt("end" * show);		catch
		end
		IamModern = false;		IamSubplot = false
	end
	return nothing
end

# --------------------------------------------------------------------------
function helper_sub_F(arg, dumb=nothing)
	# dims=(1,2)
	# dims=(size=(1,2), frac=((2,3),(3,4,5)))
	# dims=(width=xx, height=yy, fwidth=(), fheight=())
	out = ""
	if (isa(arg, String))
		out = arg2str(arg)
	elseif (isa(arg, NamedTuple) || isa(arg, Dict))
		d = (isa(arg, NamedTuple)) ? nt2dict(arg) : arg
		val, symb = find_in_dict(d, [:size :sizes :frac :fractions])	# ex: dims=(size=(1,2), frac=((2,3),(3,4,5)))
		if (val !== nothing)
			if (isa(val, Tuple{Tuple, Tuple}))		# ex: dims=(panels=true, sizes=((2,4),(2.5,5,1.25)))
				out *= arg2str(val[1], ',') * '/' * arg2str(val[2], ',')
			else
				if (symb == :size || symb == :sizes)  out = arg2str(val)
				else                                  error("'frac' option must be a tuple(tuple, tuple)")
				end
			end
		else
			if (!haskey(d, :width))  error("SUBPLOT: 'width' is a mandatory parameter")  end
			out = string(d[:width], '/')
			if (!haskey(d, :height))  out *= '0'	# Allow this for geog cases
			else                      out *= string(d[:height])
			end
			if (haskey(d, :fwidth))
				out *= "+f" * arg2str(d[:fwidth], ',')
				if (!haskey(d, :fheight))  error("SUBPLOT: when using 'fwidth' must also set 'fheight'")  end
				out *= '/' * arg2str(d[:fheight], ',')
			end
		end
	end
	if (out == "")  error("SUBPLOT: garbage in DIMS option")  end
	return out
end