text:
module Modeling

using Dates
using UnPack
using SciMLBase, OrdinaryDiffEq, DiffEqCallbacks
using ComponentArrays, RecursiveArrayTools

export f_cont!, f_disc!
export SystemDescriptor, SystemGroupDescriptor, NullSystemDescriptor, System, Model

############################# SystemDescriptor ############################

abstract type SystemDescriptor end #anything from which we can build a System

#every SystemDescriptor's init_x0 must return an AbstractVector{<:Real}, even if
#its inherently discrete and its f_cont! does nothing. this is ugly but ensures
#composability of Systems without the hassle of dealing automatically with empty
#state vector blocks, which is exacerbated by the need to assign views from the
#root System's state vector to its children System's state vectors
init_x0(::SystemDescriptor) = [0.0]
init_y0(::SystemDescriptor) = nothing #sytems are not required to have outputs
init_u0(::SystemDescriptor) = nothing #sytems are not required to have control inputs
init_d0(::SystemDescriptor) = nothing #systems are not required to have discrete states

############################# System ############################

#need the T type parameter for dispatch, the rest for type stability
#making System mutable does not hurt performance, because System instances are
#only instantiated upon initialization, so no runtime heap allocations
mutable struct System{T, X <: AbstractVector{<:Float64},
                    Y, U, D, P, S}
    ẋ::X #continuous state vector derivative
    x::X #continuous state vector (to be used as a buffer for f_cont! evaluation)
    y::Y #output state
    u::U #control inputs
    d::D #discrete state
    t::Base.RefValue{Float64} #this allows propagation of t updates down the subsystem hierarchy
    params::P
    subsystems::S
end

function System(c::SystemDescriptor, ẋ = init_x0(c), x = init_x0(c),
                    y = init_y0(c), u = init_u0(c), d = init_d0(c), t = Ref(0.0))

    params = c #assign the system descriptor itself as a System parameter
    subsystems = nothing
    System{map(typeof, (c, x, y, u, d, params, subsystems))...}(
                                    ẋ, x, y, u, d, t, params, subsystems)
end

#f_disc! is free to modify a Hybrid system's discrete state, control inputs and
#continuous state. if it modifies the latter, it must return true, false
#otherwise. no fallbacks are provided for safety reasons: if the intended
#f_cont! or f_disc! implementations for the System have the wrong interface, the
#dispatch will silently revert to the fallback, which does nothing and may not
#be obvious at all.
f_cont!(args...) = MethodError(f_cont!, args) |> throw
(f_disc!(args...)::Bool) = MethodError(f_disc!, args) |> throw

######################### NullSystem ############################

struct NullSystemDescriptor <: SystemDescriptor end

@inline f_cont!(::System{NullSystemDescriptor}, args...) = nothing
@inline (f_disc!(::System{NullSystemDescriptor}, args...)::Bool) = false

######################## SystemGroup #############################

#abstract type providing convenience methods for composite Systems

abstract type SystemGroupDescriptor <: SystemDescriptor end

Base.keys(g::SystemGroupDescriptor) = propertynames(g)
Base.values(g::SystemGroupDescriptor) = map(λ -> getproperty(g, λ), keys(g))

init_x0(g::SystemGroupDescriptor) = NamedTuple{keys(g)}(init_x0.(values(g))) |> ComponentVector
init_y0(g::SystemGroupDescriptor) = NamedTuple{keys(g)}(init_y0.(values(g)))
init_u0(g::SystemGroupDescriptor) = NamedTuple{keys(g)}(init_u0.(values(g)))
init_d0(g::SystemGroupDescriptor) = NamedTuple{keys(g)}(init_d0.(values(g)))

function System(g::SystemGroupDescriptor, ẋ = init_x0(g), x = init_x0(g),
                    y = init_y0(g), u = init_u0(g), d = init_d0(g), t = Ref(0.0))

    ss_list = Vector{System}()
    ss_labels = keys(g)

    for label in ss_labels
        push!(ss_list, System(map((λ)->getproperty(λ, label), (g, ẋ, x, y, u, d))..., t))
    end

    params = nothing
    subsystems = NamedTuple{ss_labels}(ss_list)

    System{map(typeof, (g, x, y, u, d, params, subsystems))...}(
                         ẋ, x, y, u, d, t, params, subsystems)

end

#default implementation calls f_cont! on all Node subsystems with the same
#arguments provided to the parent Node's System, then builds the NamedTuple
#can be overridden for specific SystemGroupDescriptor subtypes if needed
@inline @generated function f_cont!(sys::System{T}, args...) where {T<:SystemGroupDescriptor}

    # Core.print("Generated function called")
    ex_main = Expr(:block)

    #call f_cont! on each subsystem
    ex_calls = Expr(:block)
    ss_labels = fieldnames(T)
    for label in ss_labels
        push!(ex_calls.args,
            :(f_cont!(sys.subsystems[$(QuoteNode(label))], args...)))
    end

    #retrieve the y from each subsystem and build a tuple from them
    ex_tuple = Expr(:tuple)
    for label in ss_labels
        push!(ex_tuple.args,
            :(sys.subsystems[$(QuoteNode(label))].y))
    end

    #build a NamedTuple from the subsystem's labels and the constructed tuple
    ex_y = Expr(:call, Expr(:curly, NamedTuple, ss_labels), ex_tuple)

    #assign the result to the parent system's y
    ex_assign_y = Expr(:(=), :(sys.y), ex_y)

    #push everything into the main block expression
    push!(ex_main.args, ex_calls)
    push!(ex_main.args, ex_assign_y)
    push!(ex_main.args, :(return nothing))

    return ex_main

end

#default implementation calls f_disc! on all Node subsystems with the same
#arguments provided to the parent Node's System, then ORs their outputs.
#can be overridden for specific SystemGroupDescriptor subtypes if needed
@inline @generated function (f_disc!(sys::System{T}, args...)::Bool
    ) where {T<:SystemGroupDescriptor}

    # Core.print("Generated function called")
    ex = Expr(:block)
    push!(ex.args, :(x_mod = false))
    for label in fieldnames(T)
        #we need all f_disc! calls executed, so | must be used instead of ||
        push!(ex.args,
            :(x_mod = x_mod | f_disc!(sys.subsystems[$(QuoteNode(label))], args...)))
    end
    return ex

end

############################# Model ############################

#in this design, the t and x fields of m.sys behave only as temporary
#storage for f_cont! and f_disc! calls, so we have no guarantees about their
#status after a certain step. the only valid sources for t and x at any
#given moment is the integrator's t and u
struct Model{S <: System,
                   I <: OrdinaryDiffEq.ODEIntegrator,
                   L <: SavedValues}

    sys::S
    integrator::I
    log::L

    function Model(sys, args_c::Tuple = (), args_d::Tuple = ();
        solver = Tsit5(), t_start = 0.0, t_end = 10.0, y_saveat = Float64[],
        save_on = false, int_kwargs...)

        #save_on is set to false because we are not usually interested in saving
        #the naked state vector. the output saved by the SavingCallback is all
        #we need for insight
        saveat_arr = (y_saveat isa Real ? (t_start:y_saveat:t_end) : y_saveat)

        params = (sys = sys, args_c = args_c, args_d = args_d)

        log = SavedValues(Float64, typeof(sys.y))

        dcb = DiscreteCallback((u, t, integrator)->true, f_dcb!)
        scb = SavingCallback(f_scb, log, saveat = saveat_arr)
        cb_set = CallbackSet(dcb, scb)

        # x0 = copy(sys.x) #not needed, the integrator creates its own copy
        x0 = sys.x
        problem = ODEProblem{true}(f_update!, x0, (t_start, t_end), params)
        integrator = init(problem, solver; callback = cb_set, save_on, int_kwargs...)
        new{typeof(sys), typeof(integrator), typeof(log)}(sys, integrator, log)
    end
end

#these functions are better defined outside the constructor; closures seem to
#have some overhead (?)

#function barriers: the System is first extracted from integrator.p,
#then used as an argument in the call to the actual update & callback functions,
#forcing the compiler to specialize for the specific System subtype;
#accesing sys.x and sys.ẋ directly instead causes type instability
f_update!(ẋ, x, p, t) = f_update!(ẋ, x, t, p.sys, p.args_c)
f_scb(x, t, integrator) = f_scb(x, t, integrator.p.sys, integrator.p.args_c)
function f_dcb!(integrator)
    x = integrator.u; t = integrator.t; p = integrator.p
    x_modified = f_dcb!(x, t, p.sys, p.args_c, p.args_d)
    u_modified!(integrator, x_modified)
end

#in-place integrator update function
function f_update!(ẋ::X, x::X, t::Real, sys::System{T,X}, args_c) where {T, X}

    sys.x .= x
    sys.t[] = t
    f_cont!(sys, args_c...) #updates sys.ẋ and sys.y
    ẋ .= sys.ẋ

    return nothing
end

#DiscreteCallback function (called on every integration step). in addition to
#specific functionality implemented by f_disc!, this callback ensures that
#the System's internal x and y are up to date with the last integrator's
#solution step
function f_dcb!(x::X, t::Real, sys::System{T,X}, args_c, args_d) where {T,X}

    sys.x .= x #assign the updated integrator's state to the system's local continuous state
    sys.t[] = t #ditto for time

    #at this point sys.y holds the output from the last solver evaluation of
    #f_cont!, not the one corresponding to the updated x. with x up to date, we
    #can now compute the correct sys.y for this epoch
    f_cont!(sys, args_c...) #updates sys.y, but leaves sys.x unmodified
    #this could be commented for additional performance, at the cost of
    #a slight output inaccuracy at integration epochs

    #with the system's outputs up to date, call the discrete update function
    x_modified = f_disc!(sys, args_d...) #this may modify sys.x
    x .= sys.x #assign the (potentially modified) sys.x back to the integrator

    #note that if x is modified by f_disc!, it will not be reflected in sys.y
    #until the next call to f_dcb!

    return x_modified
end

#SavingCallback function, this gets called at the end of each step after f_disc!
function f_scb(::X, ::Real, sys::System{T,X}, args_c) where {T,X}
    return deepcopy(sys.y)
end

function Base.getproperty(m::Model, s::Symbol)
    if s === :t
        return m.integrator.t
    elseif s === :x
        return m.integrator.u
    elseif s === :y
        return m.sys.y
    elseif s === :u
        return m.sys.u
    elseif s ∈ (:sys, :integrator, :log)
        return getfield(m, s)
    else
        return getproperty(m.integrator, s)
    end
end

SciMLBase.step!(m::Model, args...) = step!(m.integrator, args...)

SciMLBase.solve!(m::Model) = solve!(m.integrator)

SciMLBase.get_proposed_dt(m::Model) = get_proposed_dt(m.integrator)

function SciMLBase.reinit!(m::Model, args...; kwargs...)

    #for an ODEIntegrator, the optional args... is simply a new initial
    #condition. if not specified, the original initial condition is used
    reinit!(m.integrator, args...; kwargs...)

    #grab the updated t and x from the integrator (in case they were reset by
    #the input arguments). this is not strictly necessary, since they are merely
    #buffers. just for consistency.
    m.sys.t[] = m.integrator.t
    m.sys.x .= m.integrator.u

    resize!(m.log.t, 1)
    resize!(m.log.saveval, 1)
    return nothing
end

