text:
struct TierrianPosition <: Position
    start::UInt16
    length::UInt16
end

struct TierrianEnvironment <: Environment end

get_time(tierrian_model::TierraModel) = tierrian_model.time
get_organisms(tierrian_model::TierraModel) = [o for o in values(tierrian_model.organisms)]

get_logger(tierrian_model::TierraModel) = tierrian_model.logger
function set_logger!(tierrian_model::TierraModel, logger::Logger)
    tierrian_model.logger = logger
end

get_id(::TierraModel, tierrian_organism::TierrianOrganism) = tierrian_organism.key
function get_genotype(tierrian_model::TierraModel, tierrian_organism::TierrianOrganism)
    read_memory(tierrian_model, tierrian_organism.start_address, tierrian_organism.length)
end
function get_position(::TierraModel, tierrian_organism::TierrianOrganism)
    TierrianPosition(tierrian_organism.start_address, tierrian_organism.length)
end
function get_environment(::TierraModel, ::TierrianOrganism)
    TierrianEnvironment()
end
function get_time_birth(::TierraModel, tierrian_organism::TierrianOrganism)
    tierrian_organism.time_birth
end
function get_time_death(::TierraModel, tierrian_organism::TierrianOrganism)
    tierrian_organism.time_birth + tierrian_organism.age
end
function get_age(::TierraModel, tierrian_organism::TierrianOrganism)
    tierrian_organism.age
end

get_daughters(::TierraModel, tierrian_organism::TierrianOrganism) = tierrian_organism.daughters

function get_abstracted_organism(model::TierraModel, model_organism::TierrianOrganism, parent_ids, parent_genotypes)
    return Organism{TierrianPosition, TierrianEnvironment, Vector{UInt8}}(
                    get_id(model, model_organism),
                    get_genotype(model, model_organism),
                    parent_ids, parent_genotypes,
                    get_position(model, model_organism),
                    get_environment(model, model_organism),
                    get_time_birth(model, model_organism),
                    get_time_death(model, model_organism),
                    -1
    )
end

function replace_organism!(model::TierraModel, key::UInt64, new_program::Vector{UInt8})
    # remove old organism first

    remove_organism!(model, key)

    # add new organism
    # TODO: same place

    new_key = add_organism!(model, new_program)

    if new_key === nothing
        throw(SimulationExpection("New organism cannot be created in Tierra -- not enough memory."))
    end

    new_organism = model.organisms[new_key]

    return new_organism
end

function run_until!(termination_predicate, model::TierraModel)
    while termination_predicate(model) == false
        execute_slice!(model)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
function run_Cylinder(parts,order::Int,dt::Real,tf::Real,mesh_file::String)

  # Parameters
  D = 0.1
  Re = 100
  ν = D/Re
  Uₘ = 1.5
  H = 0.41

  # TGV initial solution
  uin((x,y,z),t) = VectorValue(4*Uₘ*y*(H-y)/(H^2),0,0)
  u₀((x,y,z),t) = VectorValue(0,0,0)
  uin(t) = x->uin(x,t)
  u₀(t) = x->u₀(x,t)

  # Discretization
  models_path=ENV["DELFTBLUE_MODELS"]
  meshfile = joinpath(models_path,mesh_file)
  𝒯 = GmshDiscreteModel(parts,meshfile)

  # Triangulation and Integration measure
  Ω = Interior(𝒯)
  dΩ = Measure(Ω,2*order)
  Γ = Boundary(𝒯,tags=["sides"])
  dΓ = Measure(Γ,2*order)
  nΓ = get_normal_vector(Γ)

  # FE spaces
  refFEᵤ = ReferenceFE(lagrangian,VectorValue{3,Float64},order)
  refFEₚ = ReferenceFE(lagrangian,Float64,order-1)
  V = TestFESpace(𝒯,refFEᵤ,dirichlet_tags=["inlet","cylinder"])
  Q = TestFESpace(𝒯,refFEₚ)
  U = TransientTrialFESpace(V,[uin,u₀])
  P = TrialFESpace(Q)
  Y = MultiFieldFESpace([V,Q])
  X = TransientMultiFieldFESpace([U,P])

  # Explicit FE functions
  global ηₙₕ = interpolate_everywhere(VectorValue(0.0,0.0,0.0),U(0.0))
  global uₙₕ = interpolate_everywhere(u₀(0.0),U(0.0))
  global fv_u = zero_free_values(U(0.0))

  # Stabilization Parameters
  c₁ = 12.0
  c₂ = 2.0
  cc = 4.0

  # h = map_parts(Ω.trians) do trian
  #   lazy_map(h->h^(1/3),get_cell_measure(trian))
  # end
  h = 0.01
  τₘ = 1/(c₁*ν/h^2 + c₂*(meas∘uₙₕ)/h)
  τc = cc *(h^2/(c₁*τₘ))
  κ = 10.0*order*(order-1)/h

  # Weak form
  c(a,u,v) = 0.5*((∇(u)'⋅a)⋅v - u⋅(∇(v)'⋅a))
  res(t,(u,p),(v,q)) = ∫( ∂t(u)⋅v  + c(u,u,v) + ν*(∇(u)⊙∇(v)) - p*(∇⋅v) + (∇⋅u)*q +
                          τₘ*((∇(u)'⋅u - ηₙₕ)⋅(∇(v)'⋅u)) + τc*((∇⋅u)*(∇⋅v)) )dΩ +
                       ∫( - ((ν*(∇(u)⋅nΓ) - p*nΓ)⋅nΓ) ⋅ (v⋅nΓ) - ((ν*(∇(v)⋅nΓ) - q*nΓ)⋅nΓ) ⋅ (u⋅nΓ) + κ*(u⋅nΓ)*(v⋅nΓ) )dΓ
  jac(t,(u,p),(du,dp),(v,q)) = ∫( c(du,u,v) + c(u,du,v) + ν*(∇(du)⊙∇(v)) - dp*(∇⋅v) + (∇⋅du)*q +
                                  τₘ*((∇(u)'⋅u - ηₙₕ)⋅(∇(v)'⋅du) + (∇(du)'⋅u + ∇(u)'⋅du)⋅(∇(v)'⋅u)) +
                                  τc*((∇⋅du)*(∇⋅v)) )dΩ+
                               ∫( - ((ν*(∇(du)⋅nΓ) - dp*nΓ)⋅nΓ) ⋅ (v⋅nΓ) - ((ν*(∇(v)⋅nΓ) - q*nΓ)⋅nΓ) ⋅ (du⋅nΓ) + κ*(du⋅nΓ)*(v⋅nΓ) )dΓ
  jac_t(t,(u,p),(dut,dpt),(v,q)) = ∫( dut⋅v )dΩ
  op = TransientFEOperator(res,jac,jac_t,X,Y)

  # Stokes Weake form
  a₀((u,p),(v,q)) = ∫( ν*(∇(u)⊙∇(v)) - p*(∇⋅v) + (∇⋅u)*q )dΩ +
                    ∫( - ((ν*(∇(u)⋅nΓ) - p*nΓ)⋅nΓ) ⋅ (v⋅nΓ) - ((ν*(∇(v)⋅nΓ) - q*nΓ)⋅nΓ) ⋅ (u⋅nΓ) + κ*(u⋅nΓ)*(v⋅nΓ) )dΓ
  l₀((v,q)) =  ∫( 0.0*q )dΩ
  op₀ = AffineFEOperator(a₀,l₀,X(0.0),Y)

  # Orthogonal projection
  a(η,κ) = ∫( τₘ*(η⋅κ) )dΩ
  b(κ) = ∫( τₘ*((∇(uₙₕ)'⋅uₙₕ)⋅κ) )dΩ
  op_proj(t) = AffineFEOperator(a,b,U(t),V)

  # Setup solver via low level PETSC API calls
  function mykspsetup(ksp)
    pc       = Ref{GridapPETSc.PETSC.PC}()
    mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
    @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY)
    @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc)
    @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU)
    @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS)
    @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[])
    @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat)
    @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  4, 2)
    @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2)
    @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2)
    @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6)
    @check_error_code GridapPETSc.PETSC.KSPSetFromOptions(ksp[])
  end

   # Linear Solver
   ls₀ = PETScLinearSolver(mykspsetup)
   ls = PETScLinearSolver()

  # Nonlinear Solver
  nls = NLSolver(ls,show_trace=true,method=:newton,iterations=10)
  #nls = PETScNonlinearSolver()

  # Transient solver
  ode_solver = GeneralizedAlpha(nls,dt,0.0)

  # Initial solution
  xₕ₀ = solve(ls₀,op₀)
  #vₕ₀ = interpolate_everywhere([VectorValue(0.0,0.0,0.0),0.0],X(0.0))
  du₀ = interpolate_everywhere(VectorValue(0.0,0.0,0.0),U(0.0))
  dp₀ = interpolate_everywhere(0.0,P)
  vₕ₀ = interpolate_everywhere([du₀,dp₀],X(0.0))

  # Solution (lazy)
  xₕₜ = solve(ode_solver,op,(xₕ₀,vₕ₀),0,tf)

  # Iterate over steps
  createpvd(parts,"results/Cylinter") do pvd
    for (xₕ,t) in xₕₜ
      println("--------------------")
      println("Start of time: $t")
      uₕ,pₕ=xₕ
      println("updating global variables")
      uₙₕ = interpolate!(uₕ,fv_u,U(t))
      ηₙₕ = solve(ls,op_proj(t))
      pvd[t] = createvtk(Ω,"results/Cylinder_$t";cellfields=["u"=>uₕ,"p"=>pₕ,"eta"=>ηₙₕ],nsubcells=10)
      println("End of time: $t")
    end
 end

end

ext: jl


text:
using SignalAnalysis: signal
using FFTW: ifft!
using DSP: nextfastfft

export SoundSpeedProfile, soundspeed
export Bathymetry, depth, maxdepth
export Altimetry, altitude
export ReflectionModel, reflectioncoef
export UnderwaterEnvironment, altimetry, bathymetry, ssp, salinity, seasurface, seabed, noise
export AcousticSource, AcousticReceiver, location, nominalfrequency, phasor, record, recorder
export PropagationModel, arrivals, transfercoef, transmissionloss, eigenrays, rays
export NoiseModel
export impulseresponse

### interface: SoundSpeedProfile

abstract type SoundSpeedProfile end

"""
    soundspeed(ssp::SoundSpeedProfile, x, y, z)

Get sound speed at location (`x`, `y`, `z`). If a sound speed profile is range
independent, `x` and `y` may be ignored. `z` is generally negative, since the
sea surface is the datum and z-axis points upwards.
"""
function soundspeed end

### interface: Bathymetry

abstract type Bathymetry end

"""
    depth(bathy:Bathymetry, x, y)

Get water depth at location (`x`, `y`).
"""
function depth end

"""
    maxdepth(bathy::Bathymetry)

Get the maximum water depth.
"""
function maxdepth end

### interface: Altimetry

abstract type Altimetry end

"""
    altitude(alt:Altimetry, x, y)

Get water surface altitude at location (`x`, `y`). The nominal water surface
is considered to have an altitude of zero. However, the water surface may
not be flat, and the Altimetry provisions for variations of altitude around the
nominal altitutde of zero.
"""
function altitude end

### interface: ReflectionModel

abstract type ReflectionModel end

"""
    reflectioncoef(rm::ReflectionModel, f, θ)

Get complex reflection coefficient at frequency `f` Hz and incidence angle `θ`
(w.r.t. the surface normal).
"""
function reflectioncoef end

### interface: UnderwaterEnvironment

abstract type UnderwaterEnvironment end

"""
    altimetry(env::UnderwaterEnvironment)::Altimetry

Get the altimetry for the underwater environment.
"""
function altimetry end

"""
    bathymetry(env::UnderwaterEnvironment)::Bathymetry

Get the bathymetry for the underwater environment.
"""
function bathymetry end

"""
    ssp(env::UnderwaterEnvironment)::SoundSpeedProfile

Get the sound speed profile for the underwater environment.
"""
function ssp end

