text:
# A modified next reaction method
struct MNRM <: AbstractAggregatorAlgorithm end
mutable struct MNRMJumpAggregation{T,S,F1,F2,RNG,DEPGR} <: AbstractSSAJumpAggregator
  next_jump::Int
  prev_jump::Int
  next_jump_time::T
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  ma_jumps::S
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
  rng::RNG
  internal_waitingtimes::Vector{T}
  dep_gr::DEPGR
end
function MNRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; num_specs, internal_waitingtimes, dep_graph=nothing, kwargs...) where {T,S,F1,F2,RNG}

    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the Next Reaction Method (NRM) algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph

        # make sure each jump depends on itself
        add_self_dependencies!(dg)
    end
    # pq = MutableBinaryMinHeap{T}()
    MNRMJumpAggregation{T,S,F1,F2,RNG,typeof(dg)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng, internal_waitingtimes, dg)
end

function aggregate(aggregator::MNRM, u, p, t, end_time, constant_jumps, ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using tuples
    rates, affects! = get_jump_info_tuples(constant_jumps)
    num_reactions = get_num_majumps(ma_jumps) + length(constant_jumps)
    internal_waitingtimes = zeros(typeof(t),num_reactions)
    build_jump_aggregation(MNRMJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, save_positions, rng; num_specs = length(u), internal_waitingtimes = internal_waitingtimes, kwargs...)
end

function initialize!(p::MNRMJumpAggregation, integrator, u, params, t)
    fill_rates_and_get_times!(p, u, params, t)
    nothing
end
@inline function execute_jumps!(p::MNRMJumpAggregation, integrator, u, params, t)
    update_state!(p, integrator, u)
    nothing
end
# calculate the next jump / jump time
function generate_jumps!(p::MNRMJumpAggregation, integrator, u, params, t)
    update_internal_times!(p, u, params, t)
    nothing
end

@inline function update_internal_times!(p::MNRMJumpAggregation, u, params, t)
    @unpack internal_waitingtimes, cur_rates, rates, ma_jumps, rng = p
    dep_rxs = p.dep_gr[copy(p.next_jump)]
    num_majumps = get_num_majumps(ma_jumps)
    @inbounds for rx in dep_rxs
        cur_rates[rx] = calculate_jump_rate(ma_jumps,num_majumps,rates,u,partial_path,t,rx)
    end
    pqdata = [internal_waitingtimes[i]/cur_rates[i] for i in eachindex(cur_rates)]
    # ttnj, next_jump = top_with_handle(MutableBinaryMinHeap(pqdata))
    ttnj, p.next_jump = findmin(pqdata)
    @fastmath p.next_jump_time = t + ttnj
    # update internal time
    @inbounds for rx in eachindex(cur_rates)
        internal_waitingtimes[rx] -= ttnj*cur_rates[rx]
    end
    internal_waitingtimes[p.next_jump] = randexp(rng)
    nothing
end

# fill the propensity rates
function fill_rates_and_get_times!(p::MNRMJumpAggregation,  u, params,  t)
    @unpack internal_waitingtimes, cur_rates, rates, ma_jumps, rng = p
    pqdata = Vector{typeof(t)}(undef,length(cur_rates))
    @inbounds for i in 1:get_num_majumps(ma_jumps)
        cur_rates[i] = evalrxrate(u, i, ma_jumps)
        pqdata[i] = randexp(rng)/cur_rates[i]
    end
    # constant rates
    rates = p.rates
    idx   = get_num_majumps(ma_jumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        pqdata[idx] = randexp(rng)/cur_rates[idx]
        idx += 1
    end
    ttnj, p.next_jump = findmin(pqdata)
    p.next_jump_time = t +ttnj
    @inbounds for rx in eachindex(cur_rates)
        internal_waitingtimes[rx] -= ttnj*cur_rates[rx]
    end
    internal_waitingtimes[p.next_jump] = randexp(rng)
end

# Change below


text:
"""
    set_geometry(cascade, ::Type{T}; kwargs...) where T<:Geometry
Given a cascade and geometry, return a cascade with an updated type.
"""
function set_geometry(cascade, ::Type{Violin}, args...; kwargs...)
    return _set_geometry(cascade, Violin, Poly, Coloring, args...; style=:fill, kwargs...)
end

function set_geometry(cascade, ::Type{Vertical}, args...; kwargs...)
    return _set_geometry(cascade, Vertical, Box, Blending, args...;
        style=:fill,
        subdivide=true,
        space=true,
        kwargs...,
    )
end

function set_geometry(cascade, ::Type{Horizontal}, args...; usegradient=missing, kwargs...)
    # !!!! want to show gradient if there's only one cascade element. Is this working??
    usegradient = coalesce(usegradient, length(cascade)==1)

    return _set_geometry(cascade, Horizontal, Box, usegradient ? Blending : Coloring, args...;
        alpha=length(cascade),
        style=:fill,
        subdivide=false,
        space=true,
        kwargs...,
    )
end

function set_geometry(cascade, ::Type{Parallel}, args...; slope::Bool=true, kwargs...)
    return _set_geometry(cascade, Parallel, Line, Coloring, args...;
        quantile = convert(Float64, slope),
        alpha = length(cascade),
        factor = 0.5,
        style = :stroke,
        subdivide = false,
        space = false,
        kwargs...,
    )
end

"""
    _set_geometry(cascade::Cascade{Data}, Geometry, Shape, Color; kwargs...)
"""
function _set_geometry(cascade::Cascade{Data}, Geometry, Shape, Color, args...;
    colorcycle::Bool = false,
    kwargs...,
)
    pos = scale_for(cascade, Geometry, args...; kwargs...)
    sgn = sign(cascade)

    label = get_label.(collect_data(cascade))
    shape = vectorize.(_define_from(Shape, Color, pos, sgn; kwargs...))

    data = Geometry.(label, shape, length(cascade))

    return Cascade(
        start = set_hue!(first(data), "black"),
        stop = set_hue!(last(data), "black"),
        steps = colorcycle ? set_hue!.(data[2:end-1], cascade.permutation) : data[2:end-1],
        permutation = cascade.permutation,
        correlation = cascade.correlation,
        ispermuted = cascade.ispermuted,
        iscorrelated = cascade.iscorrelated,
    )
end

"""
    set_hue!(x, h)
"""
function set_hue!(x::Coloring, h; kwargs...)
    x.hue = define_from(Luxor.RGB, h; kwargs...)
    return x
end

set_hue!(x::Handle, h; kwargs...) = begin set_hue!(x.shape, h; kwargs...); return x end
set_hue!(x::Blending, h) = begin x.hue = _define_gradient(h); return x end
set_hue!(x::T, h) where T<:Geometry = begin set_hue!(x.shape, h); return x end
set_hue!(x::T, h) where T<:Shape = begin set_hue!(x.color, h); return x end
set_hue!(x::Vector{T}, h) where T<:Shape = begin set_hue!.(x, h); return x end

"""
    set_position!(handle::Handle)
This method calculates the position of a Handle's shape and label such that the shape is
centered at (x0,y0) and fits within the bounds of `(2*dx, 2*dy)`. The label is aligned
middle/left at `(x0+dx+space, y0)`
"""
function set_position!(handle::Handle;
    x0 = WIDTH-100,
    y0 = 3*SEP,
    dx = 2*SEP,
    dy = SEP,
    space = SEP,
    idx = 1,
    kwargs...,
)
    y0 = y0 + (2*dy+space)*(idx-1)

    setproperty!(handle.shape, :position, (
        Luxor.Point(x0-dx, y0+dy),
        Luxor.Point(x0+dx, y0-dy),
    ))

    setproperty!(handle.label, :position, Luxor.Point(x0+dx+space, y0))

    return handle
end

"""
    get_shape(cascade; kwargs...)
This function returns the shape used in `cascade`, with its position updated, as defined by
keyword arguments.
"""
get_shape(cascade::Cascade{T}) where T<:Geometry = copy(first(cascade.start.shape))
get_shape(cascade::Cascade{Violin}) = copy(cascade.start.shape)

ext: jl


text:
export Transform

struct Transform{T<:Number}
    m::SMatrix{4,4,T}
    m_inv::SMatrix{4,4,T}

    function Transform{T}(args::NTuple{16, T}) where T
        m = SMatrix{4, 4}(args)
        new{T}(m, inv(m))
    end

    Transform{T}(m) where T = new{T}(m, inv(m))

    Transform{T}() where T =
        new{T}(one(SMatrix{4,4}), one(SMatrix{4,4}))

    Transform{T}(m, m_inv) where T = new{T}(m, m_inv)
end

Transform(m::SMatrix{4,4,T}) where T = Transform{T}(m)
Transform(m::SMatrix{4,4,T}, m_inv::SMatrix{4,4,U}) where {T,U} = Transform{promote_type(T,U)}(m, m_inv)


text:
"""
    SeawaterBuoyancy{FT, EOS, T, S} <: AbstractBuoyancy{EOS}

Buoyancy model for seawater. `T` and `S` are either `nothing` if both
temperature and salinity are active, or of type `FT` if temperature
or salinity are constant, respectively.
"""
struct SeawaterBuoyancy{FT, EOS, T, S} <: AbstractBuoyancy{EOS}
    gravitational_acceleration :: FT
             equation_of_state :: EOS
          constant_temperature :: T
             constant_salinity :: S
end

required_tracers(::SeawaterBuoyancy) = (:T, :S)
required_tracers(::SeawaterBuoyancy{FT, EOS, <:Nothing, <:Number}) where {FT, EOS} = (:T,) # active temperature only
required_tracers(::SeawaterBuoyancy{FT, EOS, <:Number, <:Nothing}) where {FT, EOS} = (:S,) # active salinity only

"""
    SeawaterBuoyancy([FT=Float64;] gravitational_acceleration = g_Earth,
                                  equation_of_state = LinearEquationOfState(FT),
                                  constant_temperature = false, constant_salinity = false)

Returns parameters for a temperature- and salt-stratified seawater buoyancy model
with a `gravitational_acceleration` constant (typically called 'g'), and an
`equation_of_state` that related temperature and salinity (or conservative temperature
and absolute salinity) to density anomalies and buoyancy. If either `temperature` or `salinity`
are specified, buoyancy is calculated
"""
function SeawaterBuoyancy(                        FT = Float64;
                          gravitational_acceleration = g_Earth,
                                   equation_of_state = LinearEquationOfState(FT),
                                constant_temperature = nothing,
                                   constant_salinity = nothing)

    # Input validation: convert constant_temperature or constant_salinity = true to zero(FT).
    # This method of specifying constant temperature or salinity in a SeawaterBuoyancy model
    # should only be used with a LinearEquationOfState where the constant value of either temperature
    # or sailnity is irrelevant.
    constant_temperature = constant_temperature === true ? zero(FT) : constant_temperature
    constant_salinity = constant_salinity === true ? zero(FT) : constant_salinity

    return SeawaterBuoyancy{FT, typeof(equation_of_state), typeof(constant_temperature), typeof(constant_salinity)}(
                            gravitational_acceleration, equation_of_state, constant_temperature, constant_salinity)
end

const TemperatureSeawaterBuoyancy = SeawaterBuoyancy{FT, EOS, <:Nothing, <:Number} where {FT, EOS}
const SalinitySeawaterBuoyancy = SeawaterBuoyancy{FT, EOS, <:Number, <:Nothing} where {FT, EOS}

@inline get_temperature_and_salinity(::SeawaterBuoyancy, C) = C.T, C.S
@inline get_temperature_and_salinity(b::TemperatureSeawaterBuoyancy, C) = C.T, b.constant_salinity
@inline get_temperature_and_salinity(b::SalinitySeawaterBuoyancy, C) = b.constant_temperature, C.S

"""
    ∂x_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the x-derivative of buoyancy for temperature and salt-stratified water,

```math
∂_x b = g ( α ∂_x Θ - β ∂_x sᴬ ) ,
```

where `g` is gravitational acceleration, `α` is the thermal expansion
coefficient, `β` is the haline contraction coefficient, `Θ` is
conservative temperature, and `sᴬ` is absolute salinity.

Note: In Oceananigans, `model.tracers.T` is conservative temperature and
`model.tracers.S` is absolute salinity.

Note that `∂x_Θ`, `∂x_sᴬ`, `α`, and `β` are all evaluated at cell interfaces in `x`
and cell centers in `y` and `z`.
"""
@inline function ∂x_b(i, j, k, grid, b::SeawaterBuoyancy, C)
    Θ, sᴬ = get_temperature_and_salinity(b, C)
    return b.gravitational_acceleration * (
           thermal_expansionᶠᶜᶜ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂xᶠᵃᵃ(i, j, k, grid, Θ)
        - haline_contractionᶠᶜᶜ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂xᶠᵃᵃ(i, j, k, grid, sᴬ) )
end

"""
    ∂y_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the y-derivative of buoyancy for temperature and salt-stratified water,

```math
∂_y b = g ( α ∂_y Θ - β ∂_y sᴬ ) ,
```

where `g` is gravitational acceleration, `α` is the thermal expansion
coefficient, `β` is the haline contraction coefficient, `Θ` is
conservative temperature, and `sᴬ` is absolute salinity.

Note: In Oceananigans, `model.tracers.T` is conservative temperature and
`model.tracers.S` is absolute salinity.

Note that `∂y_Θ`, `∂y_sᴬ`, `α`, and `β` are all evaluated at cell interfaces in `y`
and cell centers in `x` and `z`.
"""
@inline function ∂y_b(i, j, k, grid, b::SeawaterBuoyancy, C)
    Θ, sᴬ = get_temperature_and_salinity(b, C)
    return b.gravitational_acceleration * (
           thermal_expansionᶜᶠᶜ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂yᵃᶠᵃ(i, j, k, grid, Θ)
        - haline_contractionᶜᶠᶜ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂yᵃᶠᵃ(i, j, k, grid, sᴬ) )
end

"""
    ∂z_b(i, j, k, grid, b::SeawaterBuoyancy, C)

Returns the vertical derivative of buoyancy for temperature and salt-stratified water,

```math
∂_z b = N^2 = g ( α ∂_z Θ - β ∂_z sᴬ ) ,
```

where `g` is gravitational acceleration, `α` is the thermal expansion
coefficient, `β` is the haline contraction coefficient, `Θ` is
conservative temperature, and `sᴬ` is absolute salinity.

Note: In Oceananigans, `model.tracers.T` is conservative temperature and
`model.tracers.S` is absolute salinity.

Note that `∂z_Θ`, `∂z_sᴬ`, `α`, and `β` are all evaluated at cell interfaces in `z`
and cell centers in `x` and `y`.
"""
@inline function ∂z_b(i, j, k, grid, b::SeawaterBuoyancy, C)
    Θ, sᴬ = get_temperature_and_salinity(b, C)
    return b.gravitational_acceleration * (
           thermal_expansionᶜᶜᶠ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂zᵃᵃᶠ(i, j, k, grid, Θ)
        - haline_contractionᶜᶜᶠ(i, j, k, grid, b.equation_of_state, Θ, sᴬ) * ∂zᵃᵃᶠ(i, j, k, grid, sᴬ) )
end

@inline function buoyancy_perturbation(i, j, k, grid, b::SeawaterBuoyancy{FT, <:AbstractNonlinearEquationOfState}, C) where FT
    Θ, sᴬ = get_temperature_and_salinity(b, C)
    return - b.gravitational_acceleration * ρ′(i, j, k, grid, b.equation_of_state, Θ, sᴬ) / b.equation_of_state.ρ₀
end


text:
"""
Tools for exploratory data analysis and debugging.
"""
module CelesteEDA

using DataFrames

import ..DeterministicVI:
    load_source_brightnesses, accumulate_source_pixel_brightness!, ElboArgs
import ..Model:
    linear_world_to_pix, load_bvn_mixtures, ids_names, ids, CatalogEntry,
    populate_fsm_vecs!
import ..SensitiveFloats.clear!
import ..DeterministicVIImagePSF:
    FSMSensitiveFloatMatrices
using ..DeterministicVIImagePSF:
    clear_brightness!, populate_star_fsm_image!, populate_gal_fsm_image!,
    accumulate_source_image_brightness!, load_gal_bvn_mixtures

import Base.print
function print(ce::CatalogEntry)
    for field in fieldnames(ce)
      println(field, ": ", getfield(ce, field))
    end
end

function print_vp(vp::Vector{Float64})
    df = DataFrame(ids=ids_names)
    s = 1
    df[Symbol(string("v", s))] = vp
    println(df)
end

function source_pixel_location(ea::ElboArgs, s::Int, n::Int)
    p = ea.patches[s, n]
    pix_loc = linear_world_to_pix(
        p.wcs_jacobian,
        p.center,
        p.pixel_center,
        ea.vp[s][ids.u])
    return pix_loc - p.bitmap_offset
end

function render_source(ea::ElboArgs, s::Int, n::Int;
                       include_epsilon=true, field=:E_G,
                       include_iota=true)
    local p = ea.patches[s, n]
    local image = fill(NaN, size(p.active_pixel_bitmap))
    local sbs = load_source_brightnesses(
        ea, calculate_gradient=false, calculate_hessian=false)
    star_mcs, gal_mcs = load_bvn_mixtures(ea.S, ea.patches,
                                ea.vp, ea.active_sources,
                                ea.psf_K, n,
                                calculate_gradient=false,
                                calculate_hessian=false)

    p = ea.patches[s,n]

    H2, W2 = size(p.active_pixel_bitmap)
    for w2 in 1:W2, h2 in 1:H2
        # (h2, w2) index the local patch, while (h, w) index the image
        h = p.bitmap_offset[1] + h2
        w = p.bitmap_offset[2] + w2

        if !p.active_pixel_bitmap[h2, w2]
            continue
        end

        clear!(ea.elbo_vars.E_G)
        clear!(ea.elbo_vars.var_G)

        populate_fsm_vecs!(ea.elbo_vars.bvn_derivs,
                           ea.elbo_vars.fs0m_vec,
                           ea.elbo_vars.fs1m_vec,
                           ea.patches,
                           ea.active_sources,
                           ea.num_allowed_sd,
                           n, h, w,
                           gal_mcs, star_mcs)

        accumulate_source_pixel_brightness!(
            ea.elbo_vars, ea, ea.elbo_vars.E_G, ea.elbo_vars.var_G,
            ea.elbo_vars.fs0m_vec[s], ea.elbo_vars.fs1m_vec[s],
            sbs[s], ea.images[n].b, s, false)

        if field == :E_G
            image[h2, w2] = ea.elbo_vars.E_G.v[]
        elseif field == :fs0m
            image[h2, w2] = ea.elbo_vars.fs0m_vec[s].v[]
        elseif field == :fs1m
            image[h2, w2] = ea.elbo_vars.fs1m_vec[s].v[]
        else
            error("Unknown field ", field)
        end
        if include_iota
            image[h2, w2] *= ea.images[n].iota_vec[h]
        end
        if include_epsilon
            image[h2, w2] += ea.images[n].epsilon_mat[h, w]
        end
    end

    return image
end

function render_source_fft(
    ea::ElboArgs,
    fsm_mat::Matrix{FSMSensitiveFloatMatrices},
    s::Int, n::Int;
    include_epsilon=true,
    field=:E_G, include_iota=true)

    local p = ea.patches[s, n]
    local image = fill(NaN, size(p.active_pixel_bitmap))
    local sbs = load_source_brightnesses(
        ea, calculate_gradient=false, calculate_hessian=false)
    local fsms = fsm_mat[s, n]

    local gal_mcs = load_gal_bvn_mixtures(
            ea.S, ea.patches, ea.vp, ea.active_sources, n,
            calculate_gradient=false,
            calculate_hessian=false);

    clear_brightness!(fsms)
    populate_star_fsm_image!(
        ea, s, n, fsms.psf, fsms.fs0m_conv,
        fsms.h_lower, fsms.w_lower, fsms.kernel_fun, fsms.kernel_width)
    populate_gal_fsm_image!(ea, s, n, gal_mcs, fsms)
    accumulate_source_image_brightness!(ea, s, n, fsms, sbs[s])

    H2, W2 = size(p.active_pixel_bitmap)
    for w2 in 1:W2, h2 in 1:H2
        # (h2, w2) index the local patch, while (h, w) index the image
        h = p.bitmap_offset[1] + h2
        w = p.bitmap_offset[2] + w2

        if !p.active_pixel_bitmap[h2, w2]
            continue
        end

        h_fsm = h - fsms.h_lower + 1
        w_fsm = w - fsms.w_lower + 1

        image[h2, w2] = getfield(fsms, field)[h_fsm, w_fsm].v[]
        if include_iota
            image[h2, w2] *= ea.images[n].iota_vec[h]
        end
        if include_epsilon
            image[h2, w2] += ea.images[n].epsilon_mat[h, w]
        end

    end

    return deepcopy(image)
end

function show_source_image(ea::ElboArgs, s::Int, n::Int)
    p = ea.patches[s, n]
    H2, W2 = size(p.active_pixel_bitmap)
    image = fill(NaN, H2, W2);
    for w2 in 1:W2, h2 in 1:H2
        # (h2, w2) index the local patch, while (h, w) index the image
        h = p.bitmap_offset[1] + h2
        w = p.bitmap_offset[2] + w2

        if !p.active_pixel_bitmap[h2, w2]
            continue
        end

        image[h2, w2] = ea.images[n].pixels[h, w]
    end
    return image
end

end


text:
function DoSmooth(a, l)
	if l % 2 == 1
		l = l - 1
	else
		l = l - 2
	end
	b = cumsum(a)
	Z = zeros(size(a))
	tmpLen = length(a)
	for i=1:tmpLen
		len = min(l/2,tmpLen-i,i-1)
		Z[i]=(b[i+len]-b[i-len]+a[i-len])/(2*len+1)
	end
	return Z
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
@doc raw"""
    airydisk([T=Float64], point; x, y, fwhm, ratio=0, amp=1, theta=0, bkg=0)
    airydisk([T=Float64], px, py; x, y, fwhm, ratio=0, amp=1, theta=0, bkg=0)

An unnormalized Airy disk. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. If `theta` is given, the PSF will be rotated by `theta` degrees counter-clockwise from the x-axis. If `bkg` is given it will be added as a scalar to the PSF.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). Keep in mind that `theta` has no effect for isotropic distributions and is degenerate with the `fwhm` parameters (i.e., theta=90 is the same as reversing the `fwhm` tuple)

If `ratio` is supplied, this will be the Airy pattern for a centrally-obscured aperture (e.g., a Newtonian telescope). This has a slightly expanded functional form, and in general the central Airy disk will be smaller and the first Airy ring will be brighter.

# Functional form

The Airy disk is a distribution over the radius `r` (the square-Euclidean distance)

```
f(x | x̂, FWHM) = [ 2J₁(q) / q ]^2
```
where `J₁` is the first-order Bessel function of the first kind and
```
q ≈ π * r * D/ λ ≈ π * r / (0.973 * FWHM)
```

If user a non-zero central obscuration via `ratio`, the functional form becomes

```
f(x | x̂, FWHM, ϵ) = [ 2J₁(q) / q - 2ϵJ₁(ϵq) / q ]^2 / (1 - ϵ^2)^2
```
where `ϵ` is the ratio (0 ≤ `ϵ` < 1).
"""
airydisk(T, px, py; x, y, fwhm, ratio=0, amp=one(T), theta=0, bkg=0) = convert(T, _airydisk(px, py, x, y, fwhm, ratio, amp, theta, bkg))

# factor for scaling radius in terms of the fwhm
const AIRY_PRE = π / 0.973

# isotropic
function _airydisk(px, py, x, y, fwhm, ratio, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    !iszero(theta) && @warn "isotropic airydisk is not affected by non-zero rotation angle $theta"
    # unnormalized airydisk
    r = sqrt(dx^2 + dy^2)
    # short-circuit
    iszero(r) && return amp / (1 - ratio^2)^2 + background
    q = AIRY_PRE * r / fwhm
    I1 = 2 * besselj1(q) / q
    if iszero(ratio)
        return amp * I1^2 + background
    end
    I2 = 2 * ratio * besselj1(q * ratio) / q
    return amp * ((I1 - I2) / (1 - ratio^2))^2 + background
end

# bivariate
function _airydisk(px, py, x, y, fwhm::BivariateLike, ratio, amp, theta, background)
    # find offset from center
    dx = px - x
    dy = py - y
    # rotate
    if !iszero(theta)
        dx, dy = rotate_point(dx, dy, theta)
    end
    # unnormalized airy disk
    fwhmx, fwhmy = fwhm
    q = AIRY_PRE * sqrt((dx / fwhmx)^2 + (dy / fwhmy)^2)
    # short-circuit
    iszero(q) && return amp / (1 - ratio^2)^2 + background
    I1 = 2 * besselj1(q) / q
    if iszero(ratio)
        return amp * I1^2 + background
    end
    I2 = 2 * ratio * besselj1(q * ratio) / q
    return amp * ((I1 - I2) / (1 - ratio^2))^2 + background
end


text:
# Coordinate Reference System Objects
# (has keys "type" and "properties")
# TODO: Handle full CRS spec
const CRS = Dict{String,Any}

# Bounding Boxes
# The value of the bbox member must be a 2*n array,
# where n is the number of dimensions represented in the contained geometries,
# with the lowest values for all axes followed by the highest values.

# The axes order of a bbox follows the axes order of geometries.
# In addition, the coordinate reference system for the bbox is assumed to match
# the coordinate reference system of the GeoJSON object of which it is a member.
const BBox = Vector{Float64}

const Position = Vector{Float64}
# (x, y, [z, ...]) - meaning of additional elements undefined.
# In an object's contained geometries, Positions must have uniform dimensions.
geotype(::Position) = :Position
xcoord(p::Position) = p[1]
ycoord(p::Position) = p[2]
zcoord(p::Position) = hasz(p) ? p[3] : zero(T)
hasz(p::Position) = length(p) >= 3
coordinates(obj::Position) = obj

coordinates(obj::Vector{Position}) = obj
coordinates(obj::Vector{T}) where {T <: AbstractPosition} =                    Position[map(coordinates, obj)...]
coordinates(obj::Vector{T}) where {T <: AbstractPoint} =                       Position[map(coordinates, obj)...]

coordinates(obj::Vector{Vector{Position}}) = obj
coordinates(obj::Vector{Vector{T}}) where {T <: AbstractPosition} =            Vector{Position}[map(coordinates, obj)...]
coordinates(obj::Vector{Vector{T}}) where {T <: AbstractPoint} =               Vector{Position}[map(coordinates, obj)...]
coordinates(obj::Vector{T}) where {T <: AbstractLineString} =                  Vector{Position}[map(coordinates, obj)...]

coordinates(obj::Vector{Vector{Vector{Position}}}) = obj
coordinates(obj::Vector{Vector{Vector{T}}}) where {T <: AbstractPosition} =    Vector{Vector{Position}}[map(coordinates, obj)...]
coordinates(obj::Vector{Vector{Vector{T}}}) where {T <: AbstractPoint} =       Vector{Vector{Position}}[map(coordinates, obj)...]
coordinates(obj::Vector{Vector{T}}) where {T <: AbstractLineString} =          Vector{Vector{Position}}[map(coordinates, obj)...]
coordinates(obj::Vector{T}) where {T <: AbstractPolygon} =                     Vector{Vector{Position}}[map(coordinates, obj)...]

mutable struct Point <: AbstractPoint
    coordinates::Position
end
Point(x::Float64,y::Float64) = Point([x,y])
Point(x::Float64,y::Float64,z::Float64) = Point([x,y,z])
Point(point::AbstractPosition) = Point(coordinates(point))
Point(point::AbstractPoint) = Point(coordinates(point))

mutable struct MultiPoint <: AbstractMultiPoint
    coordinates::Vector{Position}
end
MultiPoint(point::Position) = MultiPoint(Position[point])
MultiPoint(point::AbstractPosition) = MultiPoint(Position[coordinates(point)])
MultiPoint(point::AbstractPoint) = MultiPoint(Position[coordinates(point)])

MultiPoint(points::Vector{T}) where {T <: AbstractPosition} = MultiPoint(coordinates(points))
MultiPoint(points::Vector{T}) where {T <: AbstractPoint} = MultiPoint(coordinates(points))
MultiPoint(points::AbstractMultiPoint) = MultiPoint(coordinates(points))
MultiPoint(line::AbstractLineString) = MultiPoint(coordinates(line))

mutable struct LineString <: AbstractLineString
    coordinates::Vector{Position}
end
LineString(points::Vector{T}) where {T <: AbstractPosition} = LineString(coordinates(points))
LineString(points::Vector{T}) where {T <: AbstractPoint} = LineString(coordinates(points))
LineString(points::AbstractMultiPoint) = LineString(coordinates(points))
LineString(line::AbstractLineString) = LineString(coordinates(line))

mutable struct MultiLineString <: AbstractMultiLineString
    coordinates::Vector{Vector{Position}}
end
MultiLineString(line::Vector{Position}) = MultiLineString(Vector{Position}[line])
MultiLineString(line::Vector{T}) where {T <: AbstractPosition} = MultiLineString(Vector{Position}[coordinates(line)])
MultiLineString(line::Vector{T}) where {T <: AbstractPoint} = MultiLineString(Vector{Position}[coordinates(line)])
MultiLineString(line::AbstractLineString) = MultiLineString(Vector{Position}[coordinates(line)])

MultiLineString(lines::Vector{Vector{T}}) where {T <: AbstractPosition} = MultiLineString(coordinates(lines))
MultiLineString(lines::Vector{Vector{T}}) where {T <: AbstractPoint} = MultiLineString(coordinates(lines))
MultiLineString(lines::Vector{T}) where {T <: AbstractLineString} = MultiLineString(Vector{Position}[map(coordinates,lines)])
MultiLineString(lines::AbstractMultiLineString) = MultiLineString(coordinates(line))
MultiLineString(poly::AbstractPolygon) = MultiLineString(coordinates(poly))

mutable struct Polygon <: AbstractPolygon
    coordinates::Vector{Vector{Position}}
end
Polygon(line::Vector{Position}) = Polygon(Vector{Position}[line])
Polygon(line::Vector{T}) where {T <: AbstractPosition} = Polygon(Vector{Position}[coordinates(line)])
Polygon(line::Vector{T}) where {T <: AbstractPoint} = Polygon(Vector{Position}[coordinates(line)])
Polygon(line::AbstractLineString) = Polygon(Vector{Position}[coordinates(line)])

Polygon(lines::Vector{Vector{T}}) where {T <: AbstractPosition} = Polygon(coordinates(lines))
Polygon(lines::Vector{Vector{T}}) where {T <: AbstractPoint} = Polygon(coordinates(lines))
Polygon(lines::Vector{T}) where {T <: AbstractLineString} = Polygon(coordinates(lines))
Polygon(lines::AbstractMultiLineString) = Polygon(coordinates(lines))
Polygon(poly::AbstractPolygon) = Polygon(coordinates(poly))

mutable struct MultiPolygon <: AbstractMultiPolygon
    coordinates::Vector{Vector{Vector{Position}}}
end
MultiPolygon(line::Vector{Position}) = MultiPolygon(Vector{Vector{Position}}[Vector{Position}[line]])
MultiPolygon(line::Vector{T}) where {T <: AbstractPosition} = MultiPolygon(Vector{Vector{Position}}[Vector{Position}[coordinates(line)]])
MultiPolygon(line::Vector{T}) where {T <: AbstractPoint} = MultiPolygon(Vector{Vector{Position}}[Vector{Position}[coordinates(line)]])
MultiPolygon(line::AbstractLineString) = MultiPolygon(Vector{Vector{Position}}[Vector{Position}[coordinates(line)]])

MultiPolygon(poly::Vector{Vector{T}}) where {T <: AbstractPosition} = MultiPolygon(Vector{Vector{Position}}[coordinates(poly)])
MultiPolygon(poly::Vector{Vector{T}}) where {T <: AbstractPoint} = MultiPolygon(Vector{Vector{Position}}[coordinates(poly)])
MultiPolygon(poly::Vector{T}) where {T <: AbstractLineString} = MultiPolygon(Vector{Vector{Position}}[coordinates(poly)])
MultiPolygon(poly::AbstractMultiLineString) = MultiPolygon(Vector{Vector{Position}}[coordinates(poly)])
MultiPolygon(poly::AbstractPolygon) = MultiPolygon(Vector{Vector{Position}}[coordinates(poly)])

MultiPolygon(polys::Vector{Vector{Vector{T}}}) where {T <: AbstractPosition} = MultiPolygon(coordinates(polys))
MultiPolygon(polys::Vector{Vector{Vector{T}}}) where {T <: AbstractPoint} = MultiPolygon(coordinates(polys))
MultiPolygon(polys::Vector{Vector{T}}) where {T <: AbstractLineString} = MultiPolygon(coordinates(polys))
MultiPolygon(polys::Vector{T}) where {T <: AbstractPolygon} = MultiPolygon(coordinates(polys))
MultiPolygon(polys::AbstractMultiPolygon) = MultiPolygon(coordinates(polys))

for geom in (:MultiPolygon, :Polygon, :MultiLineString, :LineString, :MultiPoint, :Point)
    @eval coordinates(obj::$geom) = obj.coordinates
end

mutable struct GeometryCollection <: AbstractGeometryCollection
    geometries::Vector
end
geometries(collection::GeometryCollection) = collection.geometries

mutable struct Feature <: AbstractFeature
    geometry::Union{Nothing, AbstractGeometry}
    properties::Union{Nothing, Dict{String,Any}}
end
Feature(geometry::Union{Nothing,GeoInterface.AbstractGeometry}) = Feature(geometry, Dict{String,Any}())
Feature(properties::Dict{String,Any}) = Feature(nothing, properties)
geometry(feature::Feature) = feature.geometry
properties(feature::Feature) = feature.properties
bbox(feature::Feature) = get(feature.properties, "bbox", nothing)
crs(feature::Feature) = get(feature.properties, "crs", nothing)

mutable struct FeatureCollection{T <: AbstractFeature} <: AbstractFeatureCollection
    features::Vector{T}
    bbox::Union{Nothing, BBox}
    crs::Union{Nothing, CRS}
end
FeatureCollection(fc::Vector{T}) where {T <: AbstractFeature} = FeatureCollection(fc, nothing, nothing)
features(fc::FeatureCollection) = fc.features
bbox(fc::FeatureCollection) = fc.bbox
crs(fc::FeatureCollection) = fc.crs


text:
module LabelPropagation

using Graphs
export label_propagation, modularity, erwkpath, sigmaij, nmi, nvi, consensus, rankedge, rankedge1, nei_label_propagation, rankedgebynei
include("community_label_propagation.jl")
include("modularity.jl")
include("ERWKpath.jl")
include("sigmaij.jl")
include("nmi.jl")
include("nvi.jl")
include("consensus.jl")
include("rankedge.jl")
end # module


text:
"""
f-method ODE solver
"""
struct ForwardEuler <: ODESolver
  nls::NonlinearSolver
  dt::Float64
end

function solve_step!(uf::AbstractVector,
                     solver::ForwardEuler,
                     op::ODEOperator,
                     u0::AbstractVector,
                     t0::Real,
                     cache) # -> (uF,tF)

  if cache === nothing
    ode_cache = allocate_cache(op)
    vf = similar(u0)
    nl_cache = nothing
  else
    ode_cache, vf, nl_cache = cache
  end

  dt = solver.dt
  tf = t0+dt
  # The space should have the boundary conditions at tf
  ode_cache = update_cache!(ode_cache,op,t0)

  nlop = ForwardEulerNonlinearOperator(op,t0,dt,u0,ode_cache,vf)

  nl_cache = solve!(uf,solver.nls,nlop,nl_cache)

  cache = (ode_cache, vf, nl_cache)

  return (uf,tf,cache)

end

"""
Nonlinear operator that represents the Forward Euler nonlinear operator at a
given time step, i.e., A(t,u_n,(u_n+1-u_n)/dt)
"""
struct ForwardEulerNonlinearOperator <: NonlinearOperator
  odeop::ODEOperator
  tf::Float64
  dt::Float64
  u0::AbstractVector
  ode_cache
  vf::AbstractVector
end

function residual!(b::AbstractVector,op::ForwardEulerNonlinearOperator,x::AbstractVector)
  vf = op.vf
  vf = (x-op.u0)/op.dt
  residual!(b,op.odeop,op.tf,(op.u0,vf),op.ode_cache)
end

function jacobian!(A::AbstractMatrix,op::ForwardEulerNonlinearOperator,x::AbstractVector)
  uF = x
  vf = op.vf
  vf = (x-op.u0)/op.dt
  z = zero(eltype(A))
  fill_entries!(A,z)
  jacobian!(A,op.odeop,op.tf,(op.u0,vf),1,(1/op.dt),op.ode_cache)
end

function allocate_residual(op::ForwardEulerNonlinearOperator,x::AbstractVector)
  allocate_residual(op.odeop,x,op.ode_cache)
end

function allocate_jacobian(op::ForwardEulerNonlinearOperator,x::AbstractVector)
  allocate_jacobian(op.odeop,x,op.ode_cache)
end

function zero_initial_guess(op::ForwardEulerNonlinearOperator)
  x0 = similar(op.u0)
  fill!(x0,zero(eltype(x0)))
  x0
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

const TAG = :ReverseAD

"""
    _eval_hessian(
        d::NLPEvaluator,
        f::_FunctionStorage,
        H::AbstractVector{Float64},
        λ::Float64,
        offset::Int,
    )::Int

Evaluate the hessian matrix of the function `f` and store the result, scaled by
`λ`, in `H`, beginning at element `offset+1`. This function assumes that
`_reverse_mode(d, x)` has already been called.

Returns the number of non-zeros in the computed Hessian, which will be used to
update the offset for the next call.
"""
function _eval_hessian(
    d::NLPEvaluator,
    f::_FunctionStorage,
    H::AbstractVector{Float64},
    λ::Float64,
    offset::Int,
)::Int
    chunk = min(size(f.seed_matrix, 2), d.max_chunk)
    # As a performance optimization, skip dynamic dispatch if the chunk is 1.
    if chunk == 1
        return _eval_hessian_inner(d, f, H, λ, offset, Val(1))
    else
        return _eval_hessian_inner(d, f, H, λ, offset, Val(chunk))
    end
end

function _eval_hessian_inner(
    d::NLPEvaluator,
    ex::_FunctionStorage,
    H::AbstractVector{Float64},
    scale::Float64,
    nzcount::Int,
    ::Val{CHUNK},
) where {CHUNK}
    if ex.linearity == LINEAR
        @assert length(ex.hess_I) == 0
        return 0
    end
    T = ForwardDiff.Partials{CHUNK,Float64}  # This is our element type.
    Coloring.prepare_seed_matrix!(ex.seed_matrix, ex.rinfo)
    local_to_global_idx = ex.rinfo.local_indices
    input_ϵ_raw, output_ϵ_raw = d.input_ϵ, d.output_ϵ
    input_ϵ = _reinterpret_unsafe(T, input_ϵ_raw)
    output_ϵ = _reinterpret_unsafe(T, output_ϵ_raw)
    # Compute hessian-vector products
    num_products = size(ex.seed_matrix, 2) # number of hessian-vector products
    num_chunks = div(num_products, CHUNK)
    @assert size(ex.seed_matrix, 1) == length(local_to_global_idx)
    for k in 1:CHUNK:CHUNK*num_chunks
        for r in 1:length(local_to_global_idx)
            # set up directional derivatives
            @inbounds idx = local_to_global_idx[r]
            # load up ex.seed_matrix[r,k,k+1,...,k+CHUNK-1] into input_ϵ
            for s in 1:CHUNK
                input_ϵ_raw[(idx-1)*CHUNK+s] = ex.seed_matrix[r, k+s-1]
            end
            @inbounds output_ϵ[idx] = zero(T)
        end
        _hessian_slice_inner(d, ex, input_ϵ, output_ϵ, T)
        # collect directional derivatives
        for r in 1:length(local_to_global_idx)
            idx = local_to_global_idx[r]
            # load output_ϵ into ex.seed_matrix[r,k,k+1,...,k+CHUNK-1]
            for s in 1:CHUNK
                ex.seed_matrix[r, k+s-1] = output_ϵ_raw[(idx-1)*CHUNK+s]
            end
            @inbounds input_ϵ[idx] = zero(T)
        end
    end
    # leftover chunk
    remaining = num_products - CHUNK * num_chunks
    if remaining > 0
        k = CHUNK * num_chunks + 1
        for r in 1:length(local_to_global_idx)
            # set up directional derivatives
            @inbounds idx = local_to_global_idx[r]
            # load up ex.seed_matrix[r,k,k+1,...,k+remaining-1] into input_ϵ
            for s in 1:remaining
                # leave junk in the unused components
                input_ϵ_raw[(idx-1)*CHUNK+s] = ex.seed_matrix[r, k+s-1]
            end
            @inbounds output_ϵ[idx] = zero(T)
        end
        _hessian_slice_inner(d, ex, input_ϵ, output_ϵ, T)
        # collect directional derivatives
        for r in 1:length(local_to_global_idx)
            idx = local_to_global_idx[r]
            # load output_ϵ into ex.seed_matrix[r,k,k+1,...,k+remaining-1]
            for s in 1:remaining
                ex.seed_matrix[r, k+s-1] = output_ϵ_raw[(idx-1)*CHUNK+s]
            end
            @inbounds input_ϵ[idx] = zero(T)
        end
    end
    # TODO(odow): consider reverting to a view.
    output_slice = _UnsafeVectorView(nzcount, length(ex.hess_I), pointer(H))
    Coloring.recover_from_matmat!(
        output_slice,
        ex.seed_matrix,
        ex.rinfo,
        d.output_ϵ,
    )
    for i in 1:length(output_slice)
        output_slice[i] *= scale
    end
    return length(ex.hess_I)
end

function _hessian_slice_inner(d, ex, input_ϵ, output_ϵ, ::Type{T}) where {T}
    subexpr_forward_values_ϵ =
        _reinterpret_unsafe(T, d.subexpression_forward_values_ϵ)
    for i in ex.dependent_subexpressions
        subexpr = d.subexpressions[i]
        subexpr_forward_values_ϵ[i] = _forward_eval_ϵ(
            d,
            subexpr,
            _reinterpret_unsafe(T, subexpr.forward_storage_ϵ),
            _reinterpret_unsafe(T, subexpr.partials_storage_ϵ),
            input_ϵ,
            subexpr_forward_values_ϵ,
            d.data.operators,
        )
    end
    _forward_eval_ϵ(
        d,
        ex,
        _reinterpret_unsafe(T, d.forward_storage_ϵ),
        _reinterpret_unsafe(T, d.partials_storage_ϵ),
        input_ϵ,
        subexpr_forward_values_ϵ,
        d.data.operators,
    )
    # do a reverse pass
    subexpr_reverse_values_ϵ =
        _reinterpret_unsafe(T, d.subexpression_reverse_values_ϵ)
    for i in ex.dependent_subexpressions
        subexpr_reverse_values_ϵ[i] = zero(T)
        d.subexpression_reverse_values[i] = 0.0
    end
    _reverse_eval_ϵ(
        output_ϵ,
        ex,
        _reinterpret_unsafe(T, d.reverse_storage_ϵ),
        _reinterpret_unsafe(T, d.partials_storage_ϵ),
        d.subexpression_reverse_values,
        subexpr_reverse_values_ϵ,
        1.0,
        zero(T),
    )
    for i in length(ex.dependent_subexpressions):-1:1
        j = ex.dependent_subexpressions[i]
        subexpr = d.subexpressions[j]
        _reverse_eval_ϵ(
            output_ϵ,
            subexpr,
            _reinterpret_unsafe(T, subexpr.reverse_storage_ϵ),
            _reinterpret_unsafe(T, subexpr.partials_storage_ϵ),
            d.subexpression_reverse_values,
            subexpr_reverse_values_ϵ,
            d.subexpression_reverse_values[j],
            subexpr_reverse_values_ϵ[j],
        )
    end
    return
end

"""
    _forward_eval_ϵ(
        d::NLPEvaluator,
        ex::Union{_FunctionStorage,_SubexpressionStorage},
        storage_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
        partials_storage_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
        x_values_ϵ,
        subexpression_values_ϵ,
        user_operators::Nonlinear.OperatorRegistry,
    ) where {N,T}

Evaluate the directional derivatives of the expression tree in `ex`.

This is equivalent to evaluating the expression tree using DualNumbers or
ForwardDiff, but instead we keep the `ex.forward_storage` for the epsilon
components separate so that we don't need to recompute the real components.

This assumes that `_reverse_model(d, x)` has already been called.
"""
function _forward_eval_ϵ(
    d::NLPEvaluator,
    ex::Union{_FunctionStorage,_SubexpressionStorage},
    storage_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
    partials_storage_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
    x_values_ϵ,
    subexpression_values_ϵ,
    user_operators::Nonlinear.OperatorRegistry,
) where {N,T}
    @assert length(storage_ϵ) >= length(ex.nodes)
    @assert length(partials_storage_ϵ) >= length(ex.nodes)
    zero_ϵ = zero(ForwardDiff.Partials{N,T})
    # ex.nodes is already in order such that parents always appear before children
    # so a backwards pass through ex.nodes is a forward pass through the tree
    children_arr = SparseArrays.rowvals(ex.adj)
    for k in length(ex.nodes):-1:1
        node = ex.nodes[k]
        partials_storage_ϵ[k] = zero_ϵ
        if node.type == Nonlinear.NODE_VARIABLE
            storage_ϵ[k] = x_values_ϵ[node.index]
        elseif node.type == Nonlinear.NODE_VALUE
            storage_ϵ[k] = zero_ϵ
        elseif node.type == Nonlinear.NODE_SUBEXPRESSION
            storage_ϵ[k] = subexpression_values_ϵ[node.index]
        elseif node.type == Nonlinear.NODE_PARAMETER
            storage_ϵ[k] = zero_ϵ
        else
            @assert node.type != Nonlinear.NODE_MOI_VARIABLE
            ϵtmp = zero_ϵ
            @inbounds children_idx = SparseArrays.nzrange(ex.adj, k)
            for c_idx in children_idx
                @inbounds ix = children_arr[c_idx]
                @inbounds partial = ex.partials_storage[ix]
                @inbounds storage_val = storage_ϵ[ix]
                # TODO: This "if" statement can take 8% of the hessian
                # evaluation time! Find a more efficient way.
                if !isfinite(partial) && storage_val == zero_ϵ
                    continue
                end
                ϵtmp += storage_val * ex.partials_storage[ix]
            end
            storage_ϵ[k] = ϵtmp
            if node.type == Nonlinear.NODE_CALL_MULTIVARIATE
                # TODO(odow): consider how to refactor this into Nonlinear.
                op = node.index
                n_children = length(children_idx)
                if op == 3 # :*
                    # Lazy approach for now.
                    anyzero = false
                    tmp_prod = one(ForwardDiff.Dual{TAG,T,N})
                    for c_idx in children_idx
                        ix = children_arr[c_idx]
                        sval = ex.forward_storage[ix]
                        gnum = ForwardDiff.Dual{TAG}(sval, storage_ϵ[ix])
                        tmp_prod *= gnum
                        anyzero = ifelse(sval * sval == zero(T), true, anyzero)
                    end
                    # By a quirk of floating-point numbers, we can have
                    # anyzero == true && ForwardDiff.value(tmp_prod) != zero(T)
                    if anyzero || n_children <= 2
                        for c_idx in children_idx
                            prod_others = one(ForwardDiff.Dual{TAG,T,N})
                            for c_idx2 in children_idx
                                (c_idx == c_idx2) && continue
                                ix = children_arr[c_idx2]
                                gnum = ForwardDiff.Dual{TAG}(
                                    ex.forward_storage[ix],
                                    storage_ϵ[ix],
                                )
                                prod_others *= gnum
                            end
                            partials_storage_ϵ[children_arr[c_idx]] =
                                ForwardDiff.partials(prod_others)
                        end
                    else
                        for c_idx in children_idx
                            ix = children_arr[c_idx]
                            prod_others =
                                tmp_prod / ForwardDiff.Dual{TAG}(
                                    ex.forward_storage[ix],
                                    storage_ϵ[ix],
                                )
                            partials_storage_ϵ[ix] =
                                ForwardDiff.partials(prod_others)
                        end
                    end
                elseif op == 4 # :^
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds base = ex.forward_storage[ix1]
                    @inbounds base_ϵ = storage_ϵ[ix1]
                    @inbounds exponent = ex.forward_storage[ix2]
                    @inbounds exponent_ϵ = storage_ϵ[ix2]
                    base_gnum = ForwardDiff.Dual{TAG}(base, base_ϵ)
                    exponent_gnum = ForwardDiff.Dual{TAG}(exponent, exponent_ϵ)
                    if exponent == 2
                        partials_storage_ϵ[ix1] = 2 * base_ϵ
                    elseif exponent == 1
                        partials_storage_ϵ[ix1] = zero_ϵ
                    else
                        partials_storage_ϵ[ix1] = ForwardDiff.partials(
                            exponent_gnum * base_gnum^(exponent_gnum - 1),
                        )
                    end
                    result_gnum = ForwardDiff.Dual{TAG}(
                        ex.forward_storage[k],
                        storage_ϵ[k],
                    )
                    # TODO(odow): fix me to use NaNMath.jl instead
                    log_base_gnum = base_gnum < 0 ? NaN : log(base_gnum)
                    partials_storage_ϵ[ix2] =
                        ForwardDiff.partials(result_gnum * log_base_gnum)
                elseif op == 5 # :/
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds numerator = ex.forward_storage[ix1]
                    @inbounds numerator_ϵ = storage_ϵ[ix1]
                    @inbounds denominator = ex.forward_storage[ix2]
                    @inbounds denominator_ϵ = storage_ϵ[ix2]
                    recip_denominator =
                        1 / ForwardDiff.Dual{TAG}(denominator, denominator_ϵ)
                    partials_storage_ϵ[ix1] =
                        ForwardDiff.partials(recip_denominator)
                    partials_storage_ϵ[ix2] = ForwardDiff.partials(
                        -ForwardDiff.Dual{TAG}(numerator, numerator_ϵ) *
                        recip_denominator *
                        recip_denominator,
                    )
                elseif op > 6
                    f_input = _UnsafeVectorView(d.jac_storage, n_children)
                    for (i, c) in enumerate(children_idx)
                        f_input[i] = ex.forward_storage[children_arr[c]]
                    end
                    H = _UnsafeLowerTriangularMatrixView(
                        d.user_output_buffer,
                        n_children,
                    )
                    has_hessian = Nonlinear.eval_multivariate_hessian(
                        user_operators,
                        user_operators.multivariate_operators[node.index],
                        H,
                        f_input,
                    )
                    # This might be `false` if we extend this code to all
                    # multivariate functions.
                    @assert has_hessian
                    for col in 1:n_children
                        dual = zero(ForwardDiff.Partials{N,T})
                        for row in 1:n_children
                            # Make sure we get the lower-triangular component.
                            h = row >= col ? H[row, col] : H[col, row]
                            # Performance optimization: hessians can be quite
                            # sparse
                            if !iszero(h)
                                i = children_arr[children_idx[row]]
                                dual += h * storage_ϵ[i]
                            end
                        end
                        i = children_arr[children_idx[col]]
                        partials_storage_ϵ[i] = dual
                    end
                end
            elseif node.type == Nonlinear.NODE_CALL_UNIVARIATE
                @inbounds child_idx = children_arr[ex.adj.colptr[k]]
                f′′ = Nonlinear.eval_univariate_hessian(
                    user_operators,
                    user_operators.univariate_operators[node.index],
                    ex.forward_storage[child_idx],
                )
                partials_storage_ϵ[child_idx] = f′′ * storage_ϵ[child_idx]
            end
        end
    end
    return storage_ϵ[1]
end

# Compute directional derivatives of the reverse pass, goes with _forward_eval_ϵ
# to compute hessian-vector products.
function _reverse_eval_ϵ(
    output_ϵ::AbstractVector{ForwardDiff.Partials{N,T}},
    ex::Union{_FunctionStorage,_SubexpressionStorage},
    reverse_storage_ϵ,
    partials_storage_ϵ,
    subexpression_output,
    subexpression_output_ϵ,
    scale::T,
    scale_ϵ::ForwardDiff.Partials{N,T},
) where {N,T}
    @assert length(reverse_storage_ϵ) >= length(ex.nodes)
    @assert length(partials_storage_ϵ) >= length(ex.nodes)
    if ex.nodes[1].type == Nonlinear.NODE_VARIABLE
        @inbounds output_ϵ[ex.nodes[1].index] += scale_ϵ
        return
    elseif ex.nodes[1].type == Nonlinear.NODE_SUBEXPRESSION
        @inbounds subexpression_output[ex.nodes[1].index] +=
            scale * ex.reverse_storage[1]
        @inbounds subexpression_output_ϵ[ex.nodes[1].index] += scale_ϵ
        return
    end
    reverse_storage_ϵ[1] = scale_ϵ
    for k in 2:length(ex.nodes)
        @inbounds node = ex.nodes[k]
        if node.type == Nonlinear.NODE_VALUE ||
           node.type == Nonlinear.NODE_LOGIC ||
           node.type == Nonlinear.NODE_COMPARISON ||
           node.type == Nonlinear.NODE_PARAMETER
            continue
        end
        parent_value = scale * ex.reverse_storage[node.parent]
        if !isfinite(ex.partials_storage[k]) && iszero(parent_value)
            reverse_storage_ϵ[k] = zero(ForwardDiff.Partials{N,T})
        else
            reverse_storage_ϵ[k] = ForwardDiff._mul_partials(
                partials_storage_ϵ[k],
                reverse_storage_ϵ[node.parent],
                parent_value,
                ex.partials_storage[k],
            )
        end
        if node.type == Nonlinear.NODE_VARIABLE
            @inbounds output_ϵ[node.index] += reverse_storage_ϵ[k]
        elseif node.type == Nonlinear.NODE_SUBEXPRESSION
            @inbounds subexpression_output[node.index] +=
                scale * ex.reverse_storage[k]
            @inbounds subexpression_output_ϵ[node.index] += reverse_storage_ϵ[k]
        end
    end
    return
end


text:
function update_opinion!(
    state::Tuple{AbstractGraph, AbstractArray}, agent_idx::Integer,
    config::Config
)
    agent_list = state[2]
    this_agent = agent_list[agent_idx]
    # weighted mean of own opinion and perceived public opinion
    if (abs(this_agent.opinion - this_agent.perceiv_publ_opinion) < config.opinion_treshs.backfire)
        this_agent.opinion = (
            config.agent_props.own_opinion_weight * this_agent.opinion
            + (1 - config.agent_props.own_opinion_weight) * this_agent.perceiv_publ_opinion
        )
    else
        if ((this_agent.opinion * this_agent.perceiv_publ_opinion > 0)
            && (abs(this_agent.opinion) - abs(this_agent.perceiv_publ_opinion) < 0))
            this_agent.opinion = config.agent_props.own_opinion_weight * this_agent.opinion
        else
            this_agent.opinion = (2 - config.agent_props.own_opinion_weight) * this_agent.opinion
        end
        if this_agent.opinion > 1
            this_agent.opinion = 1
        elseif this_agent.opinion < -1
            this_agent.opinion = -1
        end
    end
    return state
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using Test

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

for typ in Base.uniontypes(MPI.MPIDatatype)

    # Allocating version
    a = fill(convert(typ, rank),size)
    b = MPI.Alltoall(a, 1, comm)
    @test b == convert.(typ, collect(0:(size-1)))

    # Non Allocating version
    a = fill(convert(typ, rank),size)
    b = Array{typ}(undef, size*1)
    MPI.Alltoall!(a, b, 1, comm)
    @test b == convert.(typ, collect(0:(size-1)))

    # IN_PLACE version
    a = fill(convert(typ, rank),size)
    MPI.Alltoall!(MPI.IN_PLACE, a, 1, comm)
    @test a == convert.(typ, collect(0:(size-1)))

end

MPI.Finalize()


text:
module Kepler

function ecc_anom(M::Float64, e::Float64)
    M = mod2pi(M)
    if M < zero(M)
        M = M + 2.0*pi
    end

    Emin = zero(M)
    fmin = zero(M)
    Emax = 2.0*pi*one(M)
    fmax = zero(M)
    Eguess = M
    while true
        f = M - Eguess + e*sin(Eguess)

        if abs(f) < convert(typeof(f), 1e-10)
            return Eguess
        end
        df = -one(Eguess) + e*cos(Eguess)

        Enewton = Eguess - f/df

        if Enewton < Emin || Enewton > Emax
            fmin = M - Emin - e*sin(Emin)
            fmax = M - Emax - e*sin(Emax)

            if fmin == zero(fmin)
                return Emin
            elseif fmax == zero(fmin)
                return Emax
            end

            if fmin*f < zero(fmin)
                Emax = Eguess
            elseif fmax*f < zero(fmax)
                Emin = Eguess
            else
                error("Not bracketing a root in bisection.")
            end
            Eguess = (Emax-Emin)/2
        else
            Eguess = Enewton
        end
    end
end

function true_anom(ea::Float64, e::Float64)
    ta = 2.0*atan(sqrt((one(e)+e)/(one(e)-e))*tan(ea/2.0))
    if ta < zero(ta)
        ta + 2.0*pi
    else
        ta
    end
end

function rv(t::Float64, K::Float64, P::Float64, e::Float64, omega::Float64, chi::Float64)
    fp, ip = modf(t/P - chi)

    if fp < zero(fp)
        fp = fp + one(fp)
    end

    m = 2*pi*fp

    ta = true_anom(ecc_anom(m, e), e)

    K*(cos(ta + omega) + e*cos(omega))
end

end


text:
function trunc(x::Float64)
    xu = reinterpret(Uint64,x)
    k = int(xu >> 52) & 0x07ff - 1023
    if k > 1023 # Inf or NaN
        x+x
    elseif k < 0
        copysign(0.0,x)
    else
        reinterpret(Float64, xu & ~(0x000f_ffff_ffff_ffff >> k))
    end
end

function trunc(x::Float32)
    xu = reinterpret(Uint32,x)
    k = int(xu >> 23) & 0x00ff - 127
    if k > 127 # Inf or NaN
        x+x
    elseif k < 0
        copysign(0f0,x)
    else
        reinterpret(Float32, xu & ~(0x007f_ffff >> k))
    end
end


text:
using MAT

@testset "simulateSNR" begin
    # MATLAB output
    test = matread("resources/simulateSNR.mat")

    # test vectors
    osr  = [32,    64,  8,  8]
    f0   = [ 0, 0.125,  0,  0]
    nlev = [ 2,     2, 17, 17]
    ntf = [
        synthesizeNTF(5, osr[1], 1, 1.5, f0[1]),
        synthesizeNTF(8, osr[2], 1, 1.5, f0[2]),
        synthesizeNTF(7, osr[3], 1,   2, f0[3]),
        synthesizeNTF(7, osr[4], 1,   8, f0[4])
    ]
    amp_test = [-120:10:-20; -15; -10:0]

    for i = 1:length(ntf)
        snr, amp = simulateSNR(ntf[i], osr[i], amp_test, f0[i], nlev[i])
        @test isapprox(vec(test["snr_$i"]), snr)
        @test amp_test == amp   # passed through without changes
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
###############################################################################
#
# Calculate gas exchange from flow rate only for required parameters
#
###############################################################################
"""
    leaf_gas_exchange_nonopt!(
                node::SPACSimple{FT},
                flow::FT,
                par::FT,
                rad::FT,
                la::FT,
                container::SPACContainer1L{FT}
    ) where {FT<:AbstractFloat}
    leaf_gas_exchange_nonopt!(
                node::SPACSimple{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    leaf_gas_exchange_nonopt!(
                node::SPACSimple{FT},
                f_sl::FT,
                f_sh::FT
    ) where {FT<:AbstractFloat}

Simulate leaf level gas exchange and fill it into the `container` for 1-layer
    or 2-layer canopy, given
- `node` [`SPACSimple`] type struct
- `flow` Flow rate per basal area into the leaves (e.g., for sunlit leaves)
- `f_sl` Flow rate per basal area into the sunlit leaves
- `f_sh` Flow rate per basal area into the shaded leaves
- `par` Leaf-level photosynthetic active radiation
- `rad` Leaf-level absorbed radiative energy
- `la` Leaf area of the leaves (total or each layer)
- `container` [`SPACContainer1L`] type container
"""
function leaf_gas_exchange_nonopt!(
            node::SPACSimple{FT},
            flow::FT,
            par::FT,
            rad::FT,
            la::FT,
            container::SPACContainer1L{FT}
) where {FT<:AbstractFloat}
    # if flow < 0
    if flow < 0
        container.an = FT(-Inf);

    # if flow >= 0
    else
        # 0. unpack required variables
        @unpack envir, g_max, width = node;
        @unpack P_AIR, p_H₂O, wind = envir;

        # 1. calculate leaf temperature from the flow rate
        t_leaf = max(200, leaf_temperature(node, rad, flow));

        # 2. update leaf photosynthetic variables and leaf-to-air VPD
        node.ps.apar = par;
        photosystem_temperature_dependence!(node.ps.PSM, envir, t_leaf);
        node.ps.p_H₂O_sat = saturation_vapor_pressure(t_leaf);
        d_leaf = node.ps.p_H₂O_sat - p_H₂O;

        # 3. update f_vis and f_st in leaf and calculate water potentials
        # TODO do not account for temperature effects for flow now
        # More reasonable functions need to be added
        # vc_temperature_effects!(node.hs.leaf, t_leaf);

        # if flow == 0
        if flow == 0
            g_lw = FT(0);
            g_lc = FT(1e-6);
            leaf_photosynthesis!(node.ps, envir, GCO₂Mode(), g_lc);
            container.an = node.ps.PSM.a_net;

        # if flow > 0 and reasonable
        elseif (d_leaf > 0) && (t_leaf > T_0(FT))
            t_cor = relative_diffusive_coefficient( (t_leaf+node.envir.t)/2 );
            g_lw  = flow / la / d_leaf * P_AIR;
            g_bw  = boundary_layer_conductance(wind, width);
            g_bc  = g_bw / FT(1.35);
            g_sw  = 1 / (1/g_lw - 1/(g_bw*t_cor));
            g_sc  = g_sw / FT(1.6);
            g_lc  = max(FT(1e-6), 1 / (1/g_bc + 1/g_sc));
            g_lim = 1 / (1/g_bw + 1/g_max);
            if g_lw < g_lim * t_cor
                leaf_photosynthesis!(node.ps, envir, GCO₂Mode(), g_lc);
                container.an = node.ps.PSM.a_net;
            else
                container.an = FT(-Inf);
            end

        # if flow > 0 and unrealistic
        else
            container.an = FT(-Inf);
        end
    end

    return nothing
end

function leaf_gas_exchange_nonopt!(
            node::SPACSimple{FT},
            flow::FT
) where {FT<:AbstractFloat}
    # unpack the data
    @unpack frac_sh, frac_sl, par_sh, par_sl, rad_sh, rad_sl = node.container2L;

    # calculate mean par and rad per leaf area, then gas exchange rate
    par_mean = par_sl * frac_sl + par_sh * frac_sh;
    rad_mean = rad_sl * frac_sl + rad_sh * frac_sh;
    leaf_gas_exchange_nonopt!(node, flow, par_mean, rad_mean, node.laba, node.container1L);

    node.containerOP = (node.ec - flow) * (node.container1L).an;

    return nothing
end

function leaf_gas_exchange_nonopt!(
            node::SPACSimple{FT},
            f_sl::FT,
            f_sh::FT
) where {FT<:AbstractFloat}
    # unpack the data
    @unpack frac_sh, frac_sl, la_sh, la_sl, par_sh, par_sl, rad_sh,
            rad_sl = node.container2L;

    # calculate gas exchangr for sunlit and shaded layers
    leaf_gas_exchange_nonopt!(node, f_sl, par_sl, rad_sl, la_sl, node.container2L.cont_sl);
    leaf_gas_exchange_nonopt!(node, f_sh, par_sh, rad_sh, la_sh, node.container2L.cont_sh);

    a_sum = frac_sl * node.container2L.cont_sl.an + frac_sh * node.container2L.cont_sh.an;
    e_sum = f_sl + f_sh;
    node.containerOP = (node.ec - e_sum) * a_sum;

    return nothing
end

###############################################################################
#
# Calculate gas exchange from flow rate considering temperature dynamics
#
###############################################################################
"""
    leaf_gas_exchange!(
                node::SPACSimple{FT},
                flow::FT,
                par::FT,
                rad::FT,
                la::FT,
                container::SPACContainer1L{FT}
    ) where {FT<:AbstractFloat}
    leaf_gas_exchange!(
                node::SPACSimple{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    leaf_gas_exchange!(
                node::SPACSimple{FT},
                f_sl::FT,
                f_sh::FT
    ) where {FT<:AbstractFloat}

Simulate leaf level gas exchange and fill it into the `container` for 1-layer
    or 2-layer canopy, given
- `node` [`SPACSimple`] type struct
- `flow` Flow rate per basal area into the leaves (e.g., for sunlit leaves)
- `f_sl` Flow rate per basal area into the sunlit leaves
- `f_sh` Flow rate per basal area into the shaded leaves
- `par` Leaf-level photosynthetic active radiation
- `rad` Leaf-level absorbed radiative energy
- `la` Leaf area of the leaves (total or each layer)
- `container` [`SPACContainer1L`] type container
"""
function leaf_gas_exchange!(
            node::SPACSimple{FT},
            flow::FT,
            par::FT,
            rad::FT,
            la::FT,
            container::SPACContainer1L{FT}
) where {FT<:AbstractFloat}
    # 0. unpack required variables
    @unpack envir = node;
    @unpack P_AIR, p_H₂O = envir;

    # 1. calculate leaf temperature from the flow rate
    t_leaf = max(200, leaf_temperature(node, rad, flow));

    # 2. update leaf photosynthetic variables and leaf-to-air VPD
    node.ps.apar = par;
    photosystem_temperature_dependence!(node.ps.PSM, envir, t_leaf);
    node.ps.p_H₂O_sat = saturation_vapor_pressure(t_leaf);
    d_leaf = node.ps.p_H₂O_sat - p_H₂O;

    # 3. update f_vis and f_st in leaf and calculate water potentials
    # TODO do not account for temperature effects for flow now
    # More reasonable functions need to be added
    # vc_temperature_effects!(node.hs.leaf, t_leaf);
    p_leaf = end_pressure(node.hs, flow);

    # 4. calculate photosynthesis
    g_lw = flow / la / d_leaf * P_AIR;
    g_lc = max(FT(1e-6), g_lw / FT(1.6));
    leaf_photosynthesis!(node.ps, envir, GCO₂Mode(), g_lc);
    container.ag = node.ps.PSM.a_gross;
    container.an = node.ps.PSM.a_net;
    container.c  = node.ps.p_CO₂_i;
    container.e  = flow;
    container.gh = g_lw;
    container.p  = p_leaf;
    container.t  = t_leaf;

    return nothing
end

function leaf_gas_exchange!(
            node::SPACSimple{FT},
            flow::FT
) where {FT<:AbstractFloat}
    # unpack the data
    @unpack container2L, laba = node;
    @unpack frac_sh, frac_sl, par_sh, par_sl, rad_sh, rad_sl = container2L;

    # calculate mean par and rad per leaf area, then gas exchange rate
    par_mean = par_sl * frac_sl + par_sh * frac_sh;
    rad_mean = rad_sl * frac_sl + rad_sh * frac_sh;
    leaf_gas_exchange!(node, flow, par_mean, rad_mean, laba, node.container1L);

    return nothing
end

function leaf_gas_exchange!(
            node::SPACSimple{FT},
            f_sl::FT,
            f_sh::FT
) where {FT<:AbstractFloat}
    # unpack the data
    @unpack frac_sh, frac_sl, la_sh, la_sl, par_sh, par_sl, rad_sh,
            rad_sl = node.container2L;

    # calculate gas exchangr for sunlit and shaded layers
    leaf_gas_exchange!(node, f_sl, par_sl, rad_sl, la_sl, node.container2L.cont_sl);
    leaf_gas_exchange!(node, f_sh, par_sh, rad_sh, la_sh, node.container2L.cont_sh);

    return nothing
end


text:
using Test
using OpenQuantumSystems
using Random, SparseArrays, LinearAlgebra, StableRNGs
import QuantumOpticsBase

@testset "memory kernel" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    # TODO: change to macro
    HR = 0.01
    shift = (2.0 * HR)^0.5
    mode1 = Mode(300.0, shift)
    Energy = [12500.0, 12800.0]
    mol1 = Molecule([mode1], 2, Energy)
    mol2 = Molecule([mode1], 2, Energy)
    agg = Aggregate([mol1, mol2])
    agg.coupling[1, 2] = 100
    agg.coupling[2, 1] = 100
    aggInds = getIndices(agg)
    vibindices = getVibIndices(agg, aggInds)
    aggIndsLen = length(aggInds)
    basis = GenericBasis([aggIndsLen])
    FCFact = getFranckCondonFactors(agg, aggInds)
    FCProd = getFCProd(agg, FCFact, aggInds, vibindices)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    Ham_I = getAggHamInteraction(agg, aggInds, FCFact)
    Ham_0 = getAggHamSystemBath(agg, aggInds, FCFact)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    # T = 300
    # W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], Ham, aggInds; diagonalize=false, diagonal=true)
    data = Matrix(Hermitian(rand(ComplexF64, aggIndsLen, aggIndsLen)))
    W0 = DenseOperator(basis, basis, data)
    normalize!(W0)
    W0_bath = get_rho_bath(
        W0,
        agg,
        FCProd,
        aggInds,
        vibindices;
        justCopy = true,
    )
    rho0 = trace_bath(W0.data, agg, FCProd, aggInds, vibindices)

    W_ab = take_el_part(W0.data, 1, 1, vibindices)
    W_ab_len = size(W_ab, 1)
    @test 1e-15 > D(W_ab, W0.data[1:W_ab_len, 1:W_ab_len])

    tspan = [0.0:0.1:0.3;]
    for t_i = 1:length(tspan)
        for s_i = 1:t_i
            t = tspan[t_i]
            s = tspan[s_i]
            # println(t, " ", s)
            H_II_t = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, t)
            H_II_s = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, s)

            MemoryKernel_1 = MemoryKernel_1_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_1[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_1_ref = H_II_t * H_II_s * W0.data
            rho_ref = trace_bath(
                MemoryKernel_1_ref,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_2 = MemoryKernel_2_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_2[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_2_ref = H_II_t * W0.data * H_II_s
            rho_ref = trace_bath(
                MemoryKernel_2_ref,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_3 = MemoryKernel_3_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_3[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_3_ref = H_II_s * W0.data * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_3_ref,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))

            MemoryKernel_4 = MemoryKernel_4_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel_4[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_4_ref = W0.data * H_II_s * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_4_ref,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            @test 1e-6 > D(rho, rho_ref)
            # println(D(rho, rho_ref))
        end
    end

    tspan = [0.0:0.1:0.3;]
    for t_i = 1:length(tspan)
        for s_i = 1:t_i
            t = tspan[t_i]
            s = tspan[s_i]
            H_II_t = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, t)
            H_II_s = getInteractionHamIPictureA(Ham_I.data, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, s)

            MemoryKernel = MemoryKernel_traced(
                H_II_t,
                H_II_s,
                W0_bath.data,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            rho = zero(complex(rho0))
            for a = 1:size(rho, 1),
                b = 1:size(rho, 2),
                c = 1:size(rho, 1),
                d = 1:size(rho, 2)

                rho[a, b] += MemoryKernel[a, b, c, d] * rho0[c, d]
            end
            MemoryKernel_ref = H_II_s * W0.data - W0.data * H_II_s
            MemoryKernel_ref = H_II_t * MemoryKernel_ref - MemoryKernel_ref * H_II_t
            rho_ref = trace_bath(
                MemoryKernel_ref,
                agg,
                FCProd,
                aggInds,
                vibindices
            )
            @test 1e-6 > D(rho, rho_ref)
        end
    end

end # testset


text:
## history.jl : a history object for bookkeeping
#
# A history object contains bookkeeping data for multilevel estimators.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

## History ##
"""
    History

A History is the result of a call to [`run`](@ref) and contains useful diagnostic information about the state of the `Estimator`.

The information is stored as a `Vector` of `Dict`'s, where each entry in the vector is the result of a single run. See below for a list of keys added to these `Dict`s. Access to the key corresponding to the last entry in the `History` (that is, the one with the highest accuracy) is provided as `history[key]`.

See also the package [`Reporter`](https://github.com/PieterjanRobbe/Reporter.jl) for automatic generation of reports based on a History object.

# Keys

- `type::AbstractIndexSet`: the index set type of the `Estimator`.
- `ndims::Integer`: the number of dimensions of the `Estimator`.
- `name::String`: the name of this `Estimator`.
- `folder::String`: the directory where the results are stored.
- `elapsed::Real`: the elapsed time (in seconds).
- `tol::Real`: the request accuracy.
- `current_index_set::Vector{::Index}`: the levels (or indices) currently in use by the `Estimator`.
- `index_set::Vector{::Index}`: all levels (or indices) that are added to the `Estimator`.
- `mse::Real`: the measured mean square error of the `Estimator`.
- `rmse::Real`: the measured root mean square error of the `Estimator`.
- `mean::Real`: the estimated expected value of the quantity of interest, up to an absolute error of `tol`.
- `var::Real`: the measured variance of the quantity of interest.
- `varest::Real`: the measured variance of the `Estimator` for the expected value of the quantity of interest.
- `bias::Real`: the measured bias of the `Estimator`.
- `E::Dict{::Index, ::Real}`: a `Dict` with the measured expected value of the quantity of interest on each level (or index).
- `V::Dict{::Index, ::Real}`: a `Dict` with the measured variance of the quantity of interest on each level (or index).
- `dE::Dict{::Index, ::Real}`: a `Dict` with the measured expected value of the **difference** of the quantity of interest on each level (or index).
- `dV::Dict{::Index, ::Real}`: a `Dict` with the measured variance of the **difference** of the quantity of interest on each level (or index).
- `T::Dict{::Index, ::Real}`: a `Dict` with the measured computation time for a sample of the quantity of interest on each level (or index).
- `W::Dict{::Index, ::Real}`: a `Dict` with the computational cost for a sample of the quantity of interest on each level (or index). *Only available when a cost model was provided!*
- `α::Tuple`: the estimated rates of decay of the expected value of the difference of the quantity of interest.
- `β::Tuple`: the estimated rates of decay of the variance of the difference of the quantity of interest.
- `γ::Tuple`: the estimated rates of increase of the cost of the difference of the quantity of interest.
- `nb_of_samples::Array`: an `Array` with, on each index, the number of samples taken on that index.
- `logbook`: a logbook with information about the adaptive algorithm *Only available for `Estimator`s with index sets of type `AD`*
- `samples::Array`: an `Array` with on each index, the collection of samples of the quantity of interest taken on that index. *Only available when the optional key `save_samples = true` was used*
- `samples_diff::Array`: an `Array` with on each index, the collection of samples of the **difference** of the quantity of interest taken on that index. *Only available when the optional key `save_samples = true` was used*
See also: [`Estimator`](@ref), [`run`](@ref)
"""
mutable struct History{T}
    data::T
end

History() = History(Vector{Dict{Symbol, Any}}(undef, 0))

function update_history!(history::History, estimator::Estimator, tol::Real, elapsed::Real)
    h = Dict{Symbol, Any}()
    h[:type]              = typeof(estimator)
    h[:ndims]             = ndims(estimator)
    h[:name]              = estimator[:name]
    h[:folder]            = estimator[:folder]
    h[:elapsed]           = elapsed
    h[:tol]               = tol
    h[:current_index_set] = copy(collect(keys(estimator)))
    h[:index_set]         = copy(collect(all_keys(estimator)))
    h[:mse]               = mse(estimator)
    h[:rmse]              = rmse(estimator)
    h[:mean]              = mean(estimator)
    h[:var]               = var(estimator)
    h[:varest]            = varest(estimator)
    h[:bias]              = bias(estimator)
    h[:E]                 = apply(mean0, estimator)
    h[:V]                 = apply(var0, estimator)
    h[:dE]                = apply(mean, estimator)
    h[:dV]                = apply(var, estimator)
    h[:T]                 = apply(time, estimator)
    h[:α]                 = α(estimator)
    h[:β]                 = β(estimator)
    h[:γ]                 = γ(estimator)
    h[:nb_of_samples]     = copy(nb_of_samples(estimator))
    h[:qoi_with_max_var]  = qoi_with_max_var(estimator)

    # add cost model is it was provided
    if !(estimator[:cost_model] isa EmptyFunction)
        h[:W]             = apply(cost, estimator)
    end

    # add logbook when using adpative index sets
    if estimator isa Estimator{<:AD}
        h[:logbook]       = copy(logbook(estimator))
    end

    # save samples if required
    if estimator[:save]
        if estimator[:save_samples]
            h[:samples] = samples(estimator)
            h[:samples_diff] = samples_diff(estimator)
        end

        push!(history.data, h)

        # save history
        save(history)
    end
end

apply(f::Function, estimator::Estimator) = Dict(i => f(estimator, i) for i in all_keys(estimator))

save(history::History) = @save joinpath(history[:folder], history[:name]) history

getindex(history::History, s::Symbol) = history.data[end][s]

getindex(history::History, i::Integer) = history.data[i]

haskey(history::History, s::Symbol) = haskey(history.data[end], s)

length(history::History) = length(history.data)

show(io::IO, history::History) = print(io, "MultilevelEstimators.jl history file")


text:
# ---
# title: 518. Coin Change 2
# id: problem518
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/coin-change-2/description/>
# hidden: true
# ---
#
# You are given coins of different denominations and a total amount of money.
# Write a function to compute the number of combinations that make up that
# amount. You may assume that you have infinite number of each kind of coin.
#
#
#
# **Example 1:**
#
#
#
#     Input: amount = 5, coins = [1, 2, 5]
#     Output: 4
#     Explanation: there are four ways to make up the amount:
#     5=5
#     5=2+2+1
#     5=2+1+1+1
#     5=1+1+1+1+1
#
#
# **Example 2:**
#
#
#
#     Input: amount = 3, coins = [2]
#     Output: 0
#     Explanation: the amount of 3 cannot be made up just with coins of 2.
#
#
# **Example 3:**
#
#
#
#     Input: amount = 10, coins = [10]
#     Output: 1
#
#
#
#
# **Note:**
#
# You can assume that
#
#   * 0 <= amount <= 5000
#   * 1 <= coin <= 5000
#   * the number of coins is less than 500
#   * the answer is guaranteed to fit into signed 32-bit integer
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Compat.String

#
# This form maps JuliaBox ids to MIT email ids
#
function show_mit_form()
    display(MIME("text/html"), """Loading...""")
end
#
# Save answer in a closure for the submit button event,
# make the button and display it. return the answer as is.
#
function attempt_prompt(metadata_json, answer)

    metadata = JSON.parse(metadata_json)

    metadata_channel = Signal(Any, Dict())

    question = metadata["question"]
    map(metadata_channel, init=script("")) do x
        set_metadata(question, x)
    end |> display

    if !get(metadata, "finished", false)
        b = button("Submit »")
        display(b)

        map(signal(b), init=nothing) do _
            evaluate(metadata, answer, metadata_channel)
        end
    end
    answer
end

#
# Evaluate an answer.
#

function evaluate(metadata, answer, meta)
    @assert haskey(metadata, "question")

    if !haskey(global_config, "host")
        global_config["host"] = "https://juliabox.com"
    end
    question_no = string(metadata["question"])

    @async push!(meta, alert("info", "Evaluating your answer..."))

    # The HTTP requests to evaluate answer goes here...
    # After the request, you can push the
    query_params = Dict(
        "mode" => "submit",
        "params" => JSON.json(Dict(
            "course" => global_config["course"],
            "problemset" => global_config["problemset"],
            "question" => question_no,
            "answer" => JSON.json(encode(metadata, answer))
        ))
    )

    response = nothing
    try
        response = get(string(strip(global_config["host"], ['/']), "/jboxplugin/hw/");
            query = query_params,
            headers = Dict("Cookie" => global_config["cookie"]))

    catch err
        open("homework_log.log", "a") do f
            bt = catch_backtrace()
            Base.show_backtrace(f, bt)
            println(f, err)
        end
    end

    try
        if response != nothing && statuscode(response) == 200
            result = Requests.json(response)
            if result["code"] != 0
                open("homework_log.log", "a") do f
                    println(f, "Code Non-zero")
                end
                @async push!(meta, alert("danger", "Something went wrong while verifying your code!"))
            else
                @async report_evaluation(result, metadata, meta)
            end
        else
            open("homework_log.log", "a") do f
                write(f, string(response))
            end
            @async push!(meta, alert("danger",
                "There was an unexpected error while accessing the homework server."))
        end

    catch err
        open("homework_log.log", "a") do f
            bt = catch_backtrace()
            Base.show_backtrace(f, bt)
            println(f, err)
        end
        @async push!(meta, alert("danger",
            "Something went wrong in Homework.jl"))
    end
    # return the answer itself for consistency
    answer
end

function report_evaluation(result, metadata, meta_channel)

    data = result["data"]
    status = data["status"]
    score = data["score"]
    attempts = data["attempts"]
    explanation = get(data, "explanation", nothing)
    data["max_attempts"] = data["max_attempts"] == 0 ? get(metadata, "max_attempts", 0) : data["max_attempts"]
    max_attempts = data["max_attempts"]
    max_score = data["max_score"]

    if status == 1
        msg = " Your last attempt was correct!"

        if isa(explanation, String) && explanation != ""
            msg = msg*"

Explanation: $explanation
"
        end

        merge!(data,
            alert("success", msg))

        data["finished"] = true
    else
        if max_attempts != 0 && attempts > max_attempts
            msg = " You exceeded the maximum number of failed attempts allowed for this question."
            data["finished"] = true
        else
            msg = " Wrong answer. Try again."
        end

        merge!(data,
            alert("warning", msg,))
    end

    push!(meta_channel, data)
end


text:
using Distributions
import Random: rand!, rand
import Distributions: logpdf
import Base: eltype, length

"""
    RandomWalk(ϵ::T, pos::S)

Defines a random walk on `|ϵ|`-dimensional space. ``ϵ` defines the maximal
one-sided range of a single step and `pos` is a vector of indicators for whether
a respective index is restricted to take only positive values. For elements
restricted to take positive values, the update is done via: x⁽ⁿᵉʷ⁾ <- x⁽ᵒˡᵈ⁾eᵁ,
where U ∼ Unif(-ϵ,ϵ). For unrestricted: x⁽ⁿᵉʷ⁾ <- x⁽ᵒˡᵈ⁾ + U,
where U ∼ Unif(-ϵ,ϵ).
"""
struct RandomWalk{T,N}
    ϵ::NTuple{N,T}
    pos::NTuple{N,Bool}

    function RandomWalk(ϵ::S, pos=nothing
                        ) where S<:Union{T,Vector{T},NTuple{N,T}} where {N,T}
        _RandomWalk(T, ϵ, pos)
    end

    function _RandomWalk(::Type{T}, ϵ, pos::Nothing) where T
        ϵ, N = Tuple(ϵ), length(ϵ)
        new{T,N}(ϵ, Tuple([false for i in 1:N]))
    end

    function _RandomWalk(::Type{T}, ϵ, pos::Array{Any,1}) where T
        ϵ, N = Tuple(ϵ), length(ϵ)
        @assert length(pos) == 0
        new{T,N}(ϵ, Tuple([false for i in 1:N]))
    end

    function _RandomWalk(::Type{T}, ϵ, pos::R) where {T,R<:Union{S,Vector{S},NTuple{N1,S}}} where {N1,S<:Integer}
        ϵ, N = Tuple(ϵ), length(ϵ)
        @assert minimum(pos) >= 1 && maximum(pos) <= N
        new{T,N}(ϵ, Tuple([i in pos for i in 1:N]))
    end

    function _RandomWalk(::Type{T}, ϵ, pos::R) where {T,R<:Union{S,Vector{S},NTuple{N1,S}}} where {N1,S<:Bool}
        ϵ, N = Tuple(ϵ), length(ϵ)
        pos, N_temp = Tuple(pos), length(pos)
        @assert N == N_temp
        new{eltype(ϵ),N}(ϵ, pos)
    end

end

additiveStep = (x,pos)->pos ? 0.0 : x
multipStep = (x,pos)->pos ? exp(x) : 1.0
eltype(::RandomWalk{T}) where T = T
length(::RandomWalk{T,N}) where {T,N} = N

"""
    rand!(rw::RandomWalk, θ)

Update all elements of a random walker in-place.
"""
function rand!(rw::RandomWalk, θ)
    θ .+= additiveStep.(rand.(map(Uniform,-rw.ϵ, rw.ϵ)), rw.pos)
    θ .*= multipStep.(rand.(map(Uniform,-rw.ϵ, rw.ϵ)), rw.pos)
    θ
end

"""
    rand(rw::RandomWalk, θ)

Return a newly sampled state of a random walker, with all element updated.
"""
function rand(rw::RandomWalk, θ)
    θc = copy(θ)
    rand!(rw, θc)
end

"""
    rand!(rw::RandomWalk, θ, ::UpdtIdx)

Update elements of a random walker on indices specified by the object `UpdtIdx`
in-place.
"""
function rand!(rw::RandomWalk, θ, ::UpdtIdx) where UpdtIdx
    for i in idx(UpdtIdx())
        θ[i] += additiveStep(rand(Uniform(-rw.ϵ[i], rw.ϵ[i])), rw.pos[i])
        θ[i] *= multipStep(rand(Uniform(-rw.ϵ[i], rw.ϵ[i])), rw.pos[i])
    end
    θ
end

"""
    rand(rw::RandomWalk, θ, ::UpdtIdx)

Return a newly sampled state of a random walker, with updated elements only on
the indices specified by the object `UpdtIdx`.
"""
function rand(rw::RandomWalk, θ, ::UpdtIdx) where UpdtIdx
    θc = copy(θ)
    rand!(rw, θc, UpdtIdx())
end

"""
    logpdf(rw::RandomWalk, θ, θᵒ)
Log-transition density of a random walker `rw` for going from `θ` to `θᵒ`.
"""
logpdf(rw::RandomWalk, θ, θᵒ) = sum( map((x,ϵ,pos)->pos ? -log(2.0*ϵ)-log(x) :
                                                        0.0, θᵒ, rw.ϵ, rw.pos) )


text:
# by Lazaro Alonso
using CairoMakie, LaTeXStrings, SpecialFunctions
CairoMakie.activate!() # HIDE
let
    cycle = Cycle([:color, :linestyle, :marker], covary = true)
    set_theme!(Lines = (cycle = cycle,), Scatter = (cycle = cycle,))
    #set_theme!()
    fig = Figure(resolution = (600,400), font = "CMU Serif") # probably you need to install this font in your system
    ax = Axis(fig, xlabel = L"x", ylabel = L"Y_{\nu}(x)", ylabelsize = 22,
        xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
        xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10)
    x = 0.1:0.1:15
    for ν in 0:4
        lines!(ax, x, bessely.(ν, x), label = latexstring("Y_{$(ν)}(x)"), linewidth = 2)
        #scatter!(ax, x, bessely.(ν, x), label = latexstring("Y_{$(ν)}(x)"), linewidth = 2)
    end
    axislegend(; position = :rb, nbanks = 2, framecolor = (:grey,0.5))
    ylims!(-1.8,0.7)
    #xlims!(0.1,nothing)
    fig[1,1] = ax
    #save("filename.png", fig, px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "besselSecondKind.png"), fig, px_per_unit = 2.0) # HIDE
    display(fig)
end
using Pkg # HIDE
Pkg.status() # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text: isarmstrong(n) = sum(k^length(digits(n)) for k in digits(n)) == n


text:
include(Pkg.dir("ReinforcementLearning", "examples", "maze","maze.jl"))

function main()
    maze_env = MazeEnv((5,5))
    discount=0.99

    print_maze(maze_env.maze)
    println("Start: $(maze_env.start)")
    println("Goal: $(maze_env.goal)")

    policy = Dict{MazeState, MazeAction}()

    # Initialize the policy with random actions
    for state in getAllStates(maze_env)
        policy[state] = rand(getActions(state, maze_env))
    end

    println("Initial Policy: ")
    print_policy(policy)

    V = policy_evaluation(maze_env, policy, discount)
    println("Initial State Values: ")
    print_state_values(V)

    opt_policy = policy_iteration(maze_env, policy, discount)
    println("Final Policy: ")
    print_policy(opt_policy)

    V = policy_evaluation(maze_env, opt_policy, discount)
    println("Final State Values: ")
    print_state_values(V)
end

function print_policy(policy::Dict{MazeState, MazeAction})
    # Print the policy by iterating over keys
    # MazeState -> MazeAction
    # e.g. MazeState(1,1) -> MazeAction(RIGHT)

    for (k,v) in policy
        println("$k -> $v")
    end
end

function print_state_values(V::Dict{MazeState, Float64})
    # Print the state values by iterating over keys
    # MazeState -> Float64
    # e.g. MazeState(1,1) -> 2.0

    for (k,v) in V
        println("$k -> $v")
    end
end

function policy_evaluation(env::MazeEnv, policy::Dict{MazeState, MazeAction}, discount=0.99)
    V = Dict{MazeState, Float64}() # State values
    # Initialize all state values with 0

    for state in getAllStates(env)
        V[state] = 0.0
    end

    eps = 1e-10

    # Evaluate the policy until there is no change between old values and new values
    # The maximum difference between values for a state must be less than eps
    while true
        delta = 0
        for state in getAllStates(env)
            temp = V[state]
            action = policy[state]
            successors = getSuccessors(state,action,env)
            V[state] = sum(map(i->i[3] * (i[2] + discount * V[i[1]]), successors))
            delta = max(delta, abs(temp-V[state]))
        end

        if delta < eps
            break
        end
    end

    return V
end

function policy_iteration(env::MazeEnv, policy::Dict{MazeState, MazeAction}, discount=0.99)
    # evaluate the policy
    # improve the policy
    # Implement the policy iteration algorithm

    while true
        V = policy_evaluation(env, policy, discount)
        policy_stable = true
        for state in getAllStates(env)
            temp = policy[state]
            x0 = map(a->(a,getSuccessors(state,a,env)), getActions(state,env))
            x1 = map(i->(i[1], map(j->j[3] * (j[2] + discount*V[j[1]]), i[2])), x0)
            x2 = map(i->(i[1],sum(i[2])), x1)
            x3 = sort(x2, by=i->i[2], rev=true)
            policy[state] = x3[1][1]
            if temp != policy[state]
                policy_stable = false
            end
        end

        if policy_stable
            break
        end
    end

    return policy
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/rl

max_issues_repo_name: ilkerkesen/rl

max_forks_repo_name: ilkerkesen/rl


text:
# authors: Georgios Bakirtzis (bakirtzis.net)
#          Raul Gonzalez Garcia (raulg@iastate.edu)
#
# the following example is a mechanization from
# 1. compositional cyber-physical systems modeling (http://dx.doi.org/10.4204/EPTCS.333.9)
# 2. categorical semantics of cyber-physical systems theory (https://doi.org/10.1145/3461669)
# please make sure to consult those to understand what is going on below

using AlgebraicDynamics
using AlgebraicDynamics.DWDDynam

using Catlab.WiringDiagrams
using Catlab.Graphics
using Catlab.Graphics.Graphviz

using DifferentialEquations

using Plots

using LabelledArrays

# we use the notion of functorial semantics, in the wiring diagram formalism
# that means that we define empty boxes first to define the architecture
# and then assign a behavior that will inhabit those boxes compositionally

# we first have to define our boxes, with what the inputs and outputs are,
# for example, the sensor box has two input :e and :s and one output s_prime.

s = Box(:sensor, [:s, :e], [:s′])
c = Box(:controller, [:d, :s′], [:c])
d = Box(:dynamics, [:c], [:s])

# a wiring diagram has (inputs, outputs) which define the ports of entire diagram
UAV = WiringDiagram([:e,:d], [:s])

sensor     = add_box!(UAV, s)
controller = add_box!(UAV, c)
dynamics   = add_box!(UAV, d)

add_wires!(UAV, [
    # net inputs
    (input_id(UAV), 1) => (sensor, 2),
    (input_id(UAV), 2) => (controller, 2),

    # connections
    (sensor, 1) => (controller, 1),
    (controller, 1) => (dynamics, 1),
    (dynamics, 1) => (sensor, 1),

    # net output
    (dynamics, 1) => (output_id(UAV), 1)
])

function 𝗟(𝐖)
    𝐿(u, x, p, t) = LVector( sc = -p.𝓐l * (u[1] - x[1] - x[2]) );
    𝐶(u, x, p, t) = LVector( sl = -p.𝓐c * (u[1] + p.𝓑c*x[1] - x[2]) );
    𝐷(u, x, p, t) = LVector( α = -0.313*u[1] +  56.7*u[2] +  0.232*x[1],
                             q = -0.0139*u[1] - 0.426*u[2] + 0.0203*x[1],
                             θ =  56.7*u[2]);

    u_𝐿(u) = [ u[1] ];  # outputs sl
    u_𝐶(u) = [ u[1] ];  # outputs sc
    u_𝐷(u) = [ u[3] ];  # outputs θ

    return oapply(𝐖,
                  Dict(:sensor     => ContinuousMachine{Float64}(2, 1, 1, 𝐿, u_𝐿),
                       :controller => ContinuousMachine{Float64}(2, 1, 1, 𝐶, u_𝐶),
                       :dynamics   => ContinuousMachine{Float64}(1, 3, 1, 𝐷, u_𝐷)));
end

# we then have to assign the behavior that ought to inhabit the boxes
𝑢ᵤₐᵥ = 𝗟(UAV)

# initial values
xₒ = LVector( e = 0.01,  # [e, d] -> [θ offset, 𝛿 control input]
              d = 0.05);

uₒ = LVector( sl = 0,
              sc = 0,
              α = 0,
              q = 0,
              θ = 0);

# integration interval
t = (0, 20);

parameters = (𝓐l = 100,  # decay constant of sensor
              𝓐c = 100,  # decay constant of controller
              𝓑c = 0);   # ratio of velocity to reference velocity

# then compute the total dynamics
solution = solve(ODEProblem(𝑢ᵤₐᵥ, uₒ, xₒ, t, parameters), alg_hints=[:stiff]);

# plot the behavior of the system based on the wiring diagram definition
fsize = 28
simulation = plot( solution.t,        # x values
                   [solution[1,:],    # y values [q is amplified]
                    solution[2,:],
                    solution[3,:],
                    solution[4,:] * 1e2,
                    solution[5,:]],

                   # labels
                   label  = ["sl" "sc" "α" "q" "θ"],
                   xlabel = "Time parameter",
                   ylabel = "Response",
                   title  = "Aircraft pitch behaviour",
                   size   = (2000,1400),

                   # font size
                   legendfontsize = fsize,
                   xtickfontsize  = fsize,
                   ytickfontsize  = fsize,
                   xguidefontsize = fsize,
                   yguidefontsize = fsize,
                   titlefontsize  = fsize,
                   lw = 4,

                   # margins
                   left_margin   = 15Plots.mm,
                   bottom_margin = 15Plots.mm);

display(simulation) # this can be removed if not running from the command line