text:
function assemble(op::Identity,
        testnfs::AbstractTimeBasisFunction,
        trialfns::AbstractTimeBasisFunction)

    tbf = convolve(testnfs, trialfns)
    has_zero_tail = all(tbf.polys[end].data .== 0)
    # @show has_zero_tail

    T = scalartype(tbf)
    if has_zero_tail
        z = zeros(T, numintervals(tbf)-1)
    else
        z = zeros(T, numfunctions(tbf))
    end

    Δt = timestep(tbf)
    #for i in eachindex(z)
    for i in 1:numintervals(tbf)-1
        p = tbf.polys[i]
        t = (i-1)*Δt
        z[i] = evaluate(p,t)
    end

    for i in numintervals(tbf):length(z)
        p = tbf.polys[end]
        t = (i-1)*Δt
        z[i] = evaluate(p,t)
    end

    return z
end

mutable struct TensorOperator <: Operator
    spatial_factor
    temporal_factor
end

⊗(A::AbstractOperator, B::AbstractOperator) = TensorOperator(A, B)
function scalartype(A::TensorOperator)
    promote_type(
        scalartype(A.spatial_factor),
        scalartype(A.temporal_factor))
end

function allocatestorage(op::TensorOperator, test_functions, trial_functions,
    ::Type{Val{:bandedstorage}}, ::Type{LongDelays{:ignore}},)

    M = numfunctions(spatialbasis(test_functions))
    N = numfunctions(spatialbasis(trial_functions))

    time_basis_function = BEAST.convolve(
        temporalbasis(test_functions),
        temporalbasis(trial_functions))

    space_operator = op.spatial_factor
    A = assemble(space_operator, spatialbasis(test_functions), spatialbasis(trial_functions))

    K0 = ones(M,N)
    bandwidth = numintervals(time_basis_function) - 1
    data = zeros(scalartype(op), bandwidth, M, N)
    maxk1 = bandwidth
    Z = SparseND.Banded3D(K0, data, maxk1)
    return ()->Z, (v,m,n,k)->(Z[m,n,k] += v)
end

function allocatestorage(op::TensorOperator, test_functions, trial_functions,
    ::Type{Val{:bandedstorage}}, ::Type{LongDelays{:compress}},)

    M = numfunctions(spatialbasis(test_functions))
    N = numfunctions(spatialbasis(trial_functions))

    time_basis_function = BEAST.convolve(
        temporalbasis(test_functions),
        temporalbasis(trial_functions))

    space_operator = op.spatial_factor
    A = assemble(space_operator, spatialbasis(test_functions), spatialbasis(trial_functions))

    K0 = ones(Int, M, N)
    bandwidth = numintervals(time_basis_function) - 1
    K1 = ones(Int,M,N) .+ (bandwidth - 1)
    T  = scalartype(op)
    data = zeros(T, bandwidth, M, N)
    tail = zeros(T, M, N)

    Nt = numfunctions(temporalbasis(trial_functions))
    Z = ConvOp(data, K0, K1, tail, Nt)
    function store1(v,m,n,k)
		if Z.k0[m,n] ≤ k ≤ Z.k1[m,n]
			Z.data[k - Z.k0[m,n] + 1,m,n] += v
		elseif k == Z.k1[m,n]+1
			Z.tail[m,n] += v
		end
	end
    return Z, store1
end

# function allocatestorage(op::TensorOperator, test_functions, trial_functions)

#     M = numfunctions(spatialbasis(test_functions))
#     N = numfunctions(spatialbasis(trial_functions))

#     time_basis_function = BEAST.convolve(
#         temporalbasis(test_functions),
#         temporalbasis(trial_functions))

#     tbf = time_basis_function
#     has_zero_tail = all(tbf.polys[end].data .== 0)
#     @show has_zero_tail

#     if has_zero_tail
#         K = numintervals(time_basis_function)-1
#     else
#         speedoflight = 1.0
#         @warn "Assuming speed of light to be equal to 1!"
#         Δt = timestep(tbf)
#         ct, hs = boundingbox(geometry(spatialbasis(trial_functions)).vertices)
#         diam = 2 * sqrt(3) * hs
#         K = ceil(Int, (numintervals(tbf)-1) + diam/speedoflight/Δt)+1
#     end
#     @assert K > 0

#     Z = zeros(M, N, K)
#     return MatrixConvolution(Z), (v,m,n,k)->(Z[m,n,k] += v)
# end

function assemble!(operator::TensorOperator, testfns, trialfns, store,
    threading = Threading{:multi})

    space_operator = operator.spatial_factor
    time_operator  = operator.temporal_factor

    space_testfns = spatialbasis(testfns)
    space_trialfns = spatialbasis(trialfns)

    time_testfns = temporalbasis(testfns)
    time_trialfns = temporalbasis(trialfns)

    zt = assemble(time_operator, time_testfns, time_trialfns)

    # Truncate to a reasonable size in case the tbf has a semi-infinite support
    tbf = convolve(time_testfns, time_trialfns)
    has_zero_tail = all(tbf.polys[end].data .== 0)
    if !has_zero_tail
        speedoflight = 1.0
        @warn "Assuming speed of light to be equal to 1!"
        Δt = timestep(tbf)
        ct, hs = boundingbox(geometry(space_trialfns).vertices)
        diam = 2 * sqrt(3) * hs
        kmax = ceil(Int, (numintervals(tbf)-1) + diam/speedoflight/Δt)+1
        zt = zt[1:kmax]
    end

    function store1(v,m,n)
        for (k,w) in enumerate(zt)
            store(w*v,m,n,k)
        end
    end
    assemble!(space_operator, space_testfns, space_trialfns, store1)

end

mutable struct TemporalDifferentiation <: Operator
    operator
end

derive(op::AbstractOperator) = TemporalDifferentiation(op)
scalartype(op::TemporalDifferentiation) = scalartype(op.operator)
Base.:*(a::Number, op::TemporalDifferentiation) = TemporalDifferentiation(a * op.operator)

function allocatestorage(op::TemporalDifferentiation, testfns, trialfns,
	storage_trait, longdelays_trait)

	trial_time_fns  = temporalbasis(trialfns)
	trial_space_fns = spatialbasis(trialfns)

	trialfns = SpaceTimeBasis(
		trial_space_fns,
		derive(trial_time_fns)
	)

    return allocatestorage(op.operator, testfns, trialfns, storage_trait, longdelays_trait)
end

function assemble!(operator::TemporalDifferentiation, testfns, trialfns, store, threading = Threading{:multi})

    trial_time_fns  = temporalbasis(trialfns)
    trial_space_fns = spatialbasis(trialfns)

    trialfns = SpaceTimeBasis(
        trial_space_fns,
        derive(trial_time_fns)
    )

    assemble!(operator.operator, testfns, trialfns, store, threading)

end

struct TemporalIntegration <: AbstractOperator
    operator::AbstractOperator
end

integrate(op::AbstractOperator) = TemporalIntegration(op)
derive(op::TemporalIntegration) = op.operator
scalartype(op::TemporalIntegration) = scalartype(op.operator)
Base.:*(a::Number, op::TemporalIntegration) = TemporalIntegration(a * op.operator)

function allocatestorage(op::TemporalIntegration, testfns, trialfns,
	storage_trait, longdelays_trait)

	trial_time_fns  = temporalbasis(trialfns)
	trial_space_fns = spatialbasis(trialfns)

	trialfns = SpaceTimeBasis(
		trial_space_fns,
		integrate(trial_time_fns)
	)

	return allocatestorage(op.operator, testfns, trialfns, storage_trait, longdelays_trait)
end

function assemble!(operator::TemporalIntegration, testfns, trialfns, store,
    threading = Threading{:multi})

    trial_time_fns  = temporalbasis(trialfns)
    trial_space_fns = spatialbasis(trialfns)

    trialfns = SpaceTimeBasis(
        trial_space_fns,
        integrate(trial_time_fns)
    )

    assemble!(operator.operator, testfns, trialfns, store)

end


text:
"""
    chord_length(array, phase[; directions,] nbins = 10)

Calculate chord length correlation function for one-, two- or
three-dimensional multiphase system.

Cord length function `p(x)` equals to probability of finding a chord
whose length is in the range `[x, x+δx]` and which lies entirely in the
phase `phase`. A chord is a line segment which touches the boundary of
a same-phase cluster with its ends.

This implementation bins chord lengths into `nbins` bins and returns
normalized histogram on collected data and the mean chord length in a
tuple.

# Examples
```jldoctest
julia> chord_length([1, 0, 0, 0, 0, 1], 0)
(StatsBase.Histogram{Float64,1,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}
edges:
  4.0:1.0:5.0
weights: [1.0]
closed: left
isdensity: false, 4.0)
```

For a list of possible dimensions, see also: [`direction1Dp`](@ref),
[`direction2Dp`](@ref), [`direction3Dp`](@ref).
"""
function chord_length(array      :: AbstractArray,
                      phase;
                      directions :: Vector{Symbol} = array |> default_directions,
                      nbins      :: Integer       = 10)
    # Select needed phase by applying indicator function to array
    ph = map(x -> x == phase, array)

    # Find a new "bogus" space for the edge
    edge_phase = maximum(array) + 1

    # Arary of chord lengths
    lengths = Vector{Int}(undef, 0)

    # Poor man's edge detection
    dist = ph |> feature_transform |> distance_transform
    edge = edge_phase * map(x -> x == 1.0, dist)

    # Combine the edge with the original picture
    combo = min.(edge + array, edge_phase)
    for direction in directions
        slicer = slice_generators(combo, false, Val(direction))
        for slice in slicer
            len = 0
            startonedge = false
            for x in slice
                # Edge found
                if x == edge_phase
                    len > 1 && startonedge && push!(lengths, len - 1)
                    len = 0
                    startonedge = true
                # Not our phase
                elseif x != phase
                    startonedge = false
                end
                len += 1
            end
        end
    end

    hist = fit(Histogram, lengths; nbins = nbins)
    mean_len = sum(lengths) ./ length(lengths)
    return normalize(hist; mode = :probability), mean_len
end

ext: jl


text:
@enum RxnType stick=1 arrhenius=2

"""
An anstract type for surface reaction mechanism parameter definitions
"""
abstract type SurfaceMechDefinition end
export SurfaceMechDefinition
abstract type ReactionState end
export ReactionState
"""
composite type for the definition of site definition of a surface reaction mechanism
"""
struct SiteInfo <: SurfaceMechDefinition
    name::String    # name of the surface reaction mechanism
    density::Float64    # surface site density in mol/cm2
    site_coordination::Array{Float64,1} # site coordination of adsorbed species by default 1 for all
    ini_covg::Array{Float64,1} # intitial surface coverage as read from the mechanism
end
export  SiteInfo
mutable struct Parameters <: SurfaceMechDefinition
    s::Float64  # stickig coefficient
    k0::Float64 # pre-exponential factor cm-mol-s
    β::Float64  # temperature exponent
    E::Float64  # Activation energy in J/mol
end
export  Parameters
struct Stoichiometry    <: SurfaceMechDefinition
    reversible::Bool    #whether the reaction is reversible or not
    reactant_ids::Array{Int64,1} # id of reactant species participating in an elementary reaction
    product_ids::Array{Int64,1} #id of product species prarticipating in an elementary reaction
    cov_dep::Dict{Int64,Float64} #coverage dependency of the reaction: species_id=>cov dependent activation energy
    order_dep::Dict{Int64,Float64} #order dependency of the reaction: species_id=>order w.r.t the species
end
export Stoichiometry
struct ElementaryReactions <: SurfaceMechDefinition
    id::Int64 #reaction number
    type::RxnType #reaction type stick or arrhenius
    params::Parameters #composite type Parameters
    mwc::Bool #whether MWC is applied or not
    rxn_stoic::Stoichiometry #composite type Stoichiometry
end
export ElementaryReactions

struct SurfaceMechanism <: SurfaceMechDefinition
    energy_factor::Float64 # Conversion factor for activation energy to SI units
    si::SiteInfo # Composite type SiteInfo
    species::Array{String,1} #list of surface species present in the mechanism
    reactions::Array{ElementaryReactions,1} #Array of composite type ElementaryReactions
