text:
macro import_axis_huge()
    :(
        using ComplexVisual:
            CV_TickLabel, cv_format_ticks, ⇒,
            CV_TickLabelAppearance, CV_Ruler, ↦,
            cv_create_2daxis_canvas, cv_ticks_labels
    )
end

import Base: show

"""
```
struct CV_TickLabel{LocT}
    location    :: LocT
    text        :: AbstractString
end
```
A location (math coordinate) and a description for a tick
(typically placed at axis).
"""
struct CV_TickLabel{LocT}
    location    :: LocT
    text        :: AbstractString
end

"""
```
⇒(label, location)
    label     AbstractString
    location  Real
```

short for `CV_TickLabel(Flaot64(location), label)`.
"""
⇒(label::AbstractString, location::Real) =
    CV_TickLabel(Float64(location), label)

show(io::IO, tl::CV_TickLabel) = cv_show_impl(io, tl)
show(io::IO, m::MIME{Symbol("text/plain")}, tl::CV_TickLabel) =
    cv_show_impl(io, m, tl)

"""
Styles and data needed to render/draw axis-ticks with labels.
"""
struct CV_TickLabelAppearance{tsT <: CV_ContextStyle, lsT <: CV_ContextStyle} # {{{
    tick_length  :: Int32
    gap          :: Int32
    tick_style   :: tsT
    label_style  :: lsT

    function CV_TickLabelAppearance(tick_length::Integer, gap::Integer,
        tick_style::CV_ContextStyle, label_style::CV_ContextStyle)
        return new{typeof(tick_style), typeof(label_style)}(
            Int32(tick_length), Int32(gap), tick_style, label_style)
    end
end

show(io::IO, app::CV_TickLabelAppearance) = cv_show_impl(io, app)

show(io::IO, m::MIME{Symbol("text/plain")}, app::CV_TickLabelAppearance) =
    cv_show_impl(io, m, app)
# }}}

"""
```
CV_TickLabelAppearance(; tick_length, gap, tick_style, label_style)

tick_length   Integer          Int32(10)
gap           Integer          cv_half(tick_length)
tick_style    CV_ContextStyle  cv_linewidth(2) → cv_black
label_style   CV_ContextStyle  cv_black → cv_fontface("serif") → cv_fontsize(15)
```
"""
function CV_TickLabelAppearance(;
        tick_length::Integer=Int32(10),
        gap::Integer=cv_half(tick_length),
        tick_style::CV_ContextStyle=cv_linewidth(2) → cv_black,
        label_style::CV_ContextStyle=cv_black →
            cv_fontface("serif") → cv_fontsize(15))
    return CV_TickLabelAppearance(tick_length, gap, tick_style, label_style)
end

"""
```
CV_TickLabelAppearance(old; tick_length, gap, tick_style, label_style)

old           CV_TickLabelAppearance
tick_length   Integer                 old.tick_length
gap           Integer                 old.gap
tick_style    CV_ContextStyle         old.tick_style
label_style   CV_ContextStyle         old.label_style
```
"""
function CV_TickLabelAppearance(old::CV_TickLabelAppearance;
        tick_length::Integer=old.tick_length, gap::Integer=old.gap,
        tick_style::CV_ContextStyle=old.tick_style,
        label_style::CV_ContextStyle=old.label_style)
    return CV_TickLabelAppearance(; tick_length, gap, tick_style, label_style)
end

"""
Ticks with their labels and appearance.

```
CV_Ruler(ticklabels[, app])

ticklabels  NTuple{N, CV_TickLabel{LocT}}
app         CV_TickLabelAppearance{tsT, lsT}
```
"""
struct CV_Ruler{N, LocT, tsT, lsT}    # {{{
    ticklabels :: NTuple{N, CV_TickLabel{LocT}}
    app        :: CV_TickLabelAppearance{tsT, lsT}
end

show(io::IO, ruler::CV_Ruler) = cv_show_impl(io, ruler)
show(io::IO, m::MIME{Symbol("text/plain")}, ruler::CV_Ruler) =
    cv_show_impl(io, m, ruler)

function CV_Ruler(ticklabels::NTuple{N, CV_TickLabel{LocT}}) where {N, LocT}
    return CV_Ruler(ticklabels, CV_TickLabelAppearance())
end

"""
```
↦(app, ticklabels)
    app          CV_TickLabelAppearance,
    ticklabels   NTuple{N, CV_TickLabel{LocT}}  where {N, LocT}
```

short for `CV_Ruler(ticklabels, app)`.
"""
function ↦(app::CV_TickLabelAppearance,
        ticklabels::NTuple{N, CV_TickLabel{LocT}})  where {N, LocT}
    return CV_Ruler(ticklabels, app)
end

"""
```
CV_Ruler(ticklabels)

ticklabels   Vararg{CV_TickLabel{LocT}, N}
```
"""
function CV_Ruler(ticklabels::Vararg{CV_TickLabel{LocT}, N}) where {N, LocT}
    return CV_Ruler(NTuple{N, CV_TickLabel{LocT}}(ticklabels))
end

# }}}

"""
```
cv_format_ticks(printf_format, locations)

printf_format   AbstractString     used for @sprintf
locations       Vararg{Real, N}
```

create for every location a `CV_TickLabel` by formatting the locations
with the given printf-format.
"""
function cv_format_ticks(printf_format::AbstractString,
        locations::Vararg{Real, N}) where {N}
    fmt_func = isempty(printf_format) ?
        x -> "" : @eval x -> @sprintf($printf_format, x)
    return NTuple{N, CV_TickLabel{Float64}}(
        map(x -> CV_TickLabel{Float64}(
                Float64(x), Base.invokelatest(fmt_func, Float64(x))),
            locations))
end

"""
```
⇒(format, locations)
  format     AbstractString
  locations  AbstractVector
```

short for `cv_format_ticks(format, locations...)`.
"""
⇒(format::AbstractString, locations::AbstractVector) =
    cv_format_ticks(format, locations...)

"""
```
cv_format_ticks(locations)

locations   Vararg{Real, N}
```

create for every location a `CV_TickLabel` by formatting the locations
with the `"%.1f"` printf-format.
"""
function cv_format_ticks(locations::Vararg{Real, N}) where {N}
    return cv_format_ticks("%.1f", locations...)
end

"""
struct used internally for constructing the ticks and labels.

Needed in advance (before allocating the axis canvas) in order to
be able to determine the size needed for the axis canvas.

```
                     ┌─────── cv_north ────────┐
                     │ label                   │
                     │                         │               ┌─────────┐
                     │   │                     │               │         │
              ┌──────┴─*─┴─────────────────────┴───────────────*         c
              │        │                                       │         v
              │        │                                       │         _
              c        │                                       ├── lab   e
              v lab1 ──┤                                       │         a
              _        │                                       │         s
              w        │                                       │         t
              e        │      math canvas                      │         │
              s        │                                       │         │
              t lab2 ──┤                                       ├─────────┘
              │        │                                       │
              │        │                                       │
              │        │                                       │
              └────────┤                                       │
                       │                                       │
                       │                                       │
                       │                                       │
                  ┌────*──┬───────────────┬───────────────┬────┘
                  │       │tick_length    │               │
                  │       │               │               │
                  │         gap                           │
                  │                                       │
                  │  tick label         label2            │
                  └─────────────── cv_south ──────────────┘

         *   is (0,0) for axis canvas units (and :default-anchor from outside)
```
"""
mutable struct CV_TickLabelData{LocT}
    text           :: AbstractString
    text_extents   :: CV_TextExtents
    math_location  :: LocT
    tick_location  :: Int32
    startpoint     :: Tuple{Int32, Int32}
    bounding_box   :: CV_Rectangle{Int32}
end

"""
stores the maximum label sizes (height, depth, width, etc.)
"""
struct CV_TicksLabelsMetric
    max_height :: Int32
    max_depth  :: Int32
    max_width  :: Int32
end

function cv_ticks_labels_draw(con::CV_2DCanvasContext,
        attach::CV_AttachType, app::CV_TickLabelAppearance,
        ticklabelsdata::NTuple{N, CV_TickLabelData{Float64}};
        debug::Union{Val{true}, Val{false}}=Val(false)) where {N}  # {{{

    tick_length, tick_style = app.tick_length, app.tick_style
    label_style = app.label_style
    canvas = con.canvas

    ctx = con.ctx
    ubox = canvas.user_box
    # ticks
    cv_prepare(con, tick_style)
    for td in ticklabelsdata
        if attach isa CV_southT
            move_to(ctx, td.tick_location, 0)
            rel_line_to(ctx, 0, tick_length)
        elseif attach isa CV_northT
            move_to(ctx, td.tick_location, 0)
            rel_line_to(ctx, 0, -tick_length)
        elseif attach isa CV_eastT
            move_to(ctx, 0, td.tick_location)
            rel_line_to(ctx, tick_length, 0)
        elseif attach isa CV_westT
            move_to(ctx, 0, td.tick_location)
            rel_line_to(ctx, -tick_length, 0)
        else
            cv_error("Unknown attach value")
        end
        stroke(ctx)
    end

    # labels
    cv_prepare(con, label_style)
    for td in ticklabelsdata
        if !td.bounding_box.empty
            move_to(ctx, td.startpoint...)
            show_text(ctx, td.text)
        end
    end

    # bb check
    if debug isa Val{true}
        set_line_width(ctx, 1)
        set_source_rgb(ctx, 1, 0, 0)
        for td in ticklabelsdata
            bb = td.bounding_box
            rectangle(ctx, bb.left, bb.bottom, cv_width(bb), cv_height(bb))
            stroke(ctx)
        end

        set_source_rgb(ctx, 0, 1, 0)
        rectangle(ctx, ubox.left, ubox.bottom, cv_width(ubox), cv_height(ubox))
        stroke(ctx)
    end
end # }}}

# Compute positions for ticks and labels and their bounding boxes, etc. {{{

function cv_init_ticklables_data(
        label_style::CV_ContextStyle,
        ticklabels::Vararg{CV_TickLabel{Float64}, N}) where {N}  # {{{
    temp_can = CV_Std2DCanvas(5, 5)
    con = cv_create_context(temp_can)
    use_mem = Vector{Float64}(undef, 6)
    try
        cv_prepare(con, label_style)
        return NTuple{N, CV_TickLabelData}(map(ticklabels) do ticklabel
            CV_TickLabelData(
                ticklabel.text,
                cv_get_text_extents(con.ctx, ticklabel.text; use_mem=use_mem),
                ticklabel.location,
                Int32(-1), (Int32(-1), Int32(-1)),
                CV_Rectangle(Int32))
        end)
    finally
        cv_destroy(con)
        cv_destroy(temp_can)
    end
end # }}}

function cv_ticks_labels_get_metric(ticklabelsdata) # {{{
    return CV_TicksLabelsMetric(
        ceil(Int32,  # max_height
            maximum(map(td -> td.text_extents.height, ticklabelsdata))),
        ceil(Int32,  # max_depth
            maximum(map(td -> td.text_extents.depth, ticklabelsdata))),
        ceil(Int32,  # max_width
            maximum(map(td -> td.text_extents.bb_width, ticklabelsdata))))
end # }}}

function cv_ticks_labels_base(tick_length, gap, tlm, attach) # {{{
    return (
        attach isa CV_southT ? tick_length + gap + tlm.max_height  :
        attach isa CV_northT ? -tick_length - gap - tlm.max_depth  :
        attach isa CV_eastT  ? tick_length + gap                   :
        attach isa CV_westT  ? -tick_length - gap                  :
        cv_error("Unknown attach value")) :: Int32
end # }}}