"""
    salinity(env::UnderwaterEnvironment)

Get the salinity of the underwater environment.
"""
function salinity end

"""
    seasurface(env::UnderwaterEnvironment)::ReflectionModel

Get the sea surface reflection model for the underwater environment.
"""
function seasurface end

"""
    seabed(env::UnderwaterEnvironment)::ReflectionModel

Get the seabed reflection model for the underwater environment.
"""
function seabed end

"""
    noise(env::UnderwaterEnvironment)::NoiseModel

Get the noise model for the underwater environment.
"""
function noise end

### interface: AcousticSource

abstract type AcousticSource end

"""
    location(src::AcousticSource)
    location(src::AcousticReceiver)

Get the location of an acoustic source or receiver as a 3-tuple (`x`, `y`, `z`).
"""
function location end

"""
    nominalfrequency(src::AcousticSource)

Get the nominal frequency of an acoustic source in Hz.
"""
function nominalfrequency end

"""
    nominalfrequency(src::AcousticSource)

Get the complex phasor representation (amplitude & phase) of a narrowband
acoustic source at the nominal frequency.
"""
function phasor end

"""
    record(src::AcousticSource, duration, fs; start=0.0)
    record(noise::NoiseModel, duration, fs; start=0.0)
    record(model::PropagationModel, tx, rx, duration, fs; start=0.0)

Make a recording of an acoustic source or ambient noise. The `start` time and
`duration` are specified in seconds, and the recording is made at a sampling
rate of `fs` Hz.

For an recording of an acoustic source, free space propagation is assumed,
and the recording is made at a nominal range of 1 meter from the acoustic
center of the source.

For a recording through a propagation model, `tx` and `rx` may be single `AcousticSource`
and `AcousticReceiver`, or an array each.
"""
function record end

### interface: AcousticReceiver

abstract type AcousticReceiver end
function location end

### interface: NoiseModel

abstract type NoiseModel end
function record end

### interface: PropagationModel

abstract type PropagationModel{T<:UnderwaterEnvironment} end

"""
    environment(pm::PropagationModel)

Get the environment associated with the propagation model.
"""
function environment end

"""
    check(pm::Type{<:PropagationModel}, env::UnderwaterEnvironment)
    check(pm::Type{<:PropagationModel}, env=missing)

Check if an propagation model is available, and can simulate the specified
environment. Returns the environment if it can be simulated, or throws an
error with a descriptive error message if it cannot be simulated.

This function is internally used by the propagation modeling toolbox to choose
a model or offer a selection of models to the user.
"""
function check end

"""
    arrivals(pm::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver)

Compute the arrivals from `tx1` to `rx1`. Returns an array of `Arrival` structs.
"""
function arrivals end

"""
    transfercoef(pm::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver; mode=:coherent)
    transfercoef(pm::PropagationModel, tx1::AcousticSource, rx::AbstractArray{<:AcousticReceiver}; mode=:coherent)

Compute the complex transfer coefficients from `tx1` to `rx1` or all receivers in `rx`.
The mode may be `:coherent` or `:incoherent`.
"""
function transfercoef end

"""
    transmissionloss(pm::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver; mode=:coherent)
    transmissionloss(pm::PropagationModel, tx1::AcousticSource, rx::AbstractArray{<:AcousticReceiver}; mode=:coherent)

Compute the transmission loss in dB from `tx1` to `rx1` or all receivers in `rx`.
The mode may be `:coherent` or `:incoherent`.
"""
function transmissionloss end

"""
    eigenrays(pm::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver)

Compute the eigenrays from `tx1` to `rx1`. Returns an array of `RayArrival` structs.
"""
function eigenrays end

"""
    rays(pm::PropagationModel, tx1::AcousticSource, θ::Real, rmax)

Compute the rays from `tx1` launched at angle `θ` (or all angles in `θ`, if it is a vector).
Returns an array of `RayArrival` datatypes. `rmax` is the maximum horizontal range in meters
to track the rays over.
"""
function rays end

function record end

"""
    recorder(model::PropagationModel, tx, rx)

Create a recorder function that may be called later to make an acoustic recording
of sources in `tx` at receviers `rx`. `tx` and `rx` may be single `AcousticSource` and
`AcousticReceiver`, or an array each.

The recorder function may be called later with `duration`, `fs`, and optionally a `start`
time. It functions in a similar way as the `record()` function.

# Examples:
```julia-repl
julia> rec = recorder(pm, tx, rx);
julia> s = rec(1.0, 44100.0; start=0.0);  # make a recording of 1 second at 44.1 kHz
```
"""
function recorder end

### arrival types

abstract type Arrival end

struct RayArrival{T1,T2} <: Arrival
  time::T1
  phasor::Complex{T1}
  surface::Int
  bottom::Int
  launchangle::T1
  arrivalangle::T1
  raypath::Union{Vector{NTuple{3,T2}},Missing}
end

function RayArrival(time::T1, phasor::Complex{T1}, surface::Int, bottom::Int, launchangle::T1, arrivalangle::T1, raypath::Vector{NTuple{3,T2}}) where {T1,T2}
  RayArrival{T1,T2}(time, phasor, surface, bottom, launchangle, arrivalangle, raypath)
end

function RayArrival(time, phasor, surface::Int, bottom::Int, launchangle, arrivalangle, raypath::Vector{NTuple{3,T}}) where T
  t, r, i, θ1, θ2 = promote(time, real(phasor), imag(phasor), launchangle, arrivalangle)
  RayArrival{typeof(t),T}(t, Complex(r, i), surface, bottom, θ1, θ2, raypath)
end

function RayArrival(time, phasor, surface::Int, bottom::Int, launchangle, arrivalangle)
  t, r, i, θ1, θ2 = promote(time, real(phasor), imag(phasor), launchangle, arrivalangle)
  RayArrival{typeof(t),Missing}(t, Complex(r, i), surface, bottom, θ1, θ2, missing)
end

phasortype(::Type{RayArrival{T1,T2}}) where {T1,T2} = T1

### fallbacks & helpers

location(x::NTuple{3,T}) where T = x
location(x::NTuple{2,T}) where T = (x[1], zero(T), x[2])

check(model::PropagationModel, env) = env
environment(model::PropagationModel) = model.env

function transfercoef(model::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver; mode=:coherent)
  arr = arrivals(model, tx1, rx1)
  length(arr) == 0 && return zero(phasortype(eltype(arr)))
  if mode === :coherent
    f = nominalfrequency(tx1)
    tc = sum(a.phasor * cis(-2π * a.time * f) for a ∈ arr)
  elseif mode === :incoherent
    tc = Complex(√sum(abs2(a.phasor) for a ∈ arr), 0)
  else
    throw(ArgumentError("Unknown mode :" * string(mode)))
  end
  tc
end

function transfercoef(model::PropagationModel, tx1::AcousticSource, rx::AbstractArray{<:AcousticReceiver}; mode=:coherent)
  # threaded version of [transfercoef(model, tx1, rx1; mode=mode) for rx1 ∈ rx], seems to be faster than tmap()
  rx1 = first(rx)
  tc1 = transfercoef(model, tx1, rx1; mode=mode)
  tc = Array{typeof(tc1)}(undef, size(rx))
  Threads.@threads for i ∈ eachindex(rx)
    tc[i] = rx1 === rx[i] ? tc1 : transfercoef(model, tx1, rx[i]; mode=mode)
  end
  tc
end

function rays(model::PropagationModel, tx1::AcousticSource, θ::AbstractArray, rmax)
  tmap(θ1 -> rays(model, tx1, θ1, rmax), θ)
end

transmissionloss(model, tx, rx; mode=:coherent) = -amp2db.(abs.(transfercoef(model, tx, rx; mode=mode)))

function recorder(model::PropagationModel, tx1::AcousticSource, rx1::AcousticReceiver)
  f = recorder(model, [tx1], [rx1])
  function rec(duration, fs; start=0.0)
    x = f(duration, fs; start=start)
    signal(dropdims(samples(x); dims=2), framerate(x))
  end
end

function recorder(model::PropagationModel, tx1::AcousticSource, rx::AbstractArray{<:AcousticReceiver})
  recorder(model, [tx1], rx)
end

function recorder(model::PropagationModel, tx::AbstractArray{<:AcousticSource}, rx1::AcousticReceiver)
  f = recorder(model, tx, [rx1])
  function rec(duration, fs; start=0.0)
    x = f(duration, fs; start=start)
    signal(dropdims(samples(x); dims=2), framerate(x))
  end
end

function recorder(model::PropagationModel, tx::AbstractArray{<:AcousticSource}, rx::AbstractArray{<:AcousticReceiver})
  arr = [arrivals(model, tx1, rx1) for tx1 ∈ tx, rx1 ∈ rx]
  mindelay, maxdelay = extrema(Iterators.flatten([[a1.time for a1 ∈ a] for a ∈ arr]))
  function rec(duration, fs; start=0.0)
    src = [record(tx1, duration + (maxdelay-mindelay) + 1/fs, fs; start=start-maxdelay) for tx1 ∈ tx]
    nsamples = round(Int, duration * fs)
    x = zeros(Base.promote_eltype(src...), nsamples, length(rx))
    for j = 1:length(tx)
      for k = 1:length(rx)
        for a ∈ arr[j,k]
          t = round(Int, (maxdelay - a.time) * fs)
          x[:,k] .+= a.phasor .* src[j][t+1:t+nsamples]
        end
      end
    end
    noisemodel = noise(environment(model))
    if noisemodel !== missing
      for k = 1:length(rx)
        x[:,k] .+= record(noisemodel, duration, fs; start=start)
      end
    end
    signal(x, fs)
  end
end

function record(model::PropagationModel, tx, rx, duration, fs; start=0.0)
  recorder(model, tx, rx)(duration, fs; start=start)
end

"""
$(SIGNATURES)
Convert a vector of arrivals to a sampled impulse response time series at a
sampling rate of `fs` Hz. If `ntaps` is zero, the number of taps of the impulse
response are chosen automatically.

If `reltime` is `true`, the impulse response start time is relative to the
first arrival, otherwise it is relative to the absolute time. If `approx`
is `true`, a fast algorithm is used to generate a sparse impulse response
that assigns an arrival to the nearest sampling time.
"""
function impulseresponse(arrivals::Vector{<:Arrival}, fs, ntaps=0; reltime=false, approx=false)
  length(arrivals) == 0 && throw(ArgumentError("No arrivals"))
  mintime, maxtime = extrema(a.time for a ∈ arrivals)
  reltime || (mintime = zero(typeof(arrivals[1].time)))
  mintaps = ceil(Int, (maxtime-mintime) * fs) + 1
  if approx
    ntaps == 0 && (ntaps = mintaps)
    ir = zeros(typeof(arrivals[1].phasor), ntaps)
    for a ∈ arrivals
      ndx = round(Int, (a.time - mintime) * fs) + 1
      ndx ≤ ntaps && (ir[ndx] = a.phasor)
    end
  else
    N = nextfastfft(4*max(256, mintaps))
    x = zeros(typeof(arrivals[1].phasor), N)
    for a ∈ arrivals
      δ = (a.time - mintime) * fs
      for i ∈ 1:N
        x[i] += a.phasor * cis(-2π * (i-1)/N * δ)
      end
    end
    ifft!(x)
    if ntaps == 0
      ndx = findfirst(abs.(x[mintaps+1:end]) .≤ abs(arrivals[1].phasor)/100)
      ndx === nothing && (ndx = argmin(abs.(x[mintaps+1:end])))
      ir = x[1:mintaps+ndx]
    elseif ntaps ≤ N/2
      ir = x[1:ntaps]
    else
      Nby2 = floor(Int, N/2)
      ir = vcat(x[1:Nby2], zeros(typeof(arrivals[1].phasor), ntaps-Nby2))
    end
  end
  ir
end

function tmap(f, itr)
  refs = [Threads.@spawn(f(i)) for i ∈ itr]
  fetch.(refs)
end