end
export  SurfaceMechanism

struct SpeciesRxnMap <: SurfaceMechDefinition
    species_id::Int64 #id of the species
    rxns::Array{Int64,1}    #id of reactions in which the species participates
    stoic_coeff::Array{Int64,1} #stoichiometric coefficient of the species in the reaction
end
export SpeciesRxnMap

struct MechanismDefinition <: SurfaceMechDefinition
    sm::SurfaceMechanism
    srm_array::Array{SpeciesRxnMap,1}
end

mutable struct State <: ReactionState
    mole_frac::Array{Float64,1}
    covg::Array{Float64,1}
    T::Float64
    p::Float64
    s_rate::Array{Float64,1}
end
export State

"""
Separate the reactants and products from the stoichiometric equation
#   Usage:
    parase_rxn_string(rxn_str)
-   rxn_str::String:    String representing the reaction equation
"""
function parase_rxn_string(rxn_str::String)
    #check whether the reaction is reversible
    reversible = false
    m = match(r"<=>|=>|=|>",rxn_str)
    if m.match == "=" || m.match == "<=>"
        reversible = true
    end
    species = collect(split(rxn_str,m.match))
    r_species = [String(strip(uppercase(i))) for i in split(species[1],"+")]
    p_species = [String(strip(uppercase(i))) for i in split(species[2],"+")]
    return (reversible,r_species,p_species)
end

function parse_rxn_params(params::String)
    k = β = E = 0.0
    all_params = split(params)
    k = parse(Float64,all_params[1])
    if length(all_params) == 3
        β = parse(Float64,all_params[2])
        E = parse(Float64,all_params[3])
    elseif length(all_params)==2
        E = parse(Float64,all_params[3])
    end
    return k, β, E
end

"""
Function to determine the coverage dependency, order dependency and MWC
for a reaction
#   Usage
    get_dependencies(cov_dep_rxns,order_dep_rxns,mwc_rxns,rxn_id)
-   cov_dep_rxns::Dict{Int64=>Dict{Int64=>Float64}}
-   order_dep_rxns::Dict{Int64=>Dict{Int64=>Float64}}
-   mwc_rxns::Array{Int64}
-   rxn_id::Int64
The function returns a tuple of cov,order and mwc
"""
function get_dependencies(cov_dep_rxns,order_dep_rxns,mwc_rxns,rxn_id)
    #Check for coverage dependency
    cov = Dict{Int64,Float64}()
    ord = Dict{Int64,Float64}()
    if haskey(cov_dep_rxns,rxn_id)
        cov = cov_dep_rxns[rxn_id]
    end
    #Check for order dependency
    if haskey(order_dep_rxns,rxn_id)
        ord = order_dep_rxns[rxn_id]
    end
    #Check for mwc
    mwc = in(rxn_id,mwc_rxns)
    return cov,ord,mwc

end


text:
mutable struct RhoVsFixedthreshold <: AbstractModel
    # params::Vector{Float64}
    params::NamedTuple{(:threshold,:patch_aplha,:beta,:reset,:bias), NTuple{5,Float64}}
end

function init(::Type{WaltonForaging.RhoVsFixedthreshold})
    return RhoVsFixedthreshold((
        threshold = rand(0.0:0.0001:0.1), #=alpha_environment = [0,0.1] learning rate of environment=#
        patch_aplha = rand(0.0:0.0001:0.2), #=alpha_patch = [0, 0.2] learning rate of the patch=#
        beta = rand(2:0.0001:10), #=beta = [2, 10] noise parameter or inverse temperature=#
        reset = rand(0.0:0.0001:0.5), #=reset = [0, 0.5] patch env rew rate at the entrance of a trials, because it has decayed in the last trial=#
        bias = 0.1 #=bias = [2, 10] bias towards staying in the softmax=#
        ))
end

Distributions.params(m::RhoVsFixedthreshold) = collect(m.params)

function updateparams!(m::RhoVsFixedthreshold,v::Vector{Float64})
    (length(v) == 5) || error("Incorrect number of params for $(typeof(m))")
    threshold,patch_aplha,beta,reset,bias = tuple(v...)
    m.params = (;threshold,patch_aplha,beta,reset,bias)
end

function Pstay(m::RhoVsFixedthreshold, patch_rew_rate::Float64)
    threshold, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    return exp(bias + beta*patch_rew_rate)/(exp(bias + beta*patch_rew_rate) + exp(beta*threshold))
end

function Pleave(m::RhoVsFixedthreshold, patch_rew_rate::Float64)
    threshold, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    return exp(beta*threshold)/(exp(beta*threshold) + exp(bias + beta*patch_rew_rate))
end

function Poutcome(m::RhoVsFixedthreshold, threshold::Float64, patch_rew_rate::Float64, leave::Int64)
    if leave == 1
        Pleave(m, threshold, patch_rew_rate)
    elseif leave == 0
        Pstay(m, threshold, patch_rew_rate)
    else
        error("Wrong argument for Poutcome calculation: $leave")
    end
end

function Likelihood(m::RhoVsFixedthreshold,df)
    # if I start using subdf for patches conditions I will need:
    # an initial env reward rate before splitting per session (maybe a separate variable to feed to loglikelihood)
    # how do I minimize the same patch reset value for the whole session?

    threshold, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    final_env_rho = sum(df.REWARD/df[end,:TIME])
    transform!(df, [:REWARD,:BIN,:NEWTRIAL] => ((r,b,n) -> patch_get_rew_rate(r, b, reset_val, n, patch_alpha)) => :PatchRewRate)
    transform!(df, [:PatchRewRate,:LEAVE] => ByRow((p,l) -> Poutcome(m,p,l)) => :Probability)
    return filter(r -> r.SIDE == "travel", df).Probability
end


text:
using LinearAlgebra: norm, dot

export nsoli

"""
function nsoli(x,f; kwargs...) -> (sol, it_hist, ierr, x_hist)

Newton-Krylov solver, globally convergent solver for f(x) = 0, using
Inexact-Newton-Armijo iteration, Eisenstat-Walker forcing term
and Parabolic line search via three point interpolation.

This code has been transcribed with permission from
http://www4.ncsu.edu/~ctk/newton/SOLVERS/nsoli.m
C. T. Kelley, April 27, 2001

### Required Input Parameters

* x: initial iterate
* f: function for which we are seeking a root

### Output Parameters

* sol = solution
* it_hist(maxit,3) = l2 norms of nonlinear residuals
           for the iteration, number of function evaluations,
           and number of steplength reductions
* ierr: 0 upon successful termination; 1 if after maxit iterations
           the termination criterion is not satsified; 2 failure in the line search. The iteration
           is terminated if too many steplength reductions
           are taken.
* x_hist: matrix of the entire interation history.
           The columns are the nonlinear iterates. This
           is useful for making movies, for example, but
           can consume way too much storage. This is an
           OPTIONAL argument. Storage is only allocated
           if x_hist is in the output argument list.

### Keyword Parameters

* atol: absolute error tolerance for the nonlinear iteration
* rtol: relative error tolerance for the nonlinear iteration
* lmeth: choice of linear iterative method
        1 (GMRES, no restarts), 2 GMRES(m), 3 (BICGSTAB), 4 (TFQMR);
     default = 1
* restart_limit: max number of restarts for GMRES if lmeth = 2, default = 20
* maxit: maxmium number of nonlinear iterations, default = 40
* maxitl: maximum number of inner iterations before restart in GMRES(m), m = maxitl, default = 40. For iterative methods other than GMRES(m) maxitl is the upper bound on linear iterations.
* |etamax| : Maximum error tolerance for residual in inner
     iteration. The inner iteration terminates
     when the relative linear residual is
     smaller than eta*| F(x_c) |. eta is determined
     by the modified Eisenstat-Walker formula if etamax > 0.
     If etamax < 0, then eta = |etamax| for the entire
     iteration. default: etamax = .9
* debug = turns on/off iteration statistics display as
               the iteration progresses
* alpha = 1.d-4, parameter to measure sufficient decrease
* sigma0 = .1, sigma1 = .5, safeguarding bounds for the linesearch
* maxarm = 20, maximum number of steplength reductions before
                    failure is reported
* debug = 0:  Set the debug parameter; 1 turns display on, otherwise off.

SOLVING NONLINEAR EQUATIONS WITH NEWTON'S METHOD
 C. T. KELLEY
 2003 / xiv + 104 pages / Softcover / ISBN: 978-0-898715-46-0
"""
function nsoli(x::Vector{T}, f;
   atol = 1e-5, rtol = 1e-5, dkrylov = dgmres,
   maxit = 40, lmaxit = 40, etamax = 0.9, reorth = 1,
   alpha = 1e-4, sigma0 = 0.1, sigma1 = 0.5, maxarm = 20, gamma = 0.9,
   debug = 0 ) where {T}

   # there was a parameter restart_limit = 20,  >>>> move into defn of dkrylov

   # Initialize some variables
   ierr = 0
   it_histx = zeros(maxit,3)
   x_hist = Vector{T}[x]

   # Initialize parameters for the iterative methods, Check for optional inputs.
   it_histx = zeros(maxit+1, 3)
   # gmparms = [abs(etamax), lmaxit, parms(5), 1]; TODO: remove this
   eta = abs(etamax)
   n = length(x)
   itc = 0          # iteration counter

   # Evaluate f at the initial iterate, and compute the stop tolerance.
   f0 = f(x)                             # TODO: count f evaluations instead of iterations?
   fnrm = norm(f0)                       # TODO: switch to Inf-norm???
   it_histx[itc+1, :] = [fnrm, 0, 0]

   fnrmo = 1.0      # old fnrm
   stop_tol = atol + rtol * fnrm

   # TODO: replace this with a nice type collecting iteration information
   outstat = zeros(1, 5)
   outstat[itc+1, :] = [itc, fnrm, 0, 0, 0]

   # main iteration loop
   while fnrm > stop_tol && itc < maxit
      # Keep track of the ratio (rat = fnrm/frnmo) of successive residual norms
      # and the iteration counter (itc).
      rat = fnrm / fnrmo
      fnrmo = fnrm
      itc = itc + 1

      # compute the Newton direction # TODO: is zero the best initialisation?
      step, errstep, inner_it_count, inner_f_evals =
            dkrylov(f0, f, x, eta, lmaxit, reorth)

      # ~~~~~~~~~~~~~~ BEGINNING OF LINESEARCH ~~~~~~~~~~~~~~~~
      xold = x;
      lambda = 1
      lamm = 1
      lamc = lambda
      iarm = 0
      xt = x + lambda * step
      ft = f(xt)
      nft = norm(ft)
      nf0 = norm(f0)
      ff0 = nf0 * nf0
      ffc = nft * nft
      ffm = nft * nft
      while nft >= (1 - alpha * lambda) * nf0
         # Apply the three point parabolic model.
         if iarm == 0
            lambda = sigma1 * lambda
         else
            lambda = parab3p( lamc, lamm, ff0, ffc, ffm;
                              sigma0 = sigma0, sigma1 = sigma1 )
         end

         # Update x; keep the books on lambda.
         xt = x + lambda * step
         lamm = lamc
         lamc = lambda

         # Keep the books on the function norms.
         ft = f(xt)
         nft = norm(ft)
         ffm = ffc
         ffc = nft*nft
         iarm = iarm+1
         if iarm > maxarm
            @warn(" Armijo failure, too many reductions ")
            ierr = 2
            display(outstat)
            it_hist = it_histx[1:itc+1, :]
            push!(x_hist, x)
            sol = xold
            return sol, it_hist, ierr, x_hist
         end
      end
      x = xt
      f0 = ft
      # ~~~~~~~~~~~~~~~ END OF LINESEARCH ~~~~~~~~~~~~~~~~

      push!(x_hist, x)
      fnrm = norm(f0)
      it_histx[itc+1, 1] = fnrm     # TODO: there is a bug here when we reach too many iterations

      # How many function evaluations did this iteration require?
      it_histx[itc+1, 2] = it_histx[itc, 2] + inner_f_evals + iarm + 1
      if itc == 1
         it_histx[itc+1, 2] = it_histx[itc+1, 2] + 1
      end
      it_histx[itc+1, 3] = iarm

      rat = fnrm/fnrmo
      # Adjust eta as per Eisenstat-Walker.   # TODO: make this a flag!
      if etamax > 0
         etaold = eta
         etanew = gamma * rat^2
         if gamma * etaold^2 > 0.1
            etanew = max(etanew, gamma * etaold^2)
         end
         eta = min(etanew, etamax)
         eta = max(eta, 0.5 * stop_tol / fnrm)
      end

      outstat = [outstat; [itc fnrm inner_it_count rat iarm]]
      #  outstat[itc+1, :] = [itc fnrm inner_it_count rat iarm]
   end
   sol = x
   it_hist = it_histx[1:itc+1, :]
   if debug == 1
      @info("Debugging information:")
      display(outstat)
      it_hist = it_histx[1:itc+1, :]
   end

   # on failure, set the error flag
   if fnrm > stop_tol
      ierr = 1
   end

   return sol, it_hist, ierr, x_hist
