text:
immutable Exponential <: ContinuousUnivariateDistribution
    scale::Float64 # note: scale not rate
    function Exponential(sc::Real)
        sc > zero(sc) || error("scale must be positive")
        new(float64(sc))
    end
end

Exponential() = Exponential(1.0)

scale(d::Exponential) = d.scale
rate(d::Exponential) = 1.0 / d.scale

cdf(d::Exponential, q::Real) = q <= zero(q) ? 0.0 : -expm1(-q / d.scale)

logcdf(d::Exponential, q::Real) = q > zero(q) ? log1mexp(-q / d.scale) : -Inf

ccdf(d::Exponential, q::Real) = q <= zero(q) ? 1.0 : exp(-q / d.scale)

logccdf(d::Exponential, q::Real) = q <= zero(q) ? 0.0 : -q / d.scale

invlogcdf(d::Exponential, lp::Real) = lp > zero(lp) ? NaN : -d.scale * log1mexp(lp)

invlogccdf(d::Exponential, lp::Real) = lp <= zero(lp) ? -d.scale * lp : NaN

entropy(d::Exponential) = 1.0 - log(1.0 / d.scale)

insupport(::Exponential, x::Real) = zero(x) <= x < Inf
insupport(::Type{Exponential}, x::Real) = zero(x) <= x < Inf

kurtosis(d::Exponential) = 6.0

mean(d::Exponential) = d.scale

median(d::Exponential) = d.scale * log(2.0)

mgf(d::Exponential, t::Real) = 1.0/(1.0 - t * d.scale)

cf(d::Exponential, t::Real) = (1.0 - t * im * d.scale)^(-1)

mode(d::Exponential) = 0.0
modes(d::Exponential) = [0.0]

pdf(d::Exponential, x::Real) = x < zero(x) ? 0.0 : exp(-x / d.scale) / d.scale

logpdf(d::Exponential, x::Real) =  x < zero(x) ? -Inf : -x / d.scale - log(d.scale)

quantile(d::Exponential, p::Real) = zero(p) <= p <= one(p) ? -d.scale * log1p(-p) : NaN

cquantile(d::Exponential, p::Real) = zero(p) <= p <= one(p) ? -d.scale * log(p) : NaN

rand(d::Exponential) = d.scale * Base.Random.randmtzig_exprnd()

function rand!(d::Exponential, A::Array{Float64})
    Base.Random.randmtzig_fill_exprnd!(A)
    for i in 1:length(A)
        A[i] *= d.scale
    end
    A
end

skewness(d::Exponential) = 2.0

var(d::Exponential) = d.scale * d.scale

## Fit model

immutable ExponentialStats <: SufficientStats
    sx::Float64   # (weighted) sum of x
    sw::Float64   # sum of sample weights

    ExponentialStats(sx::Real, sw::Real) = new(float64(sx), float64(sw))
end

suffstats(::Type{Exponential}, x::Array) = ExponentialStats(sum(x), length(x))

suffstats(::Type{Exponential}, x::Array, w::Array) = ExponentialStats(dot(x, w), sum(w))

fit_mle(::Type{Exponential}, ss::ExponentialStats) = Exponential(ss.sx / ss.sw)


text:
function std_checkbox()
    map(x->map(RGBA{U8}, x),
        (loadasset("checked.png"), loadasset("unchecked.png"))
    )
end

function widget{T<:Union{Bool, UInt8}}(tick::Signal{T}, window;
        checkbox=std_checkbox(), area=(50, 50),
        kw_args...
    )
    icon = map(tick) do t
        checkbox[Bool(t) ? 1 : 2]
    end
    robj = visualize(icon; primitive=SimpleRectangle(0,0,area[2],area[2]), kw_args...)
    robj, togl = toggle_button(tick, robj, window)
    togl, robj
end

function choice_visual(x::Union{Function, DataType, Enum}; kw_args...)
    visualize(string(x); kw_args...)
end
function choice_visual(x; kw_args...)
    visualize(x; kw_args...)
end

function widget{T<:Enum}(enum::Signal{T}, window;
        text_scale=Vec2f0(1), kw_args...
    )
    all_enums = collect(instances(T))
    i0 = findfirst(x->x==value(enum), all_enums)
    choice_widget(all_enums, window; start_idx=i0, text_scale=text_scale, kw_args...)
end

function myscale!(robj, target)
    bb = value(boundingbox(robj))
    w1, w2, w3 = widths(bb)
    w = Vec2f0(w1, w2)
    wt1, wt2, wt3 = widths(target)
    wt = Vec2f0(wt1, wt2)
    t = translationmatrix(-minimum(bb))
    t2 = translationmatrix(minimum(target))
    sm = minimum((1f0./w) .* wt)
    s = scalematrix(Vec3f0(sm, sm, 1))
    m = t2*s*t
    set_arg!(robj, :model, m)
    set_arg!(robj, :boundingbox, Signal(m*bb))
end

function choice_widget(choices::AbstractVector, window;
        text_scale=Vec2f0(1), start_idx=1, area=(150, 30), kw_args...
    )
    show_area = SimpleRectangle(0,0, area...)
    vizzes = map(enumerate(choices)) do i_c
        i, c = i_c
        vis = choice_visual(c;
            relative_scale=text_scale,
            preferred_camera=:fixed_pixel,
            kw_args...
        )
        set_arg!(vis, :visible, i==start_idx)
        myscale!(vis, AABB{Float32}(Vec3f0(2,2,0), Vec3f0(area..., 1)-Vec3f0(2,2,0)))
        vis
    end
    vis = visualize(show_area, color=RGBA{Float32}(0.95, 0.95, 0.95, 0.4)).children[]
    visual = Context(vis, vizzes...)
    selected = foldp(start_idx, toggle(visual, window)) do i0, _
        idx = mod1(i0+1, length(vizzes))
        for i=1:length(vizzes)
            set_arg!(vizzes[i], :visible, i==idx)
        end
        idx
    end
    map(getindex, Signal(choices), selected, typ=eltype(choices)), visual
end

ext: jl


text:
## Scanner

type Scanner
    string::AbstractString
    tail::AbstractString
    pos::Integer
end
Scanner(string::AbstractString) = Scanner("", string, 0)

## Returns `true` if the tail is empty (end of string).
function eos(s::Scanner)
    s.tail == ""
end

## Tries to match the given regular expression at the current position.
## Returns the matched text if it can match, the empty string otherwise.
function scan(s::Scanner, re::Regex)
    if !ismatch(re, s.tail)
        return ""
    end

    m = match(re, s.tail)
    if m.offset >= 1
        ## move past match
        no_chars = endof(m.match) + m.offset - 1
        s.pos += no_chars
        s.tail = s.tail[(no_chars + 1):end]
    end

    m.match
end
scan(s::Scanner, re::AbstractString) = scan(s, Regex(re))
scan(s::Scanner, re::Char) = scan(s, string(re))

## Skips all text until the given regular expression can be matched. Returns
## the skipped string, which is the entire tail if no match can be made.
function scanUntil!(s::Scanner, re::Regex)
    m = match(re, s.tail)

    if m == nothing
        ourmatch = s.tail
        s.pos += endof(s.tail)
        s.tail = ""
    else
        pos = m.offset
        ourmatch = s.tail[1:(pos-1)]
        s.tail = s.tail[pos:(endof(s.tail))]
        s.pos += pos
    end

    return(ourmatch)
end

scanUntil!(s::Scanner, re::AbstractString) = scanUntil!(s, Regex(re))
scanUntil!(s::Scanner, re::Char) = scanUntil!(s, string(re))


text:
### Signal Processing

import Base: fft, ifft, fft!, ifft!, conv, conv2

# Export Methods

export  fftC2R,
        fftR2C,
        conv3,
        convolve,
        fir,
        iir,
        approx1,
        approx2

# Export Constants

export  AF_CONV_DEFAULT,
        AF_CONV_EXPAND,
        AF_CONV_AUTO,
        AF_CONV_SPATIAL,
        AF_CONV_FREQ

# Constants

const AF_CONV_DEFAULT = 0
const AF_CONV_EXPAND = 1

const AF_CONV_AUTO = 0
const AF_CONV_SPATIAL = 1
const AF_CONV_FREQ = 2

# FFTs

for (op, fn) in ((:fft, :af_fft), (:ifft, :af_ifft))

    @eval function ($op){T<:Real}(a::AFVector{T}; norm_factor = 1., dim1 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1)
        AFArray{Complex{T}}(out[])
    end

    @eval function ($op){T<:Complex}(a::AFVector{T}; norm_factor = 1., dim1 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1)
        AFArray{T}(out[])
    end

end

for (op, fn) in ((:fft, :af_fft2), (:ifft, :af_ifft2))

    @eval function ($op){T<:Real}(a::AFMatrix{T}; norm_factor = 1., dim1 = 0, dim2 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1, dim2)
        AFArray{Complex{T}}(out[])
    end

    @eval function ($op){T<:Complex}(a::AFMatrix{T}; norm_factor = 1., dim1 = 0, dim2 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1, dim2)
        AFArray{T}(out[])
    end

end

for (op, fn) in ((:fft, :af_fft3), (:ifft, :af_ifft3))

    @eval function ($op){T<:Real}(a::AFArray{T,3}; norm_factor = 1., dim1 = 0, dim2 = 0, dim3 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1, dim2, dim3)
        AFArray{Complex{T}}(out[])
    end

    @eval function ($op){T<:Complex}(a::AFArray{T,3}; norm_factor = 1., dim1 = 0, dim2 = 0, dim3 = 0)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, dim1, dim2, dim3)
        AFArray{T}(out[])
    end

end

for (op, fn) in ((:fft, :af_fft), (:ifft, :af_ifft))

    @eval function ($op){T<:Real}(a::AFArray{T}, dim::Integer)
        out = new_ptr()
        if dim == 1
            eval($fn)(out, a, 1., 0)
        else
            throw("This dimension is not currently supported")
        end
        AFArray{Complex{T}}(out[])
    end

    @eval function ($op){T<:Complex}(a::AFArray{T}, dim::Integer)
        out = new_ptr()
        if dim == 1
            eval($fn)(out, a, 1., 0)
        else
            throw("This dimension is not currently supported")
        end
        AFArray{T}(out[])
    end

end