function envrealtype(env::UnderwaterEnvironment)
  a = altitude(altimetry(env), 0.0, 0.0)
  b = depth(bathymetry(env), 0.0, 0.0)
  c = soundspeed(ssp(env), 0.0, 0.0, 0.0)
  s = salinity(env)
  r1 = real(reflectioncoef(seasurface(env), 1000.0, 0.0))
  r2 = real(reflectioncoef(seabed(env), 1000.0, 0.0))
  promote_type(typeof(a), typeof(b), typeof(c), typeof(s), typeof(r1), typeof(r2))
end

### pretty printing

function Base.show(io::IO, env::UnderwaterEnvironment)
  println(io, split(string(typeof(env).name), ".")[end], ":")
  println(io, "  altimetry = ", altimetry(env))
  println(io, "  bathymetry = ", bathymetry(env))
  println(io, "  ssp = ", ssp(env))
  println(io, "  salinity = ", salinity(env))
  println(io, "  seasurface = ", seasurface(env))
  println(io, "  seabed = ", seabed(env))
  println(io, "  noise = ", noise(env))
end

function Base.show(io::IO, model::PropagationModel)
  print(io, split(string(typeof(model).name), ".")[end], " with ")
  show(io, environment(model))
end

function Base.show(io::IO, a::Arrival)
  if isnan(a.time) || isnan(a.phasor)
    @printf(io, "∠%5.1f° %2d↑ %2d↓%s",
      rad2deg(a.launchangle), a.surface, a.bottom, a.raypath === missing || length(a.raypath) == 0 ? "" : " ⤷")
  else
    @printf(io, "∠%5.1f° %2d↑ %2d↓ ∠%5.1f° | %6.2f ms | %5.1f dB ϕ%6.1f°%s",
      rad2deg(a.launchangle), a.surface, a.bottom, rad2deg(a.arrivalangle),
      1000*a.time, amp2db(abs(a.phasor)), rad2deg(angle(a.phasor)),
      a.raypath === missing || length(a.raypath) == 0 ? "" : " ⤷")
  end
end


text:
using BEASTSimulation
using Test

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


text:
function readbin(filename)

    a = reinterpret(Float64,read(filename));
    return a

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
"""
  Returns the a_j,k and c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )
"""
function get_AC_coefficients(K::Int64
                            , times::Vector
                            , z_p::Vector{Complex{T}}
                            , alpha::Float64
                            , m::Float64
                            ; method = "nfft") where T

  if method == "leastsquare"
    A,C = get_AC_coefficients_lsqr(K,times,z_p)
  elseif method == "nfft"
    A,C = get_AC_coefficients_nfft(K,times,z_p,alpha,m)
  elseif method == "hist"
    A,C = get_AC_coefficients_hist_lsqr(K,times,z_p)
  else
    A = zeros(Complex{T},length(times),K)
    C = zeros(Complex{T},K,length(z_p))
    error("Not implemented yet")
  end
  return A,C
end

function _precompute_coeffs_nfft(K::Int64
                                , times::Vector     # readout times
                                , z_p::Vector{Complex{T}}      # Correction map
                                , alpha::Float64   # oversampling factor of NFFT
                                , m::Float64       # kernel size of NFFT
                                ; adaptK::Bool=true
                                , centering::Bool=true) where T

  # Centering time and correctionterm to lower computational effort
  t_hat = (times[1] + times[end])/2
  z_hat = minimum(real(z_p)) + maximum(real(z_p)) + 1im*(minimum(imag(z_p)) + maximum(imag(z_p)))
  z_hat *= 0.5
  if !centering
    t_hat=0.0
    z_hat=0.0
  end
  # Centering readouttimes and correctionmap
  times = times .- t_hat
  z_p = z_p .- z_hat

  # Calculating timescaling factor so it satisfies: t_j / T in [0.5,0.5)
  t_scale = maximum(abs.(times)) / (0.5-1e-7)

  # determining N and K,
  # K can be choosen freely with tradeoff comp. complexity vs precision
  N = Int64(4*ceil( maximum(abs.(z_p)) * maximum(abs.(times)) / (2*pi) ))
  if adaptK
    K = ceil(Int64, alpha*N + 2*m)
    # println("Calculated K for NFFT: ",K)
  end

  win, win_hat = NFFT.getWindow(:kaiser_bessel)

  return win, win_hat, K, N, t_scale, times, z_p
end

function get_A_coefficients_nfft(K::Int64
                                , times::Vector     # readout times
                                , z_p::Vector{Complex{T}}       # Correction map
                                , alpha::Float64    # oversampling factor(NFFT)
                                , m::Float64        # kernel size(NFFT)
                                ; adaptK::Bool=true
                                , centering::Bool=true) where T

  win, win_hat, K, N, t_scale, times, z_p = _precompute_coeffs_nfft(K, times, z_p, alpha, m; adaptK=adaptK,centering=centering)

  A = zeros(Complex{T}, length(times), K)

  for kappa=1:K
    for j=1:length(times)
        argument = times[j]/t_scale - (kappa - K/2)/(K - 2*m)
        A[j,kappa] = win(argument,N*alpha,m,alpha)
    end
  end
  return A
end

function get_C_coefficients_nfft(K::Int64
                                , times::Vector     # readout times
                                , z_p::Vector{Complex{T}}       # Correction map
                                , alpha::Float64    # oversampling factor(NFFT)
                                , m::Float64       # kernel size(NFFT)
                                ; adaptK::Bool=true
                                , centering::Bool=true) where T

  win, win_hat, K, N, t_scale, times, z_p = _precompute_coeffs_nfft(K, times, z_p, alpha, m; adaptK=adaptK,centering=centering)

  C = zeros(Complex{T}, K, length(z_p))

  for p=1:length(z_p)
    for kappa=1:K
        fourierarg = 1im * z_p[p]*t_scale/(2pi)
        exparg = -z_p[p]*t_scale * ( (kappa - K/2)/(K-2*m) )
        C[kappa,p] = (1/ ( (K-2*m)* 1/(N*alpha)*win_hat(fourierarg,N*alpha,m,alpha))) * exp(exparg)
    end
  end
  return C
end

"""
  Returns the a_j,k and c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the nfft method
"""
function get_AC_coefficients_nfft(K::Int64
                                , times::Vector     # readout times
                                , z_p::Vector{Complex{T}}       # Correction map
                                , alpha::Float64   # oversampling factor(NFFT)
                                , m::Float64        # kernel size(NFFT)
                                ; adaptK::Bool=true
                                , centering::Bool=true) where T

  A = get_A_coefficients_nfft(K, times, z_p, alpha, m; adaptK=adaptK, centering=centering)
  C = get_C_coefficients_nfft(K, times, z_p, alpha, m; adaptK=adaptK, centering=centering)
  return A, C
end

"""
  Returns the a_j,k and c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the least squares method and time-segmentation
"""
function get_AC_coefficients_lsqr(K::Int64, times::Vector{T}, z_p::Vector{Complex{T}}) where T
  A = zeros(Complex{T},length(times),K)
  C = zeros(Complex{T},K,length(z_p))
  G = zeros(Complex{T},length(z_p),K)

  t_hat = zeros(K)
  t_hat = collect(range(times[1], stop=times[length(times)], length=K))

  z_p = vec(z_p)
  for kappa=1:K
    for p=1:length(z_p)
      G[p,kappa] = exp(-t_hat[kappa]*z_p[p])
      C[kappa,p] = exp(-t_hat[kappa]*z_p[p])
    end
  end

  systemMat = G'*G
  b = [exp(-times[j]*z_p[p]) for p=1:length(z_p), j=1:length(times)]
  At = systemMat \ (G'*b)
  A .= transpose(At)

  return A,C
end

"""
  Returns the a_j,k Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the least squares method with histogramms and time-segmentation
"""
function get_A_coefficients_hist_lsqr(K::Int64, times::Vector, z_p::Vector{Complex{T}}) where T

  numBins = 10*K
  if real.(z_p)== zeros(length(z_p))
    fmax_im = maximum(imag.(z_p))
    fmin_im = minimum(imag.(z_p))
    δf = (fmax_im-fmin_im)/(numBins-1)
    edges_im = [fmin_im-δf/2 + i*δf for i = 0:numBins]
    h = fit(Histogram, imag.(z_p),edges_im,closed=:left)
    z_center = 0.5*1im*vec( [edges_im[i]+edges_im[i+1] for i=1:numBins] )
    z_weights = h.weights
  elseif imag.(z_p)== zeros(length(z_p))
    fmax_re = maximum(real.(z_p))
    fmin_re = minimum(real.(z_p))
    δf = (fmax_re-fmin_re)/(numBins-1)
    edges_re = collect(fmin_re-δf/2.:δf:fmax_re+δf/2)
    h = fit(Histogram, real.(z_p),edges_re,closed=:left)
    z_center = 0.5*vec( [edges_re[i]+edges_re[i+1] for i=1:numBins] )
    z_weights = h.weights
  else
    nbins = ceil(Int, sqrt(numBins))
    fmax_re = maximum(real.(z_p))
    fmin_re = minimum(real.(z_p))
    δf_re = (fmax_re-fmin_re)/(nbins-1)
    edges_re = collect(fmin_re-δf/2.:δf:fmax_re+δf/2)
    fmax_im = maximum(imag.(z_p))
    fmin_im = minimum(imag.(z_p))
    δf = (fmax_im-fmin_im)/(nbins-1)
    edges_im = collect(fmin_im-δf/2.:δf:fmax_im+δf/2)
    h = fit(Histogram, (real.(z_p),imag.(z_p)),(edges_re,edges_im),closed=:left)
    z_center = 0.5*vec( [edges_re[i]+edges_re[i+1]+1im*(edges_im[j]+edges_im[j+1]) for i=1:numBins[1], j=1:numBins[2] ] )
    z_weights = vec(h.weights)
    numBins = nbins^2
  end

  A = zeros(Complex{T},length(times),K)
  G = zeros(Complex{T},numBins,K)

  t_hat = zeros(K)
  t_hat = collect(range(times[1], stop=times[end], length=K))

  z_p = vec(z_p)
  for p=1:numBins
    for kappa=1:K
      G[p,kappa] = exp(-t_hat[kappa]*z_center[p])
    end
  end

  systemMat = G'*Diagonal(z_weights)*G
  b = [z_weights[p]*exp(-times[j]*z_center[p]) for p=1:numBins, j=1:length(times)]
  At = systemMat \ (G'*b)
  A .= transpose(At)

  return A
end

"""
  Returns the c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the least squares method with histogramms
"""
function get_C_coefficients_hist_lsqr(K::Int64, times::Vector, z_p::Vector{Complex{T}}) where T

  C = zeros(Complex{T},K,length(z_p))
  t_hat = collect(range(times[1],stop=times[end],length=K))
  # z_p = vec(z_p)
  for p=1:length(z_p)
    for kappa=1:K
      C[kappa,p] = exp(-t_hat[kappa]*z_p[p])
    end
  end

  return C
end

"""
  Returns the a_j,k and c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the least squares method and frequency-segmentation
"""
function get_AC_coefficients_lsqr_fs(K::Int64, times::Vector, z_p::Vector{Complex{T}}) where T

  ω = imag.(z_p) # field map
  # ω_hat = collect(range(minimum(ω), stop=maximum(ω),length=K)) # translate frequencies
  ω_hat = loydMax(vec(ω),K, iterations=1000)

  A = zeros(Complex{T},length(times),K)
  G = zeros(Complex{T},length(times),K)

  for kappa=1:K
    for j=1:length(times)
      G[j,kappa] = exp(-1im*times[j]*ω_hat[kappa])
      A[j,kappa] = exp(-1im*times[j]*ω_hat[kappa])
    end
  end

  systemMat = G'*G
  b = [exp(-times[j]*z_p[p]) for j=1:length(times), p=1:length(z_p)]
  C = systemMat \ (G'*b)

  return A,C
end

"""
  Returns the a_j,k Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Using the least squares method with histogramms and frequency-segmentation
"""
function get_A_coefficients_hist_lsqr_fs(K::Int64, times::Vector, z_p::Vector{Complex{T}}) where T
  ω = imag.(z_p) # field map
  # ω_hat = collect(range(minimum(ω), stop=maximum(ω),length=K)) # translate frequencies
  ω_hat = loydMax(vec(ω),K, iterations=1000)

  A = zeros(Complex{T},length(times),K)
  for kappa=1:K
    for j=1:length(times)
      A[j,kappa] = exp(-1im*times[j]*ω_hat[kappa])
    end
  end

  return A
end

function get_A_coefficients_svd(K::Int64, times::Vector, z_p::Vector{Complex{T}}; numBins=10*K, K_tol::T=1.e-4) where T
  ω = imag.(z_p)

  # create histrogram of field map
  fmax = maximum(ω)
  fmin = minimum(ω)
  δf = (fmax-fmin)/(numBins-1)
  edges = [fmin-δf/2 + i*δf for i = 0:numBins]
  h = fit( Histogram, ω, edges, closed=:left )
  z_center = 0.5*1im * ( edges[1:numBins].+edges[2:numBins+1] )
  z_weights = h.weights

  # from reduced exponential matrix
  A = zeros(Complex{T}, length(times), numBins)
  for j=1:numBins, i=1:length(times)
    A[i,j] = sqrt(z_weights[j])*exp(-z_center[j]*times[i])
  end

  # compute SVD
  U,S,V = psvd(A, rank=2*K)
  # estimate accuracy of the approximation and determine rank threshold
  acc = sqrt.(cumsum(S.^2)./sum(S.^2))
  K0 = findfirst(x->x>=(1-K_tol-1e-14), acc)
  @debug "accuracy: $(acc)"
  @info " K_tol=$(K_tol) => chose K= $(K0)"

  return U[:,1:K0]
end

"""
  Returns the c_k,p Coefficients considering the approximation of the
  exponential term exp(-t_j * z_p) = sum( k : a_j,k * c_k,p )

  Uses the LSQR-method for a given set of a_j,k-coefficients
"""
function get_C_coefficients_lsqr(times::Vector, z_p::Vector{Complex{T}}, A::Matrix{Complex{T}}, numSamp::Int64, step=10) where T

  systemMat = adjoint(A[1:step:numSamp,:])* A[1:step:numSamp,:]
  b = [exp(-times[j]*z_p[i]) for j=1:step:numSamp, i=1:length(z_p)]
  C = systemMat \ ( adjoint(A[1:step:numSamp,:])*b )

  return C
end


text:
Base.:(+)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x + DoubleFloat{T}(y[1], y[2])
Base.:(+)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) + y

Base.:(-)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x - DoubleFloat{T}(y[1], y[2])
Base.:(-)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) - y

Base.:(*)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x * DoubleFloat{T}(y[1], y[2])
Base.:(*)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) * y