end


text:
module TestDeprecated

using Test
using Transducers

@testset "|>" begin
    @test_deprecated Map(identity) |> Map(identity)
end

@testset "reduce" begin
    @test_deprecated reduce(+, Map(identity), 1:3)
    @test_deprecated reduce(+, 1:3 |> Map(identity))
end

@testset "dreduce" begin
    @test_deprecated dreduce(+, Map(identity), 1:3)
    @test_deprecated dreduce(+, 1:3 |> Map(identity))
end

end  # module


text:
push!(LOAD_PATH,"../src/")

using Documenter, DataDrivenDiffEq
using Flux, SymbolicRegression

ENV["GKSwstype"] = "100"

makedocs(
    sitename="DataDrivenDiffEq.jl",
    authors="Julius Martensen, Christopher Rackauckas",
    modules=[DataDrivenDiffEq],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://datadriven.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Tutorials" => Any[
          "examples/linear_systems.md",
          "examples/nonlinear_systems.md",
          "examples/implicit_systems.md",
          "examples/symbolic_regression.md"
        ],
        "Unifying SINDy and DMD" => "sindy_dmd.md",
        "Problems" => "problems.md",
        "Basis" => "basis.md",
        "Solvers" => Any[
          "solvers/koopman.md",
          "solvers/optimization.md",
          "solvers/symbolic_regression.md"
        ],
        "Solutions" => "solutions.md",
        "Utilities" => "utils.md",
        "Contributing" => "contributions.md",
        "Citing" => "citations.md"
        ]
)

deploydocs(
   repo = "github.com/SciML/DataDrivenDiffEq.jl.git";
   push_preview = true
)


text:
"""
    plotcontour_old(shape, plotnormals = false, noutw=false)

Plot the contour of a `ExtractedTranslational`.
`noutw==false` means that the normals point outwards, otherwise they point where the were fitted.
"""
function plotcontour_old(shape, plotnormals = false, noutw=false)
    segs = shape.contour

    sct = lines(segs, scale_plot=false)
    if plotnormals
        tms = [midpoint(segs, i) for i in eachindex(segs)]
        if noutw
            tns = [outwardsnormal(shape, i) for i in eachindex(segs)]
        else
            tns = [contournormal(shape, i) for i in eachindex(segs)]
        end
        scatter!(sct, tms, color=:red, scale_plot=false)
        x_ = (x->x[1]).(tms)
        y_ = (x->x[2]).(tms)
        u_ = (x->x[1]).(tns)
        v_ = (x->x[2]).(tns)
        arrows!(sct, x_, y_, u_, v_, scale_plot=false)
    end
    sct
end

"""
    plotcontour(shape, plotnormals = false, noutw=false)

Plot the contour of a `ExtractedTranslational`.
`noutw==true` means that the normals point outwards, otherwise they point where the were fitted.
"""
function plotcontour(shape, plotnormals = false, noutw=false)
    function interp(a, b, n)
        r = collect(range(0, stop=1, length=n+2))
        deleteat!(r, 1)
        deleteat!(r, lastindex(r))
        return [(1-i)*a+i*b for i in r]
    end

    segs = shape.contour

    sct = lines(segs, scale_plot=false)

    if plotnormals
        tms = [midpoint(segs, i) for i in eachindex(segs)]
        midtms = [interp(segs[i], segs[i+1], 1) for i in 1:size(segs,1)-1]
        tmss = vcat(tms, midtms...)
        if noutw
            tns = [dn2shape_outw(tmss[i],shape)[2] for i in eachindex(tmss)]
        else
            tns = [dn2shape_contour(tmss[i],shape)[2] for i in eachindex(tmss)]
        end
        scatter!(sct, tms, color=:red, scale_plot=false)
        x_ = (x->x[1]).(tmss)
        y_ = (x->x[2]).(tmss)
        u_ = (x->x[1]).(tns)
        v_ = (x->x[2]).(tns)
        arrows!(sct, x_, y_, u_, v_, scale_plot=false)
    end
    sct
end

function compandcenter(shap)
    c = shap.center
    segs = shap.contour
    mi1 = []
    mi2 = []
    cent1 = []
    cent2 = []
    cn1 = []
    cn2 = []
    for i in 1:size(segs,1)
        mi = midpoint(segs, i)
        push!(mi1, mi[1])
        push!(mi2, mi[2])

        to_cent = normalize(c-midpoint(segs, i))
        push!(cent1, to_cent[1])
        push!(cent2, to_cent[2])

        c_n = segmentnormal(segs, i)
        push!(cn1, c_n[1])
        push!(cn2, c_n[2])
    end
    sc = Scene()
    arrows!(sc, mi1, mi2, cent1, cent2, arrowcolor=:red, scale_plot=false)
    arrows!(sc, mi1, mi2, cn1, cn2, arrowcolor=:blue, scale_plot=false)
    sc
end

function translheatmap!(sc, shap, sidelength, dc=6; kwargs...)
    cont = shap.contour
    aabb = findAABB(cont)
    midp = sum(aabb)/2
    sidel = aabb[2]-aabb[1]

    xr = range(aabb[1][1]-sidel[1]/4, stop=aabb[2][1]+sidel[1]/4, length=sidelength)
    yr = range(aabb[1][2]-sidel[2]/4, stop=aabb[2][2]+sidel[2]/4, length=sidelength)
    if dc == 1
        mm = [RANSAC.impldistance2segment([i,j], shap)[1] for i in xr, j in yr]
    elseif dc == 2
        mm = [RANSAC.dist2segment([i,j], cont)[1] for i in xr, j in yr]
    elseif dc == 3
        mm = [RANSAC.impldistance2segment2([i,j], shap)[1] for i in xr, j in yr]
    elseif dc == 4
        mm = [RANSAC.impldistance2segment3([i,j], shap)[1] for i in xr, j in yr]
    elseif dc == 5
        mm = [RANSAC._dn2shape_outwards([i,j], shap)[1] for i in xr, j in yr]
    elseif dc == 6
        mm = [dn2shape_outw([i,j], shap)[1] for i in xr, j in yr]
    end
    lines!(sc, cont, linewidth=1.5; kwargs...)
    heatmap!(sc, collect(xr), collect(yr), mm; kwargs...)
    leg = colorlegend(sc[end], raw=true, camera=campixel!)
    vbox(leg, sc)
end

function translheatmap(shap, sidelength, dc=1; kwargs...)
    sc = Scene()
    return translheatmap!(sc, shap, sidelength, dc; kwargs...)
end


text:
@inline nfill(v, n::Val) = ntuple(i -> v, n)

# cartesian
@inline block_index(child_dims::NTuple{N, Integer}, I::Vararg{Int, N}) where {N} = @. div(I-1, child_dims) + 1
@inline function block_local_index(child_dims::NTuple{N, Integer}, I::Vararg{Int, N}) where {N}
    blockindex = block_index(child_dims, I...)
    localindex = @. I - (blockindex-1) * child_dims
    blockindex, localindex
end

# getindex like function
for f in (:unsafe_getindex, :isactive, :allocate!, :isallocated)
    _f = Symbol(:_, f)
    @eval begin
        function $f(A::AbstractArray, I...)
            @_propagate_inbounds_meta
            Base.error_if_canonical_getindex(IndexStyle(A), A, I...)
            $_f(IndexStyle(A), A, to_indices(A, I)...)
        end
        function $_f(::IndexLinear, A::AbstractArray, I::Int...)
            @_inline_meta
            @boundscheck checkbounds(A, I...) # generally _to_linear_index requires bounds checking
            @inbounds $f(A, Base._to_linear_index(A, I...))
        end
        function $_f(::IndexCartesian, A::AbstractArray, I::Int)
            @_inline_meta
            @boundscheck checkbounds(A, I...) # generally _to_subscript_indices requires bounds checking
            @inbounds $f(A, Base._to_subscript_indices(A, I...)...)
        end
        function $_f(::IndexCartesian, A::AbstractArray{T,N}, I::Vararg{Int, N}) where {T,N}
            @_propagate_inbounds_meta
            $f(A, I...)
        end
    end
end

# setindex! like function
for f! in (:unsafe_setindex!,)
    _f! = Symbol(:_, f!)
    @eval begin
        function $f!(x::AbstractArray, v, I...)
            @_propagate_inbounds_meta
            Base.error_if_canonical_setindex(IndexStyle(A), A, I...)
            $_f!(IndexStyle(A), A, v, to_indices(A, I)...)
        end
        function $_f!(::IndexLinear, A::AbstractArray, v, I::Int...)
            @_inline_meta
            @boundscheck checkbounds(A, I...) # generally _to_linear_index requires bounds checking
            @inbounds $f!(A, v, Base._to_linear_index(A, I...))
            A
        end
        function $_f!(::IndexCartesian, A::AbstractArray, v, I::Int)
            @_inline_meta
            @boundscheck checkbounds(A, I...) # generally _to_subscript_indices requires bounds checking
            @inbounds $f!(A, v, Base._to_subscript_indices(A, I...)...)
            A
        end
        function $_f!(::IndexCartesian, A::AbstractArray{T,N}, v, I::Vararg{Int, N}) where {T,N}
            @_propagate_inbounds_meta
            $f!(A, v, I...)
            A
        end
    end
end

struct Power2{n} <: Integer
    function Power2{n}() where {n}
        new{n::UInt}()
    end
end
@pure Power2(n::Int) = (@assert n ≥ 0; Power2{unsigned(n)}())
@pure Power2(n::UInt) = Power2{n}()

function print_superscript(io::IO, i::Int)
    i == 0 && return print(io, Char(0x02070))
    i == 1 && return print(io, Char(0x000B9))
    i == 2 && return print(io, Char(0x000B2))
    i == 3 && return print(io, Char(0x000B3))
    i == 4 && return print(io, Char(0x02074))
    i == 5 && return print(io, Char(0x02075))
    i == 6 && return print(io, Char(0x02076))
    i == 7 && return print(io, Char(0x02077))
    i == 8 && return print(io, Char(0x02078))
    i == 9 && return print(io, Char(0x02079))
    error()
end
function Base.show(io::IO, p::Power2{n}) where {n}
    print(io, "2")
    for char in string(n)
        print_superscript(io, parse(Int, char))
    end
end

@inline Base.convert(::Type{T}, p::Power2{n}) where {T <: Integer, n} = one(T) << n
Base.promote_type(::Type{<: Power2}, ::Type{T}) where {T} = T
Base.promote_type(::Type{T}, ::Type{<: Power2}) where {T} = T

@inline Base.:*(a::Integer, p::Power2{n}) where {n} = a << n
@inline Base.:*(p::Power2{n}, a::Integer) where {n} = a << n
@inline Base.:*(p::Power2{m}, q::Power2{n}) where {m, n} = Power2(m + n)

@inline Base.div(a::Integer, p::Power2{n}) where {n} = a >> n
@inline Base.rem(a::Integer, p::Power2{n}) where {n} = a & (one(a) << n - 1)
@inline Base.divrem(a::Integer, p::Power2) = (div(a, p), rem(a, p))