function cv_set_tick_location!(ticklabelsdata, for_canvas, attach) # {{{
    for td in ticklabelsdata
        td.tick_location = (
            ((attach isa CV_northT) || (attach isa CV_southT))    ?
            cv_math2pixel(for_canvas, td.math_location, 0.0)[1]   :
            cv_math2pixel(for_canvas, 0.0, td.math_location)[2])
    end
    return nothing
end # }}}

function cv_get_label_xpos(td, base, attach::Union{CV_southT, CV_northT})
    return td.tick_location - round(Int32, cv_half(td.text_extents.bb_width))
end

function cv_get_label_xpos(td, base, attach::CV_eastT)
    return base
end

function cv_get_label_xpos(td, base, attach::CV_westT)
    return base - round(Int32, td.text_extents.bb_width)
end

function cv_get_label_ypos(td, base, attach::Union{CV_southT, CV_northT})
    return base
end

function cv_get_label_ypos(td, base, attach::Union{CV_eastT, CV_westT})
    return td.tick_location + round(Int32,
        td.text_extents.depth + td.text_extents.height/3)
end

function cv_get_label_bb(td, xpos, ypos)
    return CV_Rectangle(
        ceil(Int32, ypos + td.text_extents.depth),
        xpos,
        floor(Int32, ypos - td.text_extents.height),
        ceil(Int32, xpos + td.text_extents.bb_width))
end

function cv_add_rect_for_ticks(rstore, ticklabelsdata, tick_length,
        gap, attach) # {{{
    max_loc = ceil(Int32, maximum(map(td -> td.tick_location, ticklabelsdata)))
    min_loc = floor(Int32, minimum(map(td -> td.tick_location, ticklabelsdata)))
    z, o = zero(Int32), one(Int32)
    if min_loc == max_loc
        max_loc += o
    end
    if tick_length == 0
        tick_length = o
    end
    cv_add_rectangle!(rstore, (
        attach isa CV_southT ? CV_Rectangle(tick_length, min_loc, z, max_loc)  :
        attach isa CV_northT ? CV_Rectangle(z, min_loc, -tick_length, max_loc) :
        attach isa CV_eastT  ? CV_Rectangle(max_loc, z, min_loc, tick_length)  :
        attach isa CV_westT  ? CV_Rectangle(max_loc, -tick_length, min_loc, z) :
        cv_error("Unknown attach value")) :: CV_Rectangle{Int32})
    return nothing
end # }}}

function cv_add_rect_for_default_anchor(rstore, attach) # {{{
    z, o = zero(Int32), one(Int32)
    cv_add_rectangle!(rstore, (
        attach isa CV_southT ? CV_Rectangle(o, z, z, o)      :
        attach isa CV_northT ? CV_Rectangle(z, z, -o, o)     :
        attach isa CV_eastT  ? CV_Rectangle(o, z, z, o)      :
        attach isa CV_westT  ? CV_Rectangle(o, -o, z, z)     :
        cv_error("Unknown attach value")) :: CV_Rectangle{Int32})
    return nothing
end # }}}
# }}}

function cv_get_default_anchor(rstore, attach)  # {{{
    z, o = zero(Int32), one(Int32)
    bb = rstore.bounding_box
    return (
        attach isa CV_southT ? (-bb.left, z)              :
        attach isa CV_northT ? (-bb.left, cv_height(bb))  :
        attach isa CV_eastT  ? (z, -bb.bottom)            :
        attach isa CV_westT  ? (cv_width(bb), -bb.bottom) :
        cv_error("Unknown attach value")) :: Tuple{Int32, Int32}

end # }}}

function cv_layout_ruler(for_canvas::CV_Math2DCanvas, attach::CV_AttachType,
        ruler::CV_Ruler, rstore::CV_RectangleStore{Int32})  # {{{
    app = ruler.app
    tick_length, tick_style = app.tick_length, app.tick_style
    gap, label_style = app.gap, app.label_style

    ticklabelsdata = cv_init_ticklables_data(label_style, ruler.ticklabels...)
    tlm = cv_ticks_labels_get_metric(ticklabelsdata)

    base = cv_ticks_labels_base(tick_length, gap, tlm, attach)
    cv_set_tick_location!(ticklabelsdata, for_canvas, attach)

    for td in ticklabelsdata
        xpos = cv_get_label_xpos(td, base, attach)
        ypos = cv_get_label_ypos(td, base, attach)
        td.bounding_box = cv_get_label_bb(td, xpos, ypos)
        td.startpoint = (xpos - round(Int32, td.text_extents.x_bearing), ypos)
        !td.bounding_box.empty && cv_add_rectangle!(rstore, td.bounding_box)
    end
    cv_add_rect_for_ticks(rstore, ticklabelsdata, tick_length, gap, attach)
    cv_add_rect_for_default_anchor(rstore, attach)

    return ticklabelsdata
end # }}}

"""
```
cv_create_2daxis_canvas(for_canvas, attach, rulers)

for_canvas   CV_Math2DCanvas
attach       Union{cv_north, cv_south, cv_east, cv_west}
rulers       Tuple{Vararg{CV_Ruler, N}}
```

create a `CV_2DLayoutCanvas` to render the given rulers. Position the
ticks and labels such that the axis-canvas can be attached to `for_canvas`
at the given `attach` position.
"""
function cv_create_2daxis_canvas(for_canvas::CV_Math2DCanvas,
        attach::CV_AttachType,
        rulers::NTuple{N, CV_Ruler}) where {N} # {{{

    N < 1 && cv_error("Need at least one ruler")
    rstore = CV_RectangleStore(Int32)     # all bounding boxes of the labels

    ticklabelsdata = tuple(
        map(ruler -> cv_layout_ruler(
            for_canvas, attach, ruler, rstore), rulers)...)

    default_anchor = cv_get_default_anchor(rstore, attach)
    anchor_func = (can, anchor_name) -> (anchor_name === :default ?
            default_anchor : cv_anchor(can.bounding_box, anchor_name))
    canvas = CV_2DLayoutCanvas(rstore.bounding_box, anchor_func)

    con = cv_create_context(canvas; fill_with=cv_color(0,0,0,0))
    for index in 1:N
        cv_ticks_labels_draw(con,
            attach, rulers[index].app, ticklabelsdata[index])
    end
    cv_destroy(con)

    return canvas
end # }}}

"""
```
cv_create_2daxis_canvas(for_canvas, attach, ticklabels; app)

for_canvas   CV_Math2DCanvas
attach       Union{cv_north, cv_south, cv_east, cv_west}
ticklabels   Vararg{CV_TickLabel{Float64}, N}
app          CV_TickLabelAppearance
```
"""
function cv_create_2daxis_canvas(for_canvas::CV_Math2DCanvas,
        attach::CV_AttachType,
        ticklabels::Vararg{CV_TickLabel{Float64}, N};
        app::CV_TickLabelAppearance=CV_TickLabelAppearance()) where {N}
    return cv_create_2daxis_canvas(for_canvas,
        attach, (CV_Ruler(ticklabels, app),))
end

"""
```
cv_ticks_labels(layout, for_canvas_l, attach, rulers)

layout         CV_Abstract2DLayout
for_canvas_l   CV_2DLayoutPosition
attach         Union{cv_north, cv_south, cv_east, cv_west}
rulers         Tuple{Vararg{ComplexVisual.CV_Ruler, N}}
```

short version for `cv_create_2daxis_canvas` and `cv_add_canvas!`.
"""
function cv_ticks_labels(layout::CV_Abstract2DLayout,
        for_canvas_l::CV_2DLayoutPosition{CV_Math2DCanvas, dcbT, styleT},
        attach::CV_AttachType,
        rulers::NTuple{N, CV_Ruler}) where {dcbT, styleT, N}
    canvas = cv_create_2daxis_canvas(for_canvas_l.canvas, attach, rulers)

    return cv_add_canvas!(layout, canvas, cv_anchor(canvas, :default),
        cv_anchor(for_canvas_l, (
            attach isa CV_southT ? :southwest   :
            attach isa CV_northT ? :northwest   :
            attach isa CV_eastT  ? :northeast   : :northwest)))
end

"""
```
cv_ticks_labels(layout, for_canvas_l, attach, ticklabels; app)

layout         CV_Abstract2DLayout
for_canvas_l   CV_2DLayoutPosition
attach         Union{cv_north, cv_south,
                     cv_east, cv_west}
ticklabels     Vararg{CV_TickLabel{Float64}, N}
app            CV_TickLabelAppearance              CV_TickLabelAppearance()
```
"""
function cv_ticks_labels(layout::CV_Abstract2DLayout,
        for_canvas_l::CV_2DLayoutPosition, attach::CV_AttachType,
        ticklabels::Vararg{CV_TickLabel{Float64}, N};
        app::CV_TickLabelAppearance=CV_TickLabelAppearance()) where {N}

    return cv_ticks_labels(layout, for_canvas_l,
        attach, (CV_Ruler(ticklabels, app), ))
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using Documenter
#include("../src/QSWalk.jl")
using QSWalk

# same for contributing and license
cp(normpath(@__FILE__, "../../LICENSE"), normpath(@__FILE__, "../src/license.md"); force=true)