for (op, fn) in ((:fft!, :af_fft_inplace), (:ifft!, :af_ifft_inplace))

    @eval function ($op){T<:Real}(a::AFVector{T}; norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

    @eval function ($op){T<:Complex}(a::AFVector{T}; norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

end

for (op, fn) in ((:fft!, :af_fft2_inplace), (:ifft!, :af_ifft2_inplace))

    @eval function ($op){T<:Real}(a::AFMatrix{T}, norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

    @eval function ($op){T<:Complex}(a::AFMatrix{T}, norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

end

for (op, fn) in ((:fft!, :af_fft3_inplace), (:ifft!, :af_ifft3_inplace))

    @eval function ($op){T<:Real}(a::AFArray{T,3}, norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

    @eval function ($op){T<:Complex}(a::AFArray{T,3}, norm_factor = 1.)
        eval($fn)(a, norm_factor)
        a
    end

end

for (arr, fn) in ((:(AFVector{Complex{T}}), :af_fft_c2r), (:(AFMatrix{Complex{T}}), :af_fft2_c2r),
                    (:(AFArray{Complex{T},3}), :af_fft3_c2r))

    @eval function fftC2R{T<:Real}(a::($arr); is_odd = false, norm_factor = 0.)
        out = new_ptr()
        eval($fn)(out, a, norm_factor, is_odd)
        AFArray{T}(out[])
    end

end

function fftR2C{T}(a::AFVector{T}, pad1::Integer; norm_factor = 0.)
    out = new_ptr()
    af_fft_r2c(out, a, norm_factor, pad1)
    AFArray{Complex{T}}(out[])
end

function fftR2C{T}(a::AFMatrix{T}, pad1::Integer, pad2::Integer; norm_factor = 0.)
    out = new_ptr()
    af_fft2_r2c(out, a, norm_factor, pad1)
    AFArray{Complex{T}}(out[])
end

function fftR2C{T}(a::AFArray{T,3}, pad1::Integer, pad2::Integer, pad3::Integer; norm_factor = 0.)
    out = new_ptr()
    af_fft3_r2c(out, a, norm_factor, pad1)
    AFArray{Complex{T}}(out[])
end

# Convolutions

for (op, arr1, arr2, fn) in ((:conv, :(sig::AFVector{T}), :(fil::AFVector{S}), :af_fft_convolve),
                        (:conv2, :(sig::AFMatrix{T}), :(fil::AFMatrix{S}), :af_fft_convolve2),
                        (:conv3, :(sig::AFArray{T,3}), :(fil::AFArray{S,3}), :af_fft_convolve3))

    @eval function ($op){T,S}($arr1, $arr2; mode = AF_CONV_EXPAND)
        out = new_ptr()
        eval($fn)(out, sig, fil, mode)
        AFArray{af_promote(T,S)}(out[])
    end

end

for (arr1, arr2, fn) in ((:(sig::AFVector{T}), :(fil::AFVector{S}), :af_convolve1),
                            (:(sig::AFMatrix{T}), :(fil::AFMatrix{S}), :af_convolve2),
                            (:(sig::AFArray{T,3}), :(fil::AFArray{S,3}), :af_convolve3))

    @eval function convolve{T,S}($arr1, $arr2; mode = AF_CONV_EXPAND, domain = AF_CONV_AUTO)
        out = new_ptr()
        eval($fn)(out, sig, fil, mode, domain)
        AFArray{af_promote(T,S)}(out[])
    end

end

# Filters

function fir(b::AFArray, x::AFArray)
    out = new_ptr()
    af_fir(out, b, x)
    AFArray{backend_eltype(out[])}(out[])
end

function iir(b::AFArray, a::AFArray, x::AFArray)
    out = new_ptr()
    af_iir(out, b, a, x)
    AFArray{backend_eltype(out[])}(out[])
end

# Approx

function approx1(a::AFArray, pos::AFArray; method = AF_INTERP_LINEAR, offGrid = 0.0)
    out = new_ptr()
    af_approx1(out, a, pos, method, offGrid)
    AFArray{backend_eltype{out[]}}(out[])
end

function approx2(a::AFArray, pos1::AFArray, pos2::AFArray; method = AF_INTERP_LINEAR, offGrid = 0.0)
    out = new_ptr()
    af_approx2(out, a, pos1, pos1, method, offGrid)
    AFArray{backend_eltype{out[]}}(out[])
end


text:
#!/usr/bin/env julia --project
using Test
using ClimateMachine
ClimateMachine.init()
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.Ocean.ShallowWater
using ClimateMachine.Ocean.OceanProblems

using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.BalanceLaws: vars_state, Prognostic, Auxiliary
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.VTK

using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function run_hydrostatic_spindown(; refDat = ())
    mpicomm = MPI.COMM_WORLD
    ArrayType = ClimateMachine.array_type()

    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel = ll == "DEBUG" ? Logging.Debug :
        ll == "WARN" ? Logging.Warn :
        ll == "ERROR" ? Logging.Error : Logging.Info
    logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    brickrange_2D = (xrange, yrange)
    topl_2D = BrickTopology(
        mpicomm,
        brickrange_2D,
        periodicity = (true, true),
        boundary = ((0, 0), (0, 0)),
    )
    grid_2D = DiscontinuousSpectralElementGrid(
        topl_2D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    problem = SimpleBox{FT}(Lˣ, Lʸ, H)

    model_2D = ShallowWaterModel{FT}(
        param_set,
        problem,
        ShallowWater.ConstantViscosity{FT}(5e3),
        nothing;
        c = FT(1),
        fₒ = FT(0),
        β = FT(0),
    )

    dt_fast = 300
    nout = ceil(Int64, tout / dt_fast)
    dt_fast = tout / nout

    dg_2D = DGModel(
        model_2D,
        grid_2D,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q_2D = init_ode_state(dg_2D, FT(0); init_on_cpu = true)

    lsrk_2D = LSRK54CarpenterKennedy(dg_2D, Q_2D, dt = dt_fast, t0 = 0)

    odesolver = lsrk_2D

    vtkstep = [0, 0]
    cbvector = make_callbacks(
        vtkpath,
        vtkstep,
        nout,
        mpicomm,
        odesolver,
        dg_2D,
        model_2D,
        Q_2D,
    )

    eng0 = norm(Q_2D)
    @info @sprintf """Starting
    norm(Q₀) = %.16e
    ArrayType = %s""" eng0 ArrayType

    solve!(Q_2D, odesolver; timeend = timeend, callbacks = cbvector)

    Qe_2D = init_ode_state(dg_2D, timeend, init_on_cpu = true)

    error_2D = euclidean_distance(Q_2D, Qe_2D) / norm(Qe_2D)

    println("2D error = ", error_2D)
    @test isapprox(error_2D, FT(0.0); atol = 0.005)

    ## Check results against reference
    ClimateMachine.StateCheck.scprintref(cbvector[end])
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(cbvector[end], refDat)
    end

    return nothing
end

function make_callbacks(
    vtkpath,
    vtkstep,
    nout,
    mpicomm,
    odesolver,
    dg_fast,
    model_fast,
    Q_fast,
)
    if isdir(vtkpath)
        rm(vtkpath, recursive = true)
    end
    mkpath(vtkpath)
    mkpath(vtkpath * "/fast")

    function do_output(span, vtkstep, model, dg, Q)
        outprefix = @sprintf(
            "%s/%s/mpirank%04d_step%04d",
            vtkpath,
            span,
            MPI.Comm_rank(mpicomm),
            vtkstep
        )
        @info "doing VTK output" outprefix
        statenames = flattenednames(vars_state(model, Prognostic(), eltype(Q)))
        auxnames = flattenednames(vars_state(model, Auxiliary(), eltype(Q)))
        writevtk(outprefix, Q, dg, statenames, dg.state_auxiliary, auxnames)
    end

    do_output("fast", vtkstep[2], model_fast, dg_fast, Q_fast)
    cbvtk_fast = GenericCallbacks.EveryXSimulationSteps(nout) do (init = false)
        do_output("fast", vtkstep[2], model_fast, dg_fast, Q_fast)
        vtkstep[2] += 1
        nothing
    end

    starttime = Ref(now())
    cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q_fast)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(
                    convert(Dates.DateTime, Dates.now() - starttime[]),
                    Dates.dateformat"HH:MM:SS",
                ),
                energy
            )
        end
    end

    cbcs_dg = ClimateMachine.StateCheck.sccreate(
        [(Q_fast, "2D state")],
        nout;
        prec = 12,
    )

    # don't write vtk during CI testing
    # return (cbvtk_fast, cbinfo, cbcs_dg)
    return (cbinfo, cbcs_dg)
end

#################
# RUN THE TESTS #
#################
FT = Float64
vtkpath = abspath(joinpath(
    ClimateMachine.Settings.output_dir,
    "vtk_shallow_spindown",
))

const timeend = FT(24 * 3600) # s
const tout = FT(2 * 3600) # s
# const timeend = 1200 # s
# const tout = 600 # s

const N = 4
const Nˣ = 5
const Nʸ = 5
const Lˣ = 1e6  # m
const Lʸ = 1e6  # m
const H = 400  # m

xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)

@testset "$(@__FILE__)" begin
    include("../refvals/2D_hydrostatic_spindown_refvals.jl")

    run_hydrostatic_spindown(refDat = refVals.explicit) # error = 0.00011327920483879001
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct Rrda
    fm
    lev::AbstractVector
    ni::AbstractVector
end

struct KrrDa
    fm
    lev::AbstractVector
    ni::AbstractVector
end

"""
    rrda(X, y, weights = ones(size(X, 1)); lb)
Discrimination (DA) based on ridge regression (RR).
* `X` : X-data.
* `y` : y-data (class membership).
* `weights` : Weights of the observations.
* `lb` : A value of the regularization parameter "lambda".

The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a RR is implemented on the X-data and each column of Ydummy,
returning predictions of the dummy variables (= object `posterior` returned by
function `predict`).
These predictions can be considered as unbounded (i.e. eventuall outside of [0, 1])
estimates of the class membership probabilities.
For a given observation, the final prediction is the class corresponding
to the dummy variable for which the prediction is the highest.
"""
function rrda(X, y, weights = ones(size(X, 1)); lb)
    z = dummy(y)
    fm = rr(X, z.Y, weights; lb = lb)
    Rrda(fm, z.lev, z.ni)
end

"""
    krrda(X, y, weights = ones(size(X, 1)); lb)
Discrimination (DA) based on kernel ridge regression (KRR).
* `X` : X-data.
* `y` : Univariate class membership.
* `weights` : Weights of the observations.
* `lb` : A value of the regularization parameter "lambda".
* 'kern' : Type of kernel used to compute the Gram matrices.
    Possible values are "krbf" of "kpol" (see respective functions `krbf` and `kpol`.
* `kwargs` : Named arguments to pass in the kernel function.

The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a KRR is implemented on the X-data and each column of Ydummy,
returning predictions of the dummy variables (= object `posterior` returned by
function `predict`).
These predictions can be considered as unbounded (i.e. eventuall outside of [0, 1])
estimates of the class membership probabilities.
For a given observation, the final prediction is the class corresponding
to the dummy variable for which the prediction is the highest.
"""
function krrda(X, y, weights = ones(size(X, 1)); lb, kern = "krbf", kwargs...)
    z = dummy(y)
    fm = krr(X, z.Y, weights; lb = lb, kern = kern, kwargs...)
    KrrDa(fm, z.lev, z.ni)
end

"""
    predict(object::Union{Rrda, KrrDa}, X; lb = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
* `lb` : Regularization parameter, or collection of regularization parameters,
    "lambda" to consider. If nothing, it is the parameter stored in the
    fitted model.
"""
function predict(object::Union{Rrda, KrrDa}, X; lb = nothing)
    X = ensure_mat(X)
    m = size(X, 1)
    isnothing(lb) ? lb = object.fm.lb : nothing
    le_lb = length(lb)
    pred = list(le_lb, Matrix{Float64})
    posterior = list(le_lb, Matrix{Float64})
    @inbounds for i = 1:le_lb
        zp = predict(object.fm, X; lb = lb[i]).pred
        z =  mapslices(argmax, zp; dims = 2)  # if equal, argmax takes the first
        pred[i] = reshape(replacebylev(z, object.lev), m, 1)
        posterior[i] = zp
    end
    if le_lb == 1
        pred = pred[1]
        posterior = posterior[1]
    end
    (pred = pred, posterior = posterior)
end


text:
module TestUnordered
include("../preamble.jl")
using Transducers: append_unordered!, channel_unordered, reduce_commutative

@testset "reduce_commutative" begin
    @testset "Map(inc)(+)" begin
        @test reduce_commutative(+, Map(inc), 1:100) == sum((1:100) .+ 1)
    end
end

@testset "append_unordered!" begin
    input = Channel(Map(identity), 1:100)
    output = Channel{Int}(Inf)
    append_unordered!(output, Map(inc), input)
    close(output)
    ys = sort!(collect(output))
    @test ys == (1:100) .+ 1
end

@testset "append_unordered!(output, (x + 1 for x in input))" begin
    input = Channel(Map(identity), 1:100)
    output = Channel{Int}(Inf)
    append_unordered!(output, (x + 1 for x in input))
    close(output)
    ys = sort!(collect(output))
    @test ys == (1:100) .+ 1
end

@testset "channel_unordered" begin
    @testset "Map(inc)" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(Map(inc), input)))
        @test ys == (1:100) .+ 1
    end
    @testset "x + 1 for x in input" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(x + 1 for x in input)))
        @test ys == (1:100) .+ 1
    end
    @testset "eduction(x + 1 for x in input)" begin
        input = Channel(Map(identity), 1:100)
        ys = sort!(collect(channel_unordered(eduction(x + 1 for x in input))))
        @test ys == (1:100) .+ 1
    end
    @testset "ReduceIf(isodd): early termination" begin
        input_data = fill(2, 100)
        input_data[10] = 1
        input = Channel(Map(identity), input_data)
        ys = collect(channel_unordered(ReduceIf(isodd), input))
        @test count(x -> x == 1, ys) == 1
        @test sort!(unique(ys)) == [1, 2]
        @test length(ys) < length(input_data)
    end
    @testset "error handling (no dangling tasks)" begin
        input_length = 100
        input = Channel(Map(identity), 1:input_length)
        trace = Channel(input_length)
        ntasks = 10
        block = Channel{Nothing}(1)
        xs = input |> Map() do x
            put!(trace, x)
            wait(block)
            x
        end |> Map() do x
            error("Throwing error with x = $x")
        end
        bg = @async begin
            traced = Int[]
            for _ in 1:ntasks
                push!(traced, popfirst!(trace))
            end
            # Now all tasks completed `put!(trace, x)`.
            put!(block, nothing)
            for x in trace
                push!(traced, x)
            end
            traced
        end
        @test_throws Exception take!(channel_unordered(xs; ntasks=ntasks))
        close(trace)
        consumed = sort!(fetch(bg))
        @test consumed == 1:ntasks
    end
    @testset "error handling (terminate other tasks)" begin
        input_length = 100
        input = Channel(Map(identity), 1:input_length)
        trace = Channel(input_length)
        ntasks = 10
        block = Channel{Nothing}(1)
        xs = input |> Map() do x
            put!(trace, x)
            wait(block)
            x
        end |> Map() do x
            if x == 1
                error("Throwing error with x = $x")
            end
            sleep(0.01)
            x
        end
        bg = @async begin
            traced = Int[]
            for _ in 1:ntasks
                push!(traced, popfirst!(trace))
            end
            # Now all tasks completed `put!(trace, x)`.
            put!(block, nothing)
            for x in trace
                push!(traced, x)
            end
            traced
        end
        @test_throws Exception collect(channel_unordered(xs; ntasks=ntasks))
        close(trace)
        consumed = sort!(fetch(bg))
        @test length(consumed) >= ntasks
        @test length(consumed) < input_length
    end