@inline Base.zero(::Power2) = Power2(0)


text:
# TODO: remove tests when removed
@deprecate MahalanobisKernel(; P::AbstractMatrix{<:Real}) transform(
    SqExponentialKernel(), LinearTransform(sqrt(2) .* cholesky(P).U)
)

# TODO: remove keyword argument `maha` when removed
@deprecate PiecewisePolynomialKernel{V}(A::AbstractMatrix{<:Real}) where {V} transform(
    PiecewisePolynomialKernel{V}(size(A, 1)), LinearTransform(cholesky(A).U)
)

Base.@deprecate_binding TensorProduct KernelTensorProduct

# TODO: remove tests when removed
@deprecate kerneldiagmatrix(k::Kernel, x::AbstractVector) kernelmatrix_diag(k, x)
@deprecate kerneldiagmatrix(k::Kernel, x::AbstractVector, y::AbstractVector) kernelmatrix_diag(
    k, x, y
)
@deprecate kerneldiagmatrix!(K, k::Kernel, x::AbstractVector) kernelmatrix_diag!(K, k, x)
@deprecate kerneldiagmatrix!(K, k::Kernel, x::AbstractVector, y::AbstractVector) kernelmatrix_diag!(
    K, k, x, y
)
@deprecate kerneldiagmatrix(k::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs) kernelmatrix_diag(
    k, X; obsdim=obsdim
)
@deprecate kerneldiagmatrix(
    k::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=defaultobs
) kernelmatrix_diag(k, X, Y; obsdim=obsdim)
@deprecate kerneldiagmatrix!(K, k::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs) kernelmatrix_diag!(
    K, k, X; obsdim=obsdim
)
@deprecate kerneldiagmatrix!(
    K, k::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=defaultobs
) kernelmatrix_diag!(K, k, X, Y; obsdim=obsdim)


text:
using .MiniQhull
import .MiniQhull: delaunay
export delaunay

"""
    delaunay(X::LazySet)

Compute the Delaunay triangulation of the given convex set.

### Input

- `X` -- set

### Output

A union of polytopes in vertex representation.

### Notes

This function requires that you have properly installed
the Julia package [MiniQhull.jl](https://github.com/gridap/MiniQhull.jl), including
the library [Qhull](http://www.qhull.org/).

The method works in arbitrary dimension and the requirement is that the list of
vertices of `X` can be obtained.
"""
function delaunay(X::LazySet)
    n = dim(X)
    v = vertices_list(X)
    m = length(v)
    coordinates = vcat(v...)
    connectivity_matrix = delaunay(n, m, coordinates)
    nelements = size(connectivity_matrix, 2)
    elements = [VPolytope(v[connectivity_matrix[:, i]]) for i in 1:nelements]
    return UnionSetArray(elements)
end


text:
mutable struct axisObj

    active::Bool

    # width::Float64  # width of figure (in inches)
    # height::Float64 # height of figure (in inches)

    # lw::Float64                  # default line width
    # linecolor::Vector{Float64}   # default line color
    # axcolor::Vector{Float64}     # axis color

end

function axisObj( obj::tikzObj;
                  width = -1., height = -1., unit = "in",
                  lw = 0.,
                  pos = [0.,0.],
                  linecolor = [-1., -1., -1.],
                  axcolor = [0.1, 0.1, 0.1],
                  xlabel = "",
                  ylabel = "",
                  xmin = -Inf,
                  xmax =  Inf,
                  ymin = -Inf,
                  ymax =  Inf,
                  style = :Tufte )

    # open axis
    @printf( obj.fID, "  \\begin{axis}[%%\n" )

    #   height and width
    if width > 0.
        @printf( obj.fID, "    width  = %4.3f,\n", width )
    end
    if height > 0.
        @printf( obj.fID, "    height = %4.3f,\n", height )
    end
    @printf( obj.fID, "    at = {(%4.3f,%4.3f)},\n", pos[1], pos[2] )
    @printf( obj.fID, "    scale only axis,\n" )

    if style == :Tufte
        # Move axis out

    end

    # end of optional commands axis
    @printf( obj.fID, "    ]\n" )

    return axisObj( true )

end

function closeObj( objT::tikzObj, object::axisObj )

    @printf( objT.fID, "\n\\end{axis}\n")

    object.active = false

end

max_stars_repo_name: mopg/julia-tikz

max_issues_repo_name: mopg/julia-tikz

max_forks_repo_name: mopg/julia-tikz


text:
# Not using Artifacts for now
# This is a simple script to finally register the OpenCV package
# with the local package manager.

using Pkg

print(ARGS)
if size(ARGS) == 2
    Pkg.activate(ARGS[2])
end

Pkg.develop(PackageSpec(path=ARGS[1]))


text:
using Plots
using LightGraphs
using LayeredLayouts
using Test
using VisualRegressionTests

include("examples.jl")

@testset "LayeredLayouts.jl" begin
    include("demos.jl")
end


text:
using JuMP, GLPK

m = Model(GLPK.Optimizer)

@variable(m, 0 <= x <= 2 )
@variable(m, 0 <= y <= 30 )

@objective(m, Max, 5x + 3*y )
@constraint(m, 1x + 5y <= 3.0 )

println(m)
optimize!(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))

max_stars_repo_name: ricardoaf/juliafirststeps

max_issues_repo_name: ricardoaf/juliafirststeps

max_forks_repo_name: ricardoaf/juliafirststeps


text:
# Here I demo how the box counting process works
using FractionalDim

using PyCall
@pyimport matplotlib.pyplot as plt

# Curves
include("../test/curves.jl")

# Use flake
fractal = koch_line(8);

function plot(fractal::Vector{Segment{2}}, color="k")
    x = map(segment -> segment.O[1], fractal)
    y = map(segment -> segment.O[2], fractal)
    push!(x, last(fractal).B[1])
    push!(y, last(fractal).B[2])
    plt.plot(x, y, "-"*color)
end

counter = BoxCounter(fractal)
sizes, counts = [], []

state = start(counter, 8)

level, count = 0, 1
rate = Inf
while !done(counter, state)

    plt.figure()
    # Background
    plot(fractal)
    # Boxes on where the intersect are computed on current level
    for boxes in counter.boxes
        for box in boxes
            for line in surface(box)
                plt.plot([line.O[1], line.B[1]], [line.O[2], line.B[2]], "b")
            end
        end
    end

    state = next(counter, state)
    size, count = last(state)
    push!(sizes, size)
    push!(counts, count)

    if length(sizes) > 1
        rate = -log(counts[end]/counts[end-1])/log(sizes[end]/sizes[end-1])
    end
    plt.xlim((-0.05, 1.05))
    plt.ylim((-0.05, 1.05))
    plt.axis("off")
    plt.title("$(count)/$(4^level) -> $(rate)")
    # plt.show()
    plt.savefig("./koch_$(level).png", bbox_inches="tight")
    level += 1
end

max_stars_repo_name: MiroK/fractional-dim


text:
module TestVar

using Modia
using Modia.Measurements
using Modia.StaticArrays

struct nnnVar
	var::NamedTuple
	Var(; kwargs...) = new((; kwargs...) )
end

struct Model
	model::NamedTuple
	Model(; kwargs...) = new((; kwargs...) )
end

struct Map
	map::NamedTuple
	Map(; kwargs...) = new((; kwargs...) )
end

Var(;kwargs...) = (;kwargs...)
Var(value; kwargs...) = (; value=value, kwargs...)

#Var(args...; kwargs...) = begin println("1"); (;args..., kwargs...) end
#Var(value, args...; kwargs...) = typeof(value) <: NamedTuple ? begin println("2", value, args, kwargs); (;value..., args..., kwargs...) end : begin println("3"); (;value = value, args..., kwargs...) end
xar(value, args...; kwargs...) = begin
	@show value args kwargs merge(value, args[1], args[2], (;kwargs...))
	end

var(value, args...; kwargs...) = typeof(value) <: NamedTuple ?
	begin println("2", value, args, kwargs); value end : #merge(value, args, kwargs ) end :
	begin println("3"); (;value = value, args..., kwargs...) end
#Var(value, args...; kwargs...) = (;value = value, args..., kwargs...)
#Var(value::Union{Number, String, Expr, AbstractArray}, args...; kwargs...) = (;value = value, args..., kwargs...)

#NotPair = Union{Number, String, Expr, AbstractArray} # ..., ....}
#Var(value::NotPair, args...; kwargs...) = (;value = value, args..., kwargs...)

Model(; kwargs...) = (; kwargs...)

constant = Var(constant = true)
parameter = Var(parameter = true)
input = Var(input = true)
output = Var(output = true)
potential = Var(potential = true)
flow = Var(flow = true)

range(min, max) = Var(min=min, max=max)

recursiveMerge(x, ::Nothing) = x
recursiveMerge(x, y) = y
recursiveMerge(x::Expr, y::Expr) = begin dump(x); dump(y); Expr(x.head, x.args..., y.args...) end
recursiveMerge(x::Expr, y::Tuple) = begin x = copy(x); xargs = x.args; xargs[y[2]] = y[3]; Expr(x.head, xargs...) end

function recursiveMerge(nt1::NamedTuple, nt2::NamedTuple)
    all_keys = union(keys(nt1), keys(nt2))
    gen = Base.Generator(all_keys) do key
        v1 = get(nt1, key, nothing)
        v2 = get(nt2, key, nothing)
        key => recursiveMerge(v1, v2)
    end
    return Var(; gen...)
end

Base.:|(m::Model, n::Model) =  recursiveMerge(m.model, n.model)
Base.:|(m::Model, n::Map) = recursiveMerge(m.model, n.map)
Base.:|(m, n::Map) = recursiveMerge(m, n.map)
Base.:|(m, n) = if !(typeof(n) <: NamedTuple); recursiveMerge(m, (; value=n)) else recursiveMerge(n, (value=m,)) end

#Base.:|(m::Model, n::Model) =  (:MergeModel, m.model, n.model)
#Base.:|(m::Model, n::Map) =  (:MergeMap, m.model, n.map)
#Base.:|(m, n::Map) =  (:MergeMap, m, n.map)
#Base.:|(m::Var, n::Var) =  (:MergeModel, m.var, n.var)

# ---------------------------------------------------------------------

#v = Var(potential, nominal=10)
#@show v

v1 = input | output | flow | Var(min=0, nominal=10)
@show v1

v2 = parameter | 5
@show v2

v3 = 10u"m/s" | parameter
@show v3

v4 = parameter | 5.5 | Var(min=0)
@show v4

v5 = parameter | Var(5, min=0)
@show v5

Pin = Model( v = potential | Var(nominal=10), i = flow )
@show Pin

v6 = Var(2.0)
@show v6

v7 = Var(2.0 ± 0.1)
@show v7

v8 = Var(SVector(1.0,2.0,3.0))
@show v8

v9 = v4 | Var(parameter=false)
@show v9

v10 = v4 | Var(parameter=nothing)
@show v10

v11 = v4 | range(0,10)
@show v11

end


text:
"""
Copyright Marcus Greiff 2020

Tests of the attitude conversions general math functions
"""

# Include modules and recompile if running test standalone
if !isdefined(Main, :runAllTests)
    using Test, AerialVehicleControl, LinearAlgebra
    recompile();
end

# Testing fixture
function fixture()
    # Generate a random normalized Quaternion
    q = rand(4,1); q = q./norm(q);
    # Generate a random rotation matrix
    R = exp(skew(rand(3,1)));
    return q, R
end