Base.:(/)(x::DoubleFloat{T}, y::Tuple{T,T}) where {T<:AbstractFloat} = x / DoubleFloat{T}(y[1], y[2])
Base.:(/)(x::Tuple{T,T}, y::DoubleFloat{T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2]) / y

@inline function add_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    (isnan(LO(x)) | isnan(LO(y))) && return add_dbdb_db_nonfinite(x,y)
    return DoubleFloat{T}(add_dddd_dd(HILO(x), HILO(y)))
end

@inline function sub_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    (isnan(LO(x)) | isnan(LO(y))) && return sub_dbdb_db_nonfinite(x,y)
    return DoubleFloat{T}(sub_dddd_dd(HILO(x), HILO(y)))
end

@inline function mul_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    (isnan(LO(x)) | isnan(LO(y))) && return mul_dbdb_db_nonfinite(x,y)
    return DoubleFloat{T}(mul_dddd_dd(HILO(x), HILO(y)))
end

@inline function dvi_dbdb_db(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    (iszero(HI(y)) | isnan(LO(x)) | isnan(LO(y))) && return dvi_dbdb_db_nonfinite(x,y)
    return DoubleFloat{T}(dvi_dddd_dd(HILO(x), HILO(y)))
end

@inline function add_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) + HI(y)
    return DoubleFloat{T}(z, T(NaN))
end

@inline function sub_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) - HI(y)
    return DoubleFloat{T}(z, T(NaN))
end

@inline function mul_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) * HI(y)
    return DoubleFloat{T}(z, T(NaN))
end

@inline function dvi_dbdb_db_nonfinite(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:AbstractFloat}
    z = HI(x) / HI(y)
    return DoubleFloat{T}(z, T(NaN))
end


text:
module Bobby

    export Bitboard
    using Printf
    using Crayons

    include("constants.jl")
    include("converters.jl")
    include("bitboard.jl")
    include("parser.jl")
    include("print.jl")
    include("check.jl")
    include("moves.jl")
    include("unmove.jl")
    include("pawns.jl")
    include("king.jl")
    include("queen.jl")
    include("nights.jl")
    include("rooks.jl")
    include("bishops.jl")
    include("perft.jl")
    include("game.jl")
end


text:
# Here is the code to install packages.
# The list below is an example, none of the packages are actually required.
using Pkg
Pkg.add(
    ["Parameters", "RecursiveArrayTools", "Statistics", "DifferentialEquations", "PyPlot"]
    )


text:
"""
Adapted from sandreza/Learning/sandbox/learn_simple_convection.jl
https://github.com/sandreza/Learning/blob/master/sandbox/learn_simple_convection.jl

This example does NOT use the ProfileData struct.
See example02 for an example that uses ProfileData struct.
"""

using Statistics, LinearAlgebra, Plots

include("kernels.jl")
# include("gaussian_process.jl")
# include("../les/get_les_data.jl")

include("GP1.jl")

const save_figure = false

filename = "general_strat_16_profiles.jld2"
data = get_les_data(filename);

# pick variable to model
V_name = "T"
V = data.T;
# V_name = "wT"
# V = data.wT

t = data.t;
Nz, Nt = size(V);

D = 16 # gridpoints
zavg = custom_avg(data.z, D); # compress z vector to D values
vavg = [custom_avg(V[:,j], D) for j in 1:Nt]; # compress variable array to D values per time
x = vavg[1:(Nt-1)]; # (v₀, v₁, ... ,v_(Nt-1)) (Nt-1)-length array of D-length inputs
y = vavg[2:Nt];     # (v₁, v₂, ... ,v_Nt    ) (Nt-1)-length array of D-length targets

# reserve 25% of data for training, but across the entire time interval
total_set = 1:(Nt-1);
training_set = 1:4:(Nt-1);
verification_set = setdiff(total_set, training_set);

# n_train = length(training_set)
x_train = x[training_set];
y_train = y[training_set];

kernel = Matern12I(1e4,1.0); #[γ,σ]
𝒢 = construct_gpr(x_train, y_train, kernel);

# index_check = 1
# y_prediction = prediction([x_train[index_check]], 𝒢)
# norm(y_prediction - y_train[index_check])

if V_name=="T"; scaling = Tscaling(vavg[1][end]-vavg[1][1]) end # Tscaling(ΔT)
if V_name=="wT"; scaling = wTscaling( maximum(maximum, vavg) ) end # wTscaling(nc)

# get error at each time in the verification set
gpr_error = collect(verification_set)*1.0;
# greedy check
for j in eachindex(verification_set);
    test_index = verification_set[j];
    y_prediction = prediction(x[test_index], 𝒢, scaling);
    println(y_prediction)
    println(y[test_index])
    δ = norm(y_prediction .- y[test_index]);
    gpr_error[j] = δ;
end
histogram(gpr_error)

#error across all verification time steps
println("The mean error is " * string(sum(gpr_error)/length(gpr_error)))
println("The maximum error is " * string(maximum(gpr_error)))

###
# test_index = 100
# gpr_y = prediction([x[test_index]], 𝒢)
# norm(gpr_y - y[test_index])
# scatter(gpr_y,zavg)
###

# the true check
# time evolution given the same initial condition
Nt = length(data.t)
gpr_prediction = similar(y[total_set])
starting = x[1]
gpr_prediction[1] = starting
Nt = length(y[total_set])
for i in 1:(Nt-1)
    println(gpr_prediction[i])
    gpr_prediction[i+1] = prediction(gpr_prediction[i], 𝒢, scaling)
end

animation_set = 1:30:(Nt-1)
anim = @animate for i in animation_set
    exact = V[:,i+1]
    day_string = string(floor(Int, data.t[i]/86400))
    p1 = scatter(gpr_prediction[i+1], zavg, label = "GP")
    plot!(exact,data.z, legend = :topleft, label = "LES", xlabel = "temperature", ylabel = "depth", title = "day " * day_string, xlims = (19,20))
    display(p1)
end
if save_figure == true
    gif(anim, pwd() * "gp_emulator.gif", fps = 15)
    mp4(anim, pwd() * "gp_emulator.mp4", fps = 15)
end

####

"""
Adapted from sandreza/Learning/sandbox/learn_simple_convection.jl
https://github.com/sandreza/Learning/blob/master/sandbox/learn_simple_convection.jl
"""

using JLD2, NetCDF, Statistics, LinearAlgebra, Plots
include("kernels.jl")
include("gaussian_process.jl")
include("../les/get_les_data.jl")

const save_figure = false

filename = "general_strat_16_profiles.jld2"
data = get_les_data(filename);

# pick variable to model
V = data.T
# V = data.wT

t = data.t
Nz, Nt = size(V)

D = 16 # gridpoints
zavg = custom_avg(data.z, D) # compress z vector to D values
vavg = [custom_avg(V[:,j], D) for j in 1:Nt] # compress variable array to D values per time
x = vavg[1:(Nt-1)] # (v₀, v₁, ... ,v_(Nt-1)) (Nt-1)-length array of D-length inputs
y = vavg[2:Nt]     # (v₁, v₂, ... ,v_Nt    ) (Nt-1)-length array of D-length targets

# reserve 25% of data for training, but across the entire time interval
total_set = 1:(Nt-1)
training_set = 1:4:(Nt-1)
verification_set = setdiff(total_set, training_set)

# n_train = length(training_set)
x_train = x[training_set]
y_train = y[training_set]

kernel = SquaredExponentialKernelI(1000.0,1.0) #[γ,σ]
𝒢 = construct_gpr(x_train, y_train, kernel)

index_check = 1
y_prediction = prediction([x_train[index_check]], 𝒢)
norm(y_prediction - y_train[index_check])

##

# get error at each time in the verification set
gpr_error = collect(verification_set)*1.0;
# greedy check
for j in eachindex(verification_set)
    test_index = verification_set[j]
    y_prediction = prediction([x[test_index]], 𝒢)
    δ = norm(y_prediction - y[test_index])
    gpr_error[j] = δ
end
histogram(gpr_error)

#error across all verification time steps
println("The mean error is " * string(sum(gpr_error)/length(gpr_error)))
println("The maximum error is " * string(maximum(gpr_error)))
# the true check
# time evolution given the same initial condition
Nt = length(data.t)
set = 1:(Nt-2)
gpr_prediction = similar(y[total_set])
starting = x[1]
gpr_prediction[1] = starting
Nt = length(y[total_set])
for i in set
    gpr_prediction[i+1] = prediction([gpr_prediction[i]], 𝒢)
end

animation_set = 1:20:(Nt-1)
anim = @animate for i in animation_set
    exact = V[:,i+1]
    day_string = string(floor(Int, data.t[i]/86400))
    p1 = scatter(gpr_prediction[i+1], zavg, label = "GP")
    plot!(exact, data.z, legend = :topleft, label = "LES", xlabel = "temperature", ylabel = "depth", title = "day " * day_string, xlims = (19,20))
    display(p1)
end
if save_figure == true
    gif(anim, "gp_emulator2.gif", fps = 15)
    mp4(anim, "gp_emulator2.mp4", fps = 15)
end

# explore_hyp_grid(σ,γ)

# α_proxy(x) = x[end] - x[end-1]