end

end  # module


text:
# https://github.com/ericphanson/TravelingSalesmanExact.jl
#
using Base.Threads: @spawn, @threads, nthreads, threadid
#using FLoops
using Statistics
using Combinatorics
using Logging
using Printf: @sprintf
using BenchmarkTools
#
logger = ConsoleLogger(stdout, Logging.Info) # Warn or Info or Debug
global_logger(logger)
disable_logging(LogLevel(-100000))
#
@info @sprintf "running julia with %s threads." nthreads()


text:
"""
Credits:

 Tim Holy - JuliaCon 2016

"""
abstract type AbstractMappedArray{T,N} <: AbstractArray{T,N} end
abstract type AbstractMultiMappedArray{T,N} <: AbstractMappedArray{T,N} end

"""

trivial define array structur that does new things for you
"""
#--- is there overhead for this ?
using BenchmarkTools

#No method matching +(::Array{Float64,0}, ::Float64)
#myfun
#0-dimensional sum (1 scalar value)
function mysum(A)
    s = zeros(eltype(A))
    @inbounds for a in A
        s .+= a
    end
    s
end
#=should got 1 Dimensional vector=#
function mysum(A)
    s = zeros(eltype(A))
    @inbounds for a in A
        s += a
    end
    s
end

mysum(typeof(rand(5, 5)))

m = mysum(rand(5, 5))

#meanwhile, emulating large dataset

A = rand(10^4, 10^4)
B = Mappedarray(identity, A) # MappedArray not defined (yet)

@time mysum(A)
@btime mysum(A)
@Benchmarking mysum(A)

@time mysum(B)
@btime mysum(B)
@Benchmarking mysum(B)

#TODO: fix mysum()

"""
julia has the capabilities, that allows you to handle
Data, in a reasonably efficient fashion

Before Julia:

Each experiment needs a 4D (Large) Array

Useful things Julia does to support arrays
1. Memory mapping
can't hold any theese in Ram,

present it as a view in a single memory

transformation s

range compression

16 bit cameras , needs scaling - intensity-scaling:

2. intensity scling
i.e. 16 bit cameras , needs scaling - intensity-scaling:
How:
before Julia
2.1.1 find gui viewer offers scaling and other features you ned
2.1.1 give up

2.2.1 Find blank drive
2.2.2 Format it
2.2.3 Write script saves intensity=scale images to disk (wait serveral hours)
2.2.4 Load the scaled imge into your viewer

That's the Bottleneck, very irritating (in early days)
In Julia:
Compute represent
Virtual image "view" (of underlying data )

imgsc = mappedarray(sqrt, img )

the key here: "we don't actually perform the computation"
we signal that  we want (to perform this computation)
whn perform it, we  signal it into a Dumb Viewer

why Dumb? it doesn't know how to do transformations(i.e. Scaling)
[the kicker: but it might have (some) other features, we want]
why do a dumb viewer? (3 reasons)

1.Easy
2.Simple
3.Composable
(Hype of) workflow

make this operation v. fast (microseconds)

1.main point

imgsc = mappedarray(sqrt, img)
"""

"""
transformations are a
#1: vectorized functions -Heratige of programming languages operates element-wise on the array

amounts: Viewd as  not desirable
    -Explicit mapping function  (over the elements)
    --vectorized functions
"""

"""
Gregory(scotland Yard detective):
Is theree any other point to which you would wish to

which you wish to draw my attention

Holmes: To the curious incident of the dog in the night-time

Gregory: th dog did nothing in the night-time
Holmes: that was the curious incident
- Silver Blaze - conan Doyle

"""

"""

Meaning#1
the best optimization is
not Doing the Computation

Best computation (Strategy)
is by not doing the Computation

useful point to the things tim's lab does
i.e. performing element-wise computations on th array

how should we write these types of transformations?
on option: a Heratige from (serveral) other programming languages
1.the vectorized function : operates element-wise on the Array
that is
1.maybe undesirable, say matrix square root (3D)

2.Explicit mapping
"""

"""
take sqrt when we ask for values from that array

    nice thing is:

    Lazy is usd widely thru computing (in general)

    general arrays b. fit for comp
    as
    because this bracket []

    taking valus from array v. crispyly
    encapsulates now, is the time to enact that computation

    -allows us to i.e
    1  have dataset 10K "with  different nerouns plotting each one" (of them )
    on the scatterplot
    2. Click on single one , take us right to raw data
    3. Might  have tpoint 0.01 of the total data volume
    so, no need to compute Square root (on all pixel intensitiese)
    just compute on the one (I'm zoomed in(?)) , looking at a particular moment )]

    #Defining a new Array Type

    #immu
mapped array isa substype of abstract array type
-has 2 in it

1 Function
2.Original array (we wanting to Create a view of )

array: defines 'a' (which contains ...)
"""

"""
            (Dr.) Tim's lab's favorite is:

include("MappedArrays.jl")
#=Features
1.lazy Initialization
2.(promise-like async) [me: missed javascript programming]
=#
"""

"""
-1 Create a View of (Initial dataset)
-Only take Roots, when ask the values, from that array
things nice:

Click on a Single Pixel-Intensity 'Neuron' Lazily Computes
(Take us to its data)
Might only looking to the "0.01" (1%) percent of the Data
`No need to compute` on all Pixel-Intensities

"""

"""
#Self-made: EXPERIMENTAL - This Wrecks the Workplace Variables - DO NOT RUN!

struct Mappedarray{T,N,A<:AbstractArray,f::Int...} <: AbstractArray{T,N}
    f:Int...
    data::A
end
"""

"""
gentle
function + original array (want to create a view of )
    1.This defines a container
    2.Define that happens when you (might) take a value of it

Apply Function (You store in container) to each element
    (of the data ) as you take it out
   Notice: have to define this for Integers (it's Significant)

"""
#errors out: f is not defined #that' the point of laziness # there's more into that
#introduce MappedArray
struct Mappedarray{T,N,A<:AbstractArray,F} <: AbstractArray{T,N}
    f:F
    data::A
end #f not defined, Prof. Holy
###
#or using an already set
using MappedArrays, Test # ] add MappedArrays

A = rand(5, 5)
T = ones
N= ones
MappedArray{T,N}(f, data::AbstractArray{T,N}) =  #defines f as ones # but N Not defined
    MappedArray{typeof(f(one(T))),N,typeof(data),typeof(f)}(f, data)

B = Mappedarray(sqrt, A) #only realizes computation whn calling B[i,j,..] # no method matching

#---
#sqrt() & map(sqrt,A)
#map approximates sqrt correctly
#Q. how many floating points are necessay? #TODO:

A = rand()
sqrt(A)         #implicit element-wise computation # 0.9328629863877007

map(sqrt, A)   #explicit element-wise computation # 0.9328629863877007

@time sqrt(A)  # 0.000005 seconds (1 allocation: 16 bytes)

@assert sqrt(A) == map(sqrt, A)  #True

#------
@time map(sqrt, A) # 0.000007 seconds (2 allocations: 32 bytes)
@btime map(sqrt, A) #183.176 ns (2 allocations: 32 bytes) ;0.9298476269411314 # 181.304 ns #the BEST
@benchmark map(sqrt, A)
#=
BenchmarkTools.Trial: 10000 samples with 631 evaluations.
 Range (min … max):  183.043 ns …  1.835 μs  ┊ GC (min … max): 0.00% … 86.16%
 Time  (median):     195.563 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   209.635 ns ± 59.076 ns  ┊ GC (mean ± σ):  0.30% ±  1.75%

  █▂▄▇▄▁▅▃▂▂ ▁  ▃▃ ▁▂                                          ▁
  █████████████▇██▇████▇▆▇▇▇▆▆▆▆▇▇▇▇▆▅▆▇▇▆▆▅▆▆▇▅▅▆▆▅▇▆▄▅▄▅▅▅▄▄ █
  183 ns        Histogram: log(frequency) by time       419 ns <

 Memory estimate: 32 bytes, allocs estimate: 2.
 =#

B = mappedarray(sqrt, A) #  either A 's items reach end of program (as expected) #Or   # ERROR: tuple must be non empty
#B = mappedarray(func, a) #anyways, it's generalization

##
##
@test Base.size(A::AbstractMappedArray) = size(A.data)#copies size of original data  #true #ERROR: (UNEXPECTED) during test
###
Base.Base.@propagate_inbounds Base.getindex(A::Mappedarray, i::Int...) =
    A.f(A.data[i...])

"""
Notice F is representation
Julia's Compiler
Knows
There is a Fucntion ( stored in this object)
2.which Function it is
this give:
"""

"""
Base.@propagate_inbounds Base.getindex(A::Mappedarray, i::Int...) =
    A.f(A.data[i...])#function call
"""

"""
Define What happens when you Actually
Want to Take a Value of it right ?
This is The Line that Defines the Bracket Operation

When you Actually Extract values from it

Defined in this way , there are a few v. important
Details
See: Blake

Apply function you store it
(in this container)
To each element of the data
As you take it out
You notice you only (Have) define
This for integers
    I'll point that out (a more)

2nd, you need Define Another Function, which is just
what's the Size of One of these Objects

"""

"""#to Uncomment
@deprecate mappedarray(f_finv::Tuple{Any,Any}, args::AbstractArray...) mappedarray(f_finv[1], f_finv[2], args...)
Mappedarray{T,N}(f, data::AbstractArray{T,N}) = Mappedarray{typeof(f(one(T))),N,typeof(data),typeof(f)}(f, data)

Base.size(A::AbstractMappedArray) = size(A.data)

#size( Original Data ), you hav be abl to construct them
"""
"""
This is the line you need for the constructor

For a Bare-bones Implementation, this is it- that's all you need to do
    to create one of these views that
    Computes arbitrary function, on the Datatransformation
        -static transformation is alomost every where
"Matt Baumann" method: #TODO: Research "Matt Baumann"                 <---------------------------------------
        get index
"""

#--- isn't there some overhead for this

"""benchmark it
Creates function mysum
    over elemts of the array
    adds them to "accumulator"

    """
function mysum(A)
    s = zero(eltype(A))
    @inbounds for a in A
        s += a
    end
    s
end

@time A = rand(10^3, 10^3) #V  0.002637 seconds (2 allocations: 7.629 MiB) 1000×1000 Matrix{Float64}:
@time mysum(A) # 499774.5988719051
@benchmark mysum(A)
#=
  0.003270 seconds (2 allocations: 7.629 MiB)
1000×1000 Matrix{Float64}:

BenchmarkTools.Trial: 3940 samples with 1 evaluation.
 Range (min … max):  1.120 ms …   2.929 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.199 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.260 ms ± 163.289 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁  ▄█▇ ▁
  █▇██████▆▆▅▃▃▂▃▂▂▂▂▂▂▁▂▁▂▁▂▂▂▂▂▂▂▂▂▁▁▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  1.12 ms         Histogram: frequency by time        1.89 ms <

 Memory estimate: 16 bytes, allocs estimate: 1.