@verbose(1, "Testing math functions...")
@testset "General math functions" begin
    @testset "Test the dot product function" begin
        @testset "Good inputs" begin
            for i = 2:5
                # Fixture
                v1, v1c, v1m = get_matrix_and_copy(rand(i,1));
                v2, v2c, v2m = get_matrix_and_copy(rand(i,1));
                prod = v1'*v2;

                # Test
                output = ccall((:cont_dot_product, AerialVehicleControl.CONT_LIB_PATH),
                    Cdouble,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    v1m,
                    v2m)

                # Assert
                @test isapprox(output, prod[1,1]);
                @test isapprox(v1, v1c);
                @test isapprox(v2, v2c);
            end
        end
        # TODO write bad input tests
    end
    @testset "Test the cross-product function" begin
        for i = 1:5
            # Fixture
            yA, yAc, yAm = get_matrix_and_copy(rand(3,1));
            yB, yBc, yBm = get_matrix_and_copy(rand(3,1));
            o,  oc,  om  = get_matrix_and_copy(rand(3,1));

            out_J      = zeros(3,1);
            out_J[:,1] = cross(yA[:,1], yB[:,1]);

            # Test
            ccall((:cont_cross_product, AerialVehicleControl.CONT_LIB_PATH),
                Cvoid,
                (Ref{matrix_double_t}, Ref{matrix_double_t}, Ref{matrix_double_t},),
                yAm,
                yBm,
                om)

            # Test that only the output has changed
            @test isapprox(yA, yAc);
            @test isapprox(yB, yBc);
            @test !isapprox(o,  oc);

            # Test that the correct values are computed
            @test isapprox(o,  out_J);

        end
    end
    @testset "Test the sinc function" begin
        @testset "Good inputs" begin
            for i = 1:10
                # Fixture
                input = 10.0* rand() - 5.0;

                # Test
                output = ccall((:cont_sinc, AerialVehicleControl.CONT_LIB_PATH),
                    Cdouble,
                    (Ref{Cdouble},),
                    input)

                # Assert
                @test isapprox(sin(input)/input, output);
            end
        end
    end
    @testset "Test the sign function" begin
        @testset "Good inputs" begin
            for i = 1:10
                # Fixture
                input = 10.0* rand() - 5.0;
                if i == 1;
                    input = 0.0;
                end
                sign_in = (input<0.0) ? -1.0 : 1.0;

                # Test
                output = ccall((:cont_sign_func, AerialVehicleControl.CONT_LIB_PATH),
                    Cdouble,
                    (Ref{Cdouble},),
                    input)

                # Assert
                @test isapprox(sign_in, output);
            end
        end
    end
    @testset "Test the conversion or quaternion to R in SO(3)" begin
        @testset "Check that R is in SO(3)" begin
            for i = 1:10
                # Fixture
                q, _ = fixture();
                q, qc, qm = get_matrix_and_copy(q);
                R, Rc, Rm = get_matrix_and_copy(rand(3,3));

                # Test
                status = ccall((:cont_quat_2_SO3, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    qm,
                    Rm)

                # Assert
                @test  isapprox(status, 1);
                @test  isapprox(q, qc);
                @test !isapprox(R, Rc);
                @test  isapprox(R'*R, Matrix{Float64}(I, 3, 3));
                @test  isapprox(det(R), 1.0);
            end
        end
        @testset "Check that the rotation is corresponds to the quaternion definition" begin
            for i = 1:10
                # Fixture
                theta = 2*pi*rand();
                u      = rand(3,1);
                u      = u ./ norm(u);
                q      = rand(4,1);
                q[1]   = cos(theta/2.0);
                q[2:4] = sin(theta/2.0)*u;

                Rtrue = exp(skew(u*theta));

                q, qc, qm = get_matrix_and_copy(q);
                R, Rc, Rm = get_matrix_and_copy(rand(3,3));

                # Test
                status = ccall((:cont_quat_2_SO3, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t}, Ref{matrix_double_t}),
                    qm,
                    Rm)

                # Assert
                @test  isapprox(status, 1);
                @test  isapprox(q, qc);
                @test !isapprox(R, Rc);
                @test  isapprox(R, Rtrue);
            end
        end
    end
    @testset "Test the vector normalization function" begin
        @testset "Column vector" begin
            for N = 1:5
                v, vc, vm = get_matrix_and_copy(rand(N,1));

                v_normalized = v ./ norm(v);

                status = ccall((:cont_normalize, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t},),
                    vm)

                @test isapprox(status, 1)
                for i = 1:N
                    @test  isapprox(v_normalized[i,1], v[i,1])
                    @test !isapprox(vc[i,1],           v[i,1])
                end
            end
        end
        @testset "Row vector" begin
            for N = 1:5
                v, vc, vm = get_matrix_and_copy(rand(1,N));

                v_normalized = v ./ norm(v);

                status = ccall((:cont_normalize, AerialVehicleControl.CONT_LIB_PATH),
                    Cint,
                    (Ref{matrix_double_t},),
                    vm)

                @test isapprox(status, 1)
                for i = 1:N
                    @test  isapprox(v_normalized[1,i], v[1,i])
                    @test !isapprox(vc[1,i],           v[1,i])
                end
            end
        end
        @testset "Bad inputs" begin
            v, vc, vm = get_matrix_and_copy(rand(2,2));

            v_normalized = v ./ norm(v);

            status = ccall((:cont_normalize, AerialVehicleControl.CONT_LIB_PATH),
                Cint,
                (Ref{matrix_double_t},),
                vm)

            @test isapprox(status, 0)
        end
    end
end


text:
module Actors

using DocStringExtensions

# Misc Types
export Id, Scene

# Actors
export Stage, Troupe

# Functions
export stage, play!, enter!, invite!, leave!, forward!, expect, ask, say, hear, me, my
export delegate, shout, minder, @say_info, async, interact!, local_addresses

# Messages
export Genesis!, Leave!, LogInfo!, Died!, Left!, AsyncFail!

@template TYPES = """
$DOCSTRING

### Members

$TYPEDFIELDS
"""

"""The Address of an [`Actor`](@ref)

This is a safe reference to an [`Actor`](@ref). It is most commonly used to
send messages to an [`Actor`](@ref). However many accessor methods take an
`Addr` to safely get or set some `Actor`'s internals or associated data.

!!! warning

    This will grow in size as remote actors are added.

### Type Parameters

- `S` The type of the actor state.
- `M` The message types the actor accepts, usually Any.

"""
struct Id{S, M}
    "The local index within the owning stage"
    inner::UInt32

    Id(inner) = Id{Any, Any}(inner)
    Id{S}(inner) where S = Id{S, Any}(inner)
    function Id{S, M}(inner) where {S, M}
        @assert inner > 0 && inner <= typemax(UInt32)

        Id{S, M}(UInt32(inner))
    end
    Id{S, M}(inner::UInt32) where {S, M} = new{S, M}(inner)
end

Base.:(==)(a::Id, b::Id) = a.inner == b.inner
Base.show(io::IO, id::Id{S}) where S = print(io, "$S@", id.inner)

"""The star of the show (but not really)

This holds internal state for an Actor, it wraps the user defined state
value. We usually think of the Actor as being the state value rather than this
structure which is mostely hidden. It is rare that a user will need to access
this directly. Usually Actors are referenced by an [`Id`](@ref) and we get the
actor's details by calling accessor functions on either the `Id` or the
[`Scene`](@ref).

If you do access this, then you need to be careful to avoid concurrency
violations.

### Type Parameters

- `S` The type of the user defined actor state.
- `M` The message types the actor accepts, usually Any.

It is rare to set `M`, but if it is set then it should include at least
[`Leave!`](@ref).
"""
struct Actor{S, M}
    "How the Actor recieves messages, see [`listen!`](@ref)"
    inbox::Channel{M}
    "An arbitrary value which is usually thought of as the actor"
    state::S
    "The `Id` of another actor which manages and supports this actor"
    minder::Id
    "The Task this actor runs/ran in"
    task::Union{Task, Nothing}
    "The cannonical `Id` of this `Actor`"
    me::Id{S, M}

    Actor(inbox::Channel{M}, state::S, minder) where {S, M} =
        new{S, M}(inbox, state, minder, nothing)
    Actor(inbox::Channel{M}, state::S, minder, task, me) where {S, M} =
        new{S, M}(inbox, state, minder, task, me)
end

"Create an Actor with the given state and minder"
Actor{M}(state, task, minder, me) where M =
    Actor(Channel{M}(420), state, task, minder, me)

include("addressing.jl")

"""Abstract Stage

Mainly used to allow overriding of [`Stage`](@ref)'s methods.
"""
abstract type AbsStage end

"""The Root Actor

This contains the addresses for all the actors and bootstraps the
`play`. Currently all actors are created by sending messages to the `Stage`
although this will probably change in the future.

Any messages which the `Stage` can not handle, it forwards to the `play`.

The `Stage` bootstraps itself by putting the [`PreGenesis!`](@ref) message in
its [`inbox`](@ref). It processes this message first then sends
[`Genesis!`](@ref) to the actor specified by [`play!`](@ref).

It is unusual for the user to access this directly or override its
behaviour. Most things can be achieved by creating a new [`AbsMinder`](@ref)
type or allowing messages to be forwarded to the Play actor.
"""
mutable struct Stage <: AbsStage
    "All the `Actor`s in the play"
    actors::AddressBook
    "Grace period timer before force leaving"
    time_to_leave::Union{Timer, Nothing}
    "User defined play `Actor`"
    play::Id

    Stage() = new(AddressBook(), nothing)
end

abstract type AbsScene{S, M} end

"""The context of message processing

Contains common information which is used by many different methods during
message handling. You should assume that the members of this struct are likely
to change.

### Type Parameters

- `S` The type of the current [`Actor`](@ref)'s state. This is commonly
  specified when adding a method for [`hear`](@ref) (amongst much else).

- `M` The message types accepted by the current `Actor`, usually `Any`.

"""
struct Scene{S, M} <: AbsScene{S, M}
    "The current [`Actor`](@ref)"
    subject::Actor{S, M}
    "The [`Stage`](@ref)"
    stage::Actor{Stage}
end

struct AsyncScene{S, M} <: AbsScene{S, M}
    s::Scene{S, M}
    task::Task
end

Base.show(io::IO, s::Scene{S}) where S = print(io, "Scene{$S}")
Base.show(io::IO, s::AsyncScene{S}) where S = print(io, "AsyncScene{$S}")

const WRONG_TASK =
    "The given Scene or Actor is not associated with this task. You are either trying to access another Actor's state or have used @async (use @try_async)."

subject(s::Scene) = let a = s.subject
    @assert a.task !== nothing "Actor is not playing"
    @assert a.task === current_task() WRONG_TASK

    a
end

function subject(s::AsyncScene)
    @assert s.task === current_task() WRONG_TASK

    s.s.subject
end

stage_ref(s::Scene) = let a = s.subject
    @assert a.task !== nothing "Actor is not playing"
    @assert a.task === current_task() WRONG_TASK

    s.stage
end

function stage_ref(s::AsyncScene)
    @assert s.task === current_task() WRONG_TASK

    s.s.stage
end

"""Safely get the current [`Actor`](@ref)'s state

Usually the user passes the [`Scene`](@ref) to this and gets the executing
[`Actor`](@ref)'s state in return.
"""
my(s::AbsScene) = subject(s).state

"""Get the inbox of an [`Actor`](@ref)

Useful when overriding functions such as [`listen!`](@ref) or
[`leave!`](@ref). Otherwise it is quite unusual for the user to call
this. You should use [`say`](@ref) and [`ask`](@ref) to send messages.

The `inbox` is how an `Actor` recieves messages. In a local system the sender
directly places a message into the `inbox`. Currently inboxes are implemented
with a Julia `Channel` which can be safely accessed by multiple threads,
however it is generally expected that you don't `pop!` messages from another
actor's inbox unless that actor is dead.
"""
inbox(s::AbsScene, a::Id) = stage_ref(s).state.actors[a].inbox
inbox(s::AbsScene) = subject(s).inbox

"""Get the address of an [`Actor`](@ref)'s minder

When called on the [`Scene`](@ref) it will get the current Actor's minder. If
called on an [`Id`](@ref) it will get the minder of the Actor pointed to by
the address.

!!! warning

    It is generally not thread-safe to change `Actor.minder`. However you can
    re-assign a minder's [`Id`](@ref) to another actor.

See [`AbsMinder`](@ref).
"""
minder(s::AbsScene, a::Id)::Id = stage_ref(s).state.actors[a].minder
minder(s::AbsScene) = subject(s).minder