#the following causes type instability and destroys performance:
# function f_update!(ẋ, x, p, t)
    # @unpack sys, args_c = p
    # sys.x .= x
    # sys.t[] = t
    # f_cont!(sys, args_c...)
    # ẋ = sys.ẋ
# end

# the reason seems to be that having sys stored in p obfuscates type inference.
# when unpacking sys, the compiler can no longer tell its type, and therefore
# has no knowledge of the types of sys.x, sys.dx, sys.y and sys.t. since these
# are being assigned to and read from, the type instability kills performance.

# this can be fixed by storing the x, dx and y fields of sys directly as entries
# of p. this probably fixes their types during construction, so when they are
# accessed later in the closure, the type instability is no longer an issue.

# however, this is redundant! we already have x, dx, y and t inside of sys. a
# more elegant alternative is simply to use a function barrier, first extract
# sys, then call another function using it as an argument. this forces the
# compiler to infer its type, and therefore it specializes the time-critical
# assignment statements to their actual types.

###########################################################################

#to try:

#in System, define and extend f_branch!

# #individual Component
# f_branch!(y, dx, x, u, t, sys, args...) = f_branch!(Val(has_input(sys)), y, dx, x, u, t, args...)
# f_branch!(::Val{true}, y, dx, x, u, t, sys, args...) = f_cont!(y, dx, x, u, t, sys, args...)
# f_cont!(::HasInput, y, dx, x ,u, t, sys, args...) = f_cont!(y, dx, x, u, t, sys, args...)
# f_cont!(::HasNoInput, y, dx, x, u, t, sys, args...) = f_cont!(y, dx, x, t, sys, args...)

# #for a AirframeGroup
# f_cont!(MaybeInput(S), MaybeOutput(S), y, dx, x, u, t, sys, args...)
# f_cont!(::HasInput, ::HasOutput, y, dx, x ,u, t, sys, args...)
# #now, this method needs to consider the possibility for each component that it
# #may have or not Input or Output. so it must do
# for (label, component) in zip(keys(C), values(C))
#     if MaybeInput(typeof(component)) #need tocheck, because if it has no input, u[label] will not exist!
#         f_cont!(y_cmp, dx_cmp, x_cmp, u_cmp, t, cmp, args...)
#     else
#         f_cont!(y_cmp, dx_cmp, x_cmp, t, cmp, args...)
#     end
# end

end


text:
using Raylib
using Raylib: rayvector, IsKeyDown
using StaticArrays

const G = 400
const PLAYER_JUMP_SPD = 350
const PLAYER_HOR_SPD = 200

mutable struct Player
    pos_x::Float64
    pos_y::Float64
    speed::Float64
    canJump::Bool
end
Player(pos, s, j) = Player(pos..., s, j)

position(p::Player) = rayvector(p.pos_x, p.pos_y)

struct EnvItem
    rect::Raylib.RayRectangle
    blocking::Bool
    color::Raylib.RayColor
end

function isoverlap(item, player, movement)
    rect = item.rect
    pos = position(player)
    return rect.x <= pos.x <= rect.x + rect.width &&
        rect.y - movement < pos.y <= rect.y
end

function UpdatePlayer!(player::Player, envItems::Vector{EnvItem}, delta::Float32)
    hor_movement = PLAYER_HOR_SPD*delta
    ver_movement = player.speed * delta

    IsKeyDown(Raylib.KEY_LEFT)  && (player.pos_x -= hor_movement)
    IsKeyDown(Raylib.KEY_RIGHT) && (player.pos_x += hor_movement)

    if IsKeyDown(Raylib.KEY_SPACE) && player.canJump
        player.speed = -PLAYER_JUMP_SPD
        player.canJump = false
    end

    hitObstacle = false
    for ei in envItems
        if ei.blocking && isoverlap(ei, player, ver_movement)
            player.speed = 0
            player.pos_y = ei.rect.y
            hitObstacle = true
        end
    end

    if !hitObstacle
        player.pos_y += ver_movement
        player.speed += G*delta
        player.canJump = false
    else
        player.canJump = true
    end

    return player
end

function UpdateCameraCenter(
    camera::Raylib.RayCamera2D,
    player::Player, _, _,
    width::Int64,
    height::Int64,
)
    camera.offset = rayvector(width/2, height/2)
    camera.target = position(player)
    return camera
end

function UpdateCameraCenterInsideMap(
    camera::Raylib.RayCamera2D,
    player::Player,
    envItems::Vector{EnvItem},
    _,
    width::Int64,
    height::Int64,
)
    camera = UpdateCameraCenter(camera, player, envItems, nothing, width, height)

    minX = minY = 1000; maxX = maxY = -1000

    minX = minimum(item->item.rect.x, envItems; init=1000)
    minY = minimum(item->item.rect.y, envItems; init=1000)
    maxX = maximum(item->item.rect.x + item.rect.width, envItems; init=-1000)
    maxY = maximum(item->item.rect.y + item.rect.height, envItems; init=-1000)

    maxbound = Raylib.GetWorldToScreen2D(rayvector(maxX, maxY), camera)
    minbound = Raylib.GetWorldToScreen2D(rayvector(minX, minY), camera)

    (maxbound.x < width)  && (camera.offset_x = width  - (maxbound.x - width/2))
    (maxbound.y < height) && (camera.offset_y = height - (maxbound.y - height/2))
    (minbound.x > 0) && (camera.offset_x = width/2  - minbound.x)
    (minbound.y > 0) && (camera.offset_y = height/2 - minbound.y)
    return camera
end

function UpdateCameraCenterSmoothFollow(
    camera::Raylib.RayCamera2D,
    player::Player,
    _,
    delta::Float32,
    width::Int64,
    height::Int64
)
    minSpeed = 30
    minEffectLength = 10
    fractionSpeed = 0.8

    camera.offset = rayvector(width/2, height/2)

    diff = position(player) - camera.target
    length = sqrt(sum(diff.^2))

    if length > minEffectLength
        speed = max(fractionSpeed*length, minSpeed);
        camera.target += speed*delta/length .* diff
    end

    return camera
end

mutable struct UpdateCameraEvenOutOnLandingArgv
    evenOutSpeed::Float32
    eveningOut::Bool
    evenOutTarget::Float32
end

const uceool_argv = UpdateCameraEvenOutOnLandingArgv(700, false, 0)

function UpdateCameraEvenOutOnLanding(
    camera::Raylib.RayCamera2D,
    player::Player,
    _,
    delta::Float32,
    width::Int64,
    height::Int64
)
    camera.offset = rayvector(width/2, height/2)
    camera.target = rayvector(player.pos_x, camera.target[2])

    if uceool_argv.eveningOut
        if uceool_argv.evenOutTarget > camera.target[2]
            camera.target += rayvector(0, uceool_argv.evenOutSpeed*delta)

            if camera.target[2] > uceool_argv.evenOutTarget
                camera.target = rayvector(camera.target[1], uceool_argv.evenOutTarget)
                uceool_argv.eveningOut = 0
            end
        else
            camera.target -= rayvector(0, uceool_argv.evenOutSpeed*delta)

            if camera.target[2] < uceool_argv.evenOutTarget
                camera.target = rayvector(camera.target[1], uceool_argv.evenOutTarget)
                uceool_argv.eveningOut = 0
            end
        end
    else
        if player.canJump && player.speed == 0 && player.pos_y != camera.target[2]
            uceool_argv.eveningOut = 1
            uceool_argv.evenOutTarget = player.pos_y
        end
    end

    return camera
end

function UpdateCameraPlayerBoundsPush(
    camera::Raylib.RayCamera2D,
    player::Player,
    _, _,
    width::Int64,
    height::Int64
)
    bbox = rayvector(0.2, 0.2)
    worldbox = 0.5rayvector(width, height)

    bboxWorldMin = Raylib.GetScreenToWorld2D(worldbox .* (1 .- bbox), camera)
    bboxWorldMax = Raylib.GetScreenToWorld2D(worldbox .* (1 .+ bbox), camera)
    camera.offset = worldbox .* (1 .- bbox)

    (player.pos_x < bboxWorldMin.x) &&  (camera.target_x = player.pos_x)
    (player.pos_y < bboxWorldMin.y) &&  (camera.target_y = player.pos_y)
    (player.pos_x > bboxWorldMax.x) &&  (camera.target_x = bboxWorldMin.x + player.pos_x - bboxWorldMax.x)
    (player.pos_y > bboxWorldMax.y) &&  (camera.target_y = bboxWorldMin.y + player.pos_y - bboxWorldMax.y)

    return camera
end

const screenWidth = 800
const screenHeight = 450

function main()
    Raylib.InitWindow(
        screenWidth, screenHeight,
        "raylib [core] example - 2d camera"
    )

    player = Player(400, 280, 0, false)
    envItems = [
        EnvItem(Raylib.RayRectangle(0, 0, 1000, 400), false, Raylib.LIGHTGRAY),
        EnvItem(Raylib.RayRectangle(0, 400, 1000, 200), true, Raylib.GRAY),
        EnvItem(Raylib.RayRectangle(300, 200, 400, 10), true, Raylib.GRAY),
        EnvItem(Raylib.RayRectangle(250, 300, 100, 10), true, Raylib.GRAY),
        EnvItem(Raylib.RayRectangle(650, 300, 100, 10), true, Raylib.GRAY),
    ]

    camera = Raylib.RayCamera2D(
        rayvector(screenWidth/2, screenHeight/2),
        position(player),
        0,
        1
    )
    cameraUpdaters = [
        UpdateCameraCenter,
        UpdateCameraCenterInsideMap,
        UpdateCameraCenterSmoothFollow,
        UpdateCameraEvenOutOnLanding,
        UpdateCameraPlayerBoundsPush
    ]
    cameraOption = 1
    cameraUpdatersLength = length(cameraUpdaters)
    cameraDescriptions = [
        "Follow player center",
        "Follow player center, but clamp to map edges",
        "Follow player center; smoothed",
        "Follow player center horizontally; updateplayer center vertically after landing",
        "Player push camera on getting too close to screen edge"
    ]

    Raylib.SetTargetFPS(60)
    while !Raylib.WindowShouldClose()
        # Update
        deltaTime = Raylib.GetFrameTime();

        UpdatePlayer!(player, envItems, deltaTime)

        zoom = camera.zoom + 0.05Raylib.GetMouseWheelMove()
        camera.zoom = max(min(3, zoom), 0.25)

        if Raylib.IsKeyPressed(Raylib.KEY_R)
            camera.zoom = 1
            player.pos_x, player.pos_y = (400, 280)
        end

        if Raylib.IsKeyPressed(Raylib.KEY_C)
            cameraOption = mod1(cameraOption+1, cameraUpdatersLength)
        end

        # Call update camera function by its pointer
        camera = cameraUpdaters[cameraOption](
            camera,
            player,
            envItems,
            deltaTime,
            screenWidth,
            screenHeight
        )

        # Draw
        Raylib.BeginDrawing()

        Raylib.ClearBackground(Raylib.LIGHTGRAY)

        Raylib.BeginMode2D(camera)

        for ei in envItems
            Raylib.DrawRectangleRec(ei.rect, ei.color)
        end

        playerRect = Raylib.RayRectangle(position(player) .- (20, 40), 40, 40)
        Raylib.DrawRectangleRec(playerRect, Raylib.RED)

        Raylib.EndMode2D()

        Raylib.DrawText("Controls:", 20, 20, 10, Raylib.BLACK);
        Raylib.DrawText("- Right/Left to move", 40, 40, 10, Raylib.DARKGRAY);
        Raylib.DrawText("- Space to jump", 40, 60, 10, Raylib.DARKGRAY);
        Raylib.DrawText("- Mouse Wheel to Zoom in-out, R to reset zoom", 40, 80, 10, Raylib.DARKGRAY);
        Raylib.DrawText("- C to change camera mode", 40, 100, 10, Raylib.DARKGRAY);
        Raylib.DrawText("Current camera mode:", 20, 120, 10, Raylib.BLACK);
        Raylib.DrawText(cameraDescriptions[cameraOption], 40, 140, 10, Raylib.DARKGRAY);

        Raylib.EndDrawing()
    end

    Raylib.CloseWindow()