makedocs(
    modules     = [QSWalk],
    format      = Documenter.HTML(assets = ["assets/logo.ico"]),
    sitename    = "QSWalk",
    clean       = true,
    doctest     = true,
    checkdocs   = :exports,
    pages       = Any[
		"Home"				=> "index.md",
        "GKSL master equation" => "gksl.md",
    	"Demoralization"   => "demoralization.md",
        "Contributing"		=> "contributing.md",
        "Citing"	       	=> "citing.md",
		"Licence"			=> "license.md",
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/iitis/QSWalk.jl",
    target      = "build",
)

#rm(normpath(@__FILE__, "src/license.md"))

ext: jl


text:
"""
Transform x::Real to be in a [0.1] interval.
"""
@inline transform_unit_interval(x::Real) = 1.0 / (exp(-x) + 1.0)
@inline transform_from_unit_interval(x::Real) = logit(x)
@inline invlogit(x::Real) = 1.0 / (exp(-x) + 1.0)
@inline transform_unit_interval(x::Real, a::Real, b::Real) = (b - a) * invlogit(x) + a
@inline transform_from_unit_interval(x::Real, a::Real, b::Real) = logit((x - a) / (b - a))

const LOWER_ALPHA = 1e-04
const UPPER_ALPHA = 1 - 1e-04

const LOWER_BETA = 1e-04
const UPPER_BETA = 1 - 1e-04

abstract type AbstractExponentitalSmoothing end

struct FOES{VX, IDX, MX, R} <: AbstractExponentitalSmoothing
    parms::VX
    residuals::VX
    states::MX
    values::VX
    index::IDX
    sigma::R
end

struct Holt{VX, IDX, MX, R} <: AbstractExponentitalSmoothing
    parms::VX
    residuals::VX
    states::MX
    values::VX
    index::IDX
    sigma::R
end

struct HoltDumped{VX, IDX, MX, R} <: AbstractExponentitalSmoothing
    parms::VX
    residuals::VX
    states::MX
    values::VX
    index::IDX
    sigma::R
end

function FOES(y::T) where T<:AbstractVector
    ly = length(y)
    m = FOES(Array{Float64, 1}(undef, 2),
             Array{Float64, 1}(undef, ly),
             Array{Float64, 2}(undef, ly+1, 1),
             y,
             1:length(y),
             Ref{Float64}(NaN))
    startx!(m)
    m
end

function FOES(y::T) where T<:AbstractTimeSeries
    ly = length(y)
    m = FOES(Array{Float64, 1}(undef, 2),
             Array{Float64, 1}(undef, ly),
             Array{Float64, 2}(undef, ly+1, 1),
             float(values(y)),
             timestamp(y),
             Ref{Float64}(NaN))
    startx!(m)
    m
end

function Holt(y::T) where T<:AbstractVector
    ly = length(y)
    m = Holt(Array{Float64, 1}(undef, 4),
             Array{Float64, 1}(undef, ly),
             Array{Float64, 2}(undef, ly+1, 2),
             y,
             1:length(y),
             Ref{Float64}(NaN))
    startx!(m)
    m
end

function Holt(y::T) where T<:AbstractTimeSeries
    ly = length(y)
    m = Holt(Array{Float64, 1}(undef, 4),
             Array{Float64, 1}(undef, ly),
             Array{Float64, 2}(undef, ly+1, 2),
             float(values(y)),
             timestamp(y),
             Ref{Float64}(NaN))
    startx!(m)
    m
end

next(idx::StepRange{D, P}, h) where {D<:Date, P<:DatePeriod} = next(collect(idx), h)

function next(idx::Array{Date, 1}, h)
    ## Try to infer the frequency of the series
    mts = maximum(diff(idx))
    if mts >= Day(365)
        last(idx) + Year(1):Year(1):last(idx) + Year(h)
    elseif mts >= Day(1) & mts <= Day(4)
        last(idx) + Day(1):Day(1):last(idx) + Day(h)
    elseif mts >= Day(90)
        last(idx) + Month(3):Month(3):last(idx) + Month(3*h)
    end
end

next(idx::UnitRange{P}, h) where P<:Int64 = last(idx) + 1:last(idx) + h
next(idx::AbstractVector{P}, h) where P<:Int64 = last(idx) + 1:last(idx) + h

struct Forecast{MOD, A, M, V, L}
    method::MOD
    mean::A
    margin::M
    sigma::V
    level::L
end

function trendcoef(y)
    T = length(y)
    m = mean(y)
    tbar = (1 + T)/2
    s = 0.0
    @simd for j in eachindex(y)
        @inbounds s += y[j]*j
    end
    b1 = (s - T*tbar*m)/(T*(T^2-1)/12)
    (m - tbar*b1, b1)
end

#---
# FOES - First Order Exponential Smoothing
#---

objective(m::FOES) = begin
    (y, θ) -> begin
        ly = length(y)
        a  = transform_unit_interval(first(θ), 1e-04, 0.9999)
        s₀ = last(θ)
        s = a*y[1] + (1-a)*s₀
        e = abs2(s - y[2]) + abs2(y[1] - s₀)
        @inbounds for t in 2:ly-1
            s = a*y[t] + (1-a)*s
            e += abs2(s - y[t+1])
        end
        e
    end
end

function startx!(m::FOES)
    m.parms[1] = transform_from_unit_interval(.5, 1e-04, 0.9999)
    m.parms[2] = mean(m.values)
end

function recursion!(m::FOES)
    ly = length(m.values)
    a, s₀ = m.parms
    m.states[1] = s₀
    m.residuals[1] = m.values[1] - s₀
    for t in 1:ly-1
        m.states[t+1] = a*m.values[t] + (1-a)*m.states[t]
        m.residuals[t+1] = m.values[t+1] - m.states[t+1]
    end
    m.states[ly+1] = a*m.values[ly] + (1-a)*m.states[ly]
end

function var_prediction(m::FOES, h)
    σ² = m.sigma.x
    a  = m.parms[1]
    (σ²*(1 + a^2*(j-1)) for j in 1:h)
end

function mean_forecast(m::FOES, h)
    meanf = Array{Float64,1}(undef, h)
    fill!(meanf, last(m.states))
end

# ---
# HOLT - Second Order Exponential Smoothing
# ---

objective(m::Holt) = begin
    (y, θ) -> begin
        ## θ = [α, β, ℓ₀, s₀]
        ly = length(y)
        a  = transform_unit_interval(θ[1], 1e-04, 0.9999)
        b  = transform_unit_interval(θ[2], 1e-04, 0.9999)
        ℓ₋₁ = θ[3]
        s₋₁ = θ[4]
        e = abs2(y[1] - ℓ₋₁ - s₋₁)
        @inbounds for t in 1:ly-1
            ℓ = a*y[t] + (1-a)*(ℓ₋₁ + s₋₁)
            s = b*(ℓ - ℓ₋₁) + (1-b)*s₋₁
            ℓ₋₁ = ℓ
            s₋₁ = s
            e += abs2(y[t+1] - ℓ - s)
        end
        e
    end
end

function recursion!(m::Holt)
    ly = length(m.values)
    a, b, ℓ₀, s₀ = m.parms
    m.states[1, 1] = ℓ₀
    m.states[1, 2] = s₀
    m.residuals[1] = m.values[1] - ℓ₀ - s₀
    for t in 1:ly-1
        m.states[t+1, 1] = a*m.values[t] + (1-a)*(m.states[t, 1] + m.states[t, 2])
        m.states[t+1, 2] = b*(m.states[t+1, 1] - m.states[t, 1]) + (1-b)*m.states[t, 2]
        m.residuals[t+1] = m.values[t+1] - m.states[t+1, 1] - m.states[t+1, 2]
    end
    m.states[ly+1, 1] = a*m.values[ly] + (1-a)*(m.states[ly, 1] + m.states[ly, 2])
    m.states[ly+1, 2] = b*(m.states[ly+1, 1] - m.states[ly, 1]) + (1-b)*m.states[ly, 2]
end

function startx!(m::Holt)
    m.parms[1] = transform_from_unit_interval(.5, 1e-04, 0.9999)
    m.parms[2] = transform_from_unit_interval(.5, 1e-04, 0.9999)
    b0, b1 = trendcoef(m.values)
    m.parms[3] = b0
    m.parms[4] = b1
end

function var_prediction(m::Holt, h)
    σ² = m.sigma.x
    a, b, l₀, s₀ = m.parms
    (σ²*(1 + (j-1)*(a^2 + a*b*j + (1/6)*(b^2)*j*(2*j -1))) for j in 1:h)
end

function mean_forecast(m::Holt, h)
    meanx = Array{Float64,1}(undef, h)
    ℓ₋₁ = m.states[end,1]
    s₋₁ = m.states[end,2]
    a, b, l₀, s₀ = m.parms
    for j in 1:h
        meanx[j] = ℓ₋₁ + j*s₋₁
    end
    meanx
end

#=
HOLT - (Dumped) Second Order Exponential Smoothing
=#

holt_dumped_obj(y, θ) = begin
    ## θ = [α, β, ℓ₀, s₀, ϕ]
    a  = transform_unit_interval(θ[1])
    b  = transform_unit_interval(θ[2])
    ℓ₀ = θ[3]
    s₀ = θ[4]
    ℓ = a*y[1] + (1-a)*(ℓ₀ + s₀)
    s = b*(ℓ - ℓ₀) + (1-b)*s₀
    ℓ₋₁ = ℓ
    s₋₁ = s
    e = abs2(y[1] - ℓ₀ - ϕ*s₀) + abs2(y[2] - ℓ - ϕ*s)
    @inbounds for t in 2:ly-1
        ℓ = a*y[t] + (1-a)*(ℓ₋₁ + ϕ*s₋₁)
        s = b*(ℓ - ℓ₋₁) + (1-b)*ϕ*s₋₁
        ℓ₋₁ = ℓ
        s₋₁ = s
        e += abs2(y[t+1] - ℓ - s)
    end
    e
end

function prediction_intervals(m::T, h::Int, level::AbstractArray) where T<:AbstractExponentitalSmoothing
    v = var_prediction(m, h)
    meanf  = mean_forecast(m, h)
    margin = (norminvcdf(1-(1-j)/2).*sqrt.(v) for j in level)
    return meanf, v, margin, level
end

function Forecast(m::T; h = 12, level = [.95]) where T<:AbstractExponentitalSmoothing
    meanf, v, margin, level = prediction_intervals(m, h, level)
    Forecast(m, meanf, margin, v, level)
end

function foes(y::T) where T<:AbstractVector
    m = FOES(y)
    obj = objective(m)
    res = optimize(θ -> obj(m.values, θ), m.parms, BFGS(), autodiff = :forward)
    m.parms[1] = transform_unit_interval(Optim.minimizer(res)[1], 1e-04, 0.9999)
    m.parms[2] = Optim.minimizer(res)[2]
    recursion!(m)
    m.sigma.x = var(m.residuals)
    m
end

function foes(y::T) where T<:AbstractTimeSeries
    m = FOES(y)
    obj = objective(m)
    res = optimize(θ -> obj(m.values, θ), m.parms, BFGS(), autodiff = :forward)
    m.parms[1] = transform_unit_interval(Optim.minimizer(res)[1], 1e-04, 0.9999)
    m.parms[2] = Optim.minimizer(res)[2]
    recursion!(m)
    m.sigma.x = var(m.residuals)
    m
end

function holt(y::T) where T<:AbstractVector
    m = Holt(y)
    obj = objective(m)
    #res = optimize(θ -> obj(m.values, θ), m.parms, BFGS(), autodiff = :forward)
    res = optimize(θ -> obj(m.values, θ), m.parms, NelderMead())
    m.parms[1] = transform_unit_interval(Optim.minimizer(res)[1], 1e-04, 0.9999)
    m.parms[2] = transform_unit_interval(Optim.minimizer(res)[2], 1e-04, 0.9999)
    m.parms[3] = Optim.minimizer(res)[3]
    m.parms[4] = Optim.minimizer(res)[4]
    recursion!(m)
    m.sigma.x = var(m.residuals)
    m
end

function holt(y::T) where T<:AbstractTimeSeries
    m = Holt(y)
    obj = objective(m)
    res = optimize(θ -> obj(m.values, θ), m.parms, BFGS(), autodiff = :forward)
    m.parms[1] = transform_unit_interval(Optim.minimizer(res)[1], 1e-04, 0.9999)
    m.parms[2] = transform_unit_interval(Optim.minimizer(res)[2], 1e-04, 0.9999)
    m.parms[3] = Optim.minimizer(res)[3]
    m.parms[4] = Optim.minimizer(res)[4]
    recursion!(m)
    m.sigma.x = var(m.residuals)
    m
end

@recipe function fanchart(f::T; forecast_linecolor = :white, colorscheme = ColorSchemes.Blues_3, label = "") where T<:Forecast
    seriestype  --> :path
    insmpl = f.method.index
    outsmpl = next(f.method.index, length(f.mean))
    full_idx = [f.method.index; next(f.method.index, length(f.mean))]
    c1 = colorant"#005073"
    c2 = colorant"#71c7ec"
    col_margin = range(c1, stop=c2, length=length(f.margin))
    @series begin
        linecolor --> :red
        seriestype := :path
        insmpl, f.method.values
    end

    for (i, mar) in enumerate(Iterators.reverse(f.margin))
        @series begin
            linecolor := get(colorscheme, maximum(f.level))
            fillcolor := get(colorscheme, f.level[i])
            ribbon := (mar, mar)
            outsmpl, f.mean
        end
    end

    x = length(insmpl) + 1; y = 1
    h = Inf; w = length(full_idx) - x
    @series begin
        seriestype := :shape
        fillcolor --> :lightgray
        fillalpha := 0.4
        linewidth := 0.0
        x .+ [0, w, w, 0], y .+ [-Inf, -Inf, h, h]
    end

    ## Plot the anchor
    @series begin
        seriestype := :path
        linecolor --> :red
        linestyle := :dot
        [last(insmpl); first(outsmpl)], [last(f.method.values); first(f.mean)]
    end

    ## Plot

end


text:
function wasmer_last_error_length()
    @ccall libwasmer.wasmer_last_error_length()::Cint
end

function wasmer_last_error_message(error_length = get_last_error_length())
    error_length === 0 && return

    buffer = Vector{UInt8}(undef, error_length)
    buffer_ptr = Base.unsafe_convert(Ptr{UInt8}, buffer)
    res = @ccall libwasmer.wasmer_last_error_message(
        buffer_ptr::Cstring,
        error_length::Cint,
    )::Cint
    res === -1 && error("Failed to retrieve last wasmer error")

    unsafe_string(buffer_ptr, error_length)
end

function check_wasmer_error()
    error_length = wasmer_last_error_length()
    if error_length != 0
        error_msg = wasmer_last_error_message(error_length)
        error(error_msg)
    end
end

wat2wasm(str::AbstractString) = wat2wasm(WasmByteVec(collect(wasm_byte_t, str)))
function wat2wasm(wat::WasmByteVec)
    out = WasmByteVec()
    @ccall libwasmer.wat2wasm(wat::Ptr{wasm_byte_vec_t}, out::Ptr{wasm_byte_vec_t})::Cvoid
    check_wasmer_error()

    out
end

macro wat_str(wat::String)
    :(wat2wasm($wat))
end


text:
using Printf

include("./ParticleSwarmAlignment.jl")

function main()
    #GC.enable(false)
    if length(ARGS) > 0
        if length(ARGS) != 3
            println("Usage:")
            println("  julia main.jl number_of_sequences length_of_sequences swarm_iterations")
            throw(ErrorException("Command line args were provided but there were not 3 of them"))
        end
        t = parse(Int64, ARGS[1])
        N = parse(Int64, ARGS[2])
        iterations = parse(Int64, ARGS[3])
        @printf("Using provided args t=%d, N=%d, iterations=%d\n", t, N, iterations)
    else
        # number of sequences
        t = 5
        # length of sequences
        N = 38
        # iterations
        iterations = 200
        @printf("Using default args t=%d, N=%d, iterations=%d\n", t, N, iterations)
    end

    sequences = ParticleSwarmAlignment.generate_sequences(t, N)
    println(sequences)
    val, t, bytes, gctime, memallocs = @timed ParticleSwarmAlignment.PSO_MSA(sequences, iterations)

    @printf("Total bytes: %d\n", bytes)
end

main()


text:
using DimensionalData, Statistics, Test, Unitful, SparseArrays, Dates

using LinearAlgebra: Transpose

using DimensionalData: _Rot90, _Rot180, _Rot270, _Rot360, _rottype

@testset "map" begin
    a = [1 2; 3 4]
    dimz = X(143:2:145), Y(-38:2:-36; mode=Sampled(span=Explicit([-38 -36; -36 -34])))
    da = DimArray(a, dimz)
    @test map(x -> 2x, da) == [2 4; 6 8]
    @test map(x -> 2x, da) isa DimArray{Int64,2}
    @test map(*, da, da) == [1 4; 9 16]
    @test map(*, da, da) isa DimArray{Int64,2}
end

@testset "dimension reducing methods" begin

    # Test some reducing methods with Explicit spans
    a = [1 2; 3 4]
    dimz = X((143, 145); mode=Sampled()), Y(LinRange(-38, -36, 2); mode=Sampled(span=Explicit([-38 -36; -36 -34])))
    da = DimArray(a, dimz)

    # Test all dime combinations with maxium
    @test maximum(da; dims=2) == [2 4]'
    @test maximum(da; dims=Y()) == [2 4]'
    @test maximum(da; dims=Y) == [2 4]'
    @test maximum(da; dims=:Y) == [2 4]'
    @test maximum(da; dims=(1, 2)) == [4]'
    @test maximum(da; dims=(X, Y)) == [4]'
    @test maximum(da; dims=(X(), Y())) == [4]'
    @test maximum(da; dims=(:X, :Y)) == [4]'
    @test maximum(x -> 2x, da; dims=X) == [6 8]
    @test maximum(x -> 2x, da; dims=2) == [4 8]'
    @test maximum(x -> 2x, da; dims=(X, Y)) == [8]'
    @test maximum(x -> 2x, da; dims=(:X, :Y)) == [8]'
    @test maximum(x -> 2x, da; dims=(X(), Y())) == [8]'

    @test minimum(da; dims=1) == [1 2]
    @test minimum(da; dims=:) == 1
    @test minimum(da; dims=Y()) == [1 3]'
    ds = (X([144.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()),
                Y(LinRange(-38, -36, 2), Sampled(Ordered(), Explicit([-38 -36; -36 -34]), Intervals(Center())), NoMetadata()))
    @test typeof(dims(minimum(da; dims=X()))) == typeof(ds)
    @test val.(span(minimum(da; dims=X()))) == val.(span(ds))

    @test sum(da; dims=X()) == sum(a; dims=1)
    ds = (X(LinRange(143.0, 145.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
         Y([-37.0], Sampled(Ordered(), Explicit(reshape([-38, -34], 2, 1)), Intervals(Center())), NoMetadata()))
    @test typeof(dims(sum(da; dims=Y()))) == typeof(ds)
    @test index(sum(da; dims=Y())) == index.(ds)
    @test val.(span(sum(da; dims=Y()))) == val.(span(ds))
    @test sum(da; dims=:) == 10
    @test sum(x -> 2x, da; dims=:) == 20

    a = [1 2; 3 4]
    dimz = X((143, 145); mode=Sampled()), Y((-38, -36); mode=Sampled())
    da = DimArray(a, dimz)

    @test prod(da; dims=X) == [3 8]
    @test prod(da; dims=2) == [2 12]'
    resultdimz =
        (X([144.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()),
         Y(LinRange(-38.0, -36.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()))
    @test typeof(dims(prod(da; dims=X()))) == typeof(resultdimz)
    @test bounds(dims(prod(da; dims=X()))) == bounds(resultdimz)

    @test mean(da; dims=1) == [2.0 3.0]
    @test mean(da; dims=Y()) == [1.5 3.5]'
    @test mean(da; dims=(1, 2)) == [2.5]'
    @test mean(da; dims=(X, Y)) == [2.5]'
    @test mean(x -> 2x, da; dims=1) == [4.0 6.0]
    @test mean(x -> 2x, da; dims=Y) == [3.0 7.0]'
    @test mean(x -> 2x, da; dims=(1, 2)) == [5.0]'
    @test mean(x -> 2x, da; dims=(:X, :Y)) == [5.0]'
    @test dims(mean(da; dims=Y())) ==
        (X(LinRange(143.0, 145.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
         Y([-37.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()))

    @test mapreduce(x -> x > 3, +, da; dims=X) == [0 1]
    @test mapreduce(x -> x > 3, +, da; dims=(X, Y)) == [1]'
    @test mapreduce(x -> x > 3, +, da; dims=(:X, :Y)) == [1]'
    @test mapreduce(x -> x > 3, +, da; dims=(1, 2)) == [1]'
    @test mapreduce(x -> x > 3, +, da; dims=:) == 1
    @test dims(mapreduce(x-> x > 3, +, da; dims=Y())) ==
        (X(LinRange(143.0, 145.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
         Y([-37.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()))

    @test reduce(+, da) == reduce(+, a)
    @test reduce(+, da; dims=X) == [4 6]
    @test reduce(+, da; dims=(X, Y())) == [10]'
    @test dims(reduce(+, da; dims=Y())) ==
        (X(LinRange(143.0, 145.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
         Y([-37.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()))

    @test std(da) === std(a)
    @test std(da; dims=1) == [1.4142135623730951 1.4142135623730951]
    @test std(da; dims=Y()) == [0.7071067811865476 0.7071067811865476]'

    @test var(da; dims=1) == [2.0 2.0]
    @test var(da; dims=Y()) == [0.5 0.5]'
    @test dims(var(da; dims=Y())) ==
        (X(LinRange(143.0, 145.0, 2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
         Y([-37.0], Sampled(Ordered(), Regular(4.0), Points()), NoMetadata()))

    @test extrema(da; dims=Y) == permutedims([(1, 2) (3, 4)])
    @test extrema(da; dims=X) == [(1, 3) (2, 4)]
    @test extrema(da; dims=(X, Y)) == reshape([(1, 4)], 1, 1)

    a = [1 2 3; 4 5 6]
    da = DimArray(a, dimz)
    @test median(da) == 3.5
    @test median(da; dims=X()) == [2.5 3.5 4.5]
    @test median(da; dims=2) == [2.0 5.0]'

    a = Bool[0 1 1; 0 0 0]
    da = DimArray(a, dimz)
    @test any(da) === true
    @test any(da; dims=Y) == reshape([true, false], 2, 1)
    @test all(da) === false
    @test all(da; dims=Y) == reshape([false, false], 2, 1)
    @test all(da; dims=(X, Y)) == reshape([false], 1, 1)

end

@testset "dimension dropping methods" begin
    a = [1 2 3; 4 5 6]
    dimz = X((143, 145); mode=Sampled()), Y((-38, -36); mode=Sampled())
    da = DimArray(a, dimz)
    # Dimensions must have length 1 to be dropped
    @test dropdims(da[X(1:1)]; dims=X) == [1, 2, 3]
    @test dropdims(da[2:2, 1:1]; dims=(X(), Y()))[] == 4
    @test typeof(dropdims(da[2:2, 1:1]; dims=(X(), Y()))) <: DimArray{Int,0,Tuple{}}
    @test refdims(dropdims(da[X(1:1)]; dims=X)) ==
        (X(143.0; mode=Sampled(Ordered(), Regular(2.0), Points())),)
    dropped = dropdims(da[X(1:1)]; dims=X)
    @test dropped[1:2] == [1, 2]
    @test length.(dims(dropped[1:2])) == size(dropped[1:2])
end

@testset "eachslice" begin
    a = [1 2 3 4
         3 4 5 6
         5 6 7 8]
    da = DimArray(a, (Y((10, 30)), Ti(1:4)))
    @test [mean(s) for s in eachslice(da; dims=Ti)] == [3.0, 4.0, 5.0, 6.0]
    @test [mean(s) for s in eachslice(da; dims=2)] == [3.0, 4.0, 5.0, 6.0]

    slices = [s .* 2 for s in eachslice(da; dims=Y)]
    @test slices[1] == [2, 4, 6, 8]
    @test slices[2] == [6, 8, 10, 12]
    @test slices[3] == [10, 12, 14, 16]
    dims(slices[1]) == (Ti(1.0:1.0:4.0),)

    slices = [s .* 2 for s in eachslice(da; dims=Ti)]
    @test slices[1] == [2, 6, 10]
    dims(slices[1]) == (Y(10.0:10.0:30.0),)
    @test_throws ArgumentError [s .* 2 for s in eachslice(da; dims=(Y, Ti))]
end

@testset "simple dimension permuting methods" begin
    da = DimArray(zeros(5, 4), (Y((10, 20); mode=Sampled()),
                                        X(1:4; mode=Sampled())))
    tda = transpose(da)
    @test tda == transpose(parent(da))
    resultdims = (X(1:4; mode=Sampled(Ordered(), Regular(1), Points())),
                  Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())))
    @test typeof(dims(tda)) == typeof(resultdims)
    @test dims(tda) == resultdims
    @test size(tda) == (4, 5)

    tda = Transpose(da)
    @test tda == Transpose(parent(da))
    @test dims(tda) == (X(1:4; mode=Sampled(Ordered(), Regular(1), Points())),
                        Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())))
    @test size(tda) == (4, 5)
    @test typeof(tda) <: DimArray

    ada = adjoint(da)
    @test ada == adjoint(parent(da))
    @test dims(ada) == (X(1:4; mode=Sampled(Ordered(), Regular(1), Points())),
                        Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())))
    @test size(ada) == (4, 5)

    dsp = permutedims(da)
    @test permutedims(parent(da)) == parent(dsp)
    @test dims(dsp) == reverse(dims(da))
end

@testset "dimension permuting methods with specified permutation" begin
    da = DimArray(ones(5, 2, 4), (Y((10, 20); mode=Sampled()),
                                          Ti(10:11; mode=Sampled()),
                                          X(1:4; mode=Sampled())))
    dsp = permutedims(da, [3, 1, 2])
    @test permutedims(da, [X, Y, Ti]) == permutedims(da, (X, Y, Ti))
    @test permutedims(da, [X(), Y(), Ti()]) == permutedims(da, (X(), Y(), Ti()))
    dsp = permutedims(da, (X(), Y(), Ti()))
    @test dsp == permutedims(parent(da), (3, 1, 2))
    @test dims(dsp) == (X(1:4; mode=Sampled(Ordered(), Regular(1), Points())),
                        Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())),
                        Ti(10:11; mode=Sampled(Ordered(), Regular(1), Points())))

    dsp = PermutedDimsArray(da, (3, 1, 2))
    @test dsp == PermutedDimsArray(parent(da), (3, 1, 2))
    @test typeof(dsp) <: DimArray
end

@testset "dimension rotating methods" begin
    @test _rottype(-100) == _Rot360()
    @test _rottype(-4) == _Rot360()
    @test _rottype(-3) == _Rot90()
    @test _rottype(-2) == _Rot180()
    @test _rottype(-1) == _Rot270()
    @test _rottype(0) == _Rot360()
    @test _rottype(1) == _Rot90()
    @test _rottype(2) == _Rot180()
    @test _rottype(3) == _Rot270()
    @test _rottype(4) == _Rot360()
    @test _rottype(101) == _Rot90()

    da = DimArray([1 2; 3 4], (X([:a, :b]), Y([1.0, 2.0])))

    l90 = rotl90(da)
    r90 = rotr90(da)
    r180_1 = rot180(da)
    r180_2 = rotl90(da, 2)
    r180_3 = rotr90(da, 2)
    r270 = rotl90(da, 3)
    r360 = rotr90(da, 4)
    da[X(At(:a)), Y(At(2.0))]
    @test l90[X(At(:a)), Y(At(2.0))] == 2
    @test r90[X(At(:a)), Y(At(2.0))] == 2
    @test r180_1[X(At(:a)), Y(At(2.0))] == 2
    @test r180_2[X(At(:a)), Y(At(2.0))] == 2
    @test r180_3[X(At(:a)), Y(At(2.0))] == 2
    @test r270[X(At(:a)), Y(At(2.0))] == 2
    @test r360[X(At(:a)), Y(At(2.0))] == 2
end

@testset "dimension mirroring methods" begin
    a = rand(5, 4)
    da = DimArray(a, (Y((10, 20); mode=Sampled()),
                              X(1:4; mode=Sampled())))
    cvda = cov(da; dims=X)
    @test cvda == cov(a; dims=2)
    @test dims(cvda) == (Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())),
                         Y(LinRange(10.0, 20.0, 5); mode=Sampled(Ordered(), Regular(2.5), Points())))
    crda = cor(da; dims=Y)
    @test crda == cor(a; dims=1)
    @test dims(crda) == (X(1:4; mode=Sampled(Ordered(), Regular(1), Points())),
                         X(1:4; mode=Sampled(Ordered(), Regular(1), Points())))
end

@testset "mapslices" begin
    a = [1 2 3 4
         3 4 5 6
         5 6 7 8]
    da = DimArray(a, (Y(10:10:30; mode=Sampled(sampling=Intervals())),
                      Ti(1:4; mode=Sampled(sampling=Intervals()))))
    ms = mapslices(sum, da; dims=Y)
    @test ms == [9 12 15 18]
    @test dims(ms) ==
            (Y(10:10:30; mode=Sampled(Ordered(), Regular(10), Intervals(Center()))),
             Ti(1:4; mode=Sampled(Ordered(), Regular(1), Intervals(Start()))))
    @test refdims(ms) == ()
    ms = mapslices(sum, da; dims=Ti)
    @test parent(ms) == [10 18 26]'
end

@testset "array info" begin
    da = DimArray(zeros(5, 4), (Y((10, 20)), X(1:4)))
    @test size(da, Y) == 5
    @test size(da, X()) == 4
    @test axes(da, Y()) == Base.OneTo(5)
    @test axes(da, X) == Base.OneTo(4)
    @test firstindex(da, Y) == 1
    @test firstindex(da, X()) == 1
    @test lastindex(da, Y()) == 5
    @test lastindex(da, X) == 4
end

@testset "cat" begin
    a = [1 2 3; 4 5 6]
    da = DimArray(a, (X(4.0:5.0), Y(6.0:8.0)))
    b = [7 8 9; 10 11 12]
    db = DimArray(b, (X(6.0:7.0), Y(6.0:8.0)))

    @testset "Regular Sampled" begin
        @test cat(da, db; dims=X()) == [1 2 3; 4 5 6; 7 8 9; 10 11 12]
        testdims = (X([4.0, 5.0, 6.0, 7.0]; mode=Sampled(Ordered(), Regular(1.0), Points())),
                    Y(6.0:8.0; mode=Sampled(Ordered(), Regular(1.0), Points())))
        @test cat(da, db; dims=(X(),)) == cat(da, db; dims=X()) == cat(da, db; dims=X)
              cat(da, db; dims=1) == cat(da, db; dims=(1,))
        @test typeof(dims(cat(da, db; dims=X()))) == typeof(testdims)
        @test val(cat(da, db; dims=X())) == val(testdims)
        @test mode(cat(da, db; dims=X())) == mode(testdims)
        @test cat(da, db; dims=Y()) == [1 2 3 7 8 9; 4 5 6 10 11 12]
        @test cat(da, db; dims=Z(1:2)) == cat(a, b; dims=3)
        @test cat(da, db; dims=(Z(1:2), Ti(1:2))) == cat(a, b; dims=(3, 4))
        cat(a, b; dims=(3, 4))
        @test cat(da, db; dims=(X(), Ti(1:2))) == cat(a, b; dims=(1, 3))
        dx = cat(da, db; dims=(X(), Ti(1:2)))
        @test all(map(==, index(dx), index(DimensionalData.formatdims(dx, (X([4.0, 5.0, 6.0, 7.0]), Y(6:8), Ti(1:2))))))
    end

    @testset "Irregular Sampled" begin
        @testset "Intervals" begin
            iri_dim = vcat(X([1, 3, 4], Sampled(Ordered(), Irregular(1, 5), Intervals())),
                          X([7, 8], Sampled(Ordered(), Irregular(7, 9), Intervals())))
            @test span(iri_dim) == Irregular(1, 9)
            @test mode(iri_dim) == Sampled(Ordered(), Irregular(1, 9), Intervals())
            @test index(iri_dim) == [1, 3, 4, 7, 8]
            @test bounds(iri_dim) == (1, 9)
        end
        @testset "Points" begin
            irp_dim = vcat(X([1, 3, 4], Sampled(Ordered(), Irregular(1, 5), Points())),
                          X([7, 8], Sampled(Ordered(), Irregular(7, 9), Points())))
            @test span(irp_dim) == Irregular(nothing, nothing)
            @test mode(irp_dim) == Sampled(Ordered(), Irregular(nothing, nothing), Points())
            @test index(irp_dim) == [1, 3, 4, 7, 8]
            @test bounds(irp_dim) == (1, 8)
        end
    end

    @testset "NoIndex" begin
        ni_dim = vcat(X(Base.OneTo(10), NoIndex()), X(Base.OneTo(10), NoIndex()))
        @test mode(ni_dim) == NoIndex()
        @test index(ni_dim) == Base.OneTo(20)
    end

    @testset "rebuild dim index from refdims" begin
        slices = map(i -> view(da, Y(i)), 1:3)
        cat_da = cat(slices...; dims=Y)
        @test all(cat_da .== da)
        # The range is rebuilt as a Vector during `cat`
        @test index(cat_da) == (4.0:5.0, [6.0, 7.0, 8.0])
        @test index(cat_da) isa Tuple{<:StepRangeLen,<:Vector{Float64}}
    end
end

@testset "unique" begin
    a = [1 1 6; 1 1 6]
    da = DimArray(a, (X(1:2), Y(1:3)))
    @test unique(da; dims=X()) == [1 1 6]
    @test unique(da; dims=Y) == [1 6; 1 6]
    @test unique(da; dims=:) == [1, 6]
    @test unique(da[X(1)]) == [1, 6]
end

@testset "diff" begin
    @testset "Array 2D" begin
        y = Y(['a', 'b', 'c'])
        ti = Ti(DateTime(2021, 1):Month(1):DateTime(2021, 4))
        data = [-87  -49  107  -18
                24   44  -62  124
                122  -11   48   -7]
        A = DimArray(data, (y, ti))
        @test diff(A; dims=1) == diff(A; dims=Y) == diff(A; dims=:Y) == diff(A; dims=y) ==
            DimArray([111 93 -169 142; 98 -55 110 -131], (Y(['a', 'b']), ti))
        @test diff(A; dims=2) == diff(A; dims=Ti) == diff(A; dims=:Ti) == diff(A; dims=ti) ==
            DimArray([38 156 -125; 20 -106 186; -133 59 -55], (y, Ti(DateTime(2021, 1):Month(1):DateTime(2021, 3))))
        @test_throws ErrorException diff(A; dims='X')
        @test_throws ArgumentError diff(A; dims=Z)
        @test_throws ArgumentError diff(A; dims=3)
    end
    @testset "Vector" begin
        x = DimArray([56, -123, -60, -44, -64, 70, 52, -48, -74, 86], X(2:2:20))
        @test diff(x) == diff(x; dims=1) == diff(x; dims=X) == DimArray([-179, 63, 16, -20, 134, -18, -100, -26, 160], X(2:2:18))
    end
end


text:
# ================================================================
# Time-shifted reach-set
# ================================================================

struct ShiftedReachSet{N, RT<:AbstractLazyReachSet{N}} <: AbstractLazyReachSet{N}
    R::RT
    t0::N
end

# getter functions
@inline time_shift(srs::ShiftedReachSet) = srs.t0

# time domain interface
@inline tstart(srs::ShiftedReachSet) = tstart(srs.R) + time_shift(srs)
@inline tend(srs::ShiftedReachSet) = tend(srs.R) + time_shift(srs)
@inline tspan(srs::ShiftedReachSet) = TimeInterval(tstart(srs), tend(srs))


text:
#
# Distance function, considering PBC
#
function rpbc(x,y,side)
  # dx is a vector of the difference in position of y and x
  d = 0.
  for i in 1:2
    dx = (y[i]-x[i])%side
    if dx > side/2
      dx = dx - side
    elseif  dx < -side/2
      dx = dx + side
    end
    d = d + dx^2
  end
  return sqrt(d)
end


text:
import FactCheck2TestSets: fact2testcore, fc2ts

@testset "@fact" begin
    @testset "rhs bool" begin
        ex = quote @fact the(truth) --> true end
        @test fact2testcore(ex) == :(the(truth))

        ex = quote @fact alie --> false end
        @test fact2testcore(ex) == :(!(alie))
    end

    @testset "rhs not" begin
        ex = quote @fact 1 --> not(2) end
        @test fact2testcore(ex) == :(!(1 == 2))

        ex = quote @fact foo --> not(isnan) end
        @test fact2testcore(ex) == :(!(isnan(foo)))
    end

    @testset "rhs roughly" begin
        ex = quote @fact 1 --> roughly(foo(bar)) end
        @test  fact2testcore(ex) == :(1 ≈ foo(bar))

        ex = quote @fact 1 --> roughly(2, 2) end
        @test fact2testcore(ex) == :(isapprox(1,2,atol=2))

        ex = quote @fact 1 --> roughly(2, rtol=3) end
        @test fact2testcore(ex) == :(isapprox(1,2,rtol=3))
    end

    @testset "rhs equals" begin
        ex = quote @fact a --> b(x) end
        @test fact2testcore(ex) == :(a == b(x))
    end

    @testset "rhs property" begin
        ex = quote @fact f(x, y, atol=3) --> isfinite end
        @test fact2testcore(ex) == :(isfinite(f(x,y,atol=3)))
    end

    @testset "rhs comparison" begin
        ex = quote @fact a --> less_than(b) end
        @test fact2testcore(ex) == :(a < b)
    end

end

@testset "@fact_throws" begin
    ex = :(@fact_throws DimensionMismatch x[3])
    @test fc2ts(ex) == :(@test_throws DimensionMismatch x[3])
    ex = :(@fact_throws x[3])
    @test fc2ts(ex) == :(@test_throws Exception x[3])
end


text:
using SymPoly
using SymbolicUtils, Test
using DynamicPolynomials

function generate_rand_poly(x; min_deg=0, max_deg=5, sparcity=0.5)
    n = rand(min_deg:max_deg)
    p = one(x)
    r = rand(1:10) // rand(1:10)
    while true
        for i = 0:n
            if rand() > sparcity || i == n
                s = rand() < 0.5 ? +1 : -1
                p += x^i * s * r * rand(0:10)//rand(1:5)
            end
        end
        !isequal(p, 0) && return p
        return p
    end
end

function test_eq(x, f, name; n=10, min_deg=0, max_deg=10, sparcity=0.5)
    k = 0
    outcome = true
    for i = 1:n
        p = generate_rand_poly(x; max_deg, sparcity)
        q = generate_rand_poly(x; min_deg, max_deg, sparcity)

        if iszero(p) || iszero(q) continue end

        printstyled("P1 = ", p, '\n'; color=:green)
        printstyled("P2 = ", q, '\n'; color=:blue)

        try
            r = f(p, q)
            printstyled("Δ($name) = ", r, '\n'; color=:red)
            if !iszero(r)
                @warn "+/- mismatch: $r"
                outcome = false
            else
                k += 1
            end
        catch e
            println(e)
        end
    end
    outcome
end


text:
module PoissonTestsSeq
using PartitionedArrays
include("../PoissonTests.jl")
prun(PoissonTests.main,sequential,(2,2))
end # module


text:
n = rand(1:10)

A = Array(Float32, n, n)
copy!(A, rand(n,n))

r1 = polarfact(A, alg =:newton);

r2 = polarfact(A, alg =:qdwh);

r3 = polarfact(A, alg =:halley);

r4 = polarfact(A, alg =:svd);

r5 = polarfact(A, alg =:hybrid);

@test_approx_eq r1.U r2.U
@test_approx_eq r2.U r3.U
@test_approx_eq r3.U r4.U
@test_approx_eq r4.U r5.U

println("Type Float32 passed test...")


text:
################################################################################
# Izhikevich
################################################################################
@inline f1(v, u, input) = (0.04*v + 5.0)*v + 140.0 - u + input
@inline f2(v, u, a, b)  = a*(b*v - u)

@inline function izhupdate_rk4(v, u, input, h, p::SimpleNeuron)
    half_h = 0.5h
    sixth_h = 0.16666666666666666h

    K1 = f1(v, u, 0.0)
    L1 = f2(v, u, p.a, p.b)
    K2 = f1(v + half_h*K1, u + half_h*L1, 0.0)
    L2 = f2(v + half_h*K1, u + half_h*L1, p.a, p.b)
    K3 = f1(v + half_h*K2, u + half_h*L2, 0.0)
    L3 = f2(v + half_h*K2, u + half_h*L2, p.a, p.b)
    K4 = f1(v + h*K3, u + h*L3, 0.0)
    L4 = f2(v + h*K3, u + h*L3, p.a, p.b)

    (v + sixth_h * (K1 + 2K2 + 2K3 + K4) + h*input,
     u + sixth_h * (L1 + 2L2 + 2L3 + L4))
end

function izhupdate_euler(v, u, input, h, p::SimpleNeuron)
    (v + h*((0.04v + 5.0)*v + 140.0 - u + input),
     u + h*(p.a*(p.b*v - u)))
end

################################################################################
# Gerstner
################################################################################
function resetsrmneuron(n::SRMNeuron)
	n.u = -65.0
	n.u_max = 8.0
	n.u_rest = -65.0
	n.Θ = -50.0
	n.τ_refr = 7.0/1000.0
	n.τ_memb = 3.0/1000.0
	n.t_lastevent = 0.0
	n.t_lastspike = 0.0
	n.ϵ_last = 0.0
	return n
end

function srmneuronupdate(n::SRMNeuron, t, bumpsize)
	spiked = false
	ϵ_new = n.ϵ_last*exp((n.t_lastevent - t)/n.τ_memb) + bumpsize*n.u_max
	n.t_lastevent = t
	n.u = n.u_rest + ϵ_new
	if n.u > n.Θ && t - n.t_lastspike > n.τ_refr
		spiked = true
		n.t_lastspike = t
		n.u = n.u_rest
		n.ϵ_last = 0.0
	else
		n.ϵ_last = ϵ_new
	end
	return spiked
end

################################################################################
# Delay Lines
################################################################################
function clearbuffer(dl::DelayLine)
	dl.numstored = 0
	dl.counts .*= 0
end

function pushevent(dl::DelayLine)
	if dl.numstored < length(dl.counts)
		dl.numstored += 1
		dl.counts[dl.numstored] = dl.delaylen;
	else
		println("buffer overflow")
	end
end

function advance(dl::DelayLine)
	out = 0
	if dl.counts[1] == 1
		dl.counts[1] = 0
		out = 1
		for i ∈ 2:dl.numstored
			dl.counts[i-1] = dl.counts[i]
		end
		dl.numstored -= 1
	end
	for i ∈ 1:dl.numstored
		dl.counts[i] -= 1
	end
	return out
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
## Drag polar analysis case
using AeroMDAO
using DataFrames
using Plots

## Wing
wing_foils = Foil.(fill(naca4((0,0,1,2)), 2))
wing       = Wing(foils     = wing_foils,
                  chords    = [1.0, 0.6],
                  twists    = [0.0, 0.0],
                  spans     = [5.0],
                  dihedrals = [11.3],
                  sweep_LEs = [2.29]);
print_info(wing, "Wing")

# Horizontal tail
htail_foils = Foil.(fill(naca4((0,0,1,2)), 2))
htail       = Wing(foils     = htail_foils,
                   chords    = [0.7, 0.42],
                   twists    = [0.0, 0.0],
                   spans     = [1.25],
                   dihedrals = [0.],
                   sweep_LEs = [6.39])
print_info(htail, "Horizontal Tail")

# Vertical tail
vtail_foils = Foil.(fill(naca4((0,0,1,2)), 2))
vtail = HalfWing(foils     = vtail_foils,
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.0],
                 dihedrals = [0.],
                 sweep_LEs = [7.97])
print_info(vtail, "Vertical Tail")

# Assembly
wing_panels  = panel_wing(wing, [20], 10);
htail_panels = panel_wing(htail, [12], 12;
                          position	= [4., 0, 0],
                          angle 	= deg2rad(-0.),
                          axis 	  	= [0., 1., 0.]
                         )
vtail_panels = panel_wing(vtail, [12], 10;
                          position 	= [4., 0, 0],
                          angle 	= π/2,
                          axis 	 	= [1., 0., 0.]
                         )

aircraft = Dict(
                "Wing" 			  	=> wing_panels,
                # "Horizontal Tail" 	=> htail_panels,
                # "Vertical Tail"   	=> vtail_panels
                )

S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing)

function vlm_analysis(aircraft, fs, ρ, ref, S, b, c, print = false)
    # Evaluate case
    data =  solve_case(aircraft, fs;
                       rho_ref   = ρ,
                       r_ref     = ref,
                       area_ref  = S,
                       span_ref  = b,
                       chord_ref = c,
                       print     = print
                      );

    # Get data
    nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes, Γs = data["Aircraft"]

    # Filter relevant data
    [ ff_coeffs; nf_coeffs[4:6] ]
end

##
ac_name = "My Aircraft"
ρ       = 1.225
ref     = [0.25c, 0., 0.]
V, β    = 1.0, 0.0
Ω       = [0.0, 0.0, 0.0]
αs      = -5:0.5:5
fses    = map(α -> Freestream(V, α, β, Ω), αs);

## Evaluate cases
results = vlm_analysis.(Ref(aircraft), fses, ρ, Ref(ref), S, b, c, true);

## Generate DataFrame
data = DataFrame([ (α, CD, CY, CL, Cl, Cm, Cn) for (α, (CD, CY, CL, Cl, Cm, Cn)) in zip(αs, results) ][:])
rename!(data, [:α, :CD, :CY, :CL, :Cl, :Cm, :Cn])

##
plotly(dpi = 300)
plot(xlabel = "CD", ylabel = "CL", title = "Drag Polar")
plot!(data[!,"CD"], data[!,"CL"], label = :none, marker = :dot)

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
# functions related to normal distribution

function xval(μ::Real, σ::Real, z::Number)
    if isinf(z) && iszero(σ)
        μ + one(σ) * z
    else
        μ + σ * z
    end
end
zval(μ::Real, σ::Real, x::Number) = (x - μ) / σ

# pdf
normpdf(z::Number) = exp(-abs2(z)/2) * invsqrt2π
function normpdf(μ::Real, σ::Real, x::Number)
    if iszero(σ)
        if x == μ
            z = zval(μ, one(σ), x)
        else
            z = zval(μ, σ, x)
            σ = one(σ)
        end
    else
        z = zval(μ, σ, x)
    end
    normpdf(z) / σ
end

# logpdf
normlogpdf(z::Number) = -(abs2(z) + log2π)/2
function normlogpdf(μ::Real, σ::Real, x::Number)
    if iszero(σ)
        if x == μ
            z = zval(μ, one(σ), x)
        else
            z = zval(μ, σ, x)
            σ = one(σ)
        end
    else
        z = zval(μ, σ, x)
    end
    normlogpdf(z) - log(σ)
end

# cdf
normcdf(z::Number) = erfc(-z * invsqrt2)/2
function normcdf(μ::Real, σ::Real, x::Number)
    if iszero(σ) && x == μ
        z = zval(zero(μ), σ, one(x))
    else
        z = zval(μ, σ, x)
    end
    normcdf(z)
end
# ccdf
normccdf(z::Number) = erfc(z * invsqrt2)/2
function normccdf(μ::Real, σ::Real, x::Number)
    if iszero(σ) && x == μ
        z = zval(zero(μ), σ, one(x))
    else
        z = zval(μ, σ, x)
    end
    normccdf(z)
end

# logcdf
normlogcdf(z::Number) = z < -1.0 ?
    log(erfcx(-z * invsqrt2)/2) - abs2(z)/2 :
    log1p(-erfc(z * invsqrt2)/2)
function normlogcdf(μ::Real, σ::Real, x::Number)
    if iszero(σ) && x == μ
        z = zval(zero(μ), σ, one(x))
    else
        z = zval(μ, σ, x)
    end
    normlogcdf(z)
end

# logccdf
normlogccdf(z::Number) = z > 1.0 ?
    log(erfcx(z * invsqrt2)/2) - abs2(z)/2 :
    log1p(-erfc(-z * invsqrt2)/2)
function normlogccdf(μ::Real, σ::Real, x::Number)
    if iszero(σ) && x == μ
        z = zval(zero(μ), σ, one(x))
    else
        z = zval(μ, σ, x)
    end
    normlogccdf(z)
end

norminvcdf(p::Real) = -erfcinv(2*p) * sqrt2
norminvcdf(μ::Real, σ::Real, p::Real) = xval(μ, σ, norminvcdf(p))

norminvccdf(p::Real) = erfcinv(2*p) * sqrt2
norminvccdf(μ::Real, σ::Real, p::Real) = xval(μ, σ, norminvccdf(p))

# invlogcdf. Fixme! Support more precisions than Float64
norminvlogcdf(lp::Union{Float16,Float32}) = convert(typeof(lp), _norminvlogcdf_impl(Float64(lp)))
norminvlogcdf(lp::Real) = _norminvlogcdf_impl(Float64(lp))
norminvlogcdf(μ::Real, σ::Real, lp::Real) = xval(μ, σ, norminvlogcdf(lp))

# invlogccdf. Fixme! Support more precisions than Float64
norminvlogccdf(lp::Union{Float16,Float32}) = convert(typeof(lp), -_norminvlogcdf_impl(Float64(lp)))
norminvlogccdf(lp::Real) = -_norminvlogcdf_impl(Float64(lp))
norminvlogccdf(μ::Real, σ::Real, lp::Real) = xval(μ, σ, norminvlogccdf(lp))

# norminvcdf & norminvlogcdf implementation
#
#   Rational approximations for the inverse cdf and its logarithm, from:
#
#   Wichura, M.J. (1988) Algorithm AS 241: The Percentage Points of the Normal Distribution
#   Journal of the Royal Statistical Society. Series C (Applied Statistics), Vol. 37, No. 3, pp. 477-484
#

function _norminvlogcdf_impl(lp::Float64)
    if isfinite(lp) && lp < 0.0
        q = exp(lp) - 0.5
        # qnorm_kernel(lp, q, true)
        if abs(q) <= 0.425
            _qnorm_ker1(q)
        else
            r = sqrt(q < 0 ? -lp : -log1mexp(lp))
            return copysign(_qnorm_ker2(r), q)
        end
    elseif lp >= 0.0
        lp == 0.0 ? Inf : NaN
    else # lp is -Inf or NaN
        lp
    end
end

function _qnorm_ker1(q::Float64)
    # pre-condition: abs(q) <= 0.425
    r = 0.180625 - q*q
    return q * @horner(r,
                       3.38713_28727_96366_6080e0,
                       1.33141_66789_17843_7745e2,
                       1.97159_09503_06551_4427e3,
                       1.37316_93765_50946_1125e4,
                       4.59219_53931_54987_1457e4,
                       6.72657_70927_00870_0853e4,
                       3.34305_75583_58812_8105e4,
                       2.50908_09287_30122_6727e3) /
    @horner(r,
            1.0,
            4.23133_30701_60091_1252e1,
            6.87187_00749_20579_0830e2,
            5.39419_60214_24751_1077e3,
            2.12137_94301_58659_5867e4,
            3.93078_95800_09271_0610e4,
            2.87290_85735_72194_2674e4,
            5.22649_52788_52854_5610e3)
end

function _qnorm_ker2(r::Float64)
    if r < 5.0
        r -= 1.6
        @horner(r,
                1.42343_71107_49683_57734e0,
                4.63033_78461_56545_29590e0,
                5.76949_72214_60691_40550e0,
                3.64784_83247_63204_60504e0,
                1.27045_82524_52368_38258e0,
                2.41780_72517_74506_11770e-1,
                2.27238_44989_26918_45833e-2,
                7.74545_01427_83414_07640e-4) /
        @horner(r,
                1.0,
                2.05319_16266_37758_82187e0,
                1.67638_48301_83803_84940e0,
                6.89767_33498_51000_04550e-1,
                1.48103_97642_74800_74590e-1,
                1.51986_66563_61645_71966e-2,
                5.47593_80849_95344_94600e-4,
                1.05075_00716_44416_84324e-9)
    else
        r -= 5.0
        @horner(r,
                6.65790_46435_01103_77720e0,
                5.46378_49111_64114_36990e0,
                1.78482_65399_17291_33580e0,
                2.96560_57182_85048_91230e-1,
                2.65321_89526_57612_30930e-2,
                1.24266_09473_88078_43860e-3,
                2.71155_55687_43487_57815e-5,
                2.01033_43992_92288_13265e-7) /
        @horner(r,
                1.0,
                5.99832_20655_58879_37690e-1,
                1.36929_88092_27358_05310e-1,
                1.48753_61290_85061_48525e-2,
                7.86869_13114_56132_59100e-4,
                1.84631_83175_10054_68180e-5,
                1.42151_17583_16445_88870e-7,
                2.04426_31033_89939_78564e-15)
    end
end


text:
function readAttributes( mpSim::MPsim, xf::XF, catXF::XF )

    if !XLSX.hassheet( xf, "Attributes" )
        @warn "Configuration file is missing an 'Attributes' sheet, cannot configure attributes."
        return
    end  # if !XLSX.hassheet( xf, "Attributes" )

    ws = xf[ "Attributes" ]
    nAttributes = ws[ "B4" ] isa Integer ? ws[ "B4" ] : 0
    lineNrs = 3 * (1:nAttributes) .+ 4

    if XLSX.hassheet( catXF, "Attributes" )
        catalogue = catXF[ "Attributes" ]
        attributes = readAttribute.( Ref( ws ), Ref( catalogue ),
            lineNrs )
    else
        attributes = readAttribute.( Ref( ws ), lineNrs )
    end  # if XLSX.hassheet( catXF, "Attributes" )

    setSimulationAttributes!( mpSim, attributes )

end  # readAttributes( mpSim, xf, catXF )

function readAttribute( ws::WS, catalogue::WS, sLine::Int )

    # Read the attribute's initial values.
    attribute = readAttribute( ws, sLine )

    # If it's in the catalogue, read the attribute's possible values.
    cLine = ws[ sLine, 2 ]

    if !( cLine isa Integer ) || ( cLine <= 0 )
        return attribute
    end  # if !( cLine isa Integer ) || ...

    cLine += 1

    if string( catalogue[ cLine, 1 ] ) != attribute.name
        return attribute
    end  # if string( catalogue[ cLine, 1 ] != attribute.name

    nValues = catalogue[ cLine, 5 ]

    if !( nValues isa Integer ) || ( nValues <= 0 )
        return attribute
    end  # if !( nValues isa Integer ) || ...

    values = String.( strip.( string.( catalogue[ CR( cLine, 6, cLine,
        5 + nValues ) ] ) ) )
    addPossibleAttributeValue!( attribute, values... )

    return attribute

end  # readAttribute( ws, catalogue, sLine )

function readAttribute( ws::WS, sLine::Int )

    attribute = Attribute( string( ws[ sLine, 1 ] ) )
    nInitialValues = ws[ sLine + 2, 2 ]

    # No need to process any further if there are no initial values given.
    if !( nInitialValues isa Integer )
        return attribute, false
    elseif ( nInitialValues <= 0 )
        return attribute, true
    end  # if !( nInitialValues isa Integer )

    initialValues = Dict{String, Float64}()
    values = strip.( string.( ws[ CR( sLine, 3, sLine,
        nInitialValues + 2 ) ] ) )
    weights = ws[ CR( sLine + 1, 3, sLine + 1, nInitialValues + 2 ) ]

    # Only process initial values with a given weight.
    for ii in 1:nInitialValues
        if weights[ ii ] isa Real
            initialValues[ values[ ii ] ] = weights[ ii ] +
                get( initialValues, values[ ii ], 0.0 )
        end  # if weights[ ii ] isa Real
    end  # for ii in 1:nInitialValues

    setInitialAttributeValues!( attribute, initialValues )
    return attribute

end  # readAttribute( ws, sLine )

function readAttributes( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    attributePars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Attribute'" ) ) )

    if isempty( attributePars )
        return
    end  # if isempty( attributePars )

    attributes = map( eachindex( attributePars[ :, :parName ] ) ) do ii
        attribute = Attribute( attributePars[ ii, :parName ] )
        pars = split( attributePars[ ii, :parValue ], ";" )

        if length( pars[ 1 ] ) > 2
            setPossibleAttributeValues!( attribute, string.(
                split( pars[ 1 ][ 2:(end-1) ], "," ) ) )
        end  # if length( pars[ 1 ] ) > 2

        if length( pars[ 2 ] ) > 2
            initVals = split.( split( pars[ 2 ][ 2:(end-1) ], "," ), ":" )
            initValDict = Dict{String, Float64}()

            for initVal in initVals
                initValDict[ initVal[ 1 ] ] = parse( Float64, initVal[ 2 ] )
            end  # for initVal in initVals

            setInitialAttributeValues!( attribute, initValDict )
        end  # if length( pars[ 2 ] ) > 2

        return attribute
    end  # map( eachindex( attributePars[ :parName ] ) ) do ii

    setSimulationAttributes!( mpSim, attributes )

end  # readAttributes( mpSim, configDB, configName )


text:
function log_oneplusexp(E::Float64)
    return (E < 50) ? log(1+exp(E)) : E
end

function xlogy(x::Float64, y::Float64)
    if x == 0.0
        return 0.0
    else
        return x*log(y)
    end
end

function log_binom(n::Float64, k::Float64)
    @assert n >= k
    if (n == 0) || (k == 0)
       return 0
    end
    return n*log(n)-k*log(k)-(n-k)*log(n-k)
end


text:
"""

```
mddf_frame!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)
```

Computes the MDDF for a single frame, modifies the data in the `R` (type `Result`) structure.

"""
function mddf_frame!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)

  # Simplify code by assigning some shortened names
  trajectory = framedata.trajectory
  volume_frame = framedata.volume_frame
  rdf_count_random_frame = framedata.rdf_count_random_frame
  md_count_random_frame = framedata.md_count_random_frame
  dc = framedata.dc
  dmin_mol = framedata.dmin_mol
  dref_mol = framedata.dref_mol
  x_solvent_random = framedata.x_solvent_random
  lc_solvent = framedata.lc_solvent
  solute = trajectory.solute
  solvent = trajectory.solvent
  x_solute = trajectory.x_solute
  x_solvent = trajectory.x_solvent

  # Reset counters for this frame
  reset!(volume_frame)
  @. rdf_count_random_frame = 0.
  @. md_count_random_frame = 0.

  # get pbc sides in this frame
  sides = getsides(trajectory,iframe)

  volume_frame.total = sides[1]*sides[2]*sides[3]
  R.volume.total = R.volume.total + volume_frame.total

  R.density.solute = R.density.solute + (solute.nmols / volume_frame.total)
  R.density.solvent = R.density.solvent + (solvent.nmols / volume_frame.total)

  # Add the box side information to the box structure, in this frame
  box = Box(options.lcell,sides,R.cutoff)

  # Will wrap everthing relative to the reference atom of the first molecule
  # and move everything such that that center is in the origin. This is important
  # to simplify the computation of cell indexes, as the minimum coordinates are
  # automatically -side/2 at each direction. The reason for wrapping to the solute
  # center and not to the origin is that if the solute is a protein, for example,
  # the protein gets centered, and hardly ever it will be necessary to wrap the
  # coordinates afterwards, accelerating the computation.
  solute_center = x_solute[R.irefatom]
  wrap!(x_solute,sides,solute_center)
  center_to_origin!(x_solute,solute_center)
  wrap!(x_solvent,sides,solute_center)
  center_to_origin!(x_solvent,solute_center)

  # Initialize linked cells
  initcells!(x_solvent,box,lc_solvent)

  # Check if the cutoff is not too large considering the periodic cell size
  if R.cutoff > sides[1]/2. || R.cutoff > sides[2]/2. || R.cutoff > sides[3]/2.
    error("""
    in MDDF: cutoff or dbulk > periodic_dimension/2 in frame: $iframe
             max(cutoff,dbulk) = $(R.cutoff)
             sides read from file = $(sides)
             If sides are zero it is likely that the PBC information is not available
             in the trajectory file.
    """)
  end

  local n_dmin_in_bulk
  n_solvent_in_bulk = 0.
  for isolute in 1:solute.nmols

    # We need to do this one solute molecule at a time to avoid exploding the memory
    # requirements
    x_this_solute = viewmol(isolute,x_solute,solute)

    # Compute all distances between solute and solvent atoms which are smaller than the
    # cutoff (this is the most computationally expensive part), the distances are returned
    # in the dc structure
    cutoffdistances!(R.cutoff,x_this_solute,x_solvent,lc_solvent,box,dc)

    # For each solute molecule, update the counters (this is highly suboptimal, because
    # within updatecounters there are loops over solvent molecules, in such a way that
    # this will loop with cost nsolute*nsolvent. However, I cannot see an easy solution
    # at this point with acceptable memory requirements
    n_dmin_in_bulk, n_dref_in_bulk = updatecounters!(R,solute,solvent,dc,dmin_mol,dref_mol)
    n_solvent_in_bulk += n_dref_in_bulk
  end
  n_solvent_in_bulk = n_solvent_in_bulk / solute.nmols

  #
  # Computing the random-solvent distribution to compute the random minimum-distance count
  #
  bulk_range = (solvent.nmols-n_dmin_in_bulk+1):solvent.nmols
  for i in 1:options.n_random_samples

    # generate random solvent box, and store it in x_solvent_random
    for j in 1:solvent.nmols
      # Choose randomly one molecule from the bulk, if there are actually bulk molecules
      if n_dmin_in_bulk != 0
        jmol = dmin_mol[random(RNG,bulk_range)].jmol
      else
        jmol = random(RNG,1:solvent.nmols)
      end
      # Indexes of this molecule in the x_solvent array
      x_ref = viewmol(jmol,x_solvent,solvent)
      # Indexes of the random molecule in random array
      x_rnd = viewmol(j,x_solvent_random,solvent)
      # Generate new random coordinates (translation and rotation) for this molecule
      random_move!(x_ref,R.irefatom,sides,x_rnd,RNG)
    end

    # wrap random solvent coordinates to box, with the center at the origin
    wrap!(x_solvent_random,sides)

    # Initialize linked cells
    initcells!(x_solvent_random,box,lc_solvent)

    # Choose randomly one solute molecule to be the solute in this sample
    i_rand_mol = random(RNG,1:solute.nmols)
    x_this_solute = viewmol(i_rand_mol,x_solute,solute)

    # Compute all distances between solute and solvent atoms which are smaller than the
    # cutoff (this is the most computationally expensive part), the distances are returned
    # in the dc structure
    cutoffdistances!(R.cutoff,x_this_solute,x_solvent_random,lc_solvent,box,dc)

    # Update the counters of the random distribution
    updatecounters!(R,rdf_count_random_frame,md_count_random_frame,solvent,dc,dmin_mol,dref_mol)

  end # random solvent sampling

  # Update counters with the data of this frame
  update_counters_frame!(R, rdf_count_random_frame, md_count_random_frame, volume_frame,
                         solute, solvent, n_solvent_in_bulk)

  nothing
end


text:
module Day03

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day03PuzzleInfo <: PuzzleInfo end

"""
Day 3 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day03PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 3 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day03PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Get number of characters in first line of input stream `input`.
"""
function getfirstlinelength(input::IO)
    firstline = readline(input)
    length(firstline)
end

"""
    $(FUNCTIONNAME)(bitvector)

Convert binary representation of integer `bitvector` to corresponding (decimal) integer.
"""
function bitvectortodecimal(bitvector::BitVector)
    sum(2^(i-1) for (i, bit) in enumerate(reverse(bitvector)) if bit)
end

"""
    $(FUNCTIONNAME)(input)

Read 2D bit array from input stream `input` with (equi-length) string of bits per line.
"""
function readbitarray(input::IO)
    reduce(vcat, (BitArray(parse(UInt8, b) for b in l)' for l in eachline(input)))
end

function solve(info::Part1PuzzleInfo)
    bitarray = open(readbitarray, info.datapath)
    counts = vec(sum(bitarray, dims=1))
    numline = size(bitarray, 1)
    @assert all(counts .!= (numline / 2)) "Equal number of one and zero bits"
    moreones = counts .>= (numline / 2)
    gammarate = bitvectortodecimal(moreones)
    epsilonrate = bitvectortodecimal(.!moreones)
    gammarate * epsilonrate
end

function filterbycriteriaoncolumn(criteria, bitarray::BitArray, column::Integer)
    bitarray[criteria(bitarray[:, column]), :]
end

function oxygengenerator_bitcriteria(bitvector::BitVector)
    bitvector .== (sum(bitvector) .>= size(bitvector, 1) / 2)
end

function co2scrubber_bitcriteria(bitvector::BitVector)
    bitvector .== (sum(bitvector) .< size(bitvector, 1) / 2)
end

function solve(info::Part2PuzzleInfo)
    bitarray = open(readbitarray, info.datapath)
    oxygenbitarray, co2bitarray = bitarray, bitarray
    for column in 1:size(bitarray, 2)
        if size(oxygenbitarray, 1) != 1
            oxygenbitarray = filterbycriteriaoncolumn(
                oxygengenerator_bitcriteria, oxygenbitarray, column
            )
        end
        if size(co2bitarray, 1) != 1
            co2bitarray = filterbycriteriaoncolumn(
                co2scrubber_bitcriteria, co2bitarray, column
            )
        end
        size(oxygenbitarray, 1) == 1 && size(co2bitarray, 1) == 1 && break
    end
    oxygengeneratorrate = bitvectortodecimal(oxygenbitarray[1, :])
    co2scrubberrate = bitvectortodecimal(co2bitarray[1, :])
    oxygengeneratorrate * co2scrubberrate
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day03.txt"))) == 198
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day03.txt"))) == 230
end

end # module


text:
export TXN, TXT, gemv!

# Transpose Operators
const TXN = Val('N')
const TXT = Val('T')

# a * A * x + by → y, @code_warntype ✓
function gemv!(::Val{'N'}, a::Real, A::MatI{TA}, x::VecI{Tx}, b::Real, y::VecO{Ty}) where {TA<:Real,Tx<:Real,Ty<:Real}
    scal!(b, y)                 # y := b * y
    iszero(a) && return nothing # quick return if possible
    #### y := a * A * x + y
    if isone(a)
        for j in eachindex(x)
            @inbounds xj = x[j]
            if !iszero(xj)
                @simd for i in eachindex(y)
                    @inbounds y[i] += A[i, j] * xj
                end
            end
        end
        return nothing
    end
    for j in eachindex(x)
        @inbounds axj = a * x[j]
        if !iszero(axj)
            @simd for i in eachindex(y)
                @inbounds y[i] += A[i, j] * axj
            end
        end
    end
end
# @code_warntype ✓
function gemv!(::Val{'N'}, a::Real, A::MatI{TA}, x::NTuple{N,Tx}, b::Real, y::VecO{Ty}) where {N,TA<:Real,Tx<:Real,Ty<:Real}
    scal!(b, y)                 # y := b * y
    iszero(a) && return nothing # quick return if possible
    #### y := a * A * x + y
    if isone(a)
        for j in eachindex(1:N)
            @inbounds xj = x[j]
            if !iszero(xj)
                @simd for i in eachindex(y)
                    @inbounds y[i] += A[i, j] * xj
                end
            end
        end
        return nothing
    end
    for j in eachindex(1:N)
        @inbounds axj = a * x[j]
        if !iszero(axj)
            @simd for i in eachindex(y)
                @inbounds y[i] += A[i, j] * axj
            end
        end
    end
end


text:
include("decision_variable.jl")
include("decision_constraint.jl")
include("decision_objective.jl")
include("solve.jl")


text:
using JuLIP, JuLIP.Potentials, SKTB
using JuLIP.Potentials: site_energy, site_energy_d

TB=SKTB

@testset "Kwon" begin

println("Test Kwon TB Model")
at = (1,2,2) * bulk(:Si, pbc=(true, false, true), cubic=true)
@show length(at)
β, fixed_eF = 30.0, true
tbm = TB.Kwon.KwonTBModel(potential = FermiDiracSmearing(β),
                          # potential = TB.GrandPotential(β, 0.0),
                          # bzquad=TB.MPGrid(at, (4,0,0)) )
                          bzquad = TB.GammaPoint() )
print("Test setting Nel ... ")
TB.set_δNel!(tbm, at)
println("ok.")
@test true
@show TB.get_eF(tbm.potential)

# perturb positions a bit
X = positions(at)
X[2] += 0.1 * rand(JVecF)
set_positions!(at, X)

print("check that hamiltonian evaluates ... ")
H, M = hamiltonian(tbm, at)
println("ok.")
print("check that `densitymatrix` evaluates ... ")
Γ = densitymatrix(tbm, at)
println("ok.")
print("check that `energy` evaluates ... ")
E = energy(tbm, at)
println("ok : E = ", E, ".")
print("check that `forces` evaluates ... ")
frc = forces(tbm, at) |> mat
println("ok : |f|∞ = ", norm(frc, Inf), ".")
@test true

println("  Finite-difference test:  ")
@test fdtest(tbm, at, verbose=true)
println("  Finite-difference test without the repulsive part:  ")
Vrep, tbm.Vrep = tbm.Vrep, JuLIP.Potentials.ZeroSitePotential()
@test fdtest(tbm, at, verbose=true)
tbm.Vrep = Vrep

# ==================== site energy tests ==================

n0 = 1
NQUAD = (4, 8, 16, 20, 25, 30)

E = energy(tbm, at)
∑En = sum( site_energy(tbm, at, n) for n = 1:length(at) )
println("Testing that the spectral-decompositions site energy sums to total energy")
@show E - ∑En
@test abs(E - ∑En) < 1e-10

calc = TB.PEXSI(tbm, 5, [n0])
JuLIP.rattle!(at, 0.02)
# calibrate the PEXSI calculator on a mini-system
print("calibrating . . . ")
TB.calibrate!(calc, at; at_train = bulk(:Si, pbc=true), npoles = 8)
println("done."); @test true
# output some useful info if we are watching the tests...
@show length(at)
@show TB.get_EminEmax(at)

# compute the site energy the old way and compare against the PEXSI calculation
Eold = TB.site_energy(tbm, at, n0)
println("Old Site Energy (via spectral decomposition): ", Eold)
println("Testing convergence of PEXSI site energy")
Enew = 0.0
for nquad in NQUAD
   TB.set_npoles!(calc, nquad)
   Enew = site_energy(calc, at, n0)
   println("nquad = ", nquad, "; error = ", abs(Enew - Eold))
end
@show Enew, Eold
@test abs(Enew - Eold) < 1e-5

println("Test consistency of site forces")
TB.set_npoles!(calc, 8)
X = positions(at) |> mat
Es = site_energy(calc, at, n0)
dEs = site_energy_d(calc, at, n0) |> mat
dEsh = []
errors = Float64[]
for p = 2:9
   h = 0.1^p
   dEsh = zero(dEs)
   for n = 1:length(X)
      X[n] += h
      set_positions!(at, X)
      Esh = site_energy(calc, at, n0)
      dEsh[n] = (Esh - Es) / h
      X[n] -= h
   end
   println( " ", p, " | ", norm(dEs-dEsh, Inf) )
   push!(errors, norm(dEs-dEsh, Inf))
end
@test minimum(errors) < 1e-3 * maximum(errors)

end