"Get the address ([`Id`](@ref)) of the current [`Actor`](@ref)"
me(s::AbsScene) = subject(s).me

"Get the address ([`Id`](@ref)) of the [`Stage`](@ref)"
function stage(s::AbsScene)
    subject(s)

    Id{Stage}(UInt32(0))
end

"""Send a message asynchronously

This returns immediately and doesn't guarantee the message is processed by the
actor specified by `to`. Nor does it guarantee messages will be delivered in
the order `say` is called.

This may throw an error if, for example, `to` only accepts certain message
types.

This doesn't expect a response from `to`. For that see [`ask`](@ref).
"""
function say(s::AbsScene, to::Id, msg)
    @debug "$s send" to msg
    put!(inbox(s, to), msg)
end

say(s::AbsScene, ::Id{Stage}, msg) = put!(stage_ref(s).inbox, msg)

"""Like [`say`](@ref), but *less likely* to throw an exception if `to` is
dead. This is commonly useful when an unexpected error has occurred and you
want to perform some cleanup without knowing which actors are still
alive. Under normal operation you should expect other actors to be alive if
you are sending them a message.

!!! note

    If, for example, `to` exceeds the highest known ID in the local
    [`Stage`](@ref) then this will still throw a `BoundsError`. This is to
    help prevent silent errors.

"""
try_say(s::AbsScene, to::Id, msg) = try
    say(s, to, msg)
catch ex
    ex isa Union{InvalidStateException, KeyError} || rethrow()
end

"""Handle a received message

Usually called by [`listen!`](@ref) to handle a message taken from the
[`inbox`](@ref). The user defines new `hear` methods to handle messages for
different [`Actor`](@ref)-message type combinations.

It is passed a [`Scene`](@ref) object which has the [`Actor`](@ref) state type
as the first type parameter. By convention this argument is always called `s`,
if you use a different name it may break some non-essential macros. The other
parameter is the message, which is usually called `msg`, unless you have a
more appropriate name, and can be of any type (in a local system).

If `hear` returns a value it will most likely be ignored (unless you override
[`listen!`](@ref)). If you need to respond to a message then use [`say`](@ref)
or [`ask`](@ref).
"""
hear(s::AbsScene{<:AbsStage}, msg) = say(s, my(s).play, msg)

"""Take messages from the [`inbox`](@ref) and process them

By default this simply takes messages from the [`inbox`](@ref) and calls
[`hear`](@ref) on them. However you may wish to process some messages
specially or do some work inbetween messages, in which case you can override
this for a given [`Actor`](@ref) type.

As a general rule, you should make sure to call [`hear`](@ref) on a message if
you can't process it some other way. Otherwise you will prevent some standard
messages from working (such as [`Leave!`](@ref)).
"""
function listen!(s::Scene)
    @debug "$s listening"

    for msg in inbox(s)
        @debug "$s recv" msg

        hear(s, msg)
    end
end

function listen!(s::Scene{<:AbsStage})
    inb = inbox(s)
    as = my(s).actors

    @debug "$s listening"
    for msg in inb
        @debug "$s recv" msg

        hear(s, msg)

        if !isnothing(my(s).time_to_leave) && isempty(as)
            @debug "All actors left on time"
            close(my(s).time_to_leave)
            close(inb)
        end
    end

    @assert !isnothing(my(s).time_to_leave)
end

"""Signal that the current [`Actor`](@ref) should exit (leave)

Usually this closes the [`inbox`](@ref) which will cause the actor to exit
once it has processed any messages it already received. Some actors (e.g
[`Stage`](@ref)) have a grace period where they will continue to accept new
messages.

If you wish to exit immediately then throw an exception.

Also see [`Leave!`](@ref).
"""
leave!(s::AbsScene) = close(inbox(s))

function leave!(s::Scene{<:AbsStage})
    write_lock(my(s).actors) do table
        for a in keys(table.rev_entries)
            try
                put!(a.inbox, Leave!())
            catch ex
                ex isa InvalidStateException || rethrow()
            end
        end
    end

    inb = inbox(s)
    as = my(s).actors
    my(s).time_to_leave = timer = Timer(3)
    @async begin
        wait(timer)
        close(inb)
        @debug "$s Exit grace period over"

        write_lock(as) do table
            for a in keys(table.rev_entries)
                @debug "$a took too long to leave, forcibly closing inbox..."
                close(a.inbox)

                try
                    isnothing(a.task) || wait(a.task)
                catch ex
                    @debug "$a Errored" ex
                end
            end
        end
    end
end

function forward!(s::AbsScene, to::Id)
    as = stage_ref(s).state.actors
    ids = as[subject(s)]
    to_a = as[to]

    for id in ids
        as[Id(id)] = to_a
    end

    task = async(s) do s
        for msg in inbox(s)
            say(s, to, msg)
        end
    end

    yield()
    leave!(s)
    wait(task)
end

"""Used to capture variables from the parent thread/task

This is a workaround for integrating third party libraries which use Task
local storage or similar. By default it does nothing and it is probably best
to avoid using it. See [`enter!`](@ref) and [`play!`](@ref).

Note that this is not necessary for getting OS "environment variables", which
are typically shared amongst threads.
"""
capture_environment(id) = nothing

play!(play) = let id = Id{Stage}(UInt32(0))
    a = Actor(Channel{Any}(1024), Stage(), id, current_task(), id)

    put!(a.inbox, PreGenesis!(play))

    play!(Scene(a, a), capture_environment(id))
end

"""Ran before [`listen!`](@ref) and in [`async`](@ref)

By default does nothing, but can be overriden to mess with an
[`Actor`](@ref)'s internals. It is passed `env` which is taken by
[`capture_environment`](@ref).

See [`enter!`](@ref) and [`play!`](@ref).
"""
function prologue!(s::AbsScene, env) end

"""Start the 'Actor System' or a single [`Actor`](@ref)

Typically the user calls this method on some value which is used as the 'play'
[`Actor`](@ref) state. This then blocks the calling thread until the
[`Stage`](@ref) is sent [`Leave!`](@ref) or dies due to an error which was
allowed to propagate.

This is also called by [`enter!`](@ref), in a new Thread/Task, to run a single
[`Actor`](@ref). Nothing should happen before or after [`play!`](@ref) so it
contains the `Actor`'s full lifetime (excluding Julia internals).

The liftime of the `Actor` looks something like:
1. [`prologue!`](@ref)
2. [`listen!`](@ref)   (unless `prologue!` errors)
3. [`epilogue!`](@ref) (unless `listen!` errors)

Then, if there is an error

4. [`dieing_breath!`](@ref)

You can override any of the above, so if you are forced to override
`play!(s::Scene, env)` this is considered a library design error. Usually if
you feel the need to override any of these, then you should solve the problem
with [`AbsMinder`](@ref) or more message passing instead.
"""
play!(s::Scene, env) = try
    prologue!(s, env)
    listen!(s)
    epilogue!(s, env)
catch ex
    dieing_breath!(s, ex, env)
    rethrow()
finally
    close(inbox(s))
end

"""Run after [`listen!`](@ref)

By default says to [`minder`](@ref) that the [`Actor`](@ref)
[`Left!`](@ref). However it can be overriden to mess with the `Actor`'s
internals.It is passed `env` which is taken by [`capture_environment`](@ref).

See [`enter!`](@ref) and [`play!`](@ref).
"""
epilogue!(::AbsScene, env) = nothing
epilogue!(s::Scene, env) = say(s, minder(s), Left!(me(s)))
epilogue!(::Scene{<:AbsStage}, env) = nothing

"""Run if an exception is thrown in [`play!`](@ref)

By default says to [`minder`](@ref) that the [`Actor`](@ref)
[`Died!`](@ref). However it can be overriden to mess with the `Actor`'s
internals. It is passed `env` which is taken by [`capture_environment`](@ref).

See [`enter!`](@ref) and [`play!`](@ref).
"""
function dieing_breath!(s::Scene, ex, env)
    @debug "$(me(s)) Died" ex
    say(s, minder(s), Died!(me(s), s.subject))
end
dieing_breath!(s::Scene{<:AbsStage}, ex, env) = leave!(s)

"Create a new Task, Thread or similar primitive"
function fork(fn::Function)
    task = Task(fn)
    task.sticky = false
    schedule(task)
end

enter!(s::AbsScene, actor_state::S) where S = enter!(s, actor_state, Any)
enter!(s::AbsScene, actor_state::S, ::Type{M}) where {S, M} =
    enter!(s, actor_state, minder(s), M)
enter!(s::AbsScene, actor_state::S, minder::Id) where S =
    enter!(s, actor_state, minder, Any)
enter!(s::AbsScene, actor_state::S, minder::Id, ::Type{M}) where {S, M} =
    enter!(s, Actor(Channel{M}(512), actor_state, minder))

"""Add a new [`Actor`](@ref) to the [`Stage`](@ref)/Play and return its
[`Id`](@ref). This is usually called by a [`AbsMinder`](@ref)'s
[`Invite!`](@ref) handler.

!!! note

    The `Actor` struct `a` is recreated by the new `Actor` so you can not use
    it to do a reverse lookup of the `Id`.

"""
function enter!(s::AbsScene, a::Actor{S, M}) where {S, M}
    st = stage_ref(s)
    id = push!(st.state.actors, a)
    env = capture_environment(id)

    @debug "$s forking" a
    fork() do
        try
            a = Actor(a.inbox, a.state, a.minder, current_task(), id)
            st.state.actors[id] = a
            play!(Scene(a, st), env)
        catch ex
            @debug "Actor died" s a ex
            rethrow()
        end
    end

    id
end

"""Wait for a message of type `R` and return it

This will block waiting for a message of the right type. It will cause the
[`Actor`](@ref) to become 'insensitive', meaning that all other messages will
be buffered while waiting.

When a message of the right type is recieved, then the buffered messages will
be put back in the [`inbox`](@ref) and the matching message is returned.

In theory this could accept the wrong message if the type matches, but it was
from an old request or it was sent for some other reason. One way to
avoid this is to [`delegate`](@ref) the [`ask`](@ref) request, or the entire
operation, to a [`Stooge`](@ref) which will have a new address.
"""
function expect(s::AbsScene, ::Type{R}) where R
    inb = inbox(s)
    msg = take!(inb)
    msg isa R && return msg

    scratch = Any[msg]
    for outer msg in inb
        msg isa R && break

        push!(scratch, msg)
    end

    foreach(m -> put!(inb, m), scratch)

    msg
end

"""Like [`say`](@ref), but wait for a response of a given type `R`

Simply calls `say` then [`expect`](@ref).
"""
function ask(s::AbsScene, a::Id, favor, ::Type{R}) where R
    me(s) == a && error("Asking oneself results in deadlock")
    say(s, a, favor)
    expect(s, R)
end

# Messages

"Used by the [`Stage`](@ref) for bootstraping"
struct PreGenesis!{T}
    "The state of the first user defined [`Actor`](@ref)"
    play::T
end

function hear(s::Scene{<:AbsStage}, msg::PreGenesis!{S}) where S
    logger = enter!(s, Logger(stdout), me(s), LoggerMsgs)
    passive_minder = enter!(s, PassiveMinder(logger), me(s))
    tree_minder = enter!(s, TreeMinder(passive_minder), passive_minder)

    play = my(s).play = ask(s, tree_minder, Invite!(me(s), msg.play), Invited!{S}).who
    say(s, play, Genesis!())
end

"Sent to the users 'play' [`Actor`](@ref) to get things started"
struct Genesis! end

"Informs that an [`Actor`](@ref) left"
struct Left!
    "The `Actor` who left"
    who::Id
end

function hear(s::Scene{Stage}, msg::Left!)
    my(s).actors[msg.who] = nothing

    isnothing(my(s).time_to_leave) && leave!(s)
end

"Informs that an [`Actor`](@ref) died"
struct Died!
    "The `Actor` who died"
    who::Id
    "The `Actor` who died's data"
    corpse::Actor
end

hear(s::Scene{<:AbsStage}, msg::Died!) = leave!(s)

"Tell an [`Actor`](@ref) to [`leave!`](@ref)"
struct Leave! end