=#
"""benchmark output mysum(A)
B = Mappedarray(identity,A)
  0.003397 seconds (2 allocations: 7.629 MiB)
  0.010455 seconds (8.86 k allocations: 494.846 KiB,
      86.97% compilation time)500164.74558258074

"""
#
"""lesson learned
@inbounds macro oneach start of for loop
"""

"""[{
	"resource": "./DeepLearner/DeepLearner/src/DataStructures/Arrays/Array.jl",
	"owner": "_generated_diagnostic_collection_name_#1",
	"severity": 2,
	"message": "Cannot declare constant; it already has a value.",
	"source": "Julia",
	"startLineNumber": 347,
	"startColumn": 1,
	"endLineNumber": 350,
	"endColumn": 4
}]
#ERROR: MethodError: no method matching Mappedarray(::typeof(identity),::Matrix{Float64})
#--- Q.why is there no overhead?

##F = Base.identity #you can't pre-define it (because it throws error "Cannot Declare a Constant, already has a value")
#Cannot declare constant, already has a value

struct Mappedarray{T,N,A<:AbstractArray,F} <: AbstractArray{T,N}
    f::F
    data::A
end
"""
#set op: Assigning @propagate_inbounds #try-out, on every array Op (Please)
@time Base.@propagate_inbounds Base.getindex(A::Mappedarray, i::Int...) =
    A.f(A.data[i...])

#=  0.008769 seconds (47 allocations: 3.448 KiB, 61.02% compilation time)
    0.005372 seconds (47 allocations: 3.448 KiB, 50.32% compilation time)

    var seconds (const alloc = 47: 3.448 KiB var %compilation time)

Infer:allocations is the only Reliable indicator (across machines)
    [not time(seconds), not %compiliation]
=#

@time Base.getindex(A::Mappedarray, i::Int...) =
    A.f(A.data[i...]) # 0.000115 seconds (22 allocations: 1.415 KiB)

"""Lesson Learned
Base.@propagate_inbounds (is (called) an Inline Function Inline)
Tells the Compiler to Inline a Function while Retaining the caller's
  'Inbounds Context'
"""

"""Question essentially for all purposes,
Compiler knows that f(x) = x

SO It Doesn't Even Bother to make a function call #the whole Point
and returns Value of it immediately (that is how compiler Optimizes)
"""

#--- move as much analysis computation
"""
To Compile-time
So there is `Less work for the CPU` (at runtime)

    -This takes away the work, from the Cpu
    "Thee main trick of julia" ( low & Behold?!)
    """

"""Array Views in Julia
View is still Tied to original (line)data
Allow sth new values in the data too
(i.e. just find the `inverse function`) (me: if !isa N/A )
mappedarray: a-> f(a) for a in A

SubArray (has a Parentarray)

slicing out a Particular `Rectangular Region` #need-for-Visualization #makie
in work:
Variant of subarray (span?)

3.slicePlane
(span?])
pros
Great for 1 off viewers (10 lines of code)

Opportunity: Great for dimensionality reduction (i.e)
So no need to learn about (deal with) 4D arrays

Him:doing imaging in a `weird tilts
(Selcting place of 2D array )

4. ReshapedArray
-Preservees dimensionality
+ adds dimantionality to 3D objects

5. permutedDimsArray
view, traditional julia Column major
(v is transpose of this nature)

Don't write any code that Depends on `Internal Representation`

Datastructures present in julia

Mappedarray types:

1 SubArray : (Parentarray,selectedRegion)

2. ReshapedArray: a x b = n x m
(#there must be some linear non-linaer relation between domain, co-domain)
3. permutedDimsArray

What's not in julia (made )
4.SlicedPlane: (Parentarray, selectedRegion)

Why: it's very easy for you to create your own data Types

---------
Now MappedArray is a Scalar Transformation: A  = a*M ; (a isa scalar )
#--- implementing permutedDimsArray (added in 2016)

Concept behind impleementation: P (the view) reindexes as the parent A as
#To reindex Parent Array
        P[inds...] -> A[inds[iperm]...]

#Anytime it Passes a set of Indicies to the parent array (Permutated View)

Struct permutedDimsArray{T,N,AA <: AbstractArray} <: AbstractArray{T,N}
    parent::AA
    iperm:: Vector{Int} #storing inverse permutation
    dims::NTuple{N, Int}
end

"""
#compiles (ok)
struct permutedDimsArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    iperm::Vector{Int} #storing inverse permutation (as vectors of integers)
    dims::NTuple{N,Int}
end
"""
ERROR:
"""

#T not defined: me: so is it only defined on runtime? (as it's lazy?) - at runtime but in debugging, suffice to give it a try
#T = Int64
 N = 2
Base.@propagate_inbounds function Base.getindex{T,N}(P::PermutedDimsArray{T,N}, I::Vararg{T,N})
    getindex(P.parent, I[P.iperm]...)
end
#ERROR: in Type{...} expression, expected UnionAll, got a value of type typeof(getindex)

"""Lesson Learned:
Async Method on any Array ( via index)
    Invoke @propagate_inbounds onto any function (Array)

        getindex -then-> @propagate_inbounds

1.Apply Permutation to Indicies (being Supplied here)
2.Reference them to Parent Array

"""

#performance check

@time sum(A)#  0.001322 seconds (1 allocation: 16 bytes) 499814.9740729731

@time sum(P) # in demo: parent is undefined #ways of accumulating,it isa not mentioned

"""Point permutedDimsArray
(PermutedView , of that array) is 400 times Slower
(v. unhappy with performance)
How do you `go about` getting a better performance?

"""
#--- (Do) An impleementation based on tuples

""" tuples are the trick behind Almost Everything  #Tonight's Takeaway
(that you do, with an array)

"""
# v2 compiles
"""
PermutedDimsArray

Example: #Do-with-me

  julia> A = rand(3,5,4) #a Humongous Array

  julia> B = PermutedDimsArray(A, (3,1,2)) # a partof the Original Array
"""
#First of all, Define PermutedDimsArray
struct PermutedDimsArray{T,N,AA<:AbstractArray,perm,iperm} <: AbstractArray{T,N}
    parent::AA
    dims::NTuple{N,Int}

end

"""
Moved inverse permutation iperm
Into the type parameters
(me: that is an old C trick , inserting functions into parameters)

The Argument Issue:
1. you Can't put an array in there  (but at the same time )
2. you CAN Put a Tuple (which fastens the Computation Speed, drastically)

storing both `prm` , perminve (me:acutally, just passing (physically stored in Memory: is it? yeah))

so( that), one is `predictable`` from the other
but, essentially, (It) makes your life more convinient (a bit)
(in `Certain circumstances` to have both) #need-to-have
"""

"""
General block: for Generalization & abstraction only - Do not usE THIS code Block -  Compile the next one, please
"""
#General case: don't compile this block
struct structname end #always do this (in OOP Strategy - write (abstract) struct name then end ) #me: what if i have

Base.@propagate_inbounds function (base.getindex{}(A::structname, I::Vararg{typeof,N})
getindex(A.parent, ntuple(d -> I[[]], Val{N})...) #TODO:how to make that line in unicod
)#N above for (Rational?/IRrattional? #numberTheory (numberTheory req.)): Polynomial functions
    #for Real Analysis ok, but
    #if Complex Analysis required, then [ENTER YOUR METHOD HERE] IS REQUIRED - lvl 3
end

"""
Reindexing, at Level of Tuples

"""

#----works
Base.@propagate_inbounds _getindex(i, A, As...) = (A[i], _getindex(i, As...)...)
_getindex(i) = ()
#genetic function with 2 methods
Base.@propagate_inbounds function _setindex!(as::As, vals::Vs, inds::Vararg{Int,N}) where {As,Vs,N}
    a1, atail = as[1], Base.tail(as)
    v1, vtail = vals[1], Base.tail(vals)
    a1[inds...] = v1
    return _setindex!(atail, vtail, inds...)
end
_setindex!(as::Tuple{}, vals::Tuple{}, inds::Vararg{Int,N}) where {N} = nothing

"""
Testing Area
"""
#---testing

m = AA = 1;
n = N = 10;
a = ones(N, m)
b = ones(m, N)
function testvalue(data)
    if !isempty(data)
        first(data)
    else
        zero(eltype(data))
    end::eltype(data)
end
#--- Juliacon 2016

#=
by Professor David Sanders
=#

#=
extending arrays
how: by using the push (Bang)
=#
V = []
w = []

push!(w, "adamus")
append!(w, [3, 4.5, "coffee costs"])

push!(V, 1)
#=
Exclamation mark: is convention that says:
1 argument (or more ) is gonna be modified
this case, it's V
=#
#- hcat - horizontal cat
M = @inbounds hcat(V[1], V[2]) # for all elements( hardcoded here for 2 item) #true positive error @even with inbounds, alerts with error
#vector to Array: horizontal concat hcat: horizontally contatenated vertical vectors
M = hcat(V...) #...:the splat: for all values of container  ; take each element in array ; put it as an argument (concatenate it)
#splat ...: takes each element, in the Iterable, & make Another Argument (of the function)

#toget it back (as vector of vectors), transpose it

M = hcat(V...)'
#1x1 adjoint matrix eltype Int64  - example
#or we directly can write vecto as

v = Vector{Int}[[3, 4], [5, 6]] #2-element Vector{Vector{Int64}}:
v = [[3, 4], [5, 6]] #2-element Vector{Vector{Int64}}:
#=2-element Vector{Vector{Int64}}:
 [3, 4]
 [5, 6]
 Splat warning: Careful when splatting large objects
 (objects with many elements )
 e.g.
 1. concat things into 1 huge vector
2. then, do a reshape

A view: is a reshape (instead of copy, in Julia)

=#

"""higher-dimensional arrays have a fixed size!
-Instead of having a fixed size structure
we can Dynamically allocate vectors
How: by vector of vectors
if we want a dynamic array that accets

"""
w = [] # enmpty array, 0 element , of type {Any,1}
v = Vector{Int}[]
push!(v, [3, 4])
push!(v, [5, 6])

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function pendulum(numbers)
    sorted = sort(numbers)
    result = similar(sorted)
    destination = (1+lastindex(result))÷2
    step = 1
    for v in sorted
        result[destination] = v
        destination += step
        step = step < 0 ? -step+1 : -step-1
    end
    result
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Pkg; Pkg.activate(pwd())
using Revise, BenchmarkTools, Printf, ProgressMeter, Dates
using Random, StatsBase, Distributions, Distances, IterTools, DataFrames, DataFramesMeta
using StaticArrays, SparseArrays

#################################
# Build AirDefenseGame and GameSet

# Air defense game parameters
struct AirDefenseGame
    ncity::Int64 # number of cities
    X::Array{Float64,2} # city coordinates
    v::Array{Float64} # city populations (millions of people)
    r::Float64 # coverage radius of each air defense
    ndp::Int64 # number of air defense systems
    nap::Int64 # number of air attack targets
    ndpdc::Int64 # number of cyber-defendable defense systems
    ndpac::Int64 # number of cyber-attackable defense systems
    ndc::Int64 # number of cyber defense teams
    nac::Int64 # number of cyber attack teams
    pp::Float64 # physical defense effectiveness
    psc::Float64 # cyber sensor effectiveness
    pdc::Float64 # cyber defense effectiveness
    pac::Float64 # cyber attack effectiveness
    iads::Array{Int64} # sorted list of cities with IADS
end

# Generic game structure (except sensors and iadsdefenses pending refactoring)
struct GameSet
    players::UnitRange{Int64} # set of player indices
    A::Vector{Vector{Vector{Int64}}} # actions by stage
    An::Vector{UnitRange{Int64}} # action indices by stage
    na_stage::Vector{Int64} # number of actions by stage
    Iset::Vector{Any} # infosets
    ni_stage::Vector{Int64} # number of infosets by stage
    ni::Int64 # number of total infosets
    C::Array{Int64, 2} # coverage matrix
    sensors::Vector{Vector{Int64}} # combinations of sensed attacks
    iadsdefenses::Vector{Vector{Int64}} # combinations of defended iads
end