# 𝒟4 = construct_profile_data("general_strat_4_profiles.jld2", v_str, 16; N=N)
# am4 = α_proxy(𝒟4.x_train[50])
#
# 𝒟8 = construct_profile_data("general_strat_8_profiles.jld2", v_str, 16; N=N)
# am8 = α_proxy(𝒟8.x_train[1])
# am8 = α_proxy(𝒟8.x_train[50])
# am8 = α_proxy(𝒟8.x_train[133])
#
# 𝒟16 = construct_profile_data("general_strat_16_profiles.jld2", v_str, 16; N=N)
# am16 = α_proxy(𝒟16.x_train[1])
#
# 𝒟24 = construct_profile_data("general_strat_24_profiles.jld2", v_str, 16; N=N)
# am24 = α_proxy(𝒟24.x_train[50])
#
# 𝒟32 = construct_profile_data("general_strat_32_profiles.jld2", v_str, 16; N=N)
# am32 = α_proxy(𝒟32.x_train[50])
#
# p = plot(𝒟4.x_train[50], 𝒟4.zavg, legend=false)
# plot!(𝒟8.x_train[50], 𝒟8.zavg)
# plot!(𝒟16.x_train[50], 𝒟16.zavg)
# plot!(𝒟24.x_train[50], 𝒟24.zavg)
# plot!(𝒟32.x_train[50], 𝒟32.zavg)
#
# 𝒟4.x_train[50]

# gpr_prediction = get_gpr_pred(𝒢, 𝒟)

# animation_set = 1:10:(𝒟.Nt-2)
# anim = @animate for i in animation_set
#     exact = 𝒟.v[:,i]
#     day_string = string(floor(Int, 𝒟.t[i]/86400))
#     p1 = scatter(gpr_prediction[i], 𝒟.zavg, label = "GP")
#     # xlims=(minimum(data.v[:,1]),maximum(data.v[:,1]))
#     xlims=(18,20)
#
#     if i<𝒟.Nt
#         exact16 = 𝒟.v[:,i]
#     else
#         exact16 = 𝒟.v[:,data16.Nt]
#     end
#
#     plot!(exact, 𝒟.z, legend = :topleft, label = "LES", xlabel = "$(V_name["T"])", ylabel = "depth", title = "i = $(i)", xlims=xlims)
#     # plot!(exact16, data16.z, legend = :topleft, label = "LES gs 16", xlabel = "$(V_name["T"])", ylabel = "depth", xlims=xlims)
#
# end
#
# gif(anim, pwd() * "ignore.gif", fps = 20)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using MenuAdventures: MenuAdventures
using Documenter: doctest

doctest(MenuAdventures)


text:
using SparseArrays: spdiagm, sparse, spzeros, SparseMatrixCSC
using LinearAlgebra: diag, I

export differential_operator_sequence
"""
    differential_operator_sequence(m::Metric{N}, mesh::Mesh{N, K}, expr::String,
        k::Int, primal::Bool) where {N, K}

Compute the differential operators defined by the string `expr`. This string must consist
of the characters `d`, `★`, `δ`, and `Δ` indicating the exterior derivative, hodge dual,
codifferential, and Laplace-de Rham operators, respectively. If there are two `★` operators
in a row, a single `★★` operator is computed since this avoids unnecessary calculation.
The dimension of the form on which the operator acts is `k-1` and its primality or duality
is indicated by `primal`.
"""
function differential_operator_sequence(m::Metric{N}, mesh::Mesh{N, K}, expr::String,
    k::Int, primal::Bool) where {N, K}
    ops = SparseMatrixCSC{Float64,Int64}[]
    chars = collect(expr)
    i = length(chars)
    while i > 0
        char = chars[i]
        @assert char in ['d', '★', 'δ', 'Δ']
        if char == 'd'
            comp = primal ? mesh.primal.complex : mesh.dual.complex
            push!(ops, exterior_derivative(comp, k))
            k += 1
        elseif char == '★'
            if (i > 1) && chars[i-1] == '★'
                push!(ops, hodge_square(m, mesh, k, primal))
                i -= 1
            else
                push!(ops, circumcenter_hodge(m, mesh, k, primal))
                primal = !primal
                k = K-k+1
            end
        elseif char == 'δ'
            push!(ops, codifferential(m, mesh, k, primal))
            k -= 1
        elseif char == 'Δ'
            push!(ops, laplace_de_Rham(m, mesh, k, primal))
        end
        i -= 1
    end
    return reverse(ops)
end

export differential_operator
"""
    differential_operator(m::Metric{N}, mesh::Mesh{N}, expr::String, k::Int,
        primal::Bool) where N

Compute the differential operator defined by the string `expr`. This string must consist
of the characters `d`, `★`, `δ`, and `Δ` indicating the exterior derivative, hodge dual,
codifferential, and Laplace-de Rham operators, respectively. The dimension of the form on
which the operator acts is `k-1` and its primality or duality is indicated by `primal`.
"""
function differential_operator(m::Metric{N}, mesh::Mesh{N}, expr::String, k::Int,
    primal::Bool) where N
    ops = differential_operator_sequence(m, mesh, expr, k, primal)
    return reduce(*, ops)
end

"""
    differential_operator(m::Metric{N}, mesh::Mesh{N}, expr::String, k::Int,
        primal::Bool, v::AbstractVector{<:Real}) where N

Apply the differential operator defined by the string `expr` to the vector `v`.
"""
function differential_operator(m::Metric{N}, mesh::Mesh{N}, expr::String, k::Int,
    primal::Bool, v::AbstractVector{<:Real}) where N
    ops = differential_operator_sequence(m, mesh, expr, k, primal)
    return foldr(*, ops, init=v)
end

"""
    circumcenter_hodge(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}

Find the discrete hodge star operator using the circumcentric hodge star construction. If
`primal == true`, this operator takes primal `k-1` forms to dual `K-k+1` forms. Otherwise
it takes dual `k-1` forms to priaml `K-k+1` forms.
"""
function circumcenter_hodge(m::Metric{N}, mesh::Mesh{N, K}, k::Int,
    primal::Bool) where {N, K}
    @assert 1 <= k <= K+1
    if k == K+1
        return spzeros(0,0)
    end
    if primal
        ratios = Float64[]
        for (p_cell, d_cell) in zip(mesh.primal.complex.cells[k],
                mesh.dual.complex.cells[K-k+1])
            p_vol = volume(m, mesh.primal, p_cell)
            d_vol = volume(m, mesh.dual, d_cell)
            @assert p_vol != 0
            push!(ratios, d_vol/p_vol)
        end
    else
        pm = hodge_square_sign(m, K, k)
        v = diag(circumcenter_hodge(m, mesh, K-k+1, true))
        @assert !any(v .== 0)
        ratios = pm ./ v
    end
    return spdiagm(0 => ratios)
end

"""
    hodge_square_sign(m::Metric, K::Int, k::Int)

There is an identity `★★ = sign(det(metric)) * (-1)^((k-1) * (K-k)) * I` for `k-1` forms in
`K-1` dimensions. Compute the coefficient of `I` in this expression.
"""
hodge_square_sign(m::Metric, K::Int, k::Int) = sign(det(collect(m.mat))) *
    (mod((k-1) * (K-k), 2) == 0 ? 1 : -1)

"""
    hodge_square(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}

Compute `★★` without computing the `★` operators by using the identity
`★★ = sign(det(metric)) * (-1)^((k-1) * (K-k)) * I` for `k-1` forms in `K-1` dimensions.
"""
function hodge_square(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}
    comp = primal ? mesh.primal.complex : mesh.dual.complex
    n = k < K+1 ? length(comp.cells[k]) : 0
    return hodge_square_sign(m, K, k) * sparse(I, n, n)
end

"""
    exterior_derivative(comp::CellComplex{N, K}, k::Int) where {N, K}

Find the discrete exterior derivative operator.
"""
function exterior_derivative(comp::CellComplex{N, K}, k::Int) where {N, K}
    @assert 0 <= k <= K
    if k == 0
        return spzeros(length(comp.cells[k+1]), 0)
    end
    row_inds, col_inds, vals = Int[], Int[], Int[]
    for (col_ind, cell) in enumerate(comp.cells[k])
        for p in keys(cell.parents)
            o = cell.parents[p]
            row_ind = findfirst(isequal(p), comp.cells[k+1])
            push!(row_inds, row_ind); push!(col_inds, col_ind); push!(vals, 2 * o - 1)
        end
    end
    num_rows = k+1 <= K ? length(comp.cells[k+1]) : 0
    num_cols = length(comp.cells[k])
    return sparse(row_inds, col_inds, vals, num_rows, num_cols)
end

"""
    codifferential(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}

Compute the codifferential defined by
`δ = sign(det(collect(m.mat))) * (-1)^((K-1) * (k-2) + 1) * ★d★` for `k-1` forms in `K-1`
dimensions.
"""
function codifferential(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}
    ★d★ = differential_operator(m, mesh, "★d★", k, primal)
    s = sign(det(collect(m.mat)))
    return s * (mod((K-1) * (k-2) + 1, 2) == 0 ? 1 : -1) * ★d★
end

"""
    laplace_de_Rham(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}

Compute the Laplace-de Rham operator defined by `Δ = dδ + δd`.
"""
function laplace_de_Rham(m::Metric{N}, mesh::Mesh{N, K}, k::Int, primal::Bool) where {N, K}
    dδ = differential_operator(m, mesh, "dδ", k, primal)
    δd = differential_operator(m, mesh, "δd", k, primal)
    return dδ + δd
end

export sharp
"""
    sharp(m::Metric{N}, comp::CellComplex{N}, form::AbstractVector{<:Real}) where N

Given a 1-form on a cell complex, approximate a vector of length `N` at each vertex using
least squares.
"""
function sharp(m::Metric{N}, comp::CellComplex{N}, form::AbstractVector{<:Real}) where N
    field = Vector{Float64}[]
    for c in comp.cells[1]
        mat = zeros(length(c.parents), N)
        w = zeros(length(c.parents))
        for (row_ind, e) in enumerate(collect(keys(c.parents)))
            mat[row_ind, :] = sum([x.points[1].coords * (2 * x.parents[e] - 1)
                for x in e.children])
            w[row_ind] = form[findfirst(isequal(e), comp.cells[2])]
        end
        push!(field, (mat * m.mat) \ w)
    end
    return field
end


text:
# FITSHeader methods

# -----------------------------------------------------------------------------
# Helper functions
#
# Used here and in other files. Functions that operate on FITSFile
# start with `fits_`.

"""
```julia
try_parse_hdrval(T, s) -> x::Union{T,Nothing}
```

attempts to parse string `s` for a FITS card value of type `T` (`String`,
`Bool`, `Int` or `Float64`) and yields either a value of type `T` or
[`nothing`](@ref) if parsing is unsuccessful.

See also: [`tryparse`](@ref), [`parse_header_val`](@ref).

"""
function try_parse_hdrval(::Type{Bool}, s::String)
    if length(s) == 1
        if s[1] == 'T'
            return true
        elseif s[1] == 'F'
            return false
        end
    end
    return nothing
end

# Note that trailing whitespaces are not significant in FITS header
# keywords, but *leading* whitespaces are, so "'    '" parses as " " (a
# single space).  See CFITSIO manual section 4.5 for details.
#
# TODO: parse '' within the string as a single '.
function try_parse_hdrval(::Type{String}, s::String)
    if length(s) < 2 || s[1] != '\'' || s[end] != '\''
        return nothing
    end

    i = lastindex(s) - 1
    while i > 2
        if s[i] != ' '
            return s[2:i]
        end
        i -= 1
    end
    return s[2:i]
end

try_parse_hdrval(::Type{T}, s::String) where {T<:Union{Int,Float64}} =
    tryparse(T, s)

# functions for displaying header values in show(io, header)
hdrval_repr(v::Bool) = v ? "T" : "F"
hdrval_repr(v::String) = @sprintf "'%s'" v
hdrval_repr(v::Union{AbstractFloat, Integer}) = string(v)