hear(s::AbsScene, msg::Leave!) = leave!(s)
hear(s::AbsScene{<:AbsStage}, msg::Leave!) = say(s, my(s).play, msg)

# Actors (Other than Stage)

"""An [`Actor`](@ref) which prints messages to an `IO` stream

By default one of these is created which writes to stdout exlusively.
"""
struct Logger{I <: IO}
    io::I
end

struct LogDied!
    header::String
    died::Died!
end

hear(s::Scene{<:Logger}, msg::LogDied!) = try
    io = my(s).io

    printstyled(io, "Error: "; bold=true, color=Base.error_color())
    printstyled(io, msg.header; color=Base.error_color())
    println(io)

    task = msg.died.corpse.task
    if isnothing(task)
        println(io, "Actor task is nothing; actor was not started?")
    elseif isnothing(task.exception)
        println(io, "Actor task has no exception? $(task)")
    else
        showerror(io, task.exception, task.backtrace)
    end

    flush(io)
catch ex
    @debug "Arhhgg; Logger died while trying to do its basic duty" ex
    rethrow()
end

struct LogInfo!
    from::Id
    mod::String
    file::String
    line::Int64
    info::String
end

"""Print an info message

Uses the [`Logger`](@ref) provided by [`minder`](@ref) to print some
informational text.
"""
macro say_info(s, exp)
    (mod, file, line) = Base.CoreLogging.@_sourceinfo
    file = basename(file)
    mod = string(mod)

    esc(:(say($s, minder($s), LogInfo!(me($s), $mod, $file, $line, $exp))))
end

hear(s::Scene{<:Logger}, msg::LogInfo!) = try
    io = my(s).io

    printstyled(io, "Info"; bold=true, color=Base.info_color())
    printstyled(io, " $(msg.from) $(msg.mod) $(msg.file):$(msg.line): "; color=:light_cyan)
    println(io, msg.info)
    flush(io)
catch ex
    @debug "Arhhgg; Logger died while trying to do its basic duty" ex
    rethrow()
end

const LogMsgs = Union{LogInfo!, LogDied!}
const LoggerMsgs = Union{LogMsgs, Leave!}

"""Request an actor be created with the given `state`. This is usually sent to
an actor's minder which uses [`enter!`](@ref) to create the new actor.

It is expected that minders (subtypes of [`AbsMinder`](@ref)) implement this
unless they are only intended for use with a special set of actors.
"""
struct Invite!{S}
    re::Id
    state::S
end

"A response to [`Invite!`](@ref)"
struct Invited!{S}
    who::Id{S}
end

"""Request an actor be created with the given `state` and return the new
[`Id`](@ref). This is just a wrapper for sending [`Invite!`](@ref) to the
current actor's [`minder`](@ref) and waiting for a response.
"""
invite!(s::Scene, state::S) where S =
    ask(s, minder(s), Invite!(me(s), state), Invited!{S}).who

"""An [`Actor`](@ref) which looks after other `Actors`'s

Every `Actor` is assigned a [`minder`](@ref) which it can call upong to
provide general services (e.g. logging). The [`minder`](@ref) is also notified
when an `Actor` dies or decides to [`Leave!`](@ref). Therefor the `minder` can
take action to replace a failed actor.
"""
abstract type AbsMinder end

hear(s::Scene{<:AbsMinder}, msg::LogMsgs) = say(s, minder(s), msg)
hear(s::Scene{<:AbsMinder}, msg::Left!) = let as = stage_ref(s).state.actors
    wait(as[msg.who].task)
    as[msg.who] = nothing
    leave!(s)
end
hear(s::Scene{<:AbsMinder}, msg::Died!) = say(s, minder(s), msg)

struct PassiveMinder{L <: Id} <: AbsMinder
    logger::Union{L, Nothing}
end

hear(s::Scene{<:PassiveMinder}, msg::Invited!) = nothing
hear(s::Scene{<:PassiveMinder}, msg::LogMsgs) = say(s, my(s).logger, msg)
hear(s::Scene{<:PassiveMinder}, msg::Died!) = try
    Base._wait(msg.corpse.task)
    say(s, my(s).logger, LogDied!("$(me(s)): Actor $(msg.who) died!", msg))
    leave!(s)
catch ex
    @debug "Arrgg; $(typeof(my(s))) died while trying to do its basic duty" ex
    rethrow()
end

"""Organises actors into a hieracy or tree structure. This is used as the
default [`AbsMinder`](@ref) for the play.

When an actor with this minder uses [`invite!`](@ref) the returned actor will
be a child of the current actor (the current actor being the parent). An actor
may have many children, but only a single parent.

* If a child actor dies then the parent recieves [`Died!`](@ref). By default
  actors don't know how to handle this so they will also die. (In any case the
  death will be logged).

* If a parent leaves or dies, then its children are also asked to
  [`Leave!`](@ref).

If you don't wish for an actor to die when one of its children do, then define
[`hear`](@ref) for the [`Died!`](@ref) message.
"""
mutable struct TreeMinder <: AbsMinder
    root::Id
    minded::Union{Id, Nothing}
    children::Set{Id}
end

TreeMinder(root) = TreeMinder(root, nothing, Set())

hear(s::Scene{TreeMinder}, msg::Invite!{S}) where S = let m = my(s)
    if m.minded === nothing
        m.minded = enter!(s, msg.state, me(s))
        say(s, msg.re, Invited!(m.minded))
    else
        child_m = enter!(s, TreeMinder(m.root), me(s))
        say(s, child_m, msg)
        push!(m.children, child_m)
    end
end

hear(s::Scene{TreeMinder}, msg::Leave!) = if isnothing(my(s).minded)
    leave!(s)
else
    try_say(s, my(s).minded, msg)
end

stop_children(s::Scene, m::TreeMinder) = for child in m.children
    try_say(s, child, Leave!())
end

hear(s::Scene{TreeMinder}, msg::Left!) = let m = my(s)
    if msg.who == m.minded
        stop_children(s, m)

        as = stage_ref(s).state.actors
        wait(as[msg.who].task)
        as[msg.who] = nothing
        m.minded = nothing
    else
        @assert msg.who in m.children "$(msg.who) not in $(m.children)"
        delete!(m.children, msg.who)
    end

    isnothing(m.minded) && isempty(m.children) && leave!(s)
end

hear(s::Scene{TreeMinder}, msg::Died!) = let m = my(s)
    if msg.who in m.children
        say(s, m.root, msg)
    elseif msg.who == m.minded
        stop_children(s, m)
        say(s, minder(s), msg)

        Base._wait(msg.corpse.task)
        stage_ref(s).state.actors[msg.who] = nothing
        m.minded = nothing

        say(s, m.root, LogDied!("$(me(s)): Actor $(msg.who) died!", msg))
        isempty(m.children) && leave!(s)
    else
        isnothing(m.minded) || say(s, m.minded, msg)
    end
end

"""A temporary [`Actor`](@ref) which performs one `action` then leaves

Used by [`delegate`](@ref)
"""
struct Stooge
    "What the `Stooge` should do"
    action::Function
    "The arguments passed to `action`"
    args::Tuple
end

listen!(s::Scene{Stooge}) = let stooge = my(s)
    stooge.action(s, stooge.args...)

    leave!(s)
end

"""Perform some `action` in a temporary [`Actor`](@ref)

This creates a new `Actor`, called a [`Stooge`](@ref) and uses it to perform
some action in parallel.

!!! warning

    Do not use variables captured from the surrounding scope, pass them as
    `args...`. The `args` may be copied to avoid concurrency violations.
"""
delegate(action::Function, s::AbsScene, args...) =
    invite!(s, Stooge(action, args))
delegate(action::Function, s::AbsScene, minder::Id, args...) =
    enter!(s, Stooge(action, args), minder)

"""A group of [`Actor`](@ref)'s

Messages [`shout`](@ref)ed at the `Troupe` will be broadcast to the all the
`Actor`'s in it. To send a message to the `Troupe` itself, just use
[`say`](@ref).

### Example

Where `a1..a3` are `Actor` [`Id`](@ref)'s

```
g = enter!(s, Troupe(a1, a2, a3))
shout(s, g, Leave!())
say(s, g, Leave!())
```
"""
struct Troupe
    "The addresses of the `Actor`s in a `Troupe`"
    as::Vector{Id}

    Troupe(as...) = new([as...])
end

"Wraps a message to broadcast ([`shout`](@ref))"
struct Shout!{T}
    msg::T
end

"Broadcast a message to a [`Troupe`](@ref)"
shout(s::AbsScene, troupe::Id{Troupe}, msg) = say(s, troupe, Shout!(msg))

"Broadcast a message to a plain vector of actors"
shout(s::AbsScene, ids, msg) = for a in ids
    say(s, a, msg)
end

hear(s::Scene{Troupe}, msg::Shout!) = shout(s, my(s).as, msg.msg)

"Informs that an asynchronous `Task` failed"
struct AsyncFail!
    "What failed"
    async::Task
end

hear(s::Scene, msg::AsyncFail!) = wait(msg.async)

"""Similar to `@async`, but notifies the calling [`Actor`](@ref) if it fails.

This allows you to create an asynchronous `Task` without waiting on it to
check for errors. If an exception is thrown, the `Actor` will send
[`AsyncFail!`](@ref) to itself. By default this will cause
`TaskFailedException` to be thrown thus killing the `Actor`.
"""
function async(fn, sc::Scene)
    env = capture_environment(me(sc))

    @async begin
        s = AsyncScene(sc, current_task())

        try
            prologue!(s, env)
            fn(s)
            epilogue!(s, env)
        catch ex
            say(s, me(s), AsyncFail!(s.task))
            rethrow()
        end
    end
end

include("interactive.jl")

end # module


text:
"""
    ReadSegyHeader(filename;<keyword arguments>)

Read the headers of a input file in segy format

# Arguments
- `group="all"` : Options are all, some or gather
- `swap_bytes=true` : If the flag equals true, the function swaps bytes
- `key=[]`
- `minval=0`
- `maxval=0`
"""
function ReadSegyHeader(filename_in;swap_bytes="true",group="all",key=" ",minval=0,maxval=0)

#1) Stats on segy header

stream = open(filename_in)

position = 3200

seek(stream, position)

fh = SeisMain.GrabFileHeader(stream)

ntfh = swap_bytes == true ? bswap(fh.netfh) : fh.netfh
fh = 0

if ntfh == -1
	error("add instructions to deal with variable extended text header")
end
if ntfh == 0
	file_hsize = 3600
elseif ntfh > 0
	# file_hsize = 3200 * (ntfh+1) + 400
	file_hsize = 3200 * 1 + 400
else
	error("unknown data format")
end

stream = open(filename_in)

seek(stream, SeisMain.segy_count["ns"] + file_hsize)

if (swap_bytes=="true")
	nt = bswap(read(stream,Int16))
else
	nt = read(stream,Int16)
end
total = 60 + nt

nx = round(Int,(filesize(stream)-file_hsize)/4/total)

println("total number of traces: ",nx)
println("number of samples per trace: ",nt)

itrace = 1
ntrace = nx
tt = nx

if group == "gather"
# Using this option means that the traces are organized as the gather. If traces are not organized, this option won't work

       itrace = nx
       ntrace = 1
       for j=1:nx
       	   h1 = SeisMain.GrabSegyHeader(stream,swap_bytes,nt,file_hsize,j)
	   aux = getfield(h1,Symbol(key))
	   if aux >=minval && aux<=maxval
	      if itrace>j
	      	 itrace = j
	      elseif j > ntrace
	      	 ntrace = j
	      end
	   end
	end
	tt = ntrace - itrace +1
	println(" number of traces in gather: ",tt)
end

h_segy = Array{SeisMain.SegyHeader}(undef,tt)

for i=1:tt
# este position me parece que esta mal + no lo necesito
  j = i + itrace -1
  position = file_hsize + total*(j-1)*4 + SeisMain.segy_count["trace"]
    seek(stream,position)

#este solo ya sabe donde ir

    h_segy[i] = SeisMain.GrabSegyHeader(stream,swap_bytes,nt,file_hsize,j)