# data structure for optimization problem
struct GameOptSet
    ns1::Int64 # number of player 1 sequences
    ns2::Int64 # number of player 2 sequences
    ni1::Int64 # number of player 1 infosets
    ni2::Int64 # number of player 2 infosets
    ns1_stage::Vector{Int64} # number of player 1 sequences by stage
    ns2_stage::Vector{Int64} # number of player 2 sequences by stage
    reward # player 1 utility by seq1, seq2 (default is sparse array)
    na_stage::Vector{Int64} # number of actions by stage
    sensors::Vector{Vector{Int64}} # combinations of sensed attacks
    seqI1::Vector{Int64} # player 1: info set to sequence arriving at info set
    seqI2::Vector{Int64} # player 2: info set to sequence arriving at info set
    nextseqI1::Vector{UnitRange{Int64}} # player 1: info set to all sequences available from infoset
    nextseqI2::Vector{UnitRange{Int64}} # player 2: info set to all sequences available from infoset
end

AirDefenseGame() = AirDefenseGame(10,
        reshape([0.03674972063074833,0.7471235124528937,0.8006182686839263,0.1134119509104532,0.3720810993450727,0.6452465195265953,0.9927983579441313,0.8545575129265128,0.29473362663716407,0.9808676508834953,0.6242283949611096,0.6173536676330247,0.4739652533798455,0.6677615345281047,0.6338233451643072,0.8231268712815414,0.011816562466650193,0.048514560425303443,0.7974455234366102,0.22065505430418475], 10, 2),
        [0.47973468500756566,0.1596483573101528,0.7116043145363229,0.969005005082471,0.7139100481289631,0.4747849042596852,0.5066187533113924,0.7074385618615113,0.21009348198668265,0.29520436695319363],
        0.3,
        6,
        3,
        4,
        4,
        2,
        2,
        0.9,
        0.8,
        0.7,
        0.8,
        [1, 3, 4, 5, 8, 9])

AirDefenseGame(ncity, ndp, nap, ndpdc, ndpac, ndc, nac) = AirDefenseGame(ncity,
        rand(ncity, 2), rand(ncity), 0.3, ndp, nap, ndpdc, ndpac, ndc, nac,
        0.9, 0.8, 0.7, 0.8, sort!(sample(1:ncity, ndp, replace = false)))

AirDefenseGame(ncity, v, ndp, nap, ndpdc, ndpac, ndc, nac) = AirDefenseGame(ncity,
        rand(ncity, 2), v, 0.3, ndp, nap, ndpdc, ndpac, ndc, nac,
        0.9, 0.8, 0.7, 0.8, sort!(sample(1:ncity, ndp, replace = false)))

# AirDefenseGame(psc) = AirDefenseGame(10,
#         rand(10, 2), rand(10), 0.3, 6, 3, 4, 4, 2, 2,
#         0.9, psc, 0.7, 0.8, sort!(sample(1:10, 6, replace = false)))

AirDefenseGame(psc) = AirDefenseGame(10,
        reshape([0.03674972063074833,0.7471235124528937,0.8006182686839263,0.1134119509104532,0.3720810993450727,0.6452465195265953,0.9927983579441313,0.8545575129265128,0.29473362663716407,0.9808676508834953,0.6242283949611096,0.6173536676330247,0.4739652533798455,0.6677615345281047,0.6338233451643072,0.8231268712815414,0.011816562466650193,0.048514560425303443,0.7974455234366102,0.22065505430418475], 10, 2),
        [0.47973468500756566,0.1596483573101528,0.7116043145363229,0.969005005082471,0.7139100481289631,0.4747849042596852,0.5066187533113924,0.7074385618615113,0.21009348198668265,0.29520436695319363],
        0.3,
        6,
        3,
        4,
        4,
        2,
        2,
        0.9,
        psc,
        0.7,
        0.8,
        [1, 3, 4, 5, 8, 9])

function AirDefenseGame(ncity, X, v, r, ndp, nap, ndpdc, ndpac, ndc, nac,
        pp, psc, pdc, pac, iads)
    if any((pp, psc, pdc, pac) .< 0)
        @warn "Negative probabilities. Returning standard AirDefenseGame."
        return AirDefenseGame()
    elseif any((pp, psc, pdc, pac) .> 1)
        @warn "Probability > 1. Returning standard AirDefenseGame."
        return AirDefenseGame()
    elseif any((ndp, nap) .> ncity)
        @warn "Too many ndp or nap. Check ncity. Returning standard AirDefenseGame."
        return AirDefenseGame()
    elseif iads != sort(iads)
        @warn "IADS list must be sorted in increasing order. Returning standard AirDefenseGame."
        return AirDefenseGame()
    else
        return AirDefenseGame(ncity, X, v, r, ndp, nap, ndpdc, ndpac, ndc, nac,
                pp, psc, pdc, pac, iads)
    end
end