"""
```julia
parse_header_val(s) -> x::Union{String,Bool,Int,Float64,Nothing}
```

parses the FITS card value in the string `s` and returns a value of type
`String`, `Bool`, `Int`, `Float64` or `Nothing`.  The latter indicates an
empty value.  This method never throws an error; if the value cannot be
parsed it is returned as it.

See also: [`try_parse_hdrval`](@ref).
"""
function parse_header_val(s::String)
    length(s) == 0 && return nothing

    vb = try_parse_hdrval(Bool, s)
    vb === nothing || return vb

    vs = try_parse_hdrval(String, s)
    vs === nothing || return vs

    vi = try_parse_hdrval(Int, s)
    vi === nothing || return vi

    vf = try_parse_hdrval(Float64, s)
    vf === nothing || return vf

    # FIXME: error("invalid FITS keyword value")
    return s
end

"""
```julia
fits_try_read_keys(f, T, keys)
```

tries to read the raw FITS keys in given order if FITS handle `f` and
returns a value of type `T` or [`nothing`](@ref) if no key exists or if
parsing an existing key is unsuccessful.

See also: [`try_parse_hdrval`](@ref).

"""
function fits_try_read_keys(f::FITSFile, ::Type{T}, keys) where T
    status = Cint[0]
    value = Vector{UInt8}(undef, 71)
    for key in keys
        ccall((:ffgkey, libcfitsio), Cint,
              (Ptr{Cvoid},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{Cint}),
              f.ptr, key, value, C_NULL, status)

        # If the key is found, return it. If there was some other error
        # besides key not found, throw an error.
        if status[1] == 0
            return try_parse_hdrval(T, unsafe_string(pointer(value)))
        elseif status[1] != 202
            error(fits_get_errstatus(status[1]))
        end
    end
    return nothing
end

# Build a string with extension keywords, if present.
# This is a helper function for show(::HDU).
const EXTNAME_KEYS = ["EXTNAME", "HDUNAME"]
const EXTVER_KEYS = ["EXTVER", "HDUVER"]
fits_try_read_extname(f::FITSFile) :: Union{String,Nothing} =
    fits_try_read_keys(f, String, EXTNAME_KEYS)
fits_try_read_extver(f::FITSFile) :: Union{Int,Nothing} =
    fits_try_read_keys(f, Int, EXTVER_KEYS)

function fits_get_ext_info_string(f::FITSFile)
    extname = fits_try_read_extname(f)
    extver = fits_try_read_extver(f)
    if extname === nothing
        return ""
    elseif extver === nothing
        return " (name=$extname)"
    else
        return " (name=$extname, ver=$extver)"
    end
end

# Return indices of reserved keys in a header.
# This is more complex than you would think because some reserved keys
# are only reserved when other keys are present. Also, in general a key
# may appear more than once in a header.
const RESERVED_KEYS = ["SIMPLE","EXTEND","XTENSION","BITPIX","PCOUNT","GCOUNT",
                       "THEAP","EXTNAME","BSCALE","BZERO","BLANK",
                       "ZQUANTIZ","ZDITHER0","ZIMAGE","ZCMPTYPE","ZSIMPLE",
                       "ZTENSION","ZPCOUNT","ZGCOUNT","ZBITPIX","ZEXTEND",
                       "CHECKSUM","DATASUM"]
function reserved_key_indices(hdr::FITSHeader)
    nhdr = length(hdr)
    indices = Int[]
    for i=1:nhdr
        if in(hdr.keys[i], RESERVED_KEYS)
            push!(indices, i)
        end
    end

    # Note that this removes anything matching NAXIS\d regardless of # of axes.
    if in("NAXIS", hdr.keys)
        for i=1:nhdr
            if occursin(r"^NAXIS\d*$", hdr.keys[i])
                push!(indices, i)
            end
        end
    end

    if in("ZNAXIS", hdr.keys)
        for i=1:nhdr
            if (occursin(r"^ZNAXIS\d*$", hdr.keys[i]) ||
                occursin(r"^ZTILE\d*$", hdr.keys[i]) ||
                occursin(r"^ZNAME\d*$", hdr.keys[i]) ||
                occursin(r"^ZVAL\d*$", hdr.keys[i]))
                push!(indices, i)
            end
        end
    end

    if in("TFIELDS", hdr.keys)
        for i=1:nhdr
            for re in [r"^TFORM\d*$", r"^TTYPE\d*$", r"^TDIM\d*$",
                       r"^TUNIT\d*$", r"^TSCAL\d*$", r"^TZERO\d*$",
                       r"^TNULL\d*$", r"^TDISP\d*$", r"^TDMIN\d*$",
                       r"^TDMAX\d*$", r"^TDESC\d*$", r"^TROTA\d*$",
                       r"^TRPIX\d*$", r"^TRVAL\d*$", r"^TDELT\d*$",
                       r"^TCUNI\d*$", r"^TFIELDS$"]
                if occursin(re, hdr.keys[i])
                    push!(indices, i)
                end
            end
        end
    end

    return indices
end

# Write header to CHDU.
# If `clean` is true, skip writing reserved header keywords.
function fits_write_header(f::FITSFile, hdr::FITSHeader, clean::Bool=true)
    indices = clean ? reserved_key_indices(hdr) : Int[]
    for i=1:length(hdr)
        if clean && in(i, indices)
            continue
        end
        if hdr.keys[i] == "COMMENT"
            fits_write_comment(f, hdr.comments[i])
        elseif hdr.keys[i] == "HISTORY"
            fits_write_history(f, hdr.comments[i])
        elseif hdr.comments[i] == ""
            fits_update_key(f, hdr.keys[i], hdr.values[i])
        else
            fits_update_key(f, hdr.keys[i], hdr.values[i], hdr.comments[i])
        end
    end
end

# -----------------------------------------------------------------------------
# Public API

"""
```julia
read_key(hdu, key::String) -> (value, comment)
```

reads the HDU header record specified by keyword and returns a tuple where
`value` is the keyword parsed value (of type `String`, `Bool`, `Int`,
`Float64` or `Nothing`), `comment` is the keyword comment (as a string).
An error is thrown if `key` is not found.

```julia
read_key(hdu, key::Integer) -> (keyname, value, comment)
```

same as above but FITS card is specified by its position and returns a 3
element tuple where `keyname` is the keyword name (a string).

"""
function read_key(hdu::HDU, key::Integer)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    keyout, value, comment = fits_read_keyn(hdu.fitsfile, key)
    keyout, parse_header_val(value), comment
end

function read_key(hdu::HDU, key::String)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    value, comment = fits_read_keyword(hdu.fitsfile, key)
    parse_header_val(value), comment
end

"""
    write_key(hdu, key::String, value[, comment])

Write a keyword value the HDU's header. `value` can be a standard
header type (`String`, `Bool`, `Integer`, `AbstractFloat`) or
`nothing`, in which case the value part of the record will be
empty. If the keyword already exists, the value will be
overwritten. The comment will only be overwritten if given. If the
keyword does not already exist, a new record will be appended at the
end of the header.
"""
function write_key(hdu::HDU, key::String,
                   value::Union{String, Bool, Integer, AbstractFloat, Nothing},
                   comment::Union{String, Ptr{Cvoid}}=C_NULL)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    fits_update_key(hdu.fitsfile, key, value, comment)
end

"""
    read_header(hdu) -> FITSHeader

Read the entire header from the given HDU and return a `FITSHeader` object.
The value of each header record is parsed as `Int`, `Float64`, `String`,
`Bool` or `nothing` according to the FITS standard.

If the value cannot be parsed according to the FITS standard, the value is
stored as the raw unparsed `String`.
"""
function read_header(hdu::HDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)

    # Below, we use a direct call to ffgkyn so that we can keep reusing the
    # same buffers.
    key = Vector{UInt8}(undef, 81)
    value = Vector{UInt8}(undef, 81)
    comment = Vector{UInt8}(undef, 81)
    status = Cint[0]

    nkeys, morekeys = fits_get_hdrspace(hdu.fitsfile)

    # Initialize output arrays
    keys = Vector{String}(undef, nkeys)
    values = Vector{Any}(undef, nkeys)
    comments = Vector{String}(undef, nkeys)
    for i=1:nkeys
        ccall((:ffgkyn,libcfitsio), Cint,
              (Ptr{Cvoid},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{Cint}),
              hdu.fitsfile.ptr, i, key, value, comment, status)
        keys[i] = unsafe_string(pointer(key))
        values[i] = parse_header_val(unsafe_string(pointer(value)))
        comments[i] = unsafe_string(pointer(comment))
    end
    fits_assert_ok(status[1])
    FITSHeader(keys, values, comments)
end

"""
    read_header(hdu, String) -> String

Read the entire header from the given HDU as a single string.
"""
function read_header(hdu::HDU, ::Type{String})
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    fits_hdr2str(hdu.fitsfile)
end

"""
    length(hdr)

Number of records.
"""
length(hdr::FITSHeader) = length(hdr.keys)

"""
    haskey(hdr)

Header keyword exists.
"""
haskey(hdr::FITSHeader, key::String) = in(key, hdr.keys)

"""
    keys(hdr)

Array of keywords (not a copy).
"""
keys(hdr::FITSHeader) = hdr.keys

"""
    values(hdr)

Array of values (not a copy).
"""
values(hdr::FITSHeader) = hdr.values

getindex(hdr::FITSHeader, key::String) = hdr.values[hdr.map[key]]
getindex(hdr::FITSHeader, i::Integer) = hdr.values[i]

function setindex!(hdr::FITSHeader, value::Any, key::String)
    fits_assert_isascii(key)
    if in(key, hdr.keys)
        hdr.values[hdr.map[key]] = value
    else
        push!(hdr.keys, key)
        push!(hdr.values, value)
        push!(hdr.comments, "")
        hdr.map[key] = length(hdr.keys)
    end
end

function setindex!(hdr::FITSHeader, value::Any, i::Integer)
    hdr.values[i] = value
end

# Comments
"""
    get_comment(hdr, key)

Get the comment based on keyword or index.
"""
get_comment(hdr::FITSHeader, key::String) = hdr.comments[hdr.map[key]]
get_comment(hdr::FITSHeader, i::Integer) = hdr.comments[i]

"""
    set_comment!(hdr, key, comment)

Set the comment baed on keyword or index.
"""
function set_comment!(hdr::FITSHeader, key::String, comment::String)
    fits_assert_isascii(comment)
    hdr.comments[hdr.map[key]] = comment
end
function set_comment!(hdr::FITSHeader, i::Integer, comment::String)
    fits_assert_isascii(comment)
    hdr.comments[i] = comment
end

# Display the header
function show(io::IO, hdr::FITSHeader)
    n = length(hdr)
    for i=1:n
        if hdr.keys[i] == "COMMENT" || hdr.keys[i] == "HISTORY"
            @printf io "%s %s" hdr.keys[i] hdr.comments[i][1:min(71, end)]
        else
            @printf io "%-8s" hdr.keys[i]
            if hdr.values[i] === nothing
                print(io, "                      ")
                rc = 50  # remaining characters on line
            else
                val = hdrval_repr(hdr.values[i])
                @printf io "= %20s" val
                rc = length(val) <= 20 ? 50 : 70 - length(val)
            end

            if length(hdr.comments[i]) > 0
                @printf io " / %s" hdr.comments[i][1:min(rc-3, end)]
            end
        end
        i != n && println(io)
    end
end