end

#podria haber hecho h_segy = push!(h_segy, SeisMain.GrabSegyHeader(stream,swap_bytes,nt,file_hsize,j)
# me parece que no por el tema de definir h_segy pero no estoy segura

close(stream)
println("The output array has the segy data header.")
return h_segy
end


text:
module WorkshopWizard

using Pkg

const GITHUB_BASEURL = "https://github.com/carstenbauer/"
const WORKSHOPS = [
    "JuliaNRWSS21",
    "JuliaCologne21",
    "JuliaNRW21",
    "JuliaOulu20",
    "JuliaWorkshop19",
    "JuliaWorkshop18"
    ]

const WORKSHOP_JULIA_VERSION = Dict(
    "JuliaNRWSS21" => v"1.6",
    "JuliaCologne21" => v"1.5",
    "JuliaNRW21" => v"1.5",
    "JuliaOulu20" => v"1.3",
    "JuliaWorkshop19" => v"1.2",
    "JuliaWorkshop18" => v"1.0"
    )

default_workshop() = WORKSHOPS[1]
default_repo() = joinpath(GITHUB_BASEURL, default_workshop())
default_path() = begin
    if Sys.iswindows()
        return joinpath(homedir(), "Desktop")
        # homedir()
    else
        return homedir()
    end
end

include("IJulia.jl")
include("utils.jl")
include("install.jl")

end # module


text:
using Test
using ParameterEstimocean
using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity, MixingLength, SurfaceTKEFlux
using Suppressor: @suppress

using ParameterEstimocean.Parameters: closure_with_parameters, update_closure_ensemble_member!

const CAVD = ConvectiveAdjustmentVerticalDiffusivity
const CATKE = CATKEVerticalDiffusivity

@testset "Parameters tests" begin
    @info "Testing closure_with_parameters on ConvectiveAdjustmentVerticalDiffusivity"
    closures = [CAVD(background_κz=1.0) CAVD(background_κz=1.1)
                CAVD(background_κz=1.2) CAVD(background_κz=1.3)
                closure_with_parameters(CAVD(), (background_κz=1.4,)) closure_with_parameters(CAVD(), (background_κz=1.5,))]

    @test closures[3, 1].background_κz == 1.4
    @test closures[3, 2].background_κz == 1.5

    closure_tuple = closure_with_parameters((CAVD(), CAVD()), (background_κz = 14.0,))

    for closure in closure_tuple
        @test closure.background_κz == 14.0
    end

    @info "Testing closure_with_parameters on CATKEVerticalDiffusivity"
    catke_closure = @suppress CATKE()
    new_Cᴷuʳ = 1e3
    new_catke_closure = closure_with_parameters(catke_closure, (; Cᴷuʳ = new_Cᴷuʳ))
    @test new_catke_closure.mixing_length.Cᴷuʳ == new_Cᴷuʳ

    @info "Testing update_closure_ensemble_member! on ConvectiveAdjustmentVerticalDiffusivity"
    new_background_κz, new_convective_νz = 100.0, 10.0
    ensemble_member_no = 2
    update_closure_ensemble_member!(closures, ensemble_member_no, (background_κz = new_background_κz, convective_νz = new_convective_νz))

    for j in 1:2
        @test closures[ensemble_member_no, j].background_κz == new_background_κz
        @test closures[ensemble_member_no, j].convective_νz == new_convective_νz
    end

    @info "Testing update_closure_ensemble_member! on CATKEVerticalDiffusivity"
    closures = @suppress [CATKE() CATKE()
                          CATKE() CATKE()
                          CATKE() CATKE()]

    new_CᴷRiʷ = 100.0
    ensemble_member_no = 2
    update_closure_ensemble_member!(closures, ensemble_member_no, (; CᴷRiʷ = new_CᴷRiʷ))

    for j in 1:2
        @test closures[ensemble_member_no, j].mixing_length.CᴷRiʷ == new_CᴷRiʷ
    end

    #####
    ##### Tuples of closures
    #####

    @info "Testing closure_with_parameters on (CAVD, CATKE)"

    new_Cᴷuʳ = 1e3

    closures = tuple(
        closure_with_parameters(catke_closure, (; Cᴷuʳ = new_Cᴷuʳ)),
        closure_with_parameters(CAVD(), (; background_κz = 1.4))
    )

    @test closures[1].mixing_length.Cᴷuʳ == new_Cᴷuʳ
    @test closures[2].background_κz == 1.4

    @info "Testing closure_with_parameters on (CAVD, [CATKE, CATKE])"

    catkes = @suppress [CATKE(), CATKE()]
    closures = tuple(CAVD(background_κz=1.0), catkes)

    new_Cᴷuʳ = 1e3
    new_CᴷRiʷ = 100.0
    ensemble_member_no = 2
    old_CᴷRiʷ = catkes[1].mixing_length.CᴷRiʷ
    update_closure_ensemble_member!(closures, ensemble_member_no, (; CᴷRiʷ = new_CᴷRiʷ))

    @test closures[1].background_κz == 1.0
    @test closures[2][1].mixing_length.CᴷRiʷ == old_CᴷRiʷ
    @test closures[2][ensemble_member_no].mixing_length.CᴷRiʷ == new_CᴷRiʷ
end


text:
module CometLogger

using ColorTypes
using PyCall
using Requires
using FileIO

export CLogger

const comet = PyNULL()

function __init__()
    copy!(comet, pyimport("comet_ml"))

    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        include("plots.jl")
    end
    @require ImageIO="82e4d734-157c-48bb-816b-45c225c6df19" begin
        include("image.jl")
    end
end

using Base.CoreLogging: CoreLogging, AbstractLogger, LogLevel, Info, handle_message, shouldlog, min_enabled_level, catch_exceptions

# CometLogger
mutable struct CLogger <: AbstractLogger
    cexp::PyObject
    step_increment::Int
    global_step::Int
    min_level::LogLevel
end

function CLogger(; min_level=Info, step_increment=1, start_step=0, project_name=nothing, workspace=nothing, kwargs...)
    cexp = comet.Experiment(project_name=project_name, workspace=workspace, kwargs...)
    CLogger(cexp, step_increment, start_step, min_level)
end

increment_step!(lg::CLogger, Δ_Step) = lg.global_step += Δ_Step

add_tag!(lg::CLogger, tag::String) = lg.cexp.add_tag(tag)

"""
    function log_metric(lg::CLogger, name::AbstractString, value::Real; step::Int=nothing, epoch::Int=nothing)

Logs general scalar metrics.
"""
function log_metric(lg::CLogger, name::AbstractString, value::Real; step=nothing, epoch=nothing)
    lg.cexp.log_metric(name, value, step=step, epoch=epoch)
end

const log_value = log_metric

function log_parameter(lg::CLogger, name::AbstractString, value; step=nothing)
    lg.cexp.log_parameter(name, value, step=step)
end

function log_parameters(lg::CLogger, params::Dict{String,Any}; prefix=nothing, step=nothing)
    lg.cexp.log_parameters(params, prefix=prefix, step=step)
end

"""
    function log_text(lg::CLogger, name::AbstractString, value::AbstractString; step::Int=nothing, metadata::Dict{String,Any}=nothing)

Logs the text.
"""
function log_text(lg::CLogger, name::AbstractString, value; step=noting, metadata=nothing)
    lg.cexp.log_text(name, value, step=step, metadata=metadata)
end

"""
"""
function log_curve(lg::CLogger, name::AbstractString, x::AbstractVector, y::AbstractVector; overwrite=false, step=nothing)
    @assert length(x) == length(y)
    lg.cexp.log_curve(name, x, y, overwrite=overwrite, step=step)
end

"""
    function log_image(lg::CLogger, name::AbstractString, obj::AbstractString; step=nothing, kwargs...)

Logs the image from a filepath.
"""
function log_image(lg::CLogger, name::AbstractString, obj::AbstractString; step=nothing, kwargs...)
    if !ispath(obj)
        @warn "$obj not a path to png"
        return
    end
    lg.cexp.log_image(obj, name, step=step, kwargs...)
end

# AbstractLogger interface
# Kind of mimicking TensorBoardLogger.jl

CoreLogging.catch_exceptions(lg::CLogger) = false

CoreLogging.min_enabled_level(lg::CLogger) = lg.min_level

# For now, log everything that is above the lg.min_level
CoreLogging.shouldlog(lg::CLogger, level, _module, group, id) = true

function preprocess(name, val::T, data) where {T}
    if isstructtype(T)
        fn = logable_propertynames(val)
        for f=fn
            prop = getfield(val, f)
            preprocess(name*"/$f", prop, data)
        end
    else
        # If we do not know how to serialize a type, then
        # it will be simply logged as text
        push!(data, name=>val)
    end
    data
end

"""
    logable_propertynames(val::Any)

Returns a tuple with the name of the fields of the structure `val` that
should be logged to Comet.ml. This function should be overridden when
you want Comet.ml to ignore some fields in a structure when logging
it. The default behaviour is to return the  same result as `propertynames`.
See also: [`Base.propertynames`](@ref)
"""
logable_propertynames(val::Any) = propertynames(val)

## Default unpacking of key-value dictionaries
function preprocess(name, dict::AbstractDict, data)
    for (key, val) in dict
        # convert any key into a string, via interpolating it
        preprocess("$name/$key", val, data)
    end
    return data
end

# Split complex numbers into real/complex pairs
preprocess(name, val::Complex, data) = push!(data, name*"/re"=>real(val), name*"/im"=>imag(val))

function preprocess(name, (x, y)::Tuple{AbstractVector,AbstractMatrix}, data)
    # Assuming dims x T
    for i in 1:size(y, 1)
        push!(data, name*"/$i"=>(x, y[i, :]))
    end
end

process(lg::CLogger, name::AbstractString, obj::Real, step::Int) = log_metric(lg, name, obj; step=step)
process(lg::CLogger, name::AbstractString, obj::Tuple{AbstractVector,AbstractVector}, step::Int) = log_curve(lg, name, obj...; step=step)

function CoreLogging.handle_message(lg::CLogger, level, message, _module, group,
    id, file, line; kwargs...)

    i_step = lg.step_increment # :log_step_increment default value

    if !isempty(kwargs)
        data = Vector{Pair{String,Any}}()
        # ∀ (k-v) pairs, decompose values into objects that can be serialized
        for (key,val) in pairs(kwargs)
            # special key describing step increment
            if key == :log_step_increment
                i_step = val
                continue
            end
            preprocess(message*"/$key", val, data)
        end
        iter = increment_step!(lg, i_step)
        for (name, val) in data
            process(lg, name, val, iter)
        end
    end
end

Base.show(io::IO, cl::CLogger) = begin
	str  = "CLogger(\"$(cl.cexp.project_name)\", min_level=$(cl.min_level), "*
		   "current_step=$(cl.global_step))"*
           "\n@ $(cl.cexp.url)"
    Base.print(io, str)
end

Base.close(cl::CLogger) = cl.cexp.end()

end


text:
module GoogleProtoBuf
    using ..ProtoBuf
    import ..ProtoBuf: meta
    include("descriptor_pb.jl")
end
module GoogleProtoBufCompiler
    using ..ProtoBuf
    import ..ProtoBuf: meta
    include("plugin_pb.jl")
end


text:
module FMSynth

using Base: promote_op

export Oscillator,
       next!,
       setfrequency!,
       sine_wave,
       square_wave,
       saw_wave,
       triangle_wave,
       Mixer,
       Volume,
       FrequencyModulator,
       SynthesizerSource,
       PitchController

#####
##### `SAMPLE_RATE`
#####

const SAMPLE_RATE = Ref{Int}(44_100)

set_sample_rate!(sr::Int) = SAMPLE_RATE[] = sr

#####
##### `AbstractSynthComponent`
#####

abstract type AbstractSynthComponent end

function output_type end

output_type(::T) where {T<:AbstractSynthComponent} = output_type(T)

include("abstract_oscillator.jl")
include("abstract_synth_controller.jl")
include("components_util.jl")

include("sampled_signals.jl")

end