end

ext: jl


text:
"""
Provides the [`render`](@ref) methods to write the documentation as LaTeX files
(`MIME"text/latex"`).
"""
module LaTeXWriter

import ...Documenter:
    Anchors,
    Builder,
    Documents,
    Expanders,
    Formats,
    Documenter,
    Utilities

import ..Writers: Writer, render

# TODO

function render(::Writer{Formats.LaTeX}, doc::Documents.Document)
    error("LaTeX rendering is unsupported.")
end

end


text:
@doc raw"""
    GraphManifoldType

This type represents the type of data on the graph that the [`GraphManifold`](@ref)
represents.
"""
abstract type GraphManifoldType end

@doc raw"""
    EdgeManifoldManifold <: GraphManifoldType

A type for a [`GraphManifold`](@ref) where the data is given on the edges.
"""
struct EdgeManifold <: GraphManifoldType end

@doc raw"""
    VectexGraphManifold <: GraphManifoldType

A type for a [`GraphManifold`](@ref) where the data is given on the vertices.
"""
struct VertexManifold <: GraphManifoldType end

@doc raw"""
    GraphManifold{G,𝔽,M,T} <: AbstractPowerManifold{𝔽,M,NestedPowerRepresentation}

Build a manifold, that is a [`PowerManifold`](@ref) of the [`Manifold`](@ref) `M` either on
the edges or vertices of a graph `G` depending on the [`GraphManifoldType`](@ref) `T`.

# Fields
* `G` is an `AbstractSimpleGraph`
* `M` is a [`Manifold`](@ref)
"""
struct GraphManifold{G<:AbstractGraph,𝔽,TM,T<:GraphManifoldType} <:
       AbstractPowerManifold{𝔽,TM,NestedPowerRepresentation}
    graph::G
    manifold::TM
end

function GraphManifold(
    g::G,
    M::TM,
    ::VertexManifold,
) where {G<:AbstractGraph,𝔽,TM<:Manifold{<:𝔽}}
    return GraphManifold{G,𝔽,TM,VertexManifold}(g, M)
end
function GraphManifold(
    g::G,
    M::TM,
    ::EdgeManifold,
) where {G<:AbstractGraph,𝔽,TM<:Manifold{<:𝔽}}
    return GraphManifold{G,𝔽,TM,EdgeManifold}(g, M)
end

const EdgeGraphManifold{𝔽} = GraphManifold{<:AbstractGraph,𝔽,<:Manifold{𝔽},EdgeManifold}
const VertexGraphManifold{𝔽} = GraphManifold{<:AbstractGraph,𝔽,<:Manifold{𝔽},VertexManifold}