text:
module UniformlyRefinedForestOfOctreesDiscreteModelsTests
  using MPI
  using Gridap
  using Gridap.Algebra
  using Gridap.FESpaces
  using PartitionedArrays
  using GridapDistributed
  using GridapP4est
  using Test
  using ArgParse
  using P4est_wrapper

  function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--subdomains", "-s"
        help = "Tuple with the # of coarse subdomains per Cartesian direction"
        arg_type = Int64
        default=[1,1]
        nargs='+'
        "--num-uniform-refinements", "-r"
        help = "# of uniform refinements"
        arg_type = Int64
        default=1
    end
    return parse_args(s)
  end

  function run(parts,subdomains,num_uniform_refinements)
    # Manufactured solution
    u(x) = x[1] + x[2]
    f(x) = -Δ(u)(x)

    if length(subdomains)==2
      domain=(0,1,0,1)
    else
      @assert length(subdomains)==3
      domain=(0,1,0,1,0,1)
    end

    coarse_discrete_model=CartesianDiscreteModel(domain,subdomains)
    model=UniformlyRefinedForestOfOctreesDiscreteModel(parts,
                                                       coarse_discrete_model,
                                                       num_uniform_refinements)

    # FE Spaces
    order=1
    reffe = ReferenceFE(lagrangian,Float64,order)
    V = TestFESpace(model,reffe,dirichlet_tags="boundary")
    U = TrialFESpace(u,V)

    trian=Triangulation(model)
    dΩ=Measure(trian,2*(order+1))

    function a(u,v)
      ∫(∇(v)⋅∇(u))dΩ
    end
    function l(v)
      ∫(v*f)dΩ
    end
    dv = get_fe_basis(V)
    du = get_trial_fe_basis(U)
    assem = SparseMatrixAssembler(U,V)

    dof_values = PVector(0.0,V.gids)
    uh = FEFunction(U,dof_values)
    data = collect_cell_matrix_and_vector(U,V,a(du,dv),l(dv),uh)
    A,b = assemble_matrix_and_vector(assem,data)

    x = A\b
    r = A*x -b
    uh = FEFunction(U,x)

    # Error norms and print solution
    trian=Triangulation(model)
    dΩ=Measure(trian,2*order)
    e = u-uh
    e_l2 = sum(∫(e*e)dΩ)
    tol = 1.0e-9
    @test e_l2 < tol
    map_parts(parts) do part
     if (part==1)
       println("$(e_l2) < $(tol)\n")
     end
    end
  end
  if !MPI.Initialized()
    MPI.Init()
  end
  parsed_args = parse_commandline()
  subdomains = Tuple(parsed_args["subdomains"])
  num_uniform_refinements = parsed_args["num-uniform-refinements"]
  parts = get_part_ids(mpi,(prod(subdomains)))
  run(parts,subdomains,num_uniform_refinements)
  MPI.Finalize()
end # module


text:
"""
Abstract evolutionary optimizer algorithm
"""
abstract type AbstractOptimizer end

function print_header(method::AbstractOptimizer)
    println("Iter     Function value")
end
population_size(method::AbstractOptimizer) = error("`population_size` is not implemented for $(summary(method)).")

# Options
"""
Configurable options with defaults:
```
abstol::Float64 = 1e-32
reltol::Float64 = 1e-32
successive_f_tol::Integer = 10
iterations::Integer = 1000
store_trace::Bool = false
show_trace::Bool  = false
show_every::Integer = 1
callback::TCallback = nothing
```
"""
@kwdef struct Options{TCallback <: Union{Nothing, Function}}
    abstol::Float64 = 1e-32
    reltol::Float64 = 1e-32
    successive_f_tol::Integer = 10
    iterations::Integer = 1000
    store_trace::Bool = false
    show_trace::Bool  = false
    show_every::Integer = 1
    callback::TCallback = nothing
end
function show(io::IO, o::Options)
    for k in fieldnames(typeof(o))
        v = getfield(o, k)
        if v === nothing
            print(io, lpad("$(k)",24) *" = nothing\n")
        else
            print(io, lpad("$(k)",24) *" = $v\n")
        end
    end
end

# Optimizer State
"""
Abstract type for defining an optimizer state

Every algorithm have to implement a state type derived from this abstract type.
"""
abstract type AbstractOptimizerState end

"""
    value(state)

Returns a minimum value of the current `state`.
"""
value(state::AbstractOptimizerState) = error("`value` is not implemented for $(state).")

"""
    value(state)

Returns a minimizer object in the current `state`.
"""
minimizer(state::AbstractOptimizerState) = error("`minimizer` is not implemented for $(state).")

"""
    terminate(state)

Returns `true` if the `state` requires early termination.
"""
terminate(state::AbstractOptimizerState) = false

# Optimization Trace

struct OptimizationTraceRecord{T, O <: AbstractOptimizer}
    iteration::Int
    value::T
    metadata::Dict{String,Any}
end

function show(io::IO, t::OptimizationTraceRecord)
    print(io, lpad("$(t.iteration)",6))
    print(io, "   ")
    print(io, lpad("$(t.value)",14))
    for (key, value) in t.metadata
        print(io, "\n * $key: $value")
    end
    return
end

const OptimizationTrace{T,O} = Vector{OptimizationTraceRecord{T,O}}

function show(io::IO, tr::OptimizationTrace)
    print(io, "Iter     Function value\n")
    print(io, "------   --------------\n")
    for rec in tr
        show(io, rec)
        print(io, "\n")
    end
    return
end

# Miscellaneous

# Constructor for nondifferentiable function over bit vector
NonDifferentiable(f, x::AbstractVector{T}) where {T<:Real}=
    NonDifferentiable{T,typeof(x)}(f, zero(T), zeros(T, length(x)), [0,])
NonDifferentiable(f, x::AbstractVector{Bool}) = NonDifferentiable(f, BitVector(x))
function NonDifferentiable(f, x::BitArray)
    xs = BitArray(zero(eltype(x)) for i = 1:length(x))
    NonDifferentiable{Real,typeof(xs)}(f, f(xs), xs, [0,])
end
function NonDifferentiable(f, x::Expr)
    val = f(x)
    NonDifferentiable{typeof(val),typeof(x)}(f, zero(val), :(), [0,])
end

const Individual = Union{AbstractArray, Function, Nothing}

# Evolution strategies

"""Abstract evolution strategy

All evolution strategies must be derived from this type.
"""
abstract type AbstractStrategy end

"""Empty evolution strategy"""
struct NoStrategy <: AbstractStrategy end
copy(s::NoStrategy) = NoStrategy()

"""
Isotropic evolution strategy

This strategy has one mutation parameter for all object parameter components.
"""
mutable struct IsotropicStrategy{T <: Real} <: AbstractStrategy
    σ::T
    τ₀::T
    τ::T
end

"""
    IsotropicStrategy(N)

Returns an isotropic strategy object, which has an one mutation parameter for all object parameter components, with ``\\sigma = 1.0``, ``\\tau_0 = \\sqrt{2N}^{-1}``, ``\\tau = \\sqrt{2\\sqrt{N}}^{-1}``
"""
IsotropicStrategy(N::Integer) = IsotropicStrategy{Float64}(1.0, 1.0/sqrt(2N), 1.0/sqrt(2*sqrt(N)))
copy(s::IsotropicStrategy) = IsotropicStrategy{typeof(s.σ)}(s.σ, s.τ₀, s.τ)

"""
Anisotropic evolution strategy

This strategy has a mutation parameter for each object parameter component.
"""
mutable struct AnisotropicStrategy{T} <: AbstractStrategy
    σ::Vector{T}
    τ₀::T
    τ::T
end

"""
    AnisotropicStrategy(N)

Returns an anisotropic strategy object, which has an one mutation parameter for each object parameter component, with ``\\sigma = [1, \\ldots, 1]^N``, ``\\tau_0 = \\sqrt{2N}^{-1}``, ``\\tau = \\sqrt{2\\sqrt{N}}^{-1}``
"""
AnisotropicStrategy(N::Integer) = AnisotropicStrategy(ones(N), 1/sqrt(2N), 1/sqrt(2*sqrt(N)))
copy(s::AnisotropicStrategy) = AnisotropicStrategy{typeof(s.τ)}(copy(s.σ), s.τ₀, s.τ)


text:
import Distributions: Uniform
import Statistics: std
import StaticArrays: MArray, SArray, MVector, SVector

"""
    addnoise!(pts, joggle_factor)

Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`.
"""
function addnoise!(pts::AbstractArray{T, 2}; joggle_factor = 0.001) where T

    D = minimum(size(pts))
    npts = maximum(size(pts))

    if size(pts, 1) > size(pts, 2)
        dim = 1
    else
        dim = 2
    end

    # Scale standard deviation along each axis by joggle factor
    σ = joggle_factor .* std(pts, dims = dim)

    for i in 1:D
        r = [rand(Uniform(-σ[i], σ[i])) for pt in 1:npts]
        if dim == 1
            pts[:, i] .+= r
        elseif dim == 2
            pts[i, :] .+= r
        end
    end
end

"""
    addnoise!(pts, joggle_factor)

Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`.
"""
function addnoise!(pts::Vector{Vector{T}}; joggle_factor = 1e-8) where T

    D = length(pts[1])
    npts = length(pts)

    # Scale standard deviation along each axis by joggle factor
    σ = joggle_factor

    for i in 1:D
        r = rand(Uniform(-σ, σ))
        pts[i] .+= r
    end
end

import DelayEmbeddings:
	Dataset

function addnoise!(pts::Dataset; joggle_factor = 1e-8)

    D = size(pts, 2)
    npts = length(pts)

    # Scale standard deviation along each axis by joggle factor
    σ = joggle_factor

    Dataset([pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts])
end

"""
    addnoise!(pts, joggle_factor)

Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`.
"""
function addnoise!(pts::Vector{MVector{D, T}}; joggle_factor = 1e-8) where {D, T}

    npts = length(pts)

    # Scale standard deviation along each axis by joggle factor
    σ = joggle_factor

    for i in 1:D
        r = rand(Uniform(-σ, σ))
        pts[i] .+= r
    end
end

"""
    addnoise!(pts, joggle_factor)

Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`.
"""
function addnoise!(pts::Vector{SVector{D, T}}; joggle_factor = 1e-8) where {D, T}

    npts = length(pts)

    # Scale standard deviation along each axis by joggle factor
    σ = joggle_factor

    [pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts]
end


text:
# This is a collection of toys under the Yggdrasil tree for the good <s>kids</s>
# developers.  These utilities can be employed in builder files.

using BinaryBuilder, Pkg, LibGit2
# Note: if you only need `get_addable_spec`, just add the following line to your recipe
# instead of including this file.  This is reexported here for backward compatibility.
using BinaryBuilderBase: get_addable_spec

"""
    should_build_platform(platform) -> Bool

Return whether the tarballs for the given `platform` should be built.

This is useful when the builder has different platform-dependent elements
(sources, script, products, etc...) that make it hard to have a single
`build_tarballs` call.
"""
function should_build_platform(platform)
    # If you need inspiration for how to use this function, look at the builder
    # for Git:
    # https://github.com/JuliaPackaging/Yggdrasil/blob/c3e3c4a96c723306b4da23fc6d05f12995b21ed8/G/Git/build_tarballs.jl#L76-L93

    # Get the list of platforms requested from the command line.  This should be
    # the only argument not prefixed with "--".
    requested_platforms = filter(arg -> !occursin(r"^--.*", arg), ARGS)

    if isone(length(requested_platforms))
        # `requested_platforms` has only one element: the comma-separated list
        # of platform.  We'll run the platform only if it's in the list
        return any(platforms_match.(Ref(platform), split(requested_platforms[1], ",")))
    else
        # `requested_platforms` doesn't have only one element: if its length is
        # zero, no platform has been explicitely passed from the command line
        # and we we'll run all platforms, otherwise we don't know what to do, so
        # let's return false to be safe.
        return iszero(length(requested_platforms))
    end
end

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
export AbstractModel, AbstractBlock
export forward, backward

import Base: show, size

abstract type AbstractModel end
abstract type AbstractBlock <: AbstractModel end

size(x::AbstractModel) = 0

function forward end
function backward end

include("Utils.jl")

include("Device.jl")
include("NN/NN.jl")
include("TN/TN.jl")


text:
# Macros for lazy broadcasting,
# based on @dawbarton  https://discourse.julialang.org/t/19641/20
# and @tkf            https://github.com/JuliaLang/julia/issues/19198#issuecomment-457967851
# and @chethega      https://github.com/JuliaLang/julia/pull/30939

using MacroTools

lazy(::Any) = throw(ArgumentError("function `lazy` exists only for its effect on broadcasting, see the macro @~"))
struct LazyCast{T}
    value::T
end
Broadcast.broadcasted(::typeof(lazy), x) = LazyCast(x)
Broadcast.materialize(x::LazyCast) = x.value