# ```
# coverage(g::AirDefenseGame)
#
# Compute the coverage of air defenses located in certain cities. C[i,j] is 1
# if city j is covered by defender i and 0 otherwise.
# ```
function coverage(g::AirDefenseGame)
    X, r, ncity = g.X, g.r, g.ncity
    dist = pairwise(Euclidean(), X')
    C = [dist[i,j] < r ? 1 : 0 for i = 1:ncity, j = 1:ncity]
    C
end

function build_action(g::AirDefenseGame)
    ncity, ndp, nap, ndpdc, ndpac, ndc, nac = g.ncity, g.ndp, g.nap, g.ndpdc, g.ndpac, g.ndc, g.nac
    a1c = collect(subsets(1:ndp, ndpdc)) # defensive cyber defendable na1c = binomial(ndp, ndpdc)
    a2c = collect(subsets(1:ndp, ndpac)) # offensive cyber attackable na2c = binomial(ndp, ndpac)
    a3a = collect(subsets(1:ndpac, nac)) # cyber attacks na3a = binomial(ndpac, nac)
    a4c = reshape_product(collect(Iterators.product(fill(1:2, nac)...))) # defensive cyber SA na4c = 2 ^ nac
    a5d = collect(subsets(1:ndpdc, ndc)) # defensive cyber response na5d = binomial(ndpdc, ndc)
    a6a = collect(subsets(1:ncity, nap)) # physical attcks  na6a = binomial(ncity, nap)
    A = [a1c, a2c, a3a, a4c, a5d, a6a]
    na_stage = length.(A)
    An = [1:na for na in na_stage]
    A, An, na_stage
end

function build_info(g::AirDefenseGame)
    ncity, ndp, nap, ndpdc, ndpac, ndc, nac = g.ncity, g.ndp, g.nap, g.ndpdc, g.ndpac, g.ndc, g.nac
    ni3a = binomial(ndp, ndpac)
    nsensor = sum(binomial(ndp, i) for i = 0:nac)
    ni5d = binomial(ndp, ndpdc) * nsensor
    ni6a = ni3a * binomial(ndpac, nac) * binomial(ndp, ndc)
    ni_stage = [0, 0, ni3a, 0, ni5d, ni6a] # ni_stage[j] is num of infosets at stage j
    ni = sum(ni_stage) # total number of infosets
    ni1 = ni5d
    ni2 = ni3a + ni6a
    ni, ni1, ni2, ni_stage
end

# the same physical action at a different info set has a different action name (see ref in kuhnopt.jl)
function build_nseq(g::AirDefenseGame, na_stage, ni_stage)
    ns1_stage = [1, ni_stage[5] * na_stage[5]]
    ns2_stage = [1, ni_stage[3] * na_stage[3], ni_stage[6] * na_stage[6]]
    ns1, ns2 = sum(ns1_stage), sum(ns2_stage)
    ns1, ns2, ns1_stage, ns2_stage
end

function build_Σset(g::AirDefenseGame, A, ns1, ns2)
    Σ1 = Vector{String}(undef, ns1)
    Σ2 = Vector{String}(undef, ns2)
    seqn1 = Array{Int64, 2}(undef, ns1, length(A))
    seqn2 = Array{Int64, 2}(undef, ns2, length(A))
    seq1defended_iads = Array{Int64, 2}(undef, ns1, g.ndc)
    seq2attacked_iads = Array{Int64, 2}(undef, ns2, g.nac)
    seq2defended_iads = Array{Int64, 2}(undef, ns2, g.ndc)
    seq2attacked_cities = Array{Int64, 2}(undef, ns2, g.nap)
    sensors = Vector{Int64}[]
    for i = 0:g.nac
        push!(sensors, collect(subsets(1:g.ndp, i))...) # by index
    end
    iadsdefenses = collect(subsets(1:g.ndp, g.ndc))
    Σ1[1], Σ2[1] = "∅", "∅"
    seqn1[1,:], seqn2[1,:] = zeros(Int64, length(A)), zeros(Int64, length(A))
    seq1defended_iads[1, :] = zeros(Int64, g.ndc)
    seq2attacked_iads[1, :] = zeros(Int64, g.nac)
    seq2defended_iads[1, :] = zeros(Int64, g.ndc)
    seq2attacked_cities[1, :] = zeros(Int64, g.nap)
    k = 1
    for i = 1:length(A[1]), j = 1:length(sensors), l = 1:length(A[5])
        k += 1
        Σ1[k] = string("c:", A[1][i], " c:", sensors[j], " d:", A[5][l])
        seqn1[k, :] = [i, 0, 0, j, l, 0]
        iads_ind = getindex(A[1][i], A[5][l])
        seq1defended_iads[k, :] = getindex(g.iads, iads_ind)
    end
    k = 1
    for i = 1:length(A[2]), j = 1:length(A[3])
        k += 1
        Σ2[k] = string("c:", A[2][i], " a:", A[3][j])
        seqn2[k, :] = [0, i, j, 0, 0, 0]
        seq2attacked_iads[k, :] = getindex(g.iads, getindex(A[2][i], A[3][j]))
    end
    for i = 1:length(A[2]), j = 1:length(A[3]), l = 1:length(iadsdefenses), m = 1:length(A[6])
        k += 1
        Σ2[k] = string("c:", A[2][i], " a:", A[3][j], " d:", iadsdefenses[l], " a:", A[6][m])
        seqn2[k, :] = [0, i, j, 0, l, m]
        attack_ind = getindex(A[2][i], A[3][j])
        seq2attacked_iads[k, :] = getindex(g.iads, attack_ind)
        seq2attacked_cities[k, :] = A[6][m]
    end
    (Σ1, Σ2), (seqn1, seqn2), (seq1defended_iads, seq2attacked_iads, seq2attacked_cities)
end

function build_Iset(g::AirDefenseGame, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    ndp, ndc, nac, iads = g.ndp, g.ndc, g.nac, g.iads
    Iset = Vector{String}(undef, ni1 + ni2)
    I1set = Vector{String}(undef, ni1)
    I2set = Vector{String}(undef, ni2)
    seqI1 = Vector{Int64}(undef, ni1)
    seqI2 = Vector{Int64}(undef, ni2)
    nextseqI1 = Vector{UnitRange{Int64}}(undef, ni1)
    nextseqI2 = Vector{UnitRange{Int64}}(undef, ni2)
    ki, ki1, ki2, ksi1, ksi2, knsi1, knsi2 = 1, 1, 1, 1, 1, 1, 1
    lbnext = ns2_stage[1] + 1
    for a in A[2] # offensive cyber attackable
        targetable = getindex(iads, a)
        Iset[ki] = "Attacker: targetable $targetable"
        ki += 1
        I2set[ki2] = "Attacker: targetable $targetable"
        ki2 += 1
        seqI2[ksi2] = 1
        ksi2 += 1
        nextseqI2[knsi2] = lbnext:(lbnext + na_stage[3] - 1)
        knsi2 += 1
        lbnext += na_stage[3]
    end
    sensors = Vector{Int64}[]
    for i = 0:nac
        push!(sensors, collect(subsets(iads, i))...)
    end
    lbnext = ns1_stage[1] + 1
    for a in A[1] # defensive cyber defendable
        defendable = getindex(iads, a)
        for s in sensors
            Iset[ki] = "Defender: defendable $defendable, sensed $s"
            ki += 1
            I1set[ki1] = "Defender: defendable $defendable, sensed $s"
            ki1 += 1
            seqI1[ksi1] = 1
            ksi1 += 1
            nextseqI1[knsi1] = lbnext:(lbnext + na_stage[5] - 1)
            knsi1 += 1
            lbnext += na_stage[5]
        end
    end
    lbnext = sum(ns2_stage[1:2]) + 1
    for (a1i, a1) in enumerate(A[2]) # offensive cyber capes
        targetable = getindex(iads, a1)
        for (a2i, a2) in enumerate(A[3])  # cyber attcks
            attacked = getindex(iads, getindex(a1, a2))
            seq2 = 1 + (a1i - 1) * length(A[3]) + a2i
            for defended in collect(subsets(iads, ndc))
                Iset[ki] = "Attacker: targetable $targetable, cyberattack $attacked, cyberdefend $defended"
                ki += 1
                I2set[ki2] = "Attacker: targetable $targetable, cyberattack $attacked, cyberdefend $defended"
                ki2 += 1
                seqI2[ksi2] = seq2
                ksi2 += 1
                nextseqI2[knsi2] = lbnext:(lbnext + na_stage[6] - 1)
                knsi2 += 1
                lbnext += na_stage[6]
            end
        end
    end
    Iset, (I1set, I2set), (seqI1, seqI2), (nextseqI1, nextseqI2)
end

# defender's utility for a combination of player actions
function utility(g::AirDefenseGame, C, defended_iads, attacked_iads, attacked_cities)
    u = 0.0
    for city in attacked_cities
        q = 1.0
        for defender in g.iads
            if C[defender, city] == 0
                q *= (1 - 0.0)
            elseif defender ∈ attacked_iads && defender ∉ defended_iads
                q *= (1 - 0.0)
            elseif defender ∉ attacked_iads
                q *= (1 - g.pp)
            elseif defender ∈ attacked_iads && defender ∈ defended_iads
                q *= (1 - g.pp * g.pdc * (1 - g.pac)) # if a defender is cyber attacked, it has pdc * pac of being cyber-defended
            else
                @warn "Unexpected defense layout."
            end
        end
        u += -1.0 * g.v[city] * q # defender's expected utility for a city
    end
    u
end

function build_utility_hist(g::AirDefenseGame, A, An)
    nz = prod(An[i][end] for i in 1:length(An))
    leafnodes = Vector{Vector{Int64}}(undef, nz)
    U = Vector{Float64}(undef, nz)
    k = 0
    C = coverage(g)
    println("Building history-form utilities...")
    # @progress for c1 in An[1], c2 in An[2], ac in An[3], c4 in An[4], dc in An[5], ap in An[6]
    for c1 in An[1], c2 in An[2], ac in An[3], c4 in An[4], dc in An[5], ap in An[6]
        k += 1
        h = [c1, c2, ac, c4, dc, ap]
        leafnodes[k] = h
        attacked_cities = @view A[6][h[6]][:]
        defended_iads = getindex(g.iads, getindex(A[1][h[1]], A[5][h[5]]))
        attacked_iads = getindex(g.iads, getindex(A[2][h[2]], A[3][h[3]]))
        U[k] = utility(g, C, defended_iads, attacked_iads, attacked_cities) # defender's expected utility for a leaf node
    end
    U, leafnodes
end

# Expected = true, sets the utliity to the expected utility over the sequence given the chance actions
function build_utility_seq(g::AirDefenseGame, gs::GameSet, nss, seqn, seqactions;
                            expected::Bool = false, timelimit::Int = 90)
    A = gs.A
    ns1, ns2 = nss
    seqn1, seqn2 = seqn
    seq1defended_iads, seq2attacked_iads, seq2attacked_cities = seqactions
    nstage = length(A)
    iseqn1 = zeros(Int64, nstage)
    iseqn2 = zeros(Int64, nstage)
    reward = spzeros(ns1, ns2)
    C = coverage(g)
    sensors = Vector{Int64}[]
    for i = 0:g.nac
        push!(sensors, collect(subsets(1:g.ndp, i))...) # by index
    end
    iadsdefenses = collect(subsets(1:g.ndp, g.ndc)) # by index
    probs = [(1 / 15) * (1 / 15) * chance(4, c4, g, gs) for c4 in gs.An[4]] # assumes equal prob in first two rounds
    p = 1.0 # default value if expected = false
    k = 0
    println("Building sequence-form utilities, expected = $expected...")
    starttime = now()
    # @progress for (ic1, c1) in enumerate(A[1])
    for (ic1, c1) in enumerate(A[1])
        iseqn1[1] = ic1
        for (ic2, c2) in enumerate(A[2])
            iseqn2[2] = ic2
            for (iac, ac) in enumerate(A[3])
                iseqn2[3] = iac
                for (ic4, c4) in enumerate(A[4])
                    sensed = getindex(c2, getindex(ac, c4 .== 2)) # ndp id's of sensed attacks
                    isensor = findfirst(x -> x == sensed, sensors)
                    iseqn1[4] = isensor
                    expected && (p = probs[ic4])
                    for (idc, dc) in enumerate(A[5])
                        iseqn1[5] = idc
                        idefended = findfirst(x -> x == getindex(c1, dc), iadsdefenses)
                        iseqn2[5] = idefended
                        ind1 = matchrow(iseqn1, seqn1)
                        for (iap, ap) in enumerate(A[6])
                            iseqn2[6] = iap
                            ind2 = matchrow(iseqn2, seqn2)
                            reward[ind1, ind2] = p * utility(g, C,
                                seq1defended_iads[ind1, :], seq2attacked_iads[ind2, :], seq2attacked_cities[ind2, :])
                            k += 1
                            k % 1_000 == 0 && print("\rNodes (M): $(k/1_000_000)")
                            if k % 10_000 == 0
                                (tominute(now() - starttime) > timelimit) && (return reward, false)
                            end
                        end
                    end
                end
            end
        end
    end
    return reward, true
end

function GameSet(g::AirDefenseGame)
    players = 1:2
    A, An, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    Iset, (I1set, I2set), (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    C = coverage(g)
    sensors = Vector{Int64}[]
    for i = 0:g.nac
        push!(sensors, collect(subsets(1:g.ndp, i))...)
    end
    iadsdefenses = collect(subsets(1:g.ndp, g.ndc))
    GameSet(players, A, An, na_stage, Iset, ni_stage, ni, C, sensors, iadsdefenses)
end

function GameOptSet(g::AirDefenseGame, gs::GameSet, reward)
    A, An, na_stage = gs.A, gs.An, gs.na_stage
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    Σ, seqn, seqactions = build_Σset(g, A, ns1, ns2)
    Iset, (I1set, I2set), (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    return GameOptSet(ns1, ns2, ni1, ni2, ns1_stage, ns2_stage, reward,
        gs.na_stage, gs.sensors, seqI1, seqI2, nextseqI1, nextseqI2)
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
"""
Build index from input parameters, e.g.

```julia
z = [
    (:a, [1.0, 1.9], [0.1, 0.1]),
    (:b, 1.0, 0.1),
    (:c, [2.9, 5.0], [0.1, 0.05])
]
```

is turned into

```julia
Dict([
    :a => 1:2,
    :b => 3,
    :c => 4:5,
    :a_1 => 1,
    :a_2 => 2,
    :c_1 => 4,
    :c_2 => 5
])
```

indicating at which index the parameters are stored in an array
generated with `_build_var_array()`.
"""
function _build_index(params)
    i = 1
    indices = []
    for (k, v...) in params
        j = i + length(v[1]) - 1
        if i == j
            push!(indices, k => i)
        else
            push!(indices, k => i:j)
            for (m, n) in enumerate(i:j)
                # store also the location of indexed variables
                push!(indices, Symbol(k, "_$m") => n)
            end
        end
        i += (j - i + 1)
    end
    Dict(indices)
end

"""
Build array from input parameters, e.g.

    ```julia
    z = [
        (:a, [1.0, 1.9], [0.1, 0.1]),
        (:b, 1.0, 0.1),
        (:c, [2.9, 5.0], [0.1, 0.05])
    ]
    ```

    yields

    ```julia
    # read values at position 2 of tuples in `z`
    _build_var_array(z, 2) == [1.0, 1.9, 1.0, 2.9, 5.0]
    # read uncertainties at position 3 of tuples in `z`
    _build_var_array(z, 3) == [0.1, 0.1, 0.1, 0.1, 0.05]
    ```
"""
_build_var_array(var_dict, n) = reduce((acc, elm) -> append!(acc, elm[n]), values(var_dict); init=Float64[])

struct LSQInput
    β::Vector{Float64}
    z::Vector{Float64}
    Σ::Matrix{Float64}
    βindex::Dict{Symbol, Union{Int, UnitRange{Int}}}
    zindex::Dict{Symbol, Union{Int, UnitRange{Int}}}
end

"""
    function LSQInput(β, z)

Compile the measured and guessed parameters for a constrained least squares problem.

`β` is a vector of pairs, each holding a symbol and initial guess for (an) unknown parameter(s).
The guess can be a scalar or vector of floats.

`z` is a vector of 3-tuples, each holding a symbol, the measured value and standard uncertainty
of a measured parameter. Values and uncertainties can be floats or vectors of floats.

# Example

For a linear fit, one could provide:

```julia
# measured parameters
z = [
    (:x, [1.0, 1.9, 3.2, 3.8], [0.1,  0.1, 0.2, 0.3]),
    (:y, [2.9, 5.0, 6.9, 9.1], [0.1, 0.05, 0.1, 0.2])
]

# unknown parameters, a = slope, b = offset
β₀ = [
    (:a,  1.0),
    (:b,  1.0)
]

input = LSQInput(z, β₀)
```
"""
function LSQInput(β, z)
    uz = _build_var_array(z, 3) # third element in tuple = uncertainty
    return LSQInput(
        _build_var_array(β, 2),
        _build_var_array(z, 2), # second element in tuple = value
        Diagonal(uz.^2) |> Matrix,
        _build_index(β),
        _build_index(z)
    )
end

"""
    function covariance!(input::LSQInput, x::Symbol, y::Symbol, value::Float64)

Set covariance of two input parameters stored in `input`.

If parameter is part of an array, use underscore to indicate which
covariance shall be set, e.g.

```julia
function covariance!(input, :x_1, :x_2, 0.1)
```

sets the covariance of the first and second element in
the array parameter `x`.
"""
function set_covariance!(input::LSQInput, x::Symbol, y::Symbol, value::Float64)
    i = input.zindex[x]
    j = input.zindex[y]

    (!isa(i, Integer) || !isa(j, Integer)) && ArgumentError("Covariance can only be set on individual variables.")

    input.Σ[i, j] = input.Σ[j, i] = value

    return nothing
end


text:
export frame3dd,
       frame3DD_circular_tube_truss_bar_family

@kwdef struct Frame3DDBackend{K,T} <: LazyBackend{K,T}
  realized::Parameter{Bool}=Parameter(false)
  truss_nodes::Vector{<:TrussNode}=TrussNode[]
  truss_bars::Vector{<:TrussBar}=TrussBar[]
  shapes::Vector{<:Shape}=Shape[] # This contains all the rest that is not treated yet
  truss_node_data::Vector{TrussNodeData}=TrussNodeData[]
  truss_bar_data::Vector{TrussBarData}=TrussBarData[]
  view::View=default_view()
end

abstract type FR3DDKey end
const FR3DDId = Any
const FR3DDIds = Vector{FR3DDId}
const FR3DDRef = GenericRef{FR3DDKey, FR3DDId}
const FR3DDRefs = Vector{FR3DDRef}
const FR3DDEmptyRef = EmptyRef{FR3DDKey, FR3DDId}
const FR3DDUniversalRef = UniversalRef{FR3DDKey, FR3DDId}
const FR3DDNativeRef = NativeRef{FR3DDKey, FR3DDId}
const FR3DDUnionRef = UnionRef{FR3DDKey, FR3DDId}
const FR3DDSubtractionRef = SubtractionRef{FR3DDKey, FR3DDId}
const FR3DD = Frame3DDBackend{FR3DDKey, FR3DDId}

void_ref(b::FR3DD) = FR3DDNativeRef(-1)

const frame3dd = FR3DD()

KhepriBase.backend_name(b::FR3DD) = "Frame3DD"

# Frame3DD needs to merge nodes and bars
save_shape!(b::FR3DD, s::TrussNode) = maybe_merged_node(b, s)
save_shape!(b::FR3DD, s::TrussBar) = maybe_merged_bar(b, s)

# Frame3DD does not need layers
use_material_as_layer(b::FR3DD) = false
#with_material_as_layer(f::Function, backend::FR3DD, material::Material) = f()

# Frame3DD Families
abstract type Frame3DDFamily <: Family end

Base.@kwdef struct Frame3DDTrussBarGeometry
  Ax  # cross section area
  Asy # shear area y-direction
  Asz # shear area z-direction
  Jxx # torsional moment of inertia - x axis
  Iyy # bending moment of inertia - y axis
  Izz # bending moment of inertia - z axis
end

Base.@kwdef struct Frame3DDTrussBarFamily <: Frame3DDFamily
  geometry=Ref{Any}(nothing) # properties that only depend on radius and thickness
  E   # elastic modulus
  G   # shear modulus
  p   # roll angle
  d   # mass density
end

export frame3dd_truss_bar_family
frame3dd_truss_bar_family = Frame3DDTrussBarFamily

truss_bar_family_cross_section_area(f::Frame3DDTrussBarFamily) =
  f.geometry[].Ax

#=
Circular Tube (outer radius= Ro, inner radius = Ri):

Ax = π ( Ro2 - Ri2 )
Asy = Asz = Ax / ( 0.54414 + 2.97294(Ri/Ro) - 1.51899(Ri/Ro)2 ) ± 0.05%
Jxx = (1/2) π ( Ro4 - Ri4 )
Ixx = Iyy = (1/4) π ( Ro4 - Ri4 )
=#

#=
outer radius, thickness;
E   # elastic modulus
G   # shear   modulus
p   # roll angle
d   # mass density
=#
export frame3DD_circular_tube_truss_bar_geometry
frame3dd_circular_tube_truss_bar_geometry(rₒ, e) =
  let rᵢ = rₒ - e,
      Ax = annulus_area(rₒ, rᵢ),
      Asyz = Ax/(0.54414 + 2.97294*(rᵢ/rₒ) - 1.51899*(rᵢ/rₒ)^2),
      Jxx = π/2*(rₒ^4 - rᵢ^4),
      Ixxyy = Jxx/2
    Frame3DDTrussBarGeometry(
      Ax =Ax,
      Asy=Asyz,
      Asz=Asyz,
      Jxx=Jxx,
      Iyy=Ixxyy,
      Izz=Ixxyy)
 end

#=
 @rad1 = 25.0                   # Section diameter (mm)
 @thk1 = 1.6                    # Section thickness (mm)

 @section = "CHS_25x1.6"        # Tube name
 @material = "Steel C450"       # Duragal DualGrade RHS Grade C450Lo (AS 1163)
 @ftype = "Tube"                # Tube type (Tube, Beam, Channel, Angle)
 @emod = "210.0"                # Young's Modulus (GN per square metre) for steel
 @gmod = "79.0"                 # Modulus of Rigidity (GN per square metre) for steel
 @pang = "0.0"                  # section axis rotation angle (deg)
 @ax = "568.3"                  # cross-sectional area (mm^2)
 @asy = "3.2"                   # 2 x thickness (mm) shear area
 @asz = "3.2"                   # 2 x thickness (mm) shear area
 @jxx = "0.0585E+06"            # J (mm^4)
 @iyy = "0.0702E+06"            # Iyy (mm^4)
 @izz = "0.0237E+06"            # Izz (mm^4)
=#

backend_get_family_ref(b::FR3DD, f::TrussBarFamily, tbf::Frame3DDTrussBarFamily) =
  begin
    tbf.geometry[] = frame3dd_circular_tube_truss_bar_geometry(f.radius, f.radius-f.inner_radius)
    tbf
  end

#
KhepriBase.b_delete_all_refs(b::FR3DD) =
  begin
    empty!(b.truss_nodes)
    empty!(b.truss_bars)
    b.realized(false)
    b
  end

realize(b::FR3DD, s::TrussNode) =
  error("BUM")

realize(b::FR3DD, s::TrussBar) =
  error("BUM")

realize(b::FR3DD, s::Panel) =
  error("BUM")

#new_truss_analysis(v=nothing; self_weight=false, backend=frame3dd) =
displacements_from_frame3dd(b::FR3DD, filename, load::Vec, self_weight::Bool, point_loads::Dict) =
    let nodes = process_nodes(b.truss_nodes, load, point_loads),
      bars = process_bars(b.truss_bars, nodes),
      supports = unique(filter(s -> s != false, map(n -> n.family.support, nodes))),
      loaded_nodes = filter(!truss_node_is_supported, nodes),
      supported_nodes = filter(truss_node_is_supported, nodes),
      shear = 0,		     # 1: include shear deformation effects, 0: don't
      geom  = 0,		     # 1: include geometric stiffness effects, 0: don't
      exagg = 1,		     # exaggeration factor for Gnuplot output
      scale = 1.0,		   # zoom scale factor for Gnuplot plotting
      dx = 1.0		       # x-axis increment for internal force calc's
    empty!(b.truss_node_data)
    append!(b.truss_node_data, nodes)
    empty!(b.truss_bar_data)
    append!(b.truss_bar_data, bars)
    open(filename, "w") do io
      println(io, "Frame analysis for Khepri");
      println(io, "%% node data ...");
      println(io, "$(length(nodes))\t\t%% number of nodes");
      println(io, "%% J\t\tX\t\tY\t\tZ\t\tr\t\tnode coordinates");
      for n in nodes
        i = n.id
        p = n.loc
        println(io, "$(i)\t$(p.x)\t$(p.y)\t$(p.z)\t0")
      end
      println(io, "");
      println(io, "%% reaction data ...");
      #nR = count(has_support, nodes);
      println(io, "$(length(supported_nodes))    %% number of nodes with reaction forces");
      println(io, "%% j\tRx\tRy\tRz\tRxx\tRyy\tRzz");
      for n in supported_nodes
        s = n.family.support
        print(io, "$(n.id),\t$(Int(s.ux)),\t$(Int(s.uy)),\t$(Int(s.uz))\t")
        println(io, "$(Int(s.rx)),\t$(Int(s.ry)),\t$(Int(s.rz))")
      end
      println(io, "")
      println(io, "%% frame element section property data ...")
      println(io, "$(length(bars))\t\t%% number of frame elements")
      println(io, "%% m\tn1\tn2\t\tAx\t\tAsy\t\tAsz\t\tJxx\t\tIyy\t\tIzz\t\tE\t\tG\t\tp\tdensity");
      for bar in bars
        f = family_ref(b, bar.family)
        g = f.geometry[]
        print(io, "$(bar.id),\t$(bar.node1.id),\t$(bar.node2.id),")
        print(io, "\t$(g.Ax),\t$(g.Asy),\t$(g.Asz),")
        print(io, "\t$(g.Jxx),\t$(g.Iyy),\t$(g.Izz),")
        println(io, "\t$(f.E),\t$(f.G),\t$(f.p),\t$(f.d)")
      end
      println(io, "\n");
      println(io, "$(shear)\t\t%% 1: include shear deformation, 0: do not");
      println(io, "$(geom)\t\t%% 1: include geometric stiffness, 0: do not");
      println(io, "$(exagg)\t%% exagerate deformations in plotting");
      println(io, "$(scale)\t%% zoom scale factor for 3D plotting");
      println(io, "$(dx)\t%% x-axis increment for internal forces calc");
      println(io, "\n");
      println(io, "%% static load data ...");
      println(io, "1\t\t%% number of static load cases");
      println(io, "\t\t%% begin static load case 1 of 1");
      println(io, "%% gravitational acceleration for self-weight loading");
      println(io, "%% gX         gY         gZ");
      println(io, "  0.0        0.0        $(self_weight ? -9.81 : 0.0)\n");
      println(io, "$(length(loaded_nodes))\t\t%% number of loaded nodes");
      println(io, "%% j\t\tFx\t\tFy\t\tFz\t\tMxx\t\tMyy\t\tMzz");
      for n in loaded_nodes
        println(io, "$(n.id),\t$(n.load.x),\t$(n.load.y),\t$(n.load.z),\t0,\t0,\t0")
      end
      println(io, "\n");
      println(io, "0\t\t%% number of members with uniform distributed loads")
      println(io, "0\t\t%% number of members with trapezoidal loads");
      println(io, "0\t\t%% number of members with internal point loads");
      println(io, "0\t\t%% number of members with temperature loads");
      println(io, "0\t\t%% number of nodes with prescribed displacements");
      println(io, "\t\t%% end   static load case 1 of 1\n");
      println(io, "%% inertial load data ...");
      println(io, "0\t\t%% number of dynamic modes to analyze");
    end
  end

frame3dd_plugin = joinpath(dirname(abspath(@__DIR__)), "bin", "frame3dd.exe")

frame3dd_simulation_path() =
  mktempdir(tempdir(), prefix="Frame3DD_")

#=
KhepriBase.b_truss_analysis(b::FR3DD, load::Vec, self_weight::Bool, point_loads::Dict) =
  # Ensure extension is FMM to force Matlab mode
  let simulation_folder = frame3dd_simulation_path(),
       input_path = joinpath(simulation_folder, "IOdata.IN"),
       output_path = joinpath(simulation_folder, "IOdata.OUT")
    #@info input_path
    displacements_from_frame3dd(b, input_path, load, self_weight, point_loads)
    #try
      withenv("FRAME3DD_OUTDIR"=>simulation_folder) do
        run(`$frame3dd_plugin -q -i $input_path -o $output_path`)
      end
    #catch e
    #end
    lines = readlines(output_path)
    idx1 = findfirst(l->startswith(l, "N O D E   D I S P L A C E M E N T S"), lines) + 2
    idx2 = findnext(l->startswith(l, "F R A M E   E L E M E N T   E N D   F O R C E S"), lines, idx1)
    Dict([let strs = split(line)
            parse(Int, strs[1])=>parse.(Float64, strs[2:7])
          end
          for line in lines[idx1:idx2-1]])
  end

KhepriBase.b_node_displacement_function(b::FR3DD, results) =
  n -> vxyz(get(results, n.id, (0,0,0,0,0,0))[1:3]..., world_cs)
=#

#=
using Libdl
f3ddlib = dlopen(joinpath(dirname(abspath(@__DIR__)), "bin", "Frame3DDLib"))
f3ddmain = dlsym(f3ddlib, :simulate)
=#

#dlclose(f3ddlib)

const f3ddlibpath = joinpath(dirname(abspath(@__DIR__)), "bin", "Frame3DDLib")
##########################################
KhepriBase.b_truss_analysis(be::FR3DD, load::Vec, self_weight::Bool, point_loads::Dict) =
    let nodes = process_nodes(be.truss_nodes, load, point_loads),
        bars = process_bars(be.truss_bars, nodes),
        supports = unique(filter(s -> s != false, map(n -> n.family.support, nodes))),
        loaded_nodes = filter(!truss_node_is_supported, nodes),
        supported_nodes = filter(truss_node_is_supported, nodes),
        shear = 0,		     # 1: include shear deformation effects, 0: don't
        geom  = 0,		     # 1: include geometric stiffness effects, 0: don't
        exagg = 1,		     # exaggeration factor for Gnuplot output
        scale = 1.0,		   # zoom scale factor for Gnuplot plotting
        dx = 1.0,		       # x-axis increment for internal force calc's
        nN = length(nodes), # number of nodes
        nE = length(bars),  # number of bars
        nodes_coords = zeros(Float64, 3*(nN+1)), # X,Y,Z node coordinates (global)
        nodes_radius = zeros(Float32, nN + 1),   # node size radius, for finite sizes
        DoF = 6*nN,
        nR = length(supported_nodes),            # number of restrained nodes
        q = zeros(Int32, DoF+1),
        r = zeros(Int32, DoF+1),
        sumR = 0
        empty!(be.truss_node_data)
        append!(be.truss_node_data, nodes)
        empty!(be.truss_bar_data)
        append!(be.truss_bar_data, bars)
      for n in nodes
          i = n.id
          p = n.loc
        for j in 1:3
          nodes_coords[3*i+j]=p.raw[j]
        end
      end
      for n in supported_nodes
        s = n.family.support
        j = n.id
        sj = [Int(s.ux), Int(s.uy), Int(s.uz), Int(s.rx), Int(s.ry), Int(s.rz)]
        for l in 5:-1:0
          r[6*j-l+1] = sj[6-l]
        end
      end
      sumR = sum(r)
      sumR > 4 || error("Un-restrained structure: only $(sumR) imposed reactions. At least 4 reactions are required to support static loads.")
      sumR < DoF || error("Fully restrained structure: $(sumR) imposed reactions >= $(DoF) degrees of freedom")
      for i in 1:DoF
        q[i+1] = r[i+1] > 0 ? 0 : 1
      end
      L = zeros(Float64, nE+1)   # length of each element
      Le = zeros(Float64, nE+1)  # effective length of each element
      N1  = zeros(Int32, nE+1)   # node #1 of each element
      N2  = zeros(Int32, nE+1)   # node #2 of each element
    	Ax  = zeros(Float32, nE+1) # cross section area of each element	*/
    	Asy = zeros(Float32, nE+1) # shear area in local y direction 	*/
    	Asz = zeros(Float32, nE+1) # shear area in local z direction	*/
    	Jx  = zeros(Float32, nE+1) # torsional moment of inertia 		*/
    	Iy  = zeros(Float32, nE+1) # bending moment of inertia about y-axis */
    	Iz  = zeros(Float32, nE+1) # bending moment of inertia about z-axis */
    	E   = zeros(Float32, nE+1) # frame element Young's modulus	*/
    	G   = zeros(Float32, nE+1) # frame element shear modulus		*/
    	p   = zeros(Float32, nE+1) # element rotation angle about local x axis */
    	d   = zeros(Float32, nE+1) # element mass density			*/
      for bar in bars
        f = family_ref(be, bar.family)
        g = f.geometry[]
        b = bar.id
        N1[b+1] = bar.node1.id
        N2[b+1] = bar.node2.id
        Ax[b+1] = g.Ax
        Asy[b+1] = g.Asy
        Asz[b+1] = g.Asz
        Jx[b+1] = g.Jxx
        Iy[b+1] = g.Iyy
        Iz[b+1] = g.Izz
        E[b+1] = f.E
        G[b+1] = f.G
        p[b+1] = f.p
        d[b+1] = f.d
        n1 = N1[b+1]
        n2 = N2[b+1]
        L[b+1] = distance(bar.node2.loc, bar.node1.loc)
        Le[b+1] = L[b+1] - nodes_radius[n2+1] - nodes_radius[n1+1]
      end
      nF = length(loaded_nodes)
      F_mech = zeros(Float64, DoF+1)
      for n in loaded_nodes
        F_mech[6*n.id - 5 + 1] = n.load.x;
        F_mech[6*n.id - 4 + 1] = n.load.y;
        F_mech[6*n.id - 3 + 1] = n.load.z;
      end
      gX = [0f0, 0f0]
      gY = [0f0, 0f0]
      gZ = [0f0, self_weight ? -9.81f0 : 0f0]
      nFs = Int32[0, nF]
      F_mechs = [Ref(F_mech, 1), Ref(F_mech, 1)]
      D = zeros(DoF+1)	  # displacements of each node
      #ccall(f3ddmain, Int32, (
      #      Int32,        # number of Nodes
      #      Int32,        # number of frame Elements
      #      Ref{Float64}, # X,Y,Z node coordinates (global)
      #      Ref{Float32}, # node size radius, for finite sizes
      #      Int32,        # number of restrained nodes
      #      Ref{Int32},   # reaction data
      #      Ref{Int32},   # reaction data
      #      Int32,        # total no. of reactions
      #      Ref{Float64},	# length of each element
      #      Ref{Float64},	# effective length of each element
      #      Ref{Int32},   # node #1 of each element
      #      Ref{Int32},   # N2,    // node #2 of each element
      #      Ref{Float32},	# cross section area of each element
      #      Ref{Float32},	# shear area in local y direction 	*/
      #      Ref{Float32},	# shear area in local z direction	*/
      #      Ref{Float32},	# torsional moment of inertia 		*/
      #      Ref{Float32},	# bending moment of inertia about y-axis */
      #      Ref{Float32},	# bending moment of inertia about z-axis */
      #      Ref{Float32},	# frame element Young's modulus	*/
      #      Ref{Float32},	# frame element shear modulus		*/
      #      Ref{Float32},	# element rotation angle about local x axis */
      #      Ref{Float32},	# element mass density			*/
      #      Int32,        # number of Load cases
      #      Int32,        # number of Load cases
      #      Ref{Float32}, # gravitational acceleration in global X
      #      Ref{Float32}, # gravitational acceleration in global Y
      #      Ref{Float32}, # gravitational acceleration in global Z
      #      Ref{Int32},   # number of loaded nodes
      #      Ref{Ptr{Float64}},  # mechanical load vectors, all load cases
      #      # RESULTS
      #      Ref{Float64}  # displacements of each node
      #      ),
      #  nN, nE,
      #  nodes_coords, nodes_radius,
      #  nR, q, r, sumR,
      #  L, Le, N1, N2, Ax, Asy, Asz, Jx, Iy, Iz, E, G, p, d,
      #  1, 0,
      #  gX, gY, gZ,
      #  nFs,
      #  F_mechs,
      #  # RESULTS
      #  D)
      @ccall f3ddlibpath.simulate(
         nN::Int32,        # number of Nodes
         nE::Int32,        # number of frame Elements
         nodes_coords::Ref{Float64}, # X,Y,Z node coordinates (global)
         nodes_radius::Ref{Float32}, # node size radius, for finite sizes
         nR::Int32,        # number of restrained nodes
         q::Ref{Int32},   # reaction data
         r::Ref{Int32},   # reaction data
         sumR::Int32,        # total no. of reactions
         L::Ref{Float64},	# length of each element
         Le::Ref{Float64},	# effective length of each element
         N1::Ref{Int32},   # node #1 of each element
         N2::Ref{Int32},   # node #2 of each element
         Ax::Ref{Float32},	# cross section area of each element
         Asy::Ref{Float32},	# shear area in local y direction 	*/
         Asz::Ref{Float32},	# shear area in local z direction	*/
         Jx::Ref{Float32},	# torsional moment of inertia 		*/
         Iy::Ref{Float32},	# bending moment of inertia about y-axis */
         Iz::Ref{Float32},	# bending moment of inertia about z-axis */
         E::Ref{Float32},	# frame element Young's modulus	*/
         G::Ref{Float32},	# frame element shear modulus		*/
         p::Ref{Float32},	# element rotation angle about local x axis */
         d::Ref{Float32},	# element mass density			*/
         1::Int32,        # number of Load cases
         0::Int32,        # number of Load cases
         gX::Ref{Float32}, # gravitational acceleration in global X
         gY::Ref{Float32}, # gravitational acceleration in global Y
         gZ::Ref{Float32}, # gravitational acceleration in global Z
         nFs::Ref{Int32},   # number of loaded nodes
         F_mechs::Ref{Ptr{Float64}},  # mechanical load vectors, all load cases
         # RESULTS
         D::Ref{Float64}  # displacements of each node
      )::Int32
      D
    end
#  end

KhepriBase.b_node_displacement_function(b::FR3DD, results) =
  n -> vxyz(results[6*n.id-4:6*n.id-2]..., world_cs)


text:
doc"""
    LogNormal(μ,σ)

The *log normal distribution* is the distribution of the exponential of a [`Normal`](:func:`Normal`) variate: if $X \sim \operatorname{Normal}(\mu, \sigma)$ then $\exp(X) \sim \operatorname{LogNormal}(\mu,\sigma)$. The probability density function is

$f(x; \mu, \sigma) = \frac{1}{x \sqrt{2 \pi \sigma^2}}
\exp \left( - \frac{(\log(x) - \mu)^2}{2 \sigma^2} \right),
\quad x > 0$

```julia
LogNormal()          # Log-normal distribution with zero log-mean and unit scale
LogNormal(mu)        # Log-normal distribution with log-mean mu and unit scale
LogNormal(mu, sig)   # Log-normal distribution with log-mean mu and scale sig

params(d)            # Get the parameters, i.e. (mu, sig)
meanlogx(d)          # Get the mean of log(X), i.e. mu
varlogx(d)           # Get the variance of log(X), i.e. sig^2
stdlogx(d)           # Get the standard deviation of log(X), i.e. sig
```

External links

* [Log normal distribution on Wikipedia](http://en.wikipedia.org/wiki/Log-normal_distribution)

"""
immutable LogNormal <: ContinuousUnivariateDistribution
    μ::Float64
    σ::Float64

    LogNormal(μ::Real, σ::Real) = (@check_args(LogNormal, σ > zero(σ)); new(μ, σ))
    LogNormal(μ::Real) = new(μ, 1.0)
    LogNormal() = new(0.0, 1.0)
end

@distr_support LogNormal 0.0 Inf

#### Parameters

params(d::LogNormal) = (d.μ, d.σ)

#### Statistics

meanlogx(d::LogNormal) = d.μ
varlogx(d::LogNormal) = abs2(d.σ)
stdlogx(d::LogNormal) = d.σ

mean(d::LogNormal) = ((μ, σ) = params(d); exp(μ + 0.5 * σ^2))
median(d::LogNormal) = exp(d.μ)
mode(d::LogNormal) = ((μ, σ) = params(d); exp(μ - σ^2))

function var(d::LogNormal)
    (μ, σ) = params(d)
    σ2 = σ^2
    (exp(σ2) - 1.0) * exp(2.0 * μ + σ2)
end

function skewness(d::LogNormal)
    σ2 = varlogx(d)
    e = exp(σ2)
    (e + 2.0) * sqrt(e - 1.0)
end

function kurtosis(d::LogNormal)
    σ2 = varlogx(d)
    e = exp(σ2)
    e2 = e * e
    e3 = e2 * e
    e4 = e3 * e
    e4 + 2.0 * e3 + 3.0 * e2 - 6.0
end

function entropy(d::LogNormal)
    (μ, σ) = params(d)
    0.5 * (1.0 + log(twoπ * σ^2)) + μ
end

#### Evalution

pdf(d::LogNormal, x::Float64) = normpdf(d.μ, d.σ, log(x)) / x
function logpdf(d::LogNormal, x::Float64)
    if !insupport(d, x)
        return -Inf
    else
        lx = log(x)
        return normlogpdf(d.μ, d.σ, lx) - lx
    end
end

cdf(d::LogNormal, x::Float64) = x > 0.0 ? normcdf(d.μ, d.σ, log(x)) : 0.0
ccdf(d::LogNormal, x::Float64) = x > 0.0 ? normccdf(d.μ, d.σ, log(x)) : 1.0
logcdf(d::LogNormal, x::Float64) = x > 0.0 ? normlogcdf(d.μ, d.σ, log(x)) : -Inf
logccdf(d::LogNormal, x::Float64) = x > 0.0 ? normlogccdf(d.μ, d.σ, log(x)) : 0.0

quantile(d::LogNormal, q::Float64) = exp(norminvcdf(d.μ, d.σ, q))
cquantile(d::LogNormal, q::Float64) = exp(norminvccdf(d.μ, d.σ, q))
invlogcdf(d::LogNormal, lq::Float64) = exp(norminvlogcdf(d.μ, d.σ, lq))
invlogccdf(d::LogNormal, lq::Float64) = exp(norminvlogccdf(d.μ, d.σ, lq))

function gradlogpdf(d::LogNormal, x::Float64)
    (μ, σ) = params(d)
    x > 0.0 ? - ((log(x) - μ) / (σ^2) + 1.0) / x : 0.0
end

# mgf(d::LogNormal)
# cf(d::LogNormal)

#### Sampling

rand(d::LogNormal) = exp(randn() * d.σ + d.μ)

## Fitting

function fit_mle{T <: Real}(::Type{LogNormal}, x::AbstractArray{T})
    lx = log(x)
    μ, σ = mean_and_std(lx)
    LogNormal(μ, σ)
end


text:
stderr_copy = stderr

rserr, wrerr = redirect_stderr()

using Weave
# For some reason this is needed on Windows, without it we see lots of errors
@info "Ignore"

redirect_stderr(stderr_copy)
close(rserr)
close(wrerr)

input_file = readline()
output_file = readline()
doctype = readline()

if doctype == "PREVIEW"
    template_path = joinpath(dirname(@__FILE__), "preview.tpl")

    Weave.weave(input_file, out_path = output_file, template = template_path)
else
    cd(dirname(input_file))
    Weave.weave(input_file, out_path = :doc, doctype = doctype)
end