@doc raw"""
    check_manifold_point(M::GraphManifold, p)

Check whether `p` is a valid point on the [`GraphManifold`](@ref), i.e. its length equals the number of vertices
(for [`VertexManifold`](@ref)s) or the number of edges (for [`EdgeManifold`](@ref)s) and that each element of `p`
passes the [`check_manifold_point`](@ref) test for the base manifold `M.manifold`.
"""
check_manifold_point(::GraphManifold, ::Any...)
function check_manifold_point(M::VertexGraphManifold, p; kwargs...)
    if size(p) != (nv(M.graph),)
        return DomainError(
            length(p),
            "The number of points in `x` ($(length(p))) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), nv(M.graph))
    return check_manifold_point(PM, p; kwargs...)
end
function check_manifold_point(M::EdgeGraphManifold, p; kwargs...)
    if size(p) != (ne(M.graph),)
        return DomainError(
            length(p),
            "The number of points in `x` ($(size(p))) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), ne(M.graph))
    return check_manifold_point(PM, p; kwargs...)
end

@doc raw"""
    check_tangent_vector(M::GraphManifold, p, X; check_base_point = true, kwargs...)

Check whether `p` is a valid point on the [`GraphManifold`](@ref), and
`X` it from its tangent space, i.e. its
length equals the number of vertices (for [`VertexManifold`](@ref)s) or
the number of edges (for [`EdgeManifold`](@ref)s) and that each element of `X`
together with its corresponding entry of `p` passes the
[`check_tangent_vector`](@ref) test for the base manifold `M.manifold`.
The optional parameter `check_base_point` indicates, whether to call [`check_manifold_point`](@ref)  for `p`.
"""
check_tangent_vector(::GraphManifold, ::Any...)
function check_tangent_vector(
    M::VertexGraphManifold,
    p,
    X;
    check_base_point=true,
    kwargs...,
)
    if check_base_point && size(p) != (nv(M.graph),)
        return DomainError(
            length(p),
            "The number of points in `x` ($(size(p)) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    if size(X) != (nv(M.graph),)
        return DomainError(
            length(X),
            "The number of points in `v` ($(size(X)) does not match the number of nodes in the graph ($(nv(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), nv(M.graph))
    return check_tangent_vector(PM, p, X; check_base_point=check_base_point, kwargs...)
end
function check_tangent_vector(M::EdgeGraphManifold, p, X; check_base_point=true, kwargs...)
    if check_base_point && size(p) != (ne(M.graph),)
        return DomainError(
            length(p),
            "The number of elements in `x` ($(size(p)) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    if size(X) != (ne(M.graph),)
        return DomainError(
            length(X),
            "The number of elements in `v` ($(size(X)) does not match the number of edges in the graph ($(ne(M.graph))).",
        )
    end
    PM = PowerManifold(M.manifold, NestedPowerRepresentation(), ne(M.graph))
    return check_tangent_vector(PM, p, X; check_base_point=check_base_point, kwargs...)
end

get_iterator(M::EdgeGraphManifold) = 1:ne(M.graph)
get_iterator(M::VertexGraphManifold) = 1:nv(M.graph)

@doc raw"""
    incident_log(M::GraphManifold, x)

Return the tangent vector on the (vertex) [`GraphManifold`](@ref), where at
each node the sum of the [`log`](@ref)s to incident nodes is computed.
For a `SimpleGraph`, an egde is interpreted as double edge in the corresponding
SimpleDiGraph

If the internal graph is a `SimpleWeightedGraph` the weighted sum of the
tangent vectors is computed.
"""
function incident_log(M::VertexGraphManifold, p)
    v = zero_tangent_vector(M, p)
    return incident_log!(M, v, p)
end

function incident_log!(M::VertexGraphManifold, X, p)
    rep_size = representation_size(M.manifold)
    for e in edges(M.graph)
        vw = _write(M, rep_size, X, src(e))
        X[src(e)] +=
            log(M.manifold, _read(M, rep_size, p, src(e)), _read(M, rep_size, p, dst(e)))
        if !is_directed(M.graph)
            X[dst(e)] += log(
                M.manifold,
                _read(M, rep_size, p, dst(e)),
                _read(M, rep_size, p, src(e)),
            )
        end
    end
    return X
end
function incident_log!(
    M::GraphManifold{<:AbstractSimpleWeightedGraph,𝔽,<:Manifold{𝔽},VertexManifold},
    X,
    p,
) where {𝔽}
    rep_size = representation_size(M.manifold)
    for e in edges(M.graph)
        X[src(e)] += (
            get_weight(M.graph, src(e), dst(e)) * log(
                M.manifold,
                _read(M, rep_size, p, src(e)),
                _read(M, rep_size, p, dst(e)),
            )
        )
        if !is_directed(M.graph)
            X[dst(e)] += (
                get_weight(M.graph, dst(e), src(e)) * log(
                    M.manifold,
                    _read(M, rep_size, p, dst(e)),
                    _read(M, rep_size, p, src(e)),
                )
            )
        end
    end
    return X
end

@doc raw"""
    manifold_dimension(N::GraphManifold{G,𝔽,M,VertexManifold})

returns the manifold dimension of the [`GraphManifold`](@ref) `N` on the vertices of
a graph $G=(V,E)$, i.e.
````math
\dim(\mathcal N) = \lvert V \rvert \dim(\mathcal M),
````
where $\mathcal M$ is the manifold of the data on the nodes.
"""
function manifold_dimension(M::VertexGraphManifold)
    return manifold_dimension(M.manifold) * nv(M.graph)
end
@doc raw"""
    manifold_dimension(N::GraphManifold{G,𝔽,M,EdgeManifold})

returns the manifold dimension of the [`GraphManifold`](@ref) `N` on the edges of
a graph $G=(V,E)$, i.e.
````math
\dim(\mathcal N) = \lvert E \rvert \dim(\mathcal M),
````
where $\mathcal M$ is the manifold of the data on the edges.
"""
function manifold_dimension(M::EdgeGraphManifold)
    return manifold_dimension(M.manifold) * ne(M.graph)
end

function _show_graph_manifold(io::IO, M; man_desc="", pre="")
    println(io, "GraphManifold\nGraph:")
    sg = sprint(show, "text/plain", M.graph, context=io, sizehint=0)
    sg = replace(sg, '\n' => "\n$(pre)")
    println(io, pre, sg)
    println(io, "Manifold$(man_desc):")
    sm = sprint(show, "text/plain", M.manifold, context=io, sizehint=0)
    sm = replace(sm, '\n' => "\n$(pre)")
    print(io, pre, sm)
    return nothing
end

function Base.show(io::IO, mime::MIME"text/plain", M::EdgeGraphManifold)
    return _show_graph_manifold(io, M; man_desc=" on edges", pre=" ")
end
function Base.show(io::IO, mime::MIME"text/plain", M::VertexGraphManifold)
    return _show_graph_manifold(io, M; man_desc=" on vertices", pre=" ")
end


text:
abstract AbstrModel
abstract AbstrUnivarModel <: AbstrModel
abstract AbstrMultivarModel <: AbstrModel

################
## resimulate ##
################

@doc doc"""
Resimulate data from given model.
"""->
function resimulate(data::AbstractTimenum,
                    mod::AbstrUnivarModel,
                    nPaths::Int = 1)
    ## resimulate data from given model

    ## data must be univariate for univariate models
    nObs, nAss = size(data)
    if nAss != 1
        error("Original data must be univariate for univariate models.")
    end

    ## get returns without NAs
    pureRets = dropna(data.vals[1])
    nRets = size(pureRets, 1)

    simRetsDf = DataFrame()
    noNAInds = !isna(data.vals[1])
    for ii=1:nPaths
        # simulate values
        simVals = simulate(mod, nRets)

        ## fill values into DataArray
        simRetsDa = DataArray(Float64, nObs)
        simRetsDa[noNAInds] = simVals

        ## copy values to DataFrame
        simRetsDf[ii] = simRetsDa
    end

    ## fix names
    modName = string(typeof(mod))
    if nPaths == 1
        colName = symbol(string(names(data.vals)[1],
                                "_sim_",
                                modName))
        rename!(simRetsDf, :x1, colName)
    else
        nams = [symbol(string(names(data.vals)[1],
                              "_sim_",
                              modName,
                              "_",
                              "$ii")) for ii=1:nPaths]
        names!(simRetsDf, nams)
    end
    simRets = Timenum(simRetsDf, idx(data))
    return simRets
end

## function resimulateLogPrices(logPricesTn::Timenum,
##                              mod::AbstrUnivarModel,
##                              nPaths::Int = 1)
##     ## simulate historic log price path again

##     nObs, nAss = size(logPricesTn)

##     # logPricesTn should be only one asset for univariate models
##     if nAss != 1
##         error("Original data must be exactly one asset for univariate models.")
##     end

##     # get logRets
##     logRets = price2ret(logPricesTn, log=true)

##     ## get returns without NAs
##     pureRets = dropna(logRets.vals[1])
##     nRets = size(pureRets, 1)

##     ## fit model to pure returns
##     modFit = estimate(mod, pureRets)

##     simRetsDf = DataFrame()
##     noNAInds = !isna(logRets.vals[1])
##     for ii=1:nPaths
##         # simulate values
##         simVals = simulate(modFit, nRets)

##         ## fill values into DataArray
##         simRetsDa = DataArray(Float64, nObs-1)
##         simRetsDa[noNAInds] = simVals

##         ## copy values to DataFrame
##         simRetsDf[ii] = simRetsDa
##     end

##     ## fix names
##     modName = string(typeof(modFit))
##     if nPaths == 1
##         colName = symbol(string(names(logPricesTn.vals)[1],
##                                 "_sim_",
##                                 modName))
##         rename!(simRetsDf, :x1, colName)
##     else
##         nams = [symbol(string(names(logPricesTn.vals)[1],
##                               "_sim_",
##                               modName,
##                               "_",
##                               "$ii")) for ii=1:nPaths]
##         names!(simRetsDf, nams)
##     end

##     simRets = Timenum(simRetsDf, idx(logRets))
##     simPrices = ret2price(simRets)
##     idx(simPrices)[1] = idx(logPricesTn)[1]

##     return (simRets, simPrices)
## end


text:
using WaveletRegression
using Wavelets
using Statistics
using Test

macro setupvars()
    ex = quote
        n = 128
        mx = 5
        my = 3
        wt = wavelet(WT.db2)
        nlevels = 4
        L = nlevels + 1
        trans = MODWT(wt, nlevels)
    end
    return esc(ex)
end

@testset "Utilities" begin
    @setupvars
    X = randn(n, mx)
    # MODWT
    Xw = modwt_matrix(X, trans)
    X1 = imodwt_matrix(Xw, trans)
    @test size(Xw) == (n, nlevels+1, mx)
    @test size(X1) == size(X)
    @test all(X .≈ X1)

    # DWT
    trans = DWT(wt, nlevels)
    Xw = dwt_matrix(X, trans)
    X1 = idwt_matrix(Xw, trans)
    @test size(Xw) == (n, mx)
    @test size(X1) == size(X)
    @test all(X .≈ X1)
end

# vector X and Y
@testset "Univariate timeseries" begin
    @setupvars
    X = randn(n)
    Y = randn(n)
    res1 = wavelm(X, Y, trans)
    @test length(res1.coefficients) == L
    @test all(size(B) == (1, 1) for B  in res1.coefficients)
    Y_pred = predict(res1)
    Y_pred = predict(res1, randn(length(X)))
    Y_pred = predict(res1, randn(length(X) + 1))
end

# matrix X and Y
@testset "Vector timeseries" begin
    @setupvars
    X = randn(n, mx)
    Y = randn(n, my)
    res1 = wavelm(X, Y, trans)
    @test length(res1.coefficients) == L
    @test all(size(B) == (mx, my) for B  in res1.coefficients)
    @test ncoef(res1) == mx * my * (L)
    @test dof(res1) == n - mx * my * (L) - 1
    @test size(regression_stderr(res1)) == (my, L)
    Y_pred = predict(res1)
    Y_pred_shrunk = predict(res1, true)
    @test Y_pred != Y_pred_shrunk
    @test all(Y_pred .!= Y_pred_shrunk)
    @test all(regression_stderr(res1) != regression_stderr(res1, shrink=true))
    Y_pred_newdata = predict(res1, randn(size(X)))
    Y_pred_newdata = predict(res1, randn(size(X, 1) + 1, size(X, 2)))
    SE = coef_stderr(res1)
    BB = coef(res1)
    BBshrunk = coef_shrunk(res1)
    @test all([B != Bs for (B, Bs) in zip(BB, BBshrunk)])
    res1shrunk = shrink(res1)
    @test all([B != Bs for (B, Bs) in zip(coef(res1), coef(res1shrunk))])
    @test length(SE) == length(BB) == L
    @test all([size(SE[i]) == size(BB[i]) for i in 1:L])
    @test all([all(se .>= 0) for se in SE])
    trans = DWT(wt, nlevels)
    res2 = wavelm(X, Y, trans)
end


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin: issuccessful, parsescenario!, Given, When, Then, ByLineParser, ScenarioStep, ParseOptions

@testset "Scenario             " begin
    @testset "Scenario has a Given step; the parsed scenario has a Given struct" begin
        text = """
        Scenario: Some description
            Given a precondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition")]
    end

    @testset "Scenario has a When step; the parsed scenario has a When struct" begin
        text = """
        Scenario: Some description
            When some action
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[When("some action")]
    end

    @testset "Scenario has a Then step; the parsed scenario has a Then struct" begin
        text = """
        Scenario: Some description
            Then a postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Then("a postcondition")]
    end

    @testset "Scenario has an And following a Given; the And step becomes a Given" begin
        text = """
        Scenario: Some description
            Given a precondition
              And another precondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition"),
                                                 Given("another precondition")]
    end

    @testset "Scenario has an And following a When; the And step becomes a When" begin
        text = """
        Scenario: Some description
            When some action
             And another action
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[When("some action"),
                                                 When("another action")]
    end

    @testset "Scenario has an And following a Then; the And step becomes a Then" begin
        text = """
        Scenario: Some description
            Then some postcondition
             And another postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Then("some postcondition"),
                                                 Then("another postcondition")]
    end

    @testset "Blank lines" begin
        text = """
        Scenario: Some description
            Given some precondition

            When some action

            Then some postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("some precondition"),
                                             When("some action"),
                                             Then("some postcondition"),
                                             ]
    end

    @testset "But and * keywords" begin
        @testset "But follows Given/When/Then; Each But step is same as the preceding" begin
            text = """
            Scenario: Some description
                Given some precondition
                But after given
                When some action
                But after when
                Then some postcondition
                But after then
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps == ScenarioStep[Given("some precondition"),
                                                 Given("after given"),
                                                 When("some action"),
                                                 When("after when"),
                                                 Then("some postcondition"),
                                                 Then("after then")]

        end

        @testset "* follows Given/When/Then; Each * step is same as the preceding" begin
            text = """
            Scenario: Some description
                Given some precondition
                * after given
                When some action
                * after when
                Then some postcondition
                * after then
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps == ScenarioStep[Given("some precondition"),
                                                 Given("after given"),
                                                 When("some action"),
                                                 When("after when"),
                                                 Then("some postcondition"),
                                                 Then("after then")]

        end

        @testset "List items as *; Items are the same type as the preceding step" begin
            text = """
            Scenario: Some description
                Given some precondition
                * item 1
                * item 2
                * item 3
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps == ScenarioStep[Given("some precondition"),
                                                 Given("item 1"),
                                                 Given("item 2"),
                                                 Given("item 3")]

        end
    end

    @testset "Scenario is not terminated by newline; EOF is also an OK termination" begin
        text = """
        Scenario: Some description
            Then some postcondition
            And another postcondition"""

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
    end

    @testset "A given step multiple spaces before the step description; The parsed given is stripped" begin
        text = """
        Scenario: Some description
            Given      a precondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition")]
    end

    @testset "A given step multiple spaces after the step description; The parsed given is stripped" begin
        text = """
        Scenario: Some description
            Given a precondition
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition")]
    end

    @testset "Malformed scenarios" begin
        @testset "And as a first step; Expected Given, When, or Then before that" begin
            text = """
            Scenario: Some description
                And another postcondition
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason == :leading_and
            @test result.expected == :specific_step
            @test result.actual == :and_step
        end

        @testset "Given after a When; Expected When or Then" begin
            text = """
            Scenario: Some description
                 When some action
                Given some precondition
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotGiven
            @test result.actual == :Given
        end

        @testset "Given after Then; Expected Then" begin
            text = """
            Scenario: Some description
                Then some postcondition
                Given some precondition
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotGiven
            @test result.actual == :Given
        end

        @testset "When after Then; Expected Then" begin
            text = """
            Scenario: Some description
                Then some postcondition
                When some action
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotWhen
            @test result.actual == :When
        end

        @testset "A step definition without text; Expected a valid step definition" begin
            text = """
            Scenario: Some description
                Given
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason == :invalid_step
            @test result.expected == :step_definition
            @test result.actual == :invalid_step_definition
        end

        @testset "Improper scenario header; Expected valid scenario header" begin
            text = """
            Scenario malformed: This is not a proper scenario header
                Given some precondition
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test !issuccessful(result)
            @test result.reason   == :invalid_scenario_header
            @test result.expected == :scenario_or_outline
            @test result.actual   == :invalid_header
        end
    end

    @testset "Lenient parsing" begin
        @testset "Allow arbitrary step order" begin
            @testset "Given after a When; Steps are When and Given" begin
                text = """
                Scenario: Some description
                    When some action
                    Given some precondition
                """

                byline = ByLineParser(text, ParseOptions(allow_any_step_order=true))
                result = parsescenario!(byline)

                @test issuccessful(result)
                scenario = result.value
                @test scenario.steps == ScenarioStep[When("some action"), Given("some precondition")]
            end

            @testset "When after a Then; Steps are Then and When" begin
                text = """
                Scenario: Some description
                    Then some postcondition
                    When some action
                """

                byline = ByLineParser(text, ParseOptions(allow_any_step_order=true))
                result = parsescenario!(byline)

                @test issuccessful(result)
                scenario = result.value
                @test scenario.steps == ScenarioStep[Then("some postcondition"), When("some action")]
            end

            @testset "Steps are Then/When/Given; Parse result is in that order" begin
                text = """
                Scenario: Some description
                    Then some postcondition
                    When some action
                    Given some precondition
                """

                byline = ByLineParser(text, ParseOptions(allow_any_step_order=true))
                result = parsescenario!(byline)

                @test issuccessful(result)
                scenario = result.value
                @test scenario.steps == ScenarioStep[Then("some postcondition"), When("some action"), Given("some precondition")]
            end
        end
    end

    @testset "Block text" begin
        @testset "Block text in a Given; Block text is present in step" begin
            text = """
            Scenario: Some description
                Given some precondition
                \"\"\"
                This is block text.
                There are two lines.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1].block_text == """
            This is block text.
            There are two lines."""
        end

        @testset "Another block text in a Given; Block text is present in step" begin
            text = """
            Scenario: Some description
                Given some precondition
                \"\"\"
                This is another block text.
                There are three lines.
                This is the last line.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1].block_text == """
            This is another block text.
            There are three lines.
            This is the last line."""
        end

        @testset "Block text in a When step; Block text is present in step" begin
            text = """
            Scenario: Some description
                When some action
                \"\"\"
                This is block text.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1] == When("some action"; block_text="""This is block text.""")
        end

        @testset "Block text in a Then step; Block text is present in step" begin
            text = """
            Scenario: Some description
                Then some postcondition
                \"\"\"
                This is block text.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1] == Then("some postcondition"; block_text="""This is block text.""")
        end

        @testset "Block text with a blank line; Block text is present in step" begin
            text = """
            Scenario: Some description
                Then some postcondition
                \"\"\"
                This is block text.

                This is another line.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1] == Then("some postcondition"; block_text="""This is block text.\n\nThis is another line.""")
        end
    end
end


text:
using DashBootstrapComponents, DashHtmlComponents

alert = dbc_alert([
    html_h4("Well done!", className="alert-heading"),
    html_p(
        "This is a success alert with loads of extra text in it. So much " *
        "that you can see how spacing within an alert works with this " *
        "kind of content.",
    ),
    html_hr(),
    html_p(
        "Let's put some more text down here, but remove the bottom margin",
        className="mb-0",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# # Parametric Heteroscedastic Model
#
# This example is a small extension of the standard GP regression problem, in which the
# observation noise variance is a function of the input.
# It is assumed to be a simple quadratic form, with a single unknown
# scaling parameter, in addition to the usual lengthscale and variance parameters
# of the GP.
# A point estimate of all parameters is obtained using type-II maximum likelihood,
# as per usual.

using AbstractGPs
using AbstractGPsMakie
using CairoMakie
using KernelFunctions
using LinearAlgebra
using Literate
using Optim
using ParameterHandling
using Random
using Zygote

# Specify simple GP:
build_gp(θ) = GP(0, θ.s * with_lengthscale(SEKernel(), θ.l));

# Observation variance is some scaling of $x^2$:
observation_variance(θ, x::AbstractVector{<:Real}) = Diagonal(θ.σ² .* x .^ 2);

# Specify hyperparameters:
flat_init_params, unflatten = ParameterHandling.value_flatten((
    s=positive(1.0), l=positive(3.0), σ²=positive(0.1)
));
θ_init = unflatten(flat_init_params);

# Build inputs:
const x = range(0.0, 10.0; length=100);
const y = rand(Xoshiro(123456), build_gp(θ_init)(x, observation_variance(θ_init, x)));

# Specify objective function:
function objective(θ)
    f = build_gp(θ)
    Σ = observation_variance(θ, x)
    return -logpdf(f(x, Σ), y)
end;

# Optimise the hyperparameters. They've been initialised near the correct values, so
# they ought not to deviate too far.
result = optimize(
    objective ∘ unflatten,
    θ -> only(Zygote.gradient(objective ∘ unflatten, θ)),
    flat_init_params + 0.01 * randn(Xoshiro(123456), length(flat_init_params)),
    LBFGS(;
        alphaguess=Optim.LineSearches.InitialStatic(; scaled=true),
        linesearch=Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(; iterations=4_000);
    inplace=false,
);
θ_final = unflatten(result.minimizer);

# Construct the posterior GP with the optimal model parameters:
Σ_obs_final = observation_variance(θ_final, x);
fx_final = build_gp(θ_final)(x, Σ_obs_final);
f_post = posterior(fx_final, y);

# Plot the results, making use of [AbstractGPsMakie](https://github.com/JuliaGaussianProcesses/AbstractGPsMakie.jl):
using CairoMakie.Makie.ColorSchemes: Set1_4

#! format: off
set_theme!(
    palette=(color=Set1_4,),
    patchcolor=(Set1_4[2], 0.2),
    Axis=(limits=((0, 10), nothing),),
)
#! format: on

let
    fig = Figure()
    ax = Axis(fig[1, 1])
    plot!(ax, x, f_post(x, Σ_obs_final); bandscale=3, label="posterior + noise", color=(:orange, 0.3))
    plot!(ax, x, f_post(x, 1e-9); bandscale=3, label="posterior")
    gpsample!(ax, x, f_post(x, 1e-9); samples=10, color=Set1_4[3])
    scatter!(ax, x, y; label="y")
    axislegend(ax)
    fig
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
using Distributions
using Base.Test

# Mathematica code to generate this.

### entropy[dist_] :=
###  TimeConstrained[
###   Expectation[-Log[PDF[dist, x]], x \[Distributed] dist], 5,
###   None];(* Stop after a few seconds, as this computation is very very \
### slow. *)
### test[\[Mu]_, \[Sigma]_, \[Xi]_] :=
###  Module[{
###    dist = MaxStableDistribution[\[Mu], \[Sigma], \[Xi]], cdf,
###    q = {10, 25, 50, 75, 90}/100,
###    ent, max, min, mean, median, var, skew, kurt,
###    points, pointsCDF, pointsPDF, quantiles,
###    code, printv
###    },
###   cdf = CDF[dist]@x;
###
###   (* Compute reference results. *)
###   ent = entropy[dist];
###   min = If[
###     Norm[\[Xi]] < 1 E - 15 || \[Xi] <
###       0, -\[Infinity], \[Mu] - \[Sigma]/\[Xi]];
###   max = If[
###     Norm[\[Xi]] < 1 E - 15 || \[Xi] >
###       0, \[Infinity], \[Mu] - \[Sigma]/\[Xi]];
###   mean = Mean[dist]; (* Can be Indeterminate. *)
###   median = Median[dist];
###   var = Variance[dist]; (* Can be Indeterminate. *)
###   skew = Skewness[dist]; (* Can be Indeterminate. *)
###   kurt = Kurtosis[dist] - 3.0; (* Can be Indeterminate. \
###     Mathematica returns the kurtosis, \
###     Distributions.jl returns the excess kurtosis. *)
###
###   points =
###    Table[x /. FindRoot[cdf - cdfv, {x, 1}][[1]], {cdfv, .1, .9, .1}];
###   pointsCDF = CDF[dist, #] & /@ points;
###   pointsPDF = PDF[dist, #] & /@ points;
###
###   quantiles = Quantile[dist, #] & /@ q;
###
###   (* Output test code. *)
###   printv[v_] :=
###    If[v == Infinity, "Inf",
###     If[v == -Infinity, "-Inf", TextString[N[v]]], "Inf"]; (*
###   Handles Indeterminate through If's fourth argument. *)
###
###   code = "d = GeneralizedExtremeValue(" <> printv[\[Mu]] <> ", " <>
###     printv[\[Sigma]] <> ", " <> printv[\[Xi]] <> ")\n"
###     <> "@test_approx_eq location(d) " <> printv[\[Mu]] <> "\n"
###     <> "@test_approx_eq scale(d) " <> printv[\[Sigma]] <> "\n"
###     <> "@test_approx_eq shape(d) " <> printv[\[Xi]] <> "\n"
###     <> "@test_approx_eq_eps maximum(d) " <> printv[max] <> " 1e-5\n"
###     <> "@test_approx_eq_eps minimum(d) " <> printv[min] <> " 1e-5\n"
###     <> If[ent != None,
###      "@test_approx_eq_eps entropy(d) " <> printv[ent] <> " 1e-5\n",
###      "", ""]
###     <> "@test_approx_eq_eps mean(d) " <> printv[mean] <> " 1e-5\n"
###     <> "@test_approx_eq_eps median(d) " <> printv[median] <> " 1e-5\n"
###     <> "@test_approx_eq_eps var(d) " <> printv[var] <> " 1e-5\n"
###     <> "@test_approx_eq_eps skewness(d) " <> printv[skew] <> " 1e-5\n"
###     <> "@test_approx_eq_eps kurtosis(d) " <> printv[kurt] <> " 1e-5\n"
###     <> ((
###         "@test_approx_eq_eps pdf(d, " <> printv[points[[#]]] <> ") " <>
###           printv[pointsPDF[[#]]] <> " 1e-5\n"
###          <>
###          "@test_approx_eq_eps cdf(d, " <> printv[points[[#]]] <>
###          ") " <> printv[pointsCDF[[#]]] <> " 1e-5\n"
###          <>
###          "@test_approx_eq_eps logpdf(d, " <> printv[points[[#]]] <>
###          ") " <> printv[Log[pointsPDF[[#]]]] <> " 1e-5\n"
###          <>
###          "@test_approx_eq_eps logcdf(d, " <> printv[points[[#]]] <>
###          ") " <> printv[Log[pointsCDF[[#]]]] <> " 1e-5\n"
###          <>
###          "@test_approx_eq_eps ccdf(d, " <> printv[points[[#]]] <>
###          ") " <> printv[1 - pointsCDF[[#]]] <> " 1e-5\n"
###          <>
###          "@test_approx_eq_eps logccdf(d, " <> printv[points[[#]]] <>
###          ") " <> printv[Log[1 - pointsCDF[[#]]]] <> " 1e-5\n"
###         ) & /@ Table[i, {i, 1, Length@points}])
###     <> ((
###         "@test_approx_eq_eps quantile(d, " <> printv[q[[#]]] <> ") " <>
###           printv[quantiles[[#]]] <> " 1e-5\n"
###         ) & /@ Table[i, {i, 1, Length@quantiles}]);
###
###   Return[code];
###   ];

d = GeneralizedExtremeValue(1., 1., 1.)
@test GeneralizedExtremeValue(1., 1, 1) == d
@test GeneralizedExtremeValue(1, 1, 1) == d
@test typeof(convert(GeneralizedExtremeValue{Float32}, d)) == GeneralizedExtremeValue{Float32}
@test typeof(convert(GeneralizedExtremeValue{Float32}, 1, 1, 1)) == GeneralizedExtremeValue{Float32}

# test[1., 1., 1.]
d = GeneralizedExtremeValue(1., 1., 1.)
@test_approx_eq location(d) 1.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) 1.
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) 0. 1e-5
@test_approx_eq_eps mean(d) Inf 1e-5
@test_approx_eq_eps median(d) 1.4427 1e-5
@test_approx_eq_eps var(d) Inf 1e-5
@test_approx_eq_eps skewness(d) Inf 1e-5
@test_approx_eq_eps kurtosis(d) Inf 1e-5
@test_approx_eq_eps pdf(d, 0.434294) 0.53019 1e-5
@test_approx_eq_eps cdf(d, 0.434294) 0.1 1e-5
@test_approx_eq_eps logpdf(d, 0.434294) -0.63452 1e-5
@test_approx_eq_eps logcdf(d, 0.434294) -2.30259 1e-5
@test_approx_eq_eps ccdf(d, 0.434294) 0.9 1e-5
@test_approx_eq_eps logccdf(d, 0.434294) -0.105361 1e-5
@test_approx_eq_eps pdf(d, 0.621335) 0.518058 1e-5
@test_approx_eq_eps cdf(d, 0.621335) 0.2 1e-5
@test_approx_eq_eps logpdf(d, 0.621335) -0.657668 1e-5
@test_approx_eq_eps logcdf(d, 0.621335) -1.60944 1e-5
@test_approx_eq_eps ccdf(d, 0.621335) 0.8 1e-5
@test_approx_eq_eps logccdf(d, 0.621335) -0.223144 1e-5
@test_approx_eq_eps pdf(d, 0.830584) 0.434865 1e-5
@test_approx_eq_eps cdf(d, 0.830584) 0.3 1e-5
@test_approx_eq_eps logpdf(d, 0.830584) -0.832719 1e-5
@test_approx_eq_eps logcdf(d, 0.830584) -1.20397 1e-5
@test_approx_eq_eps ccdf(d, 0.830584) 0.7 1e-5
@test_approx_eq_eps logccdf(d, 0.830584) -0.356675 1e-5
@test_approx_eq_eps pdf(d, 1.09136) 0.335835 1e-5
@test_approx_eq_eps cdf(d, 1.09136) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 1.09136) -1.09113 1e-5
@test_approx_eq_eps logcdf(d, 1.09136) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 1.09136) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 1.09136) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 1.4427) 0.240227 1e-5
@test_approx_eq_eps cdf(d, 1.4427) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 1.4427) -1.42617 1e-5
@test_approx_eq_eps logcdf(d, 1.4427) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 1.4427) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 1.4427) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 1.95762) 0.156566 1e-5
@test_approx_eq_eps cdf(d, 1.95762) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 1.95762) -1.85428 1e-5
@test_approx_eq_eps logcdf(d, 1.95762) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 1.95762) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 1.95762) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 2.80367) 0.0890519 1e-5
@test_approx_eq_eps cdf(d, 2.80367) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 2.80367) -2.41854 1e-5
@test_approx_eq_eps logcdf(d, 2.80367) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 2.80367) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 2.80367) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 4.48142) 0.0398344 1e-5
@test_approx_eq_eps cdf(d, 4.48142) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 4.48142) -3.22302 1e-5
@test_approx_eq_eps logcdf(d, 4.48142) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 4.48142) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 4.48142) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 9.49122) 0.00999075 1e-5
@test_approx_eq_eps cdf(d, 9.49122) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 9.49122) -4.6061 1e-5
@test_approx_eq_eps logcdf(d, 9.49122) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 9.49122) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 9.49122) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) 0.434294 1e-5
@test_approx_eq_eps quantile(d, 0.25) 0.721348 1e-5
@test_approx_eq_eps quantile(d, 0.5) 1.4427 1e-5
@test_approx_eq_eps quantile(d, 0.75) 3.47606 1e-5
@test_approx_eq_eps quantile(d, 0.9) 9.49122 1e-5

# test[0., 1., 0.]
d = GeneralizedExtremeValue(0., 1., 0.)
@test_approx_eq location(d) 0.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) 0.
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) -Inf 1e-5
@test_approx_eq_eps mean(d) 0.577216 1e-5
@test_approx_eq_eps median(d) 0.366513 1e-5
@test_approx_eq_eps var(d) 1.64493 1e-5
@test_approx_eq_eps skewness(d) 1.13955 1e-5
@test_approx_eq_eps kurtosis(d) 2.4 1e-5
@test_approx_eq_eps pdf(d, -0.834032) 0.230259 1e-5
@test_approx_eq_eps cdf(d, -0.834032) 0.1 1e-5
@test_approx_eq_eps logpdf(d, -0.834032) -1.46855 1e-5
@test_approx_eq_eps logcdf(d, -0.834032) -2.30259 1e-5
@test_approx_eq_eps ccdf(d, -0.834032) 0.9 1e-5
@test_approx_eq_eps logccdf(d, -0.834032) -0.105361 1e-5
@test_approx_eq_eps pdf(d, -0.475885) 0.321888 1e-5
@test_approx_eq_eps cdf(d, -0.475885) 0.2 1e-5
@test_approx_eq_eps logpdf(d, -0.475885) -1.13355 1e-5
@test_approx_eq_eps logcdf(d, -0.475885) -1.60944 1e-5
@test_approx_eq_eps ccdf(d, -0.475885) 0.8 1e-5
@test_approx_eq_eps logccdf(d, -0.475885) -0.223144 1e-5
@test_approx_eq_eps pdf(d, -0.185627) 0.361192 1e-5
@test_approx_eq_eps cdf(d, -0.185627) 0.3 1e-5
@test_approx_eq_eps logpdf(d, -0.185627) -1.01835 1e-5
@test_approx_eq_eps logcdf(d, -0.185627) -1.20397 1e-5
@test_approx_eq_eps ccdf(d, -0.185627) 0.7 1e-5
@test_approx_eq_eps logccdf(d, -0.185627) -0.356675 1e-5
@test_approx_eq_eps pdf(d, 0.0874216) 0.366516 1e-5
@test_approx_eq_eps cdf(d, 0.0874216) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 0.0874216) -1.00371 1e-5
@test_approx_eq_eps logcdf(d, 0.0874216) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 0.0874216) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 0.0874216) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 0.366513) 0.346574 1e-5
@test_approx_eq_eps cdf(d, 0.366513) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 0.366513) -1.05966 1e-5
@test_approx_eq_eps logcdf(d, 0.366513) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 0.366513) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 0.366513) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 0.671727) 0.306495 1e-5
@test_approx_eq_eps cdf(d, 0.671727) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 0.671727) -1.18255 1e-5
@test_approx_eq_eps logcdf(d, 0.671727) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 0.671727) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 0.671727) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 1.03093) 0.249672 1e-5
@test_approx_eq_eps cdf(d, 1.03093) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 1.03093) -1.38761 1e-5
@test_approx_eq_eps logcdf(d, 1.03093) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 1.03093) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 1.03093) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 1.49994) 0.178515 1e-5
@test_approx_eq_eps cdf(d, 1.49994) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 1.49994) -1.72308 1e-5
@test_approx_eq_eps logcdf(d, 1.49994) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 1.49994) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 1.49994) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 2.25037) 0.0948245 1e-5
@test_approx_eq_eps cdf(d, 2.25037) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 2.25037) -2.35573 1e-5
@test_approx_eq_eps logcdf(d, 2.25037) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 2.25037) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 2.25037) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) -0.834032 1e-5
@test_approx_eq_eps quantile(d, 0.25) -0.326634 1e-5
@test_approx_eq_eps quantile(d, 0.5) 0.366513 1e-5
@test_approx_eq_eps quantile(d, 0.75) 1.2459 1e-5
@test_approx_eq_eps quantile(d, 0.9) 2.25037 1e-5

# test[0., 1., 1.1]
d = GeneralizedExtremeValue(0., 1., 1.1)
@test_approx_eq location(d) 0.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) 1.1
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) -0.909091 1e-5
@test_approx_eq_eps mean(d) Inf 1e-5
@test_approx_eq_eps median(d) 0.451411 1e-5
@test_approx_eq_eps var(d) Inf 1e-5
@test_approx_eq_eps skewness(d) Inf 1e-5
@test_approx_eq_eps kurtosis(d) Inf 1e-5
@test_approx_eq_eps pdf(d, -2.43627) 0. 1e-5
@test_approx_eq_eps cdf(d, -2.43627) 0. 1e-5
@test_approx_eq_eps logpdf(d, -2.43627) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -2.43627) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -2.43627) 1. 1e-5
@test_approx_eq_eps logccdf(d, -2.43627) 0. 1e-5
@test_approx_eq_eps pdf(d, -1.75017) 0. 1e-5
@test_approx_eq_eps cdf(d, -1.75017) 0. 1e-5
@test_approx_eq_eps logpdf(d, -1.75017) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -1.75017) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -1.75017) 1. 1e-5
@test_approx_eq_eps logccdf(d, -1.75017) 0. 1e-5
@test_approx_eq_eps pdf(d, -1.06408) 0. 1e-5
@test_approx_eq_eps cdf(d, -1.06408) 0. 1e-5
@test_approx_eq_eps logpdf(d, -1.06408) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -1.06408) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -1.06408) 1. 1e-5
@test_approx_eq_eps logccdf(d, -1.06408) 0. 1e-5
@test_approx_eq_eps pdf(d, 0.091763) 0.332912 1e-5
@test_approx_eq_eps cdf(d, 0.091763) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 0.091763) -1.09988 1e-5
@test_approx_eq_eps logcdf(d, 0.091763) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 0.091763) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 0.091763) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 0.451411) 0.231581 1e-5
@test_approx_eq_eps cdf(d, 0.451411) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 0.451411) -1.46282 1e-5
@test_approx_eq_eps logcdf(d, 0.451411) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 0.451411) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 0.451411) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 0.99421) 0.146394 1e-5
@test_approx_eq_eps cdf(d, 0.99421) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 0.99421) -1.92145 1e-5
@test_approx_eq_eps logcdf(d, 0.99421) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 0.99421) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 0.99421) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 1.91649) 0.0803287 1e-5
@test_approx_eq_eps cdf(d, 1.91649) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 1.91649) -2.52163 1e-5
@test_approx_eq_eps logcdf(d, 1.91649) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 1.91649) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 1.91649) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 3.82421) 0.034286 1e-5
@test_approx_eq_eps cdf(d, 3.82421) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 3.82421) -3.37302 1e-5
@test_approx_eq_eps logcdf(d, 3.82421) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 3.82421) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 3.82421) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 9.89683) 0.00797749 1e-5
@test_approx_eq_eps cdf(d, 9.89683) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 9.89683) -4.83113 1e-5
@test_approx_eq_eps logcdf(d, 9.89683) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 9.89683) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 9.89683) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) -0.545871 1e-5
@test_approx_eq_eps quantile(d, 0.25) -0.274394 1e-5
@test_approx_eq_eps quantile(d, 0.5) 0.451411 1e-5
@test_approx_eq_eps quantile(d, 0.75) 2.67025 1e-5
@test_approx_eq_eps quantile(d, 0.9) 9.89683 1e-5

# test[0., 1., .6]
d = GeneralizedExtremeValue(0., 1., 0.6)
@test_approx_eq location(d) 0.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) 0.6
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) -1.66667 1e-5
@test_approx_eq_eps mean(d) 2.03027 1e-5
@test_approx_eq_eps median(d) 0.409936 1e-5
@test_approx_eq_eps var(d) Inf 1e-5
@test_approx_eq_eps skewness(d) Inf 1e-5
@test_approx_eq_eps kurtosis(d) Inf 1e-5
@test_approx_eq_eps pdf(d, -1.94901) 0. 1e-5
@test_approx_eq_eps cdf(d, -1.94901) 0. 1e-5
@test_approx_eq_eps logpdf(d, -1.94901) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -1.94901) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -1.94901) 1. 1e-5
@test_approx_eq_eps logccdf(d, -1.94901) 0. 1e-5
@test_approx_eq_eps pdf(d, -0.413975) 0.428261 1e-5
@test_approx_eq_eps cdf(d, -0.413975) 0.2 1e-5
@test_approx_eq_eps logpdf(d, -0.413975) -0.848022 1e-5
@test_approx_eq_eps logcdf(d, -0.413975) -1.60944 1e-5
@test_approx_eq_eps ccdf(d, -0.413975) 0.8 1e-5
@test_approx_eq_eps logccdf(d, -0.413975) -0.223144 1e-5
@test_approx_eq_eps pdf(d, -0.175663) 0.403746 1e-5
@test_approx_eq_eps cdf(d, -0.175663) 0.3 1e-5
@test_approx_eq_eps logpdf(d, -0.175663) -0.90697 1e-5
@test_approx_eq_eps logcdf(d, -0.175663) -1.20397 1e-5
@test_approx_eq_eps ccdf(d, -0.175663) 0.7 1e-5
@test_approx_eq_eps logccdf(d, -0.175663) -0.356675 1e-5
@test_approx_eq_eps pdf(d, 0.0897549) 0.347787 1e-5
@test_approx_eq_eps cdf(d, 0.0897549) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 0.0897549) -1.05617 1e-5
@test_approx_eq_eps logcdf(d, 0.0897549) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 0.0897549) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 0.0897549) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 0.409936) 0.278157 1e-5
@test_approx_eq_eps cdf(d, 0.409936) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 0.409936) -1.27957 1e-5
@test_approx_eq_eps logcdf(d, 0.409936) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 0.409936) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 0.409936) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 0.827268) 0.204827 1e-5
@test_approx_eq_eps cdf(d, 0.827268) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 0.827268) -1.58559 1e-5
@test_approx_eq_eps logcdf(d, 0.827268) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 0.827268) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 0.827268) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 1.42708) 0.134504 1e-5
@test_approx_eq_eps cdf(d, 1.42708) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 1.42708) -2.00616 1e-5
@test_approx_eq_eps logcdf(d, 1.42708) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 1.42708) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 1.42708) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 2.43252) 0.0725813 1e-5
@test_approx_eq_eps cdf(d, 2.43252) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 2.43252) -2.62305 1e-5
@test_approx_eq_eps logcdf(d, 2.43252) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 2.43252) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 2.43252) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 4.76379) 0.0245769 1e-5
@test_approx_eq_eps cdf(d, 4.76379) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 4.76379) -3.70595 1e-5
@test_approx_eq_eps logcdf(d, 4.76379) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 4.76379) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 4.76379) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) -0.656206 1e-5
@test_approx_eq_eps quantile(d, 0.25) -0.29662 1e-5
@test_approx_eq_eps quantile(d, 0.5) 0.409936 1e-5
@test_approx_eq_eps quantile(d, 0.75) 1.853 1e-5
@test_approx_eq_eps quantile(d, 0.9) 4.76379 1e-5

# test[0., 1., .3]
d = GeneralizedExtremeValue(0., 1., 0.3)
@test_approx_eq location(d) 0.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) 0.3
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) -3.33333 1e-5
@test_approx_eq_eps mean(d) 0.993518 1e-5
@test_approx_eq_eps median(d) 0.387422 1e-5
@test_approx_eq_eps var(d) 5.92458 1e-5
@test_approx_eq_eps skewness(d) 13.4836 5e-5 # Changed precision!
@test_approx_eq_eps kurtosis(d) Inf 1e-5
@test_approx_eq_eps pdf(d, -0.737875) 0.29572 1e-5
@test_approx_eq_eps cdf(d, -0.737875) 0.1 1e-5
@test_approx_eq_eps logpdf(d, -0.737875) -1.21834 1e-5
@test_approx_eq_eps logcdf(d, -0.737875) -2.30259 1e-5
@test_approx_eq_eps ccdf(d, -0.737875) 0.9 1e-5
@test_approx_eq_eps logccdf(d, -0.737875) -0.105361 1e-5
@test_approx_eq_eps pdf(d, -0.443476) 0.371284 1e-5
@test_approx_eq_eps cdf(d, -0.443476) 0.2 1e-5
@test_approx_eq_eps logpdf(d, -0.443476) -0.990787 1e-5
@test_approx_eq_eps logcdf(d, -0.443476) -1.60944 1e-5
@test_approx_eq_eps ccdf(d, -0.443476) 0.8 1e-5
@test_approx_eq_eps logccdf(d, -0.443476) -0.223144 1e-5
@test_approx_eq_eps pdf(d, -0.180553) 0.381877 1e-5
@test_approx_eq_eps cdf(d, -0.180553) 0.3 1e-5
@test_approx_eq_eps logpdf(d, -0.180553) -0.962658 1e-5
@test_approx_eq_eps logcdf(d, -0.180553) -1.20397 1e-5
@test_approx_eq_eps ccdf(d, -0.180553) 0.7 1e-5
@test_approx_eq_eps logccdf(d, -0.180553) -0.356675 1e-5
@test_approx_eq_eps pdf(d, 0.088578) 0.357029 1e-5
@test_approx_eq_eps cdf(d, 0.088578) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 0.088578) -1.02994 1e-5
@test_approx_eq_eps logcdf(d, 0.088578) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 0.088578) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 0.088578) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 0.387422) 0.310487 1e-5
@test_approx_eq_eps cdf(d, 0.387422) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 0.387422) -1.16961 1e-5
@test_approx_eq_eps logcdf(d, 0.387422) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 0.387422) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 0.387422) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 0.744195) 0.250557 1e-5
@test_approx_eq_eps cdf(d, 0.744195) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 0.744195) -1.38407 1e-5
@test_approx_eq_eps logcdf(d, 0.744195) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 0.744195) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 0.744195) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 1.20814) 0.183254 1e-5
@test_approx_eq_eps cdf(d, 1.20814) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 1.20814) -1.69688 1e-5
@test_approx_eq_eps logcdf(d, 1.20814) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 1.20814) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 1.20814) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 1.89428) 0.113828 1e-5
@test_approx_eq_eps cdf(d, 1.89428) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 1.89428) -2.17307 1e-5
@test_approx_eq_eps logcdf(d, 1.89428) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 1.89428) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 1.89428) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 3.21416) 0.0482752 1e-5
@test_approx_eq_eps cdf(d, 3.21416) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 3.21416) -3.03084 1e-5
@test_approx_eq_eps logcdf(d, 3.21416) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 3.21416) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 3.21416) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) -0.737875 1e-5
@test_approx_eq_eps quantile(d, 0.25) -0.311141 1e-5
@test_approx_eq_eps quantile(d, 0.5) 0.387422 1e-5
@test_approx_eq_eps quantile(d, 0.75) 1.51068 1e-5
@test_approx_eq_eps quantile(d, 0.9) 3.21416 1e-5

# test[1., 1., -1.]
d = GeneralizedExtremeValue(1., 1., -1.)
@test_approx_eq location(d) 1.
@test_approx_eq scale(d) 1.
@test_approx_eq shape(d) -1.
@test_approx_eq_eps maximum(d) 2. 1e-5
@test_approx_eq_eps minimum(d) -Inf 1e-5
@test_approx_eq_eps mean(d) 1. 1e-5
@test_approx_eq_eps median(d) 1.30685 1e-5
@test_approx_eq_eps var(d) 1. 1e-5
@test_approx_eq_eps skewness(d) -2. 1e-5
@test_approx_eq_eps kurtosis(d) 6. 1e-5
@test_approx_eq_eps pdf(d, -0.302585) 0.1 1e-5
@test_approx_eq_eps cdf(d, -0.302585) 0.1 1e-5
@test_approx_eq_eps logpdf(d, -0.302585) -2.30259 1e-5
@test_approx_eq_eps logcdf(d, -0.302585) -2.30259 1e-5
@test_approx_eq_eps ccdf(d, -0.302585) 0.9 1e-5
@test_approx_eq_eps logccdf(d, -0.302585) -0.105361 1e-5
@test_approx_eq_eps pdf(d, 0.390562) 0.2 1e-5
@test_approx_eq_eps cdf(d, 0.390562) 0.2 1e-5
@test_approx_eq_eps logpdf(d, 0.390562) -1.60944 1e-5
@test_approx_eq_eps logcdf(d, 0.390562) -1.60944 1e-5
@test_approx_eq_eps ccdf(d, 0.390562) 0.8 1e-5
@test_approx_eq_eps logccdf(d, 0.390562) -0.223144 1e-5
@test_approx_eq_eps pdf(d, 0.796027) 0.3 1e-5
@test_approx_eq_eps cdf(d, 0.796027) 0.3 1e-5
@test_approx_eq_eps logpdf(d, 0.796027) -1.20397 1e-5
@test_approx_eq_eps logcdf(d, 0.796027) -1.20397 1e-5
@test_approx_eq_eps ccdf(d, 0.796027) 0.7 1e-5
@test_approx_eq_eps logccdf(d, 0.796027) -0.356675 1e-5
@test_approx_eq_eps pdf(d, 1.08371) 0.4 1e-5
@test_approx_eq_eps cdf(d, 1.08371) 0.4 1e-5
@test_approx_eq_eps logpdf(d, 1.08371) -0.916291 1e-5
@test_approx_eq_eps logcdf(d, 1.08371) -0.916291 1e-5
@test_approx_eq_eps ccdf(d, 1.08371) 0.6 1e-5
@test_approx_eq_eps logccdf(d, 1.08371) -0.510826 1e-5
@test_approx_eq_eps pdf(d, 1.30685) 0.5 1e-5
@test_approx_eq_eps cdf(d, 1.30685) 0.5 1e-5
@test_approx_eq_eps logpdf(d, 1.30685) -0.693147 1e-5
@test_approx_eq_eps logcdf(d, 1.30685) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, 1.30685) 0.5 1e-5
@test_approx_eq_eps logccdf(d, 1.30685) -0.693147 1e-5
@test_approx_eq_eps pdf(d, 1.48917) 0.6 1e-5
@test_approx_eq_eps cdf(d, 1.48917) 0.6 1e-5
@test_approx_eq_eps logpdf(d, 1.48917) -0.510826 1e-5
@test_approx_eq_eps logcdf(d, 1.48917) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, 1.48917) 0.4 1e-5
@test_approx_eq_eps logccdf(d, 1.48917) -0.916291 1e-5
@test_approx_eq_eps pdf(d, 1.64333) 0.7 1e-5
@test_approx_eq_eps cdf(d, 1.64333) 0.7 1e-5
@test_approx_eq_eps logpdf(d, 1.64333) -0.356675 1e-5
@test_approx_eq_eps logcdf(d, 1.64333) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, 1.64333) 0.3 1e-5
@test_approx_eq_eps logccdf(d, 1.64333) -1.20397 5e-5 # Changed precision!
@test_approx_eq_eps pdf(d, 2.17463) 0. 1e-5
@test_approx_eq_eps cdf(d, 2.17463) 1. 1e-5
@test_approx_eq_eps logpdf(d, 2.17463) -Inf 1e-5
@test_approx_eq_eps logcdf(d, 2.17463) 0. 1e-5
@test_approx_eq_eps ccdf(d, 2.17463) 0. 1e-5
@test_approx_eq_eps logccdf(d, 2.17463) -Inf 1e-5
@test_approx_eq_eps pdf(d, 2.44645) 0. 1e-5
@test_approx_eq_eps cdf(d, 2.44645) 1. 1e-5
@test_approx_eq_eps logpdf(d, 2.44645) -Inf 1e-5
@test_approx_eq_eps logcdf(d, 2.44645) 0. 1e-5
@test_approx_eq_eps ccdf(d, 2.44645) 0. 1e-5
@test_approx_eq_eps logccdf(d, 2.44645) -Inf 1e-5
@test_approx_eq_eps quantile(d, 0.1) -0.302585 1e-5
@test_approx_eq_eps quantile(d, 0.25) 0.613706 1e-5
@test_approx_eq_eps quantile(d, 0.5) 1.30685 1e-5
@test_approx_eq_eps quantile(d, 0.75) 1.71232 1e-5
@test_approx_eq_eps quantile(d, 0.9) 1.89464 1e-5

# test[-1, 0.5, 0.6]
d = GeneralizedExtremeValue(-1., 0.5, 0.6)
@test_approx_eq location(d) -1.
@test_approx_eq scale(d) 0.5
@test_approx_eq shape(d) 0.6
@test_approx_eq_eps maximum(d) Inf 1e-5
@test_approx_eq_eps minimum(d) -1.83333 1e-5
@test_approx_eq_eps mean(d) 0.015133 1e-5
@test_approx_eq_eps median(d) -0.795032 1e-5
@test_approx_eq_eps var(d) Inf 1e-5
@test_approx_eq_eps skewness(d) Inf 1e-5
@test_approx_eq_eps kurtosis(d) Inf 1e-5
@test_approx_eq_eps pdf(d, -10.5807) 0. 1e-5
@test_approx_eq_eps cdf(d, -10.5807) 0. 1e-5
@test_approx_eq_eps logpdf(d, -10.5807) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -10.5807) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -10.5807) 1. 1e-5
@test_approx_eq_eps logccdf(d, -10.5807) 0. 1e-5
@test_approx_eq_eps pdf(d, -9.09221) 0. 1e-5
@test_approx_eq_eps cdf(d, -9.09221) 0. 1e-5
@test_approx_eq_eps logpdf(d, -9.09221) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -9.09221) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -9.09221) 1. 1e-5
@test_approx_eq_eps logccdf(d, -9.09221) 0. 1e-5
@test_approx_eq_eps pdf(d, -7.60374) 0. 1e-5
@test_approx_eq_eps cdf(d, -7.60374) 0. 1e-5
@test_approx_eq_eps logpdf(d, -7.60374) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -7.60374) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -7.60374) 1. 1e-5
@test_approx_eq_eps logccdf(d, -7.60374) 0. 1e-5
@test_approx_eq_eps pdf(d, -6.11528) 0. 1e-5
@test_approx_eq_eps cdf(d, -6.11528) 0. 1e-5
@test_approx_eq_eps logpdf(d, -6.11528) -Inf 1e-5
@test_approx_eq_eps logcdf(d, -6.11528) -Inf 1e-5
@test_approx_eq_eps ccdf(d, -6.11528) 1. 1e-5
@test_approx_eq_eps logccdf(d, -6.11528) 0. 1e-5
@test_approx_eq_eps pdf(d, -0.795032) 0.556315 1e-5
@test_approx_eq_eps cdf(d, -0.795032) 0.5 1e-5
@test_approx_eq_eps logpdf(d, -0.795032) -0.586421 1e-5
@test_approx_eq_eps logcdf(d, -0.795032) -0.693147 1e-5
@test_approx_eq_eps ccdf(d, -0.795032) 0.5 1e-5
@test_approx_eq_eps logccdf(d, -0.795032) -0.693147 1e-5
@test_approx_eq_eps pdf(d, -0.586366) 0.409654 1e-5
@test_approx_eq_eps cdf(d, -0.586366) 0.6 1e-5
@test_approx_eq_eps logpdf(d, -0.586366) -0.892442 1e-5
@test_approx_eq_eps logcdf(d, -0.586366) -0.510826 1e-5
@test_approx_eq_eps ccdf(d, -0.586366) 0.4 1e-5
@test_approx_eq_eps logccdf(d, -0.586366) -0.916291 1e-5
@test_approx_eq_eps pdf(d, -0.286458) 0.269007 1e-5
@test_approx_eq_eps cdf(d, -0.286458) 0.7 1e-5
@test_approx_eq_eps logpdf(d, -0.286458) -1.31302 1e-5
@test_approx_eq_eps logcdf(d, -0.286458) -0.356675 1e-5
@test_approx_eq_eps ccdf(d, -0.286458) 0.3 1e-5
@test_approx_eq_eps logccdf(d, -0.286458) -1.20397 1e-5
@test_approx_eq_eps pdf(d, 0.216262) 0.145163 1e-5
@test_approx_eq_eps cdf(d, 0.216262) 0.8 1e-5
@test_approx_eq_eps logpdf(d, 0.216262) -1.9299 1e-5
@test_approx_eq_eps logcdf(d, 0.216262) -0.223144 1e-5
@test_approx_eq_eps ccdf(d, 0.216262) 0.2 1e-5
@test_approx_eq_eps logccdf(d, 0.216262) -1.60944 1e-5
@test_approx_eq_eps pdf(d, 1.3819) 0.0491538 1e-5
@test_approx_eq_eps cdf(d, 1.3819) 0.9 1e-5
@test_approx_eq_eps logpdf(d, 1.3819) -3.0128 1e-5
@test_approx_eq_eps logcdf(d, 1.3819) -0.105361 1e-5
@test_approx_eq_eps ccdf(d, 1.3819) 0.1 1e-5
@test_approx_eq_eps logccdf(d, 1.3819) -2.30259 1e-5
@test_approx_eq_eps quantile(d, 0.1) -1.3281 1e-5
@test_approx_eq_eps quantile(d, 0.25) -1.14831 1e-5
@test_approx_eq_eps quantile(d, 0.5) -0.795032 1e-5
@test_approx_eq_eps quantile(d, 0.75) -0.0735019 1e-5
@test_approx_eq_eps quantile(d, 0.9) 1.3819 1e-5


text:
@testset "Errors" begin
    err = BioviewError("sorry :(")
    iobuf = IOBuffer(25)
    show(iobuf, err)

    @test String(iobuf.data) == "\"Bioview error: sorry :(\""
end


text:
using Random

function sample_position_params!(data::BmmData, shape_prior::ShapePrior)::MvNormalF
    μ = sample_center!(data) # TODO: should we store it as non-static arrays from the beginning?
    Σ = Array(Diagonal(shuffle(sample_var(shape_prior))))

    return MvNormalF(μ, Σ)
end

function sample_center!(data::BmmData; cache_size::Int=10000)
    if length(data.center_sample_cache) == 0 # otherwise, sampling takes too long
        data.center_sample_cache = sample(1:size(data.x, 1), Weights(confidence(data)), cache_size)
    end

    return position_data(data)[:, pop!(data.center_sample_cache)]
end

function sample_composition_params(data::BmmData)
    # TODO: this doesn't work well on small samples.
    # A simple workaround could be to further sample from the smoothed probabilities after sampling an existing distribution.
    gene_counts = deepcopy(sample(data.components).composition_params.counts)
    return CategoricalSmoothed(gene_counts, smooth=data.distribution_sampler.composition_params.smooth, sum_counts=sum(gene_counts));
end

function sample_distribution!(data::BmmData)
    data.max_component_guid += 1
    shape_prior = data.distribution_sampler.shape_prior
    position_params = sample_position_params!(data, shape_prior);
    composition_params = sample_composition_params(data);

    return Component(position_params, composition_params; shape_prior=deepcopy(shape_prior), guid=data.max_component_guid);
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
using TopOpt, LinearAlgebra, StatsFuns
using TopOptMakie
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk
using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 4.0 # filter radius
xmin = 0.0001 # minimum density
problem_size = (30, 10, 10)
V = 0.5 # maximum volume fraction
p = 4.0 # penalty

x0 = fill(V, prod(problem_size)); # initial design

using Ferrite
using TopOpt
using TopOpt.TopOptProblems: RectilinearGrid, Metadata
using TopOpt.TopOptProblems: left, right, bottom, middley, middlez,
    nnodespercell, nfacespercell, find_black_and_white, find_varind
using TopOpt.Utilities: @params

@params struct NewPointLoadCantilever{dim, T, N, M} <: StiffnessTopOptProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T}
    load_dict::Dict{Int, Vector{T}}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int}
    metadata::Metadata
end

function NewPointLoadCantilever(::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim},
    E = 1.0, ν = 0.3, force = 1.0) where {dim, CellType}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) || throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    if CellType === :Linear || dim === 3
        rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes))
    else
        rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes))
    end

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    if dim == 3
        addnodeset!(rect_grid.grid, "down_force", x -> right(rect_grid, x) &&
            bottom(rect_grid, x));
            #  && middlez(rect_grid, x));
    else
        addnodeset!(rect_grid.grid, "down_force", x -> right(rect_grid, x) &&
            right(rect_grid, x) && middley(rect_grid, x));
    end

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    if CellType === :Linear || dim === 3
        push!(dh, :u, dim) # Add a displacement field
    else
        ip = Lagrange{2, RefCube, 2}()
        push!(dh, :u, dim, ip) # Add a displacement field
    end
    close!(dh)

    ch = ConstraintHandler(dh)

    dbc = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    Ferrite.update!(ch, t)

    metadata = Metadata(dh)
    load_dict = Dict{Int, Vector{T}}()
    for fnode in getnodeset(rect_grid.grid, "down_force")
    	load_dict[fnode] = [0, -force, 0]
    end

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return NewPointLoadCantilever(rect_grid, E, ν, ch, load_dict, black, white, varind, metadata)
end

# used in FEA to determine default quad order
# we don't assume the problem struct has `rect_grid` to define its grid
TopOptProblems.nnodespercell(p::NewPointLoadCantilever) = nnodespercell(p.rect_grid)

# ! important, used for specification!
function TopOptProblems.getcloaddict(p::NewPointLoadCantilever{dim, T}) where {dim, T}
    # f = T[0, -p.force, 0]
    # fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1]
    # return Dict{Int, Vector{T}}(fnode => f)
    return p.load_dict
end

problem = NewPointLoadCantilever(Val{:Linear}, problem_size, (1.0, 1.0, 1.0), E, v, f);
solver = FEASolver(Displacement, Direct, problem, xmin = xmin);
solver()
u0 = solver.u

cheqfilter = DensityFilter(solver, rmin = rmin)
stress = TopOpt.MicroVonMisesStress(solver)
comp = TopOpt.Compliance(problem, solver)

# minimize compliance
function obj(x)
    return comp(cheqfilter(x))
end

# volume bound
function constr(x)
    return sum(cheqfilter(x)) / length(x) - V
end;

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr)

## enable printout
# Nonconvex.show_residuals[] = true

options = MMAOptions(
    maxiter=300, tol = Tolerance(kkt = 1e-3, x=1e-3, f = 1e-3),
)
TopOpt.setpenalty!(solver, p)

@time r = Nonconvex.optimize(
    m, MMA87(dualoptimizer = ConjugateGradient()),
    x0, options = options,
);

# show the optimized results!
@show obj(r.minimizer)
@show constr(r.minimizer)
@show maximum(stress(cheqfilter(r.minimizer)))

topology = cheqfilter(r.minimizer);

mesh = VTKUnstructuredData(problem)
topology = round.(topology)
inds = findall(isequal(0), topology)
deleteat!(mesh.cell_connectivity, inds)
deleteat!(mesh.cell_types, inds)

pymesh = PyVTK(mesh)
dash_vtk_utils = pyimport("dash_vtk.utils")
dash_vtk_utils.to_mesh_state
state = dash_vtk_utils.to_mesh_state(pymesh)

content = vtk_view([
        vtk_geometryrepresentation([
        vtk_mesh(state=state)
    ]),
])

app = dash()

app.layout = html_div() do
    html_div(
    style = Dict(
        "width" => "100%",
        "height" => "400px"
    ),
    children=[content],
    )
end

run_server(app, "0.0.0.0", debug = true)