is_call(ex) = isexpr(ex, :call) && !is_dotcall(ex)

""" Check if `ex` is a dot-call expression like `f.(x, y, z)` or `x .+ y .+ z`. """
is_dotcall(ex) = is_dotcall_nonop(ex) || is_dotcall_op(ex)

""" Check if `ex` is an expression like `f.(x, y, z)`. """
is_dotcall_nonop(ex) =
    isexpr(ex, :.) && length(ex.args) == 2 && isexpr(ex.args[2], :tuple)

""" Check if `ex` is an expression like `x .+ y .+ z`. """
function is_dotcall_op(ex)
    ex isa Expr && !isempty(ex.args) || return false
    op = ex.args[1]
    return op isa Symbol && Base.isoperator(op) && startswith(string(op), ".")
end

lazy_expr(x) = x
function lazy_expr(ex::Expr)
    if is_dotcall(ex)
        return bc_expr(ex)
    elseif is_call(ex)
        return app_expr(ex)
    else
        # TODO: Maybe better to support `a ? b : c` etc.? But how?
        return ex
    end
end

function bc_expr(ex::Expr)
    @assert is_dotcall(ex)
    return :($lazy.($(bc_expr_impl(ex))))
end

bc_expr_impl(x) = x
function bc_expr_impl(ex::Expr)
    # walk down chain of dot calls
    if is_dotcall_nonop(ex)
        f = ex.args[1]  # function name
        args = ex.args[2].args
        return Expr(ex.head, lazy_expr(f), Expr(:tuple, bc_expr_impl.(args)...))
    elseif is_dotcall_op(ex)
        f = ex.args[1]  # function name (e.g., `.+`)
        args = ex.args[2:end]
        return Expr(ex.head, lazy_expr(f), bc_expr_impl.(args)...)
    else
        return lazy_expr(ex)
    end
end

function app_expr(ex::Expr)
    @assert is_call(ex)
    return app_expr_impl(ex)
end

app_expr_impl(x) = x
function app_expr_impl(ex::Expr)
    # walk down chain of calls and lazy-ify them
    if is_call(ex)
        return :($applied($(app_expr_impl.(ex.args)...)))
    else
        return lazy_expr(ex)
    end
end

function checkex(ex)
    if @capture(ex, (arg__,) = val_ )
        if arg[2]==:dims
            throw(ArgumentError("@~ is capturing keyword arguments, try with `; dims = $val` instead of a comma"))
        else
            throw(ArgumentError("@~ is probably capturing capturing keyword arguments, try with ; or brackets"))
        end
    end
    if @capture(ex, (arg_,rest__) )
        throw(ArgumentError("@~ is capturing more than one expression, try $name($arg) with brackets"))
    end
    ex
end

"""
    @~ expr

Macro for creating a `Broadcasted` or `Applied` object.  Regular calls
like `f(args...)` inside `expr` are replaced with `applied(f, args...)`.
Dotted-calls like `f(args...)` inside `expr` are replaced with
`broadcasted.(f, args...)`.  Use `LazyArray(@~ expr)` if you need an
array-based interface.

```
julia> @~ A .+ B ./ 2

julia> @~ @. A + B / 2

julia> @~ A * B + C
```
"""
macro ~(ex)
    checkex(ex)
    # Expanding macro here to support, e.g., `@.`
    esc(:($instantiate($(lazy_expr(macroexpand(__module__, ex))))))
end


text:
"""
    `legend(plots, labels; kwargs...)`

Plots a legend for the given plots with the given labels.
`plots` may be a single Plot or a list of Plots.

## Theme
$(ATTRIBUTES)
"""
@recipe(Legend, plots, labels) do scene
    Theme(
        outer_area = IRect(0, 0, 1, 1),
        backgroundcolor = :white,
        strokecolor = RGBA(0.3, 0.3, 0.3, 0.9),
        strokewidth = 1,
        position = (1, 1),
        gap = 20,
        textgap = 15,
        labelwidth = 20,
        padding = 10,
        outerpadding = 10,
        align = (:left, :center),
        rotation = Quaternionf0(0, 0, 0, 1),
        textcolor = :black,
        textsize = 16,
        font = theme(scene, :font),
        markersize = 5,
        linepattern = Point2f0[(0, 0), (1, 0.0)],
        scatterpattern = Point2f0[(0.5, 0.0)],
        raw = true,
        camera = campixel!
    )
end

# We need to override this particular docstring.
"""
    colorlegend(colormap, range)
    colorlegend(plot::Plot)

Creates a colorbar from the given colormap or range, or from the Attributes of the
given Plot.

## Theme
Available attributes and their defaults for ColorLegend{...} are:

  align            (:left, :hcenter)
  backgroundcolor  :white
  camera           AbstractPlotting.campixel!
  font             "Dejavu Sans"
  formatter        AbstractPlotting.Formatters.plain
  labels           AbstractPlotting.Automatic()
  outerpadding     10
  padding          10
  position         (1, 1)
  ranges           AbstractPlotting.Automatic()
  raw              true
  rotation         0.0
  strokecolor      RGBA{Float64}(0.3,0.3,0.3,0.9)
  strokewidth      0.3
  textcolor        :black
  textgap          15
  textsize         16
  width            (20, <the height of the scene> - 10)
"""
@recipe(ColorLegend, colormap, colorrange) do scene
    Theme(
        width = (20, lift(x -> x.widths[2] - 10, pixelarea(scene))),
        backgroundcolor = :white,
        strokecolor = RGBA(0.3, 0.3, 0.3, 0.9 ),
        strokewidth = 0.3,
        position = (1, 1),
        textgap = 15,
        padding = 10,
        outerpadding = 10,
        align = (:left, :hcenter),
        rotation = 0.0,
        textcolor = :black,
        textsize = 16,
        font = theme(scene, :font),
        ranges = automatic,
        labels = automatic,
        formatter = Formatters.plain,
        raw = true,
        camera = campixel!
    )
end

function make_label(scene, plot, labeltext, i, attributes)
    w, gap, msize, tsize, lpattern, mpattern, padding = getindex.(attributes, (
        :labelwidth, :gap, :markersize, :textsize,
        :linepattern, :scatterpattern, :padding
    ))

    num_of_labels = length((attributes[:labels])[])

    _scale(x, w, pad, g, t) = Point2f0(
        pad + (x[1]w),
        pad + floor(t/2) + x[2]w + ((num_of_labels - i) * g)
    )
    scale(args...) = _scale.(args...)

    if isa(plot, Union{Lines, LineSegments})
        linesegments!(
            scene, lift(scale, lpattern, w, padding, gap, tsize),
            color = plot.color, linestyle = plot.linestyle
        )
    else
        scatter!(
            scene, lift(scale, mpattern, w, padding, gap, tsize),
            markersize = msize, color = plot.color, marker = plot.marker
        )
    end
    plot = scene.plots[end]
    translate!(plot, 0, 0, 15) # zindex!
    return
end

outerbox(x::AbstractPlot) = outerbox(x.parent)
outerbox(x::Scene) = pixelarea(x)

convert_argument(::Type{<:Legend}, plots::AbstractVector, labels::AbstractVector{<: AbstractString}) = (plots, labels)

function plot!(plot::Legend)
    @extract plot (plots, labels)
    isempty(plots[]) && return
    N = length(plots[])
    position, color, stroke, strokecolor, padding, opad = getindex.(plot, (
        :position, :backgroundcolor, :strokewidth, :strokecolor, :padding,
        :outerpadding
    ))

    textbuffer = TextBuffer(plot, Point2)

    args = getindex.(plot, (
        :labelwidth, :gap, :textgap, :padding,
        :textsize, :textcolor, :rotation, :align, :font
    ))

    make_label.(plot, plots[], labels[], 1:N, plot)

    map_once(labels, args...) do labels, w, gap, tgap, padding, font...
        start!(textbuffer)
        num_of_labels = length(labels)
        for i = 1:num_of_labels
            yposition = (num_of_labels - i) * gap
            tsize = floor(font[1] / 2) # textsize at position one, half of it since we used centered align
            xy = Point2f0(w + padding + tgap, yposition + tsize + padding)
            push!(
                textbuffer,
                labels[i], xy, textsize = font[1], color = font[2],
                rotation = font[3], align = font[4],
                font = font[5]
            )
        end
        finish!(textbuffer)
        return
    end
    translate!(textbuffer, 0, 0, 15) # zindex!
    bb = boundingbox(plot)
    legendarea = map_once(position, outerbox(plot), opad, args[4:5]..., args[1:3]...) do xy, area, opad, padding, unused...
        mini = minimum(bb)
        wx, wy, _ = widths(bb) .+ mini
        xy = (Vec2f0(xy) .* widths(area))
        w, h = if isfinite(wx) && isfinite(wy)
            round.(Int, (wx, wy)) .+ padding
        else
            (0, 0)
        end
        bb_rect = IRect(xy[1], xy[2], w, h)
        rect = dont_touch(zerorect(area), bb_rect, Vec2f0(opad))
        zindex = 10 # TODO: this is a hack, to put the legend over other things
        # Of course that stops working when others start employing the same hack^^
        translate!(plot, minimum(rect)..., zindex)
        zerorect(rect)
    end
    poly!(plot, legendarea, color = :white, strokecolor = :black, strokewidth = 1)
end

function convert_arguments(::Type{<: ColorLegend}, plot::AbstractPlot)
    (plot[:colormap][], plot[:colorrange][])
end

function calculated_attributes!(::Type{<: ColorLegend}, plot)
    ranges = replace_automatic!(plot, :ranges) do
        lift(default_ticks, plot[:colorrange], automatic, identity)
    end
    replace_automatic!(plot, :labels) do
        lift(default_labels, ranges, plot[:formatter])
    end
end

function plot!(plot::ColorLegend)
    @extract(
        plot,
        (
            padding, textsize, textcolor, align, font, colormap,
            textgap, width, position, outerpadding, ranges, labels
        )
    )
    vertices = Point3f0[(0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 0, 0)]
    mesh = GLNormalUVMesh(
        vertices = copy(vertices),
        faces = GLTriangle[(1, 2, 3), (3, 4, 1)],
        texturecoordinates = UV{Float32}[(0, 1), (0, 0), (0, 0), (0, 1)]
    )

    cmap_node = lift(colormap) do cmap
        c = to_colormap(cmap)
        # TODO cover the case of a 1D colormap explicitely in the shader
        reshape(c, (length(c), 1))
    end
    tio = TextBuffer(plot, Point2)
    rect = lift(
                ranges, labels, textsize, textcolor, align, font,
                textgap, width, padding, outerpadding, position, outerbox(plot)
            ) do ranges, labels, ts, tc, a, font, tg, w, pad, opad, position, area

        start!(tio)
        N = length(ranges)
        for (i, label) in zip(1:N, labels)
            o1 = (i - 1) / (N - 1) # 0 to 1
            pos = Point2f0(w[1] + tg, (ts/2) + o1 * w[2]) .+ pad
            push!(
                tio, label, pos, textsize = ts, color = tc, rotation = 0.0,
                align = a, font = font
            )
        end
        finish!(tio)
        limits = raw_boundingbox(tio)
        bbw, bbh = widths(limits)
        rect = FRect(
            0, 0,
            w[1] + 2pad + bbw + tg,
            2pad + bbh + (ts/2)
        )
        p = (Vec2f0(position) .* widths(area))
        r2 = FRect(p, widths(rect))
        p = p .+ move_from_touch(zerorect(FRect(area)), r2, Vec2f0(opad))
        translate!(plot, p[1], p[2], 0.0)
        rect
    end
    meshnode = lift(width, padding, textsize) do w, pad, ts
        mesh.vertices .= broadcast(vertices) do v
            Point3f0(((Point2f0(v[1], v[2]) .* Point2f0(w)) .+ Point2f0(0, ts/2) .+ pad)..., 0.0)
        end
        mesh
    end
    lvis = lines!(
        plot, rect,
        linewidth = plot[:strokewidth], color = plot[:strokecolor],
    )
    mesh!(plot, meshnode, color = cmap_node, shading = false)
end