text:
module _hXMSEn
export hXMSEn
using Statistics: std, mean, median, var
using Plots
    """
        MSx, Sn, CI = hXMSEn(Sig, Mobj)

    Returns a vector of cross-entropy values (`MSx`) calculated at each node
    in the hierarchical tree, the average cross-entropy value across all
    nodes at each scale (`Sn`), and the complexity index (`CI`) of the hierarchical
    tree (i.e. sum(`Sn`)) between the data sequences contained in `Sig` using
    the parameters specified by the multiscale object (`Mobj`) over 3 temporal
    scales (default).
    The entropy values in `MSx` are ordered from the root node (S.00) to the
    Nth subnode at scale T (S.TN): i.e. S.00, S.10, S.11, S.20, S.21, S.22,
    S.23, S.30, S.31, S.32, S.33, S.34, S.35, S.36, S.37, S.40, ... , S.TN.
    The average cross-entropy values in Sn are ordered in the same way, with the
    value of the root node given first: i.e. S0, S1, S2, ..., ST

        MSx, Sn, CI = hXMSEn(Sig::AbstractArray{T,2} where T<:Real, Mobj::NamedTuple;
                                 Scales::Int=3, RadNew::Int=0, Plotx::Bool=false)

    Returns a vector of cross-entropy values (`MSx`) calculated at each node
    in the hierarchical tree, the average cross-entropy value across all
    nodes at each scale (`Sn`), and the complexity index (`CI`) of the entire
    hierarchical tree between the data sequences contained in `Sig` using
    the following name/value pair arguments:

    # Arguments:
    `Scales`   - Number of temporal scales, an integer > 1   (default: 3)
                 At each scale (T), entropy is estimated for 2^(T-1) nodes. \n
    `RadNew`   - Radius rescaling method, an integer in the range [1 4].
                 When the entropy specified by `Mobj` is `XSampEn` or `XApEn`,
                 `RadNew` allows the radius threshold to be updated at each
                 node in the tree. If a radius value is specified by `Mobj` (`r`),
                 this becomes the rescaling coefficient, otherwise it is set
                 to 0.2 (default). The value of `RadNew` specifies one of the
                 following methods: \n
                 [1]    Standard Deviation          - r*std(Xt) \n
                 [2]    Variance                    - r*var(Xt) \n
                 [3]    Mean Absolute Deviation     - r*mean_ad(Xt) \n
                 [4]    Median Absolute Deviation   - r*med_ad(Xt,1) \n
    `Plotx`    - When `Plotx` == true, returns a plot of the average cross-entropy
                 value at each time scale (i.e. the multiscale entropy curve)
                 and a hierarchical graph showing the entropy value of each node
                 in the hierarchical tree decomposition.  (default: false) \n

    # See also `MSobject`, `XMSEn`, `rXMSEn`, `cXMSEn`, `XSampEn`, `XApEn`, `hMSEn`

    # References:
        [1]   Matthew W. Flood,
            "hXMSEn - EntropyHub Project"
            2021, https://github.com/MattWillFlood/EntropyHub

        [2]   Rui Yan, Zhuo Yang, and Tao Zhang,
            "Multiscale cross entropy: a novel algorithm for analyzing two
            time series."
            5th International Conference on Natural Computation.
            Vol. 1, pp: 411-413 IEEE, 2009.

        [3] Ying Jiang, C-K. Peng and Yuesheng Xu,
            "Hierarchical entropy analysis for biological signals."
            Journal of Computational and Applied Mathematics
            236.5 (2011): 728-742.

    """
    function hXMSEn(Sig::AbstractArray{T,2} where T<:Real, Mobj::NamedTuple;
        Scales::Int=3, RadNew::Int=0, Plotx::Bool=false)

    size(Sig,1) == 2 ? Sig = Sig' : nothing

    (size(Sig,1)>10) ? nothing : error("Sig:   must be a numeric vector" )
    (length(Mobj) >= 1) ? nothing :  error("Mobj:    must be a multiscale entropy object created
            with the function EntropyHub.MSobject")
    (Scales>1) ? nothing : error("Scales:     must be an integer > 1")
    (RadNew==0 || (RadNew in 1:4 && String(Symbol(Mobj.Func)) in ("XSampEn","XApEn"))) ? nothing :
    error("RadNew:  must be 0, or an integer in range [1 4] with entropy function `XSampEn` or `XApEn`")

    Args = NamedTuple{keys(Mobj)[2:end]}(Mobj)
    if RadNew > 0
        if RadNew == 1
            Rnew = x -> std(x, corrected=false)
        elseif RadNew == 2
            Rnew = x -> var(x, corrected=false)
        elseif RadNew == 3
            Rnew = x -> mean(abs.(x .- mean(x)))
        elseif RadNew == 4
            Rnew = x -> median(abs.(x .- median(x)))
        end

        if haskey(Mobj,:r)
            Cx = Mobj.r
        else
            Cy = ("Standard Deviation","Variance","Mean Abs Deviation",
                    "Median Abs Deviation")
            @warn("No radius value provided in Mobj.
                Default set to 0.2*$(Cy[RadNew]) of each new time-series.")
            Cx = .2
        end
    end

    XX, YY, N = Hierarchy(Sig, Scales)
    MSx = zeros(size(XX,1))
    for T = 1:size(XX,1)
        print(" .")
        Temp = hcat(XX[T,1:Int(N/(2^(floor(log2(T)))))],
                YY[T,1:Int(N/(2^(floor(log2(T)))))])
        RadNew > 0 ? Args = (Args..., r=Cx*Rnew(Temp[:])) : nothing
        Temp2 = Mobj.Func(Temp; Args...)
        typeof(Temp2)<:Tuple ? MSx[T] = Temp2[1][end] : MSx[T] = Temp2[end]
    end

    Sn = zeros(Scales)
    for t = 1:Scales
        Sn[t] = mean(MSx[2^(t-1):(2^t)-1])
    end

    CI = sum(Sn)
    print("\n")
    if any(isnan.(MSx))
        println("Some entropy values may be undefined.")
    end

    if Plotx
        p1 = plot(1:Scales, Sn, c=RGB(8/255, 63/255, 77/255), lw=3)
        scatter!(1:Scales, Sn, markersize=6, c=RGB(1, 0, 1),
        xlabel = "Scale Factor", ylabel = "Entropy Value",
        guidefont = font(12, "arial", RGB(7/255, 54/255, 66/255)),
        tickfontsize = 10, tickfontfamily="arial", legend=false,
        title = "Hierarchical Multiscale $(Mobj.Func) Entropy",
        grid=false, xlim=(0.5,Scales+.5), ylim=(minimum(Sn)-.25,maximum(Sn)+.25),
        plot_titlefontsize=16, plot_titlefontcolor=RGB(7/255, 54/255, 66/255))

        N = 2^(Scales-1)
        x = zeros(2*N - 1)
        x[1] = N
        y = Scales.*(Scales .- floor.(Int, log2.(1:(2*N)-1)))
        for k = 1:2*N-1
            Q = floor(Int, log2(k))
            P = floor(Int,k/2)
            if k>1
                Bool(k%2) ? x[k] = x[P] + N/(2^Q) : x[k] = x[P] - N/(2^Q)
            end
        end

        Edges = hcat(repeat(1:N-1,inner=2),2:2*N-1)
        labx = [k for k in string.(round.(MSx,digits=3))]
        p2 = scatter(x,y,markersize=10*(MSx.-(minimum(MSx)).+1)./(abs(minimum(MSx))+1),
        c=RGB(1,0,1), legend=false, xticks=false, yticks=false)  #, txt=labx, markerfontsize=8)
        annotate!(x, y.+1, labx, 12-Scales)
        for k = 1:length(x)-1
            plot!(x[Edges[k,:]],y[Edges[k,:]],c=RGB(8/255,63/255,77/255),
            lw=2.5, ylim=(Scales-1, maximum(y)+1), xlim=(0,2*N), grid=false)
        end
        px = plot()
        pt = plot(p1, px, p2, layout = grid(3,1, heights=[0.25,.1, 0.8]))
        display(pt)

    end

    return MSx, Sn, CI
    end

    function Hierarchy(Z,sx)
        N = Int(2^floor(log2(size(Z,1))))
        if mod(log2(size(Z,1)),1) != 0
            @warn("Only first $(N) samples were used in hierarchical decomposition.
                The last $(size(Z,1)-N) samples of the data sequence were ignored.")
        end
        if N/(2^(sx-1)) < 8
           error("Data length ($(N)) is too short to estimate entropy at the lowest subtree.
           Consider reducing the number of scales.")
        end

        Z = Z[1:N,:]
        U1 = zeros((2^sx)-1,N);  U2 = zeros((2^sx)-1,N)
        U1[1,:] = Z[:,1];  U2[1,:] = Z[:,2]

        p=2
        for k = 1:sx-1
            for n = 1:2^(k-1)
                Temp = U1[2^(k-1)+n-1,:]
                U1[p,1:Int(N/2)]  = (Temp[1:2:end] + Temp[2:2:end])/2
                U1[p+1,1:Int(N/2)]= (Temp[1:2:end] - Temp[2:2:end])/2
                p +=2
            end
        end

        p=2
        for k = 1:sx-1
            for n = 1:2^(k-1)
                Temp = U2[2^(k-1)+n-1,:]
                U2[p,1:Int(N/2)]  = (Temp[1:2:end] + Temp[2:2:end])/2
                U2[p+1,1:Int(N/2)]= (Temp[1:2:end] - Temp[2:2:end])/2
                p +=2
            end
        end

        return U1, U2, N
    end
end

"""
MM = zeros(2*Scales - 1,T)
mid = Int(ceil(T/2))
for t = 1:Scales
    MM[2*t -1, mid-(2^(t-1))+1:2:mid-1+(2^(t-1))] = MSx[2^(t-1):(2^t)-1]
end

b = bar3(MM)

zlim([min(MSx(MSx>0)) max(MSx)]);
cmap = colormap('cool');  colormap([1 1 1; cmap])
for k = 1:length(b)
    zdata = b(k).ZData;
    b(k).CData = zdata;
    b(k).FaceColor = 'interp';
end
axis square, view([30 55])
yticks(1:2:(2*T -1));  yticklabels(0:T-1)
xticks(''); %     xticks(1:2:(2*p.Results.Scales - 1))
xlabel('Subtree Nodes','FontSize',12,'FontWeight','bold','Color',[7 54 66]/255)
ylabel('Scale Factor','FontSize',12,'FontWeight','bold','Color',[7 54 66]/255)
zlabel('Entropy','FontSize',12,'FontWeight','bold','Color',[7 54 66]/255)
title(sprintf('Hierarchical Multiscale (%s) Entropy',func2str(Y{1})),...
    'FontSize',16,'FontWeight','bold','Color',[7 54 66]/255)
"""

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
const Optional{T} = Union{Nothing, T}

struct AxisAspect
    aspect::Float32
end

struct DataAspect end

abstract type Ticks end

struct AutoLinearTicks{T<:Union{Int, Float32}} <: Ticks
    target::T
end

struct WilkinsonTicks <: Ticks
    k_ideal::Int
    k_min::Int
    k_max::Int
    Q::Vector{Tuple{Float64, Float64}}
    granularity_weight::Float64
    simplicity_weight::Float64
    coverage_weight::Float64
    niceness_weight::Float64
    min_px_dist::Float64
end

struct ManualTicks <: Ticks
    values::Vector{Float32}
    labels::Vector{String}
end

"""
    struct CustomTicks{F1<:Function, F2<:Function} <: Ticks

For the use of custom functions that compute tick values and labels in an `LAxis`.

- `f_tickvalues::F1`: A function that takes minimum_value, maximum_value, and pixelwidth as arguments and returns a `Float` array of tick values.
- `f_ticklabels::F2`: A function that takes in an array of `Float`s and returns an array of `String` labels.
"""
struct CustomTicks{F1<:Function, F2<:Function} <: Ticks
    f_tickvalues::F1
    f_ticklabels::F2
end

mutable struct LineAxis
    parent::Scene
    protrusion::Node{Float32}
    attributes::Attributes
    decorations::Dict{Symbol, Any}
    tickpositions::Node{Vector{Point2f0}}
    tickvalues::Node{Vector{Float32}}
    ticklabels::Node{Vector{String}}
end

abstract type LObject end

mutable struct LAxis <: AbstractPlotting.AbstractScene
    parent::Scene
    scene::Scene
    xaxislinks::Vector{LAxis}
    yaxislinks::Vector{LAxis}
    limits::Node{FRect2D}
    layoutobservables::LayoutObservables
    attributes::Attributes
    block_limit_linking::Node{Bool}
    decorations::Dict{Symbol, Any}
end

mutable struct LColorbar <: LObject
    parent::Scene
    layoutobservables::LayoutObservables
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

mutable struct LText <: LObject
    parent::Scene
    layoutobservables::LayoutObservables
    textobject::AbstractPlotting.Text
    attributes::Attributes
end

mutable struct LRect <: LObject
    parent::Scene
    layoutobservables::LayoutObservables
    rect::AbstractPlotting.Poly
    attributes::Attributes
end

struct LSlider <: LObject
    scene::Scene
    layoutobservables::LayoutObservables
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

struct LButton <: LObject
    scene::Scene
    layoutobservables::LayoutObservables
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

struct LToggle <: LObject
    scene::Scene
    layoutobservables::LayoutObservables
    attributes::Attributes
    decorations::Dict{Symbol, Any}
end

abstract type LegendElement end

struct LineElement <: LegendElement
    attributes::Attributes
end

struct MarkerElement <: LegendElement
    attributes::Attributes
end

struct PolyElement <: LegendElement
    attributes::Attributes
end

struct LegendEntry
    elements::Vector{LegendElement}
    attributes::Attributes
end

const EntryGroup = Tuple{Optional{String}, Vector{LegendEntry}}

struct LLegend <: LObject
    scene::Scene
    entrygroups::Node{Vector{EntryGroup}}
    layoutobservables::LayoutObservables
    attributes::Attributes
    decorations::Dict{Symbol, Any}
    entrytexts::Vector{LText}
    entryplots::Vector{Vector{AbstractPlot}}
end

const Indexables = Union{UnitRange, Int, Colon}

struct GridPosition
    layout::GridLayout
    rows::Indexables # this doesn't warrant type parameters I think
    cols::Indexables # as these objects will only be used briefly
end

struct LScene <: AbstractPlotting.AbstractScene
    scene::Scene
    attributes::Attributes
    layoutobservables::MakieLayout.LayoutObservables
end

ext: jl


text:
using ExpFamily
using Base.Test

srand(123)

Cd = [0.5, 1.0]
Pd = [2.0, 1.0]
m  = [1.0,-0.5]

dgNP = DiagGaussianNatParam(mean=m, cov=Cd)
dgMP = DiagGaussianMeanParam(mean=m, cov=Cd)

@test isapprox(dgNP.theta1, Pd.*m)
@test isapprox(dgNP.theta2, -Pd)
@test isapprox(dgMP.mu1, m)
@test isapprox(dgMP.mu2, (m.^2+Cd)/2)

gNP = full(dgNP)
gMP = full(dgMP)

@test isapprox(gNP.theta2, diagm(-Pd))
@test isapprox(gMP.mu2, diagm(m.^2+Cd)/2)

@test isvalid(dgNP)
@test isvalid(dgMP)

mpFromNP  = meanparam(dgNP)
npFromMP  = natparam(dgMP)

@test isapprox(mpFromNP.mu1, dgMP.mu1)
@test isapprox(mpFromNP.mu2, dgMP.mu2)
@test isapprox(npFromMP.theta1, dgNP.theta1)
@test isapprox(npFromMP.theta2, dgNP.theta2)

@test isapprox(mean(dgNP), m)
@test isapprox(mean(dgMP), m)
@test isapprox(cov(dgNP), Cd)
@test isapprox(cov(dgMP), Cd)
@test isapprox(var(dgMP), Cd)
@test isapprox(var(dgNP), Cd)
@test isapprox(std(dgMP), sqrt.(Cd))
@test isapprox(std(dgNP), sqrt.(Cd))

@test_throws DomainError DiagGaussianNatParam(mean=m,cov=-Cd,check=true)
@test_throws DomainError DiagGaussianMeanParam(mean=m,cov=-Cd,check=true)

dim = 4
dgNPones = ones(DiagGaussianNatParam, dim)
dgMPones = ones(DiagGaussianMeanParam, dim)

@test isapprox(mean(dgNPones), zeros(dim))
@test isapprox(cov(dgNPones),  ones(dim))
@test isapprox(mean(dgMPones), zeros(dim))
@test isapprox(cov(dgMPones),  ones(dim))

@test length(dgNPones) == dim
@test length(dgMPones) == dim

rNP = rand(dgNPones,5000)
rMP = rand(dgMPones,5000)

@test isapprox(mean(rNP,2), zeros(dim), atol=1e-1)
@test isapprox(mean(rMP,2), zeros(dim), atol=1e-1)
@test isapprox(cov(rNP,2),  eye(dim),  atol=1e-1)
@test isapprox(cov(rMP,2),  eye(dim),  atol=1e-1)

dgNP2 = dgNP+dgNP
dgMP2 = dgMP+dgMP

@test isapprox(dgNP2.theta1, 2dgNP.theta1)
@test isapprox(dgNP2.theta2, 2dgNP.theta2)
@test isapprox(dgMP2.mu1, 2dgMP.mu1)
@test isapprox(dgMP2.mu2, 2dgMP.mu2)

m1  = randn(dim)
m2  = randn(dim)
C1  = randn(dim,dim)
C1 *= C1'
C2  = randn(dim,dim)
C2 *= C2'

dgNPa = DiagGaussianNatParam(mean=m1, cov=diag(C1))
dgNPb = DiagGaussianNatParam(mean=m2, cov=diag(C2))
dgMPa = DiagGaussianMeanParam(mean=m1, cov=diag(C1))
dgMPb = DiagGaussianMeanParam(mean=m2, cov=diag(C2))

dgNPtest = dgNPa + dgNPb - natparam(dgMPb)
dgMPtest = dgMPa + dgMPb - meanparam(dgNPb)

@test isapprox(dgNPtest.theta1, dgNPa.theta1)
@test isapprox(dgNPtest.theta2, dgNPa.theta2)
@test isapprox(dgMPtest.mu1, dgMPa.mu1)
@test isapprox(dgMPtest.mu2, dgMPa.mu2)

@test isapprox(dgMPtest, dgMPa)

@test isapprox(DiagGaussianNatParam(vec(dgNPa),dim),dgNPa)
@test isapprox(DiagGaussianMeanParam(vec(dgMPa),dim),dgMPa)

@test isapprox(2dgNPa + dgNPa*2 - 3dgNPa, dgNPa)
@test isapprox(2dgMPa + dgMPa*2 - 3dgMPa, dgMPa)
@test isapprox(dgNPa/2 + dgNPa/2, dgNPa)
@test isapprox(dgMPa/2 + dgMPa/2, dgMPa)

nP = -[1.5;1e7]
m  = [1.0;-0.5]
lm = 1e-6

dgNPproj = project(DiagGaussianNatParam(-nP.*m, nP), minvar=lm)
dgMPproj = project(DiagGaussianMeanParam(mean=m, cov=-1.0./nP), minvar=lm)

nPthresh = -[1.5; 1./lm]

@test isapprox(dgNPproj.theta1, -nPthresh.*m)
@test isapprox(dgNPproj.theta2, -[1.5; 1./lm])
@test isapprox(dgMPproj.mu1, m)
@test isapprox(cov(dgMPproj), [1./1.5; lm])

x = randn(dim)

Ca = cov(dgNPa)
ma = mean(dgNPa)
l  = exp( -dot(x-ma, (x-ma)./Ca) / 2 )/sqrt( (2pi)^dim * prod(Ca) )
ll = log(l)

gNPa = GaussianNatParam(mean=ma, cov=diagm(Ca))

@test isapprox( loglik(dgNPa, x), ll )
@test isapprox( loglik(dgMPa, x), ll )

@test isapprox( gradloglik(dgNPa, x),  gradloglik(gNPa, x) )
@test isapprox( gradloglik(dgMPa, x), gradloglik(dgNPa, x) )


text:
mutable struct ConstantMean{T<:Real,O} <: PriorMean{T}
    C::Vector{T}
    opt::O
end

"""
    ConstantMean(c::Real = 1.0; opt=ADAM(0.01))

## Arguments
- `c::Real` : Constant value

Construct a prior mean with constant `c`
Optionally set an optimiser `opt` (`ADAM(0.01)` by default)
"""
function ConstantMean(c::T=1.0; opt=ADAM(0.01)) where {T<:Real}
    return ConstantMean{T,typeof(opt)}([c], opt)
end

function Base.show(io::IO, ::MIME"text/plain", μ₀::ConstantMean)
    return print(io, "Constant Mean Prior (c = ", only(μ₀.C), ")")
end

(μ::ConstantMean{T})(::Real) where {T<:Real} = only(μ.C)
(μ::ConstantMean{T})(x::AbstractVector) where {T<:Real} = fill(only(μ.C), length(x))

function init_priormean_state(hyperopt_state, μ₀::ConstantMean)
    μ₀_state = (; C=Optimisers.state(μ₀.opt, μ₀.C))
    return merge(hyperopt_state, (; μ₀_state))
end

function update!(μ₀::ConstantMean{T}, hyperopt_state, grad) where {T<:Real}
    μ₀_state = hyperopt_state.μ₀_state
    C, ΔC = Optimisers.apply(μ₀.opt, μ₀_state.C, μ₀.C, grad)
    μ₀.C .+= ΔC
    return merge(hyperopt_state, (; μ₀_state=(; C)))
end


text:
using Test

macro timeinclude(x)
    :(print($x); print("\t"); @time include($x))
end

@testset "Knet" begin

#                            commit 555 686 ee1 c6a 0b5 e3a 9.2 8.3 8.3 6cb 6cb 8.6 6cb 6cb
#                           machine dy3 rzr  ci rzr  sc  sc tig ai5 ai4 tr5 tr4 aws osx os4
@timeinclude("kptr.jl")          #   30  2   ?   4   3   8  16   1   1   0   0  20   0   0
@timeinclude("gpu.jl")           #   31  2   3   1   2   2   6   1   1   0   0   2   0   0
@timeinclude("distributions.jl") #   26  1   1   1   1   1   2   1   1   2   1   3   3   2
@timeinclude("dropout.jl")       #   31  3   6   1   2   8   5                   2
@timeinclude("gcnode.jl")        #   39  1
@timeinclude("jld.jl")           #   45  2  11  11   9  26
@timeinclude("statistics.jl")    #   48  7
@timeinclude("bmm.jl")           #   65  9  10   3   9
@timeinclude("serialize.jl")     #   70 20  16   2  11   1
@timeinclude("loss.jl")          #   74 17  32  13  19  20  10                   4
@timeinclude("cuarray.jl")       #   74 18  25
@timeinclude("update.jl")        #   78 64 316  61  82  60  61  29  26 100  22  72  25  23
@timeinclude("linalg.jl")        #   85 17  33  17  22  22  62  24  14  22   7  28  33  19
@timeinclude("batchnorm.jl")     #   92 19  32  14  23  22  93
@timeinclude("ops20.jl")         #   94
@timeinclude("karray.jl")        #  114 27  49  20  27  21  55  19  12   -   -  21   -   0
@timeinclude("conv.jl")          #  114 47  71  33  48  51 107  22  12  62  47  26  44  16
@timeinclude("rnn.jl")           #  126 43  95  52  54  41  81                  12
@timeinclude("reduction.jl")     #  137 42  67  37  52  55 106  40  21  29  11  57  55  29
@timeinclude("unary.jl")         #  174 76 172  86 103 103 122  42   6  36   4  56  67  11
@timeinclude("binary.jl")        #  234 89 124  64  75  80  56  34  19 491 119  51  53  25
#TODO include("data.jl")
#TODO include("hyperopt.jl")
#TODO include("progress.jl")
#TODO include("random.jl")
#TODO include("train.jl")
#TODO include("uva.jl")

end


text:
# tags used in <form></form>

# ----------------------- form block <form>
"""
    tag_form( content::String ; id::String = "", class::String = "", style::String = "" )

creates a <form>content</form> tag string.
please organize your contents before packaging them in a form block.
"""
function tag_form( content::String ; id::String = "", class::String = "", style::String = "" )
    return PairedHtmlTag( "form", id, class, style, String[], Dict(), content )
end # tag_form

# -------------------- typical tags generating functions
# input: 文本类型（用于一般输入）
"""
    tag_input_text( name::String ; id::String = "id_" * name, class::String = "", style::String = "", placeholder::String = "please type", flag_required::Bool = false, flag_readonly::Bool = false, flag_disabled::Bool = false )
"""
function tag_input_text( name::String ;
    id::String = "id_" * name, class::String = "", style::String = "",
    placeholder::String = "please type",
    flag_required::Bool = false, flag_readonly::Bool = false, flag_disabled::Bool = false )
    # prepare single attributes
    local singleattrs = [
        ( flag_required ? "required" : "" ),
        ( flag_readonly ? "readonly" : "" ),
        ( flag_disabled ? "disabled" : "" ),
    ]
    # prepare paired attributes
    local pairedattrs = Dict(
        "placeholder" => add_doublequotations(placeholder),
        "type" => add_doublequotations("text"),
    )
    # new a SingleHtmlTag
    return SingleHtmlTag("input", id, class, style, singleattrs, pairedattrs)
end # tag_input_text

# -------------
# input: 数字类型（用于年份等）
"""
    tag_input_number( name::String ; id::String = "id_" * name, class::String = "", style::String = "", num_range::Union{Vector{T},Tuple{T,T}} where T <: Real = (1920,2019), placeholder::String = "please type", flag_required::Bool = false, flag_readonly::Bool = false, flag_disabled::Bool = false )
"""
function tag_input_number( name::String ;
    id::String = "id_" * name, class::String = "", style::String = "",
    num_range::Union{Vector{T},Tuple{T,T}} where T <: Real = (1920,2019),
    placeholder::String = "type a number in $(num_range[1]) ~ $(num_range[2])",
    flag_required::Bool = false, flag_readonly::Bool = false, flag_disabled::Bool = false )
    # validation: num_range
    if isa(num_range,Vector); @assert( length(num_range) == 2, "num_range is allowed to have 2 elements" ); end;
    # prepare single attributes
    local singleattrs = [
        ( flag_required ? "required" : "" ),
        ( flag_readonly ? "readonly" : "" ),
        ( flag_disabled ? "disabled" : "" ),
    ]
    # prepare paired attributes
    local pairedattrs = Dict(
        "min" => string(num_range[1]),
        "max" => string(num_range[1]),
        "placeholder" => add_doublequotations(placeholder),
        "type" => add_doublequotations("number"),
    )
    # new a SingleHtmlTag string
    return SingleHtmlTag("input", id, class, style, singleattrs, pairedattrs)
end # tag_input_text

# -------------
# input: 单选题目(radio) NOTE: 返回的是一个<input /> string列表
"""
    tag_input_radio( name::String, options::Dict{T,String} where T <: Any ; label_id::String = "id_labelfor_" * name, label_class::String = "", label_style::String = "", id::String = "id_" * name, class::String = "", style::String = "", whichlabel_checked::Union{String,Nothing} = nothing )

1. :options should be a Dict whose keys are the contents displayed as label, and the values are the "value" attribute recorded in <input value="">
2. returns a **LIST** of PairedHtmlTag
"""
function tag_input_radio( name::String,
    options::Dict{T,String} where T <: Any ; # label text => radio (input) value
    label_id::String = "id_labelfor_" * name, label_class::String = "", label_style::String = "", # attributes of out-layer  tag
    id::String = "id_" * name, class::String = "", style::String = "", # attributes of inner <input /> tag
    whichlabel_checked::Union{String,Nothing} = nothing )  # which one (label text) (1st, 2nd ...) option should be checked by default
    # -------------
    # a tag list
    local res = PairedHtmlTag[]
    # make tags
    for labelstr in keys(options)
        local tmpSingleAttr4input::Vector{String} = whichlabel_checked == labelstr ? String["checked",] : String[]
        local in_radiotag = new_tag(SingleHtmlTag( "input",
            id, class, style, tmpSingleAttr4input,
            Dict(
                "type" => add_doublequotations("radio"),
                "value" => add_doublequotations(options[labelstr]),
            )
        ))
        local out_labtag = PairedHtmlTag( "label",
            label_id, label_class, label_style, [], Dict(),
            labelstr * in_radiotag )
        push!(res, out_labtag)
    end # labelstr
    return res::Vector{PairedHtmlTag}
end

# -------------
# input: 复选框(checkbox)
"""
    tag_input_checkbox( name::String, options::Dict{T,String} where T <: Any ; label_id::String = "id_labelfor_" * name, label_class::String = "", label_style::String = "", id::String = "id_" * name, class::String = "", style::String = "", whichlabel_checked::Union{String,Nothing,Vector{String}} = nothing )
"""
function tag_input_checkbox( name::String,
    options::Dict{T,String} where T <: Any ; # label text => radio (input) value
    label_id::String = "id_labelfor_" * name, label_class::String = "", label_style::String = "", # attributes of out-layer  tag
    id::String = "id_" * name, class::String = "", style::String = "", # attributes of inner <input /> tag
    whichlabel_checked::Union{String,Nothing,Vector{String}} = nothing )  # which (label text) (1st, 2nd ...) option(s) should be checked by default
    # -------------
    # a tag list
    local res = PairedHtmlTag[]
    # make tags
    for labelstr in keys(options)
        local tmpSingleAttr4input::Vector{String} = String[]
        if isa(whichlabel_checked, String)
            if whichlabel_checked == labelstr; push!(tmpSingleAttr4input, "checked"); end;
        elseif isa(whichlabel_checked, Vector)
            if labelstr in whichlabel_checked; push!(tmpSingleAttr4input, "checked"); end;
        end
        local in_radiotag = new_tag(SingleHtmlTag( "input",
            id, class, style, tmpSingleAttr4input,
            Dict(
                "type" => add_doublequotations("checkbox"),
                "value" => add_doublequotations(options[labelstr]),
            )
        ))
        local out_labtag = PairedHtmlTag( "label",
            label_id, label_class, label_style, [], Dict(),
            labelstr * in_radiotag )
        push!(res, out_labtag)
    end # labelstr
    return res::Vector{PairedHtmlTag}
end # tag_input_checkbox

# -------------
# select: 下拉框

# -------------
# datalist: 下拉框(用于<input>的list=属性，list=后面跟datalist的id)，即可复用的select

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
# This file includes portions from JuliaMath/Combinatorics.jl in modified form.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module IntPartitions

using GeneralizedBinomial, DataStructures

export ModuleIntPartitions
export IntegerPartitions, PartitionNumber
export PartOrder, byNumPart, byMaxPart
export PartitionCoefficientsByLength, PartitionCoefficientsByBiggestPart
export V000041, I072233, L072233, L036038, L078760, L005651, L262071, L292222
export L115621, L328917, V088887

#=
a(n) is the sum of the parts of the partition having Heinz number n. We define the Heinz number of a partition p = [p_1, p_2, ..., p_r] as Product_{j=1..r} (p_j-th prime) (concept used by Alois P. Heinz in A215366 as an "encoding" of a partition). For example, for the partition [1, 1, 2, 4, 10] we get 2*2*3*7*29 = 2436. Example: a(33) = 7 because the partition with Heinz number 33 = 3 * 11 is [2,5]. - Emeric Deutsch, May 19 2015
=#

"""

All integer partitions are listed by two orderings:
   IntegerPartitions(n, byNumPart)
   IntegerPartitions(n, byMaxPart)

Or restricted to partitions of length m:
   IntegerPartitions(n, m)

The partition coefficients, which are the multinomial coefficients applied to partitions, are given in both orderings (L036038, L078760).

The partition numbers and the number of partitions of n into k parts are given as PartitionNumber(n) and PartitionNumber(n, k), (V000041, L072233).

The sum of all partition coefficients of n is efficiently computed with L005651.

* V000041, V088887, I072233, L072233, L036038, L078760, L005651, L262071, L292222, L115621, L328917
"""
const ModuleIntPartitions = ""

"""

Conventional names for orderings of integer partitions of ``n``. byMaxPart means first order by the biggest part. byNumPart means first order by the number of parts. The secondary order is in both cases colexicographic.
"""
@enum PartOrder begin
    byNumPart
    byMaxPart
end
"Order the partitions by number of parts."
byNumPart
"Order the partitions by greatest part."
byMaxPart

"""

Return an iterator over all partitions of ``n``. List the parts of the partitions according to the ``PartitionOrder`` given as second parameter. By default the partition order is 'byNumPart'.
"""
function IntegerPartitions(n::Int, o::PartOrder = byNumPart)
    n < 0 && throw(DomainError(n, "n >= 0 required"))

    o == byNumPart && return PartitionsByLength(n)
    PartitionsByBiggestPart(n)
end

struct IntPartIterator
    n::Int
end

Base.length(p::IntPartIterator) = PartitionNumber(p.n)
Base.eltype(p::IntPartIterator) = Vector{Int}

function Base.iterate(p::IntPartIterator, xs = Int[])
    length(xs) == p.n && return
    xs = NextIntPartition(p.n, xs)
    (xs, xs)
end

function NextIntPartition(n, as)
    isempty(as) && return Int[n]
    xs = similar(as, 0)
    sizehint!(xs, length(as) + 1)

    for i = 1:length(as)-1
        if as[i+1] == 1
            x = as[i] - 1
            push!(xs, x)
            n -= x
            while n > x
                push!(xs, x)
                n -= x
            end
            push!(xs, n)
            return xs
        end
        push!(xs, as[i])
        n -= as[i]
    end
    push!(xs, as[end] - 1)
    push!(xs, 1)
    xs
end

"""

Return an iterator over all partitions of ``n``. List the parts of the partitions of ``n`` in graded colexicographic order.
"""
function PartitionsByBiggestPart(n::Int)
    n == 0 && return [Int[]]
    IntPartIterator(n)
end

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

struct IntPartIteratorByLength
    n::Int
    m::Int
end

Base.length(f::IntPartIteratorByLength) = PartitionNumber(f.n, f.m)
Base.eltype(f::IntPartIteratorByLength) = Vector{Int}

function Base.iterate(f::IntPartIteratorByLength, s::Vector{Int} = Int[])
    f.m <= f.n || return
    if !isempty(s)
        (f.m == 1 || s[1] - 1 <= s[end]) && return
    end

    xs = NextIntPartition(f.n, f.m, s)
    (xs, xs)
end

function NextIntPartition(n, m, bs)
    as = copy(bs)
    if isempty(as)
        # First iteration
        as = ones(Int, m)
        as[1] = n - m + 1
    elseif as[2] < as[1] - 1
        # Most common iteration
        as[1] -= 1
        as[2] += 1
    else
        # Iterate
        j = 0
        s = as[1] + as[2] - 1
        for jj ∈ 3:m
            j = jj
            as[jj] < as[1] - 1 && break
            s += as[jj]
        end
        x = as[j] += 1
        for k = j - 1:-1:2
            as[k] = x
            s -= x
        end
        as[1] = s
    end
    return as
end

"""

Return an iterator over all weakly descending lists with ``m`` integers that sum to ``n``. List the parts of the partitions of ``n`` in graded colexicographic order.
"""
function IntegerPartitions(n::Int, m::Int)
# Algorithm H from TAOCP 7.2.1.4.
    if m == 0
        n == 0 && return [Int[]]
        return Int[]
    end
    if (n >= 0 && m > n)
        return Int[]
    end

    if (n > 0 && m >= 0 && m <= n)
        return IntPartIteratorByLength(n, m)
    end

    throw(DomainError((n, m), "n >= 0 and 0 <= m <= n required"))
end

"""

Return an iterator over the IntegerPartitions of ``n`` ordered first by length then colexicographicaly. Chains the iterators IntegerPartitions(n, k).
"""
function PartitionsByLength(n::Int)
    n == 0 && return [Int[]]
    Iterators.flatten([IntegerPartitions(n, k) for k = 1:n])
end

const CachePartNum = Dict{Tuple{Int, Int}, Int}()

"""

Return the number of partitions of ``n`` into ``k`` parts. Cf. A072233.
"""
function PartitionNumber(n::Int, k::Int)
    (n == k || k == 1) && return 1
    (n <  k || k == 0) && return 0
    haskey(CachePartNum, (n, k)) && return CachePartNum[(n, k)]
    CachePartNum[(n, k)] =
        PartitionNumber(n - 1, k - 1) + PartitionNumber(n - k, k)
end

"""

Return the number of partitions of ``n``. Cf. A000041.
"""
function PartitionNumber(n::Int)
    haskey(CachePartNum, (n, -1)) && return CachePartNum[(n, -1)]
    CachePartNum[(n, -1)] = sum(PartitionNumber(n, k) for k ∈ 0:n)
end

"""

Return the number of partitions of ``n``.
"""
V000041(n) = PartitionNumber(n)

"""

Return an iterator over the list of partitions of ``n`` into ``k`` parts.
"""
I072233(n) = (PartitionNumber(n, k) for k ∈ 0:n)

"""

Return a list of the number of partitions of ``n`` into ``k`` parts.
"""
L072233(n) = collect(I072233(n))

"""

Return the partition coefficients of ``n``, first ordered by length.
"""
PartitionCoefficientsByLength(n) = [
    Multinomial(p) for p ∈ IntegerPartitions(n, byNumPart)]

"""

Return the partition coefficients of ``n``, first ordered by length.
"""
L036038(n) = PartitionCoefficientsByLength(n)

"""

Return the partition coefficients of ``n``, first ordered by biggest part.
"""
PartitionCoefficientsByBiggestPart(n) = [
    Multinomial(p) for p ∈ IntegerPartitions(n, byMaxPart)]

"""

Return the partition coefficients of ``n``, first ordered by biggest part.
"""
L078760(n) = PartitionCoefficientsByBiggestPart(n)

const CacheSumPartCoeff = Dict{Tuple{Int, Int}, Int}()

function H(n::Int, i::Int)
    (n == 0 || i == 1) && return factorial(n)
    if (ret = get(CacheSumPartCoeff, n, 0)) > 0
        return ret
    end

    ret = H(n, i-1) + binomial(n, i)*H(n-i, min(n-i, i))
    CacheSumPartCoeff[n, i] = ret
end

"""

Sum of all partition coefficients of n.
"""
L005651(len::Int) = [H(n, n) for n ∈ 0:len-1]

"""

Return the number of ordered partitions of an n-set with nondecreasing block sizes and maximal block size equal to k.
"""
function L262071(n::Int)
    n == 0 && return [1]
    A = zeros(Int, n + 1)
    for p ∈ IntegerPartitions(n, byMaxPart)
        A[1+p[1]] += Multinomial(p)
    end
    return A
end

"""

Return the triangle of multinomial coefficients of ``n``. [Defined as in the 'Handbook of Mathematical Functions', p. 831, as ``M1``.]
"""
function L292222(n)
    [sum(Multinomial(p) for p ∈ IntegerPartitions(n, k)) for k ∈ 1:n]
end

"""
Return the signature of partitions in Hindenburg order.
"""
function L115621(n)
    h(p) = sort(collect(values(counter(p))))
    [h(p) for p ∈ IntegerPartitions(n)]
end

"""
Return the types of signatures of partitions of n, ordered firstly by decreasing greatest parts, then decreasing sum of parts, then by increasing number of parts.
"""
function L328917(n)
    n == 0 && return [[0]]
    h(p) = sort(collect(values(counter(p))), rev=true)
    sort(unique([h(p) for p ∈ IntegerPartitions(n)]), rev=true)
end

"""
Return the number of types of signatures of partitions of n.
"""
function V088887(n)
    h(p) = sort(collect(values(counter(p))), rev=true)
    length(unique([h(p) for p ∈ IntegerPartitions(n)]))
end

#START-TEST-########################################################

using Test, SeqUtils

function test()
    @testset "IntegerPartitions" begin

        @test IntegerPartitions(0) == [[]]

        @test collect(IntegerPartitions(0, 0)) == [[]]
        @test collect(IntegerPartitions(0, 1)) == []

        @test collect(IntegerPartitions(1, 0)) == []
        @test collect(IntegerPartitions(1, 1)) == [[1]]
        @test collect(IntegerPartitions(1, 2)) == []

        @test collect(IntegerPartitions(2, 0)) == []
        @test collect(IntegerPartitions(2, 1)) == [[2]]
        @test collect(IntegerPartitions(2, 2)) == [[1, 1]]
        @test collect(IntegerPartitions(2, 3)) == []

        @test [length(b) for b ∈ IntegerPartitions(0)] == [0]
        @test [length(b) for b ∈ IntegerPartitions(1)] == [1]
        @test length(IntegerPartitions(0)) == 1
        @test length(collect(IntegerPartitions(1))) == 1
        @test_throws DomainError IntegerPartitions(-1)

        @test collect(IntegerPartitions(5)) == reverse(Any[[1, 1, 1, 1, 1],
            [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]])

        for n ∈ 0:7
            pn = PartitionNumber(n)
            @test length(collect(IntegerPartitions(n, byMaxPart))) == pn
            @test length(collect(IntegerPartitions(n, byNumPart))) == pn
        end

        for n ∈ 0:7, k ∈ 0:n
            pn = PartitionNumber(n, k)
            @test length(collect(IntegerPartitions(n, k))) == pn
        end

        @test L005651(10) == [1, 1, 3, 10, 47, 246, 1602, 11481, 95503, 871030]
        @test [V000041(n) for n ∈ 0:8] == [1, 1, 2, 3, 5, 7, 11, 15, 22]
        @test L072233(7) == [0, 1, 3, 4, 3, 2, 1, 1]
        @test L036038(6) == [1, 6, 15, 20, 30, 60, 90, 120, 180, 360, 720]
        @test L078760(6) == [1, 6, 15, 30, 20, 60, 120, 90, 180, 360, 720]
        @test L005651(10) == [1, 1, 3, 10, 47, 246, 1602, 11481, 95503, 871030]
        @test L262071(6) == [0, 720, 630, 200, 45, 6, 1]
        @test L292222(6) == [1, 41, 180, 300, 360, 720]
    end
end

function demo()
    println("\n-- all partitions of 10")
    for p ∈ IntegerPartitions(10)
        Println(p)
    end

    for n ∈ 0:7
        println("\n\n-- $n  --")

        println("\n-- by biggest part")
        for p ∈ IntegerPartitions(n, byMaxPart)
            Println(p)
        end

        println("\n-- by length (all)")
        for p ∈ IntegerPartitions(n, byNumPart)
            Println(p)
        end

        println("\n-- selected lengths")
        for k ∈ 0:2:n
            P = IntegerPartitions(n, k)
            for p ∈ P Println(p) end
        end
    end

    println("\n-- edge cases:")
    println("-- n=0, k=0 -> P(n, k) = [[]]")
    println("-- n>0, k=0 -> P(n, k) = []")
    println("-- n>0, k>n -> P(n, k) = []")

    for n ∈ [0,1,2,6]
        println("\n-- ", n)
        for k ∈ 0:n
            P = collect(IntegerPartitions(n, k))
            print("collect ", (n, k), " -> "); Println(P)
            if isempty(IntegerPartitions(n, k))
                println((n, k), " -> no partition")
            else
                for p ∈ IntegerPartitions(n, k)
                    println((n, k), " -> ", p)
                end
            end
        end

        k = n + 1
        P = collect(IntegerPartitions(n, k))
        print("collect ", (n, k), " -> "); Println(P)

        if isempty(IntegerPartitions(n, k))
            println((n, k), " -> no partition")
        else
            for p ∈ IntegerPartitions(n, k)
                println((n, k), " -> ", p)
            end
        end
    end

    println("\n-- triangle of partition numbers")
    for n ∈ 0:7
        println([PartitionNumber(n, k) for k ∈ 0:n])
    end

    println("\n-- number of partitions")
    println([PartitionNumber(n) for n ∈ 0:12])
    println()

    println("\n-- partition coefficients by biggest part")
    for n ∈ 0:6
        PartitionCoefficientsByBiggestPart(n) |> println
    end
    println()

    println("\n-- partition coefficients by length")
    for n ∈ 0:6
        PartitionCoefficientsByLength(n) |> println
    end
    println()

    println("\n-- Sum of all partition coefficients of n")
    L005651(10) |> println
    println()

    println("\n-- Signature of partitions in Hindenburg order.")
    for n ∈ 1:6
        L115621(n) |> Println
    end

    println("\n-- Types of partitions (counted by V088887).")
    for n ∈ 0:6
        print(n, " : "); Println(L328917(n))
    end
    println()
end

"""
for n ∈ 1:50, p ∈ IntegerPartitions(n, byNumPart) p end
    0.437520 seconds (3.89 M allocations: 329.296 MiB, 26.45% gc time)
for n ∈ 1:50, p ∈ IntegerPartitions(n, byMaxPart) p end
    0.509802 seconds (2.59 M allocations: 249.845 MiB, 20.57% gc time)
"""
function perf()
    println("\n-- performance test")
    @time for n ∈ 1:50, p ∈ IntegerPartitions(n, byNumPart) p end
    @time for n ∈ 1:50, p ∈ IntegerPartitions(n, byMaxPart) p end
end

function main()
    demo()
    test()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using PyCall
using Conda

# Install needed Python packages
Conda.add("xarray")
Conda.add_channel("conda-forge")
Conda.add("xgcm", channel="conda-forge")
Conda.add("netcdf4")

# Needed to import local Python modules like sose_data
# See: https://github.com/JuliaPy/PyCall.jl/issues/48
py"""
import sys
sys.path.insert(0, ".")
"""

sose = pyimport("sose_data")

function load_sose_data(sose_dir, lat, lon, day_offset, n_days, grid, buoyancy, coriolis; array_type=Array{Float64})

    @info "Opening SOSE datasets..."

    ds2 = sose.open_sose_2d_datasets(sose_dir)
    ds3 = sose.open_sose_3d_datasets(sose_dir)

    times = sose.get_times(ds2)

    zC = ds3.Z.values
    zF = ds3.Zl.values

    @info "Extracting surface forcings and profiles..."

    τx  = sose.get_scalar_time_series(ds2, "oceTAUX",  lat, lon, day_offset, n_days) |> array_type
    τy  = sose.get_scalar_time_series(ds2, "oceTAUY",  lat, lon, day_offset, n_days) |> array_type
    Qθ  = sose.get_scalar_time_series(ds2, "oceQnet",  lat, lon, day_offset, n_days) |> array_type
    Qs  = sose.get_scalar_time_series(ds2, "oceFWflx", lat, lon, day_offset, n_days) |> array_type
    mld = sose.get_scalar_time_series(ds2, "BLGMLD",   lat, lon, day_offset, n_days) |> array_type

    U = sose.get_profile_time_series(ds3, "UVEL",   lat, lon, day_offset, n_days) |> array_type
    V = sose.get_profile_time_series(ds3, "VVEL",   lat, lon, day_offset, n_days) |> array_type
    Θ = sose.get_profile_time_series(ds3, "THETA",  lat, lon, day_offset, n_days) |> array_type
    S = sose.get_profile_time_series(ds3, "SALT",   lat, lon, day_offset, n_days) |> array_type
    N = sose.get_profile_time_series(ds3, "DRHODR", lat, lon, day_offset, n_days) |> array_type

    # Nominal values for α, β to compute geostrophic velocities
    # FIXME: Linear equation of state does not apply!
    # FIXME: Diagnose buoyancy from stratification!
    α = 1.67e-4
    β = 7.80e-4

    @info "Diagnosing geostrophic velocities..."

    Ugeo, Vgeo = sose.compute_geostrophic_velocities(ds3, lat, lon, day_offset, n_days, grid.zF, α, β,
                                                     buoyancy.gravitational_acceleration, coriolis.f)

    Ugeo = array_type(Ugeo)
    Vgeo = array_type(Vgeo)

    ds2.close()
    ds3.close()

    sose_vertical_grid = (; zC, zF)
    surface_forcings = (; τx, τy, Qθ, Qs, mld)
    profiles = (; U, V, Θ, S, N, Ugeo, Vgeo)

    return times, sose_vertical_grid, surface_forcings, profiles
end


text:
using LinearAlgebra
using Distributions

"""
    mv_gaussian(μ::AbstractArray{Float32,1}
        Σ::AbstractArray{Float32,2}
        invΣ::AbstractArray{Float32,2}
        logdetΣ::Float32
        invChol::UpperTriangular)
[Multivariate Normal Distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution)
"""
struct mv_gaussian <: distribution_sample
    μ::AbstractArray{Float32,1}
    Σ::AbstractArray{Float32,2}
    invΣ::AbstractArray{Float32,2}
    logdetΣ::Float32
    invChol::UpperTriangular
end

function log_likelihood!(r::AbstractArray,x::AbstractMatrix, distribution_sample::mv_gaussian , group::Int64 = -1)
     z = x .- distribution_sample.μ
    dcolwise_dot!(r,z, distribution_sample.invΣ * z)
    r .= -((length(distribution_sample.Σ) * Float32(log(2π)) + distribution_sample.logdetΣ)/2) .-r/2
end


text:
module FuncKinds

# Many of these are written elaborately to defeat inference

function hasfoo(list)   # a test for anonymous functions
    hf = false
    hf = map(list) do item
        if isa(item, AbstractString)
            (str->occursin("foo", str))(item)
        else
            false
        end
    end
    return any(hf)
end

## A test for keyword functions
const randnums = Any[1, 1.0, UInt(14), 2.0f0]

@noinline returnsrandom() = randnums[rand(1:length(randnums))]

@noinline function haskw(x, y; a="hello", b=1, c=returnsrandom())
    if isa(b, Integer)
        return cos(rand()) + c + x + y
    elseif isa(b, AbstractFloat)
        s = 0.0
        for i = 1:rand(1:10)
            s += log(rand())
        end
        return s + c + x + y
    end
    return "string"
end

function callhaskw()
    ret = Any[]
    for i = 1:5
        push!(ret, haskw(returnsrandom(), returnsrandom()))
    end
    push!(ret, haskw(returnsrandom(), returnsrandom(); b = 2.0))
    return ret
end

@generated function gen(x::T) where T
    Tbigger = T == Float32 ? Float64 : BigFloat
    :(convert($Tbigger, x))
end

function gen2(x::Int, y)
    if @generated
        return y <: Integer ? :(x*y) : :(x+y)
    else
        return 2x+3y
    end
end

function hasinner(x, y)
    inner(z) = 2z

    s = 0
    for i = 1:10
        s += inner(returnsrandom())
    end
    return s
end

# Two kwarg generated functions; one will be called from the no-kw call, the other from a kwcall
@generated function genkw1(; b=2)
    :(string(typeof($b)))
end
@generated function genkw2(; b=2)
    :(string(typeof($b)))
end

# Function styles from JuliaInterpreter
f1(x::Int) = 1
f1(x) = 2
# where signatures
f2(x::T) where T = -1
f2(x::T) where T<:Integer = T
f2(x::T) where Unsigned<:T<:Real = 0
f2(x::V) where V<:SubArray{T} where T = 2
f2(x::V) where V<:Array{T,N} where {T,N} = 3
f2(x::V) where V<:Base.ReshapedArray{T,N} where T where N = 4
# Varargs
f3(x::Int, y...) = 1
f3(x::Int, y::Symbol...) = 2
f3(x::T, y::U...) where {T<:Integer,U} = U
f3(x::Array{Float64,K}, y::Vararg{Symbol,K}) where K = K
# Default args
f4(x, y::Int=0) = 2
f4(x::UInt, y="hello", z::Int=0) = 3
f4(x::Array{Float64,K}, y::Int=0) where K = K
# Keyword args
f5(x::Int8; y=0) = y
f5(x::Int16; y::Int=0) = 2
f5(x::Int32; y="hello", z::Int=0) = 3
f5(x::Int64;) = 4
f5(x::Array{Float64,K}; y::Int=0) where K = K
# Default and keyword args
f6(x, y="hello"; z::Int=0) = 1

end


text:
# GaussianVectors.jl

# Define and simulate normal vectors.
struct GaussianVector{R<:AbstractArray{Float64, 1}, S<:AbstractArray{Float64},
		T<:AbstractArray{Float64}}
	mean::R
	Σ::S
	var::T
	length::Int
	numberofGaussians::Int

	function GaussianVector(mean::R, Σ::S) where {R<:AbstractArray{Float64, 1}, S<:AbstractArray{Float64}}
		length = Base.length(mean)
		size(Σ,1) == length || throw(DimensionMismatch("Dimensions of Σ and μ don't match."))
		numberofGaussians = size(Σ,2)
		var = Σ * Σ'
		new{R, S, typeof(var)}(mean, Σ, var, length, numberofGaussians)
	end
end

GaussianVector(Σ::T) where T<:AbstractArray{Float64} = GaussianVector(spzeros(size(Σ,1)), Σ)

rand(d::GaussianVector) = d.mean + d.Σ * randn((d.numberofGaussians, 1))

mean(d::GaussianVector) = d.mean

var(d::GaussianVector) = d.var

length(d::GaussianVector) = d.length


text:
module SDPT3

using SparseArrays
using MATLAB

export sdpt3

#function dim(c::Char, n::Float64)
#    if c == 's'
#        return div(n * (n + 1), 2)
#    else
#        return n
#    end
#end

# See Solver/sqlparameters.m
const ALLOWED_OPTIONS = [
    "vers",
    "gam",
    "predcorr",
    "expon",
    "gaptol",
    "inftol",
    "steptol",
    "maxit",
    "printlevel",
    "stoplevel",
    "scale_data",
    "rmdepconstr",
    "smallblkdim",
    "parbarrier",
    "schurfun",
    "schurfun_par"
]

# `SparseMatrixCSC` is returned in SumOfSquares.jl test `sos_horn`
_array(x::AbstractMatrix) = x
_array(x::Vector) = x
_array(x::Float64) = [x]

# TODO log in objective, OPTION, initial iterates X0, y0, Z0
# Solve the primal/dual pair
# min c'x,      max b'y
# s.t. Ax = b,   c - A'x ∈ K
#       x ∈ K
function sdpt3(blk::Matrix,
               At::Vector{<:Union{Matrix{Float64}, SparseMatrixCSC{Float64}}},
               C::Vector{<:Union{Matrix{Float64}, Vector{Float64}}},
               b::Vector{Float64}; kws...)
               #C::Vector{<:Union{Vector{Float64}, SparseVector{Float64}}}, b::Vector{Float64})
    options = Dict{String, Any}(string(key) => value for (key, value) in kws)
    @assert all(i -> size(At[i], 2) == length(b), 1:length(At))
    @assert length(At) == size(blk, 1)
    #@assert all(i -> size(A[i], 1) == dim(blk[i, 1], blk[i, 2]), 1:length(A))
    #@assert all(i -> length(C[i], 1) == dim(blk[i, 1], blk[i, 2]), 1:length(A))
    # There are 6 output arguments so we use `6` below
    obj, X, y, Z, info, runhist = mxcall(:sdpt3, 6, blk, At, C, b, options)
    return obj, _array.(X), _array(y), _array.(Z), info, runhist
end

include("MOI_wrapper.jl")

end # module


text:
calcMidpoints(edges::AbstractVector) = Float64[0.5 * (edges[i] + edges[i+1]) for i in 1:length(edges)-1]

"Make histogram-like bins of data"
function binData(data, nbins)
  lo, hi = extrema(data)
  edges = collect(linspace(lo, hi, nbins+1))
  midpoints = calcMidpoints(edges)
  buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges)))-1 for x in data]
  counts = zeros(Int, length(midpoints))
  for b in buckets
    counts[b] += 1
  end
  edges, midpoints, buckets, counts
end

"""
A hacky replacement for a histogram when the backend doesn't support histograms directly.
Convert it into a bar chart with the appropriate x/y values.
"""
function histogramHack(; kw...)
  d = KW(kw)

  # we assume that the y kwarg is set with the data to be binned, and nbins is also defined
  edges, midpoints, buckets, counts = binData(d[:y], d[:bins])
  d[:x] = midpoints
  d[:y] = float(counts)
  d[:linetype] = :bar
  d[:fillrange] = d[:fillrange] == nothing ? 0.0 : d[:fillrange]
  d
end

"""
A hacky replacement for a bar graph when the backend doesn't support bars directly.
Convert it into a line chart with fillrange set.
"""
function barHack(; kw...)
  d = KW(kw)
  midpoints = d[:x]
  heights = d[:y]
  fillrange = d[:fillrange] == nothing ? 0.0 : d[:fillrange]

  # estimate the edges
  dists = diff(midpoints) * 0.5
  edges = zeros(length(midpoints)+1)
  for i in 1:length(edges)
    if i == 1
      edge = midpoints[1] - dists[1]
    elseif i == length(edges)
      edge = midpoints[i-1] + dists[i-2]
    else
      edge = midpoints[i-1] + dists[i-1]
    end
    edges[i] = edge
  end

  x = Float64[]
  y = Float64[]
  for i in 1:length(heights)
    e1, e2 = edges[i:i+1]
    append!(x, [e1, e1, e2, e2])
    append!(y, [fillrange, heights[i], heights[i], fillrange])
  end

  d[:x] = x
  d[:y] = y
  d[:linetype] = :path
  d[:fillrange] = fillrange
  d
end

"""
A hacky replacement for a sticks graph when the backend doesn't support sticks directly.
Convert it into a line chart that traces the sticks, and a scatter that sets markers at the points.
"""
function sticksHack(; kw...)
  dLine = KW(kw)
  dScatter = copy(dLine)

  # these are the line vertices
  x = Float64[]
  y = Float64[]
  fillrange = dLine[:fillrange] == nothing ? 0.0 : dLine[:fillrange]

  # calculate the vertices
  yScatter = dScatter[:y]
  for (i,xi) in enumerate(dScatter[:x])
    yi = yScatter[i]
    for j in 1:3 push!(x, xi) end
    append!(y, [fillrange, yScatter[i], fillrange])
  end

  # change the line args
  dLine[:x] = x
  dLine[:y] = y
  dLine[:linetype] = :path
  dLine[:markershape] = :none
  dLine[:fillrange] = nothing

  # change the scatter args
  dScatter[:linetype] = :none

  dLine, dScatter
end

function regressionXY(x, y)
  # regress
  β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)

  # make a line segment
  regx = [minimum(x), maximum(x)]
  regy = β * regx + α
  regx, regy
end

function replace_image_with_heatmap{T<:Colorant}(z::Array{T})
    @show T, size(z)
    n, m = size(z)
    # idx = 0
    colors = ColorGradient(vec(z))
    newz = reshape(linspace(0, 1, n*m), n, m)
    newz, colors
    # newz = zeros(n, m)
    # for i=1:n, j=1:m
    #     push!(colors, T(z[i,j]...))
    #     newz[i,j] = idx / (n*m-1)
    #     idx += 1
    # end
    # newz, ColorGradient(colors)
end

function imageHack(d::KW)
    :heatmap in supportedTypes() || error("Neither :image or :heatmap are supported!")
    d[:linetype] = :heatmap
    d[:z], d[:fillcolor] = replace_image_with_heatmap(d[:z].surf)
end
# ---------------------------------------------------------------
# ------------------------------------------------------------------------------------

nop() = nothing
notimpl() = error("This has not been implemented yet")

get_mod(v::AVec, idx::Int) = v[mod1(idx, length(v))]
get_mod(v::AMat, idx::Int) = size(v,1) == 1 ? v[1, mod1(idx, size(v,2))] : v[:, mod1(idx, size(v,2))]
get_mod(v, idx::Int)       = v

makevec(v::AVec) = v
makevec{T}(v::T) = T[v]

"duplicate a single value, or pass the 2-tuple through"
maketuple(x::Real)                     = (x,x)
maketuple{T,S}(x::@compat(Tuple{T,S})) = x

mapFuncOrFuncs(f::Function, u::AVec)        = map(f, u)
mapFuncOrFuncs(fs::AVec{Function}, u::AVec) = [map(f, u) for f in fs]

unzip{T,S}(xy::AVec{Tuple{T,S}})              = [x[1] for x in xy], [y[2] for y in xy]
unzip{T,S,R}(xyz::AVec{Tuple{T,S,R}})         = [x[1] for x in xyz], [y[2] for y in xyz], [z[3] for z in xyz]
unzip{T}(xy::AVec{FixedSizeArrays.Vec{2,T}})  = T[x[1] for x in xy], T[y[2] for y in xy]
unzip{T}(xy::FixedSizeArrays.Vec{2,T})        = T[xy[1]], T[xy[2]]
unzip{T}(xyz::AVec{FixedSizeArrays.Vec{3,T}}) = T[x[1] for x in xyz], T[y[2] for y in xyz], T[z[3] for z in xyz]
unzip{T}(xyz::FixedSizeArrays.Vec{3,T})       = T[xyz[1]], T[xyz[2]], T[xyz[3]]

# given 2-element lims and a vector of data x, widen lims to account for the extrema of x
function _expand_limits(lims, x)
  try
    e1, e2 = extrema(x)
    lims[1] = min(lims[1], e1)
    lims[2] = max(lims[2], e2)
  # catch err
  #   warn(err)
  end
  nothing
end

# if the type exists in a list, replace the first occurence.  otherwise add it to the end
function addOrReplace(v::AbstractVector, t::DataType, args...; kw...)
    for (i,vi) in enumerate(v)
        if isa(vi, t)
            v[i] = t(args...; kw...)
            return
        end
    end
    push!(v, t(args...; kw...))
    return
end

function replaceType(vec, val)
  filter!(x -> !isa(x, typeof(val)), vec)
  push!(vec, val)
end

function replaceAlias!(d::KW, k::Symbol, aliases::KW)
  if haskey(aliases, k)
    d[aliases[k]] = pop!(d, k)
  end
end

function replaceAliases!(d::KW, aliases::KW)
  ks = collect(keys(d))
  for k in ks
      replaceAlias!(d, k, aliases)
    # if haskey(aliases, k)
    #   d[aliases[k]] = d[k]
    #   delete!(d, k)
    # end
  end
end

createSegments(z) = collect(repmat(z',2,1))[2:end]

Base.first(c::Colorant) = c
Base.first(x::Symbol) = x

sortedkeys(d::Dict) = sort(collect(keys(d)))

"create an (n+1) list of the outsides of heatmap rectangles"
function heatmap_edges(v::AVec)
  vmin, vmax = extrema(v)
  extra = 0.5 * (vmax-vmin) / (length(v)-1)
  vcat(vmin-extra, 0.5 * (v[1:end-1] + v[2:end]), vmax+extra)
end

function fakedata(sz...)
  y = zeros(sz...)
  for r in 2:size(y,1)
    y[r,:] = 0.95 * y[r-1,:] + randn(size(y,2))'
  end
  y
end

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited
isatom() = isdefined(Main, :Atom) && Main.Atom.isconnected()

istuple(::Tuple) = true
istuple(::Any)   = false
isvector(::AVec) = true
isvector(::Any)  = false
ismatrix(::AMat) = true
ismatrix(::Any)  = false
isscalar(::Real) = true
isscalar(::Any)  = false

isvertical(d::KW) = get(d, :orientation, :vertical) in (:vertical, :v, :vert)

# ticksType{T<:Real,S<:Real}(ticks::@compat(Tuple{T,S})) = :limits
ticksType{T<:Real}(ticks::AVec{T})                      = :ticks
ticksType{T<:AbstractString}(ticks::AVec{T})            = :labels
ticksType{T<:AVec,S<:AVec}(ticks::@compat(Tuple{T,S}))  = :ticks_and_labels
ticksType(ticks)                                        = :invalid

limsType{T<:Real,S<:Real}(lims::@compat(Tuple{T,S}))    = :limits
limsType(lims::Symbol)                                  = lims == :auto ? :auto : :invalid
limsType(lims)                                          = :invalid

axis_symbol(letter, postfix) = symbol(letter * postfix)
axis_symbols(letter, postfix...) = map(s -> axis_symbol(letter, s), postfix)

Base.convert{T<:Real}(::Type{Vector{T}}, rng::Range{T})         = T[x for x in rng]
Base.convert{T<:Real,S<:Real}(::Type{Vector{T}}, rng::Range{S}) = T[x for x in rng]

Base.merge(a::AbstractVector, b::AbstractVector) = sort(unique(vcat(a,b)))

nanpush!(a::AbstractVector, b) = (push!(a, NaN); push!(a, b))
nanappend!(a::AbstractVector, b) = (push!(a, NaN); append!(a, b))

# given an array of discrete values, turn it into an array of indices of the unique values
# returns the array of indices (znew) and a vector of unique values (vals)
function indices_and_unique_values(z::AbstractArray)
    vals = sort(unique(z))
    vmap = Dict([(v,i) for (i,v) in enumerate(vals)])
    newz = map(zi -> vmap[zi], z)
    newz, vals
end

# this is a helper function to determine whether we need to transpose a surface matrix.
# it depends on whether the backend matches rows to x (transpose_on_match == true) or vice versa
# for example: PyPlot sends rows to y, so transpose_on_match should be true
function transpose_z(d::KW, z, transpose_on_match::Bool = true)
    if d[:match_dimensions] == transpose_on_match
        z'
    else
        z
    end
end

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

wraptuple(x::@compat(Tuple)) = x
wraptuple(x) = (x,)

trueOrAllTrue(f::Function, x::AbstractArray) = all(f, x)
trueOrAllTrue(f::Function, x) = f(x)

allLineTypes(arg)   = trueOrAllTrue(a -> get(_typeAliases, a, a) in _allTypes, arg)
allStyles(arg)      = trueOrAllTrue(a -> get(_styleAliases, a, a) in _allStyles, arg)
allShapes(arg)      = trueOrAllTrue(a -> get(_markerAliases, a, a) in _allMarkers, arg) ||
                        trueOrAllTrue(a -> isa(a, Shape), arg)
allAlphas(arg)      = trueOrAllTrue(a -> (typeof(a) <: Real && a > 0 && a < 1) ||
                        (typeof(a) <: AbstractFloat && (a == zero(typeof(a)) || a == one(typeof(a)))), arg)
allReals(arg)       = trueOrAllTrue(a -> typeof(a) <: Real, arg)
allFunctions(arg)   = trueOrAllTrue(a -> isa(a, Function), arg)

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

"""
Allows temporary setting of backend and defaults for Plots. Settings apply only for the `do` block.  Example:
```
with(:gadfly, size=(400,400), type=:hist) do
  plot(rand(10))
  plot(rand(10))
end
```
"""
function with(f::Function, args...; kw...)
    newdefs = KW(kw)

    if :canvas in args
        newdefs[:xticks] = nothing
        newdefs[:yticks] = nothing
        newdefs[:grid] = false
        newdefs[:legend] = false
    end

  # dict to store old and new keyword args for anything that changes
  olddefs = KW()
  for k in keys(newdefs)
    olddefs[k] = default(k)
  end

  # save the backend
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end
  oldbackend = CURRENT_BACKEND.sym

  for arg in args

    # change backend?
    if arg in backends()
      backend(arg)
    end

    # # TODO: generalize this strategy to allow args as much as possible
    # #       as in:  with(:gadfly, :scatter, :legend, :grid) do; ...; end
    # # TODO: can we generalize this enough to also do something similar in the plot commands??

    # k = :linetype
    # if arg in _allTypes
    #   olddefs[k] = default(k)
    #   newdefs[k] = arg
    # elseif haskey(_typeAliases, arg)
    #   olddefs[k] = default(k)
    #   newdefs[k] = _typeAliases[arg]
    # end

    k = :legend
    if arg in (k, :leg)
      olddefs[k] = default(k)
      newdefs[k] = true
    end

    k = :grid
    if arg == k
      olddefs[k] = default(k)
      newdefs[k] = true
    end
  end

  # display(olddefs)
  # display(newdefs)

  # now set all those defaults
  default(; newdefs...)

  # call the function
  ret = f()

  # put the defaults back
  default(; olddefs...)

  # revert the backend
  if CURRENT_BACKEND.sym != oldbackend
    backend(oldbackend)
  end

  # return the result of the function
  ret
end

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

type DebugMode
  on::Bool
end
const _debugMode = DebugMode(false)

function debugplots(on = true)
  _debugMode.on = on
end

debugshow(x) = show(x)
debugshow(x::AbstractArray) = print(summary(x))

function dumpdict(d::KW, prefix = "", alwaysshow = false)
  _debugMode.on || alwaysshow || return
  println()
  if prefix != ""
    println(prefix, ":")
  end
  for k in sort(collect(keys(d)))
    @printf("%14s: ", k)
    debugshow(d[k])
    println()
  end
  println()
end

function dumpcallstack()
  error()  # well... you wanted the stacktrace, didn't you?!?
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------
# used in updating an existing series

extendSeriesByOne(v::UnitRange{Int}, n::Int = 1) = isempty(v) ? (1:n) : (minimum(v):maximum(v)+n)
extendSeriesByOne(v::AVec, n::Integer = 1)       = isempty(v) ? (1:n) : vcat(v, (1:n) + maximum(v))
extendSeriesData{T}(v::Range{T}, z::Real)        = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::Range{T}, z::AVec)        = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::AVec{T}, z::Real)         = (push!(v, convert(T, z)); v)
extendSeriesData{T}(v::AVec{T}, z::AVec)         = (append!(v, convert(Vector{T}, z)); v)

# -------------------------------------------------------
# NOTE: backends should implement the following methods to get/set the x/y/z data objects

tovec(v::AbstractVector) = v
tovec(v::Void) = zeros(0)

function getxy(plt::Plot, i::Integer)
    d = plt.seriesargs[i]
    tovec(d[:x]), tovec(d[:y])
end
function getxyz(plt::Plot, i::Integer)
    d = plt.seriesargs[i]
    tovec(d[:x]), tovec(d[:y]), tovec(d[:z])
end

function setxy!{X,Y}(plt::Plot, xy::Tuple{X,Y}, i::Integer)
    d = plt.seriesargs[i]
    d[:x], d[:y] = xy
end
function setxyz!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer)
    d = plt.seriesargs[i]
    d[:x], d[:y], d[:z] = xyz
end

# -------------------------------------------------------
# indexing notation

Base.getindex(plt::Plot, i::Integer) = getxy(plt, i)
Base.setindex!{X,Y}(plt::Plot, xy::Tuple{X,Y}, i::Integer) = setxy!(plt, xy, i)
Base.setindex!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) = setxyz!(plt, xyz, i)

# -------------------------------------------------------
# push/append for one series

# push value to first series
Base.push!(plt::Plot, y::Real) = push!(plt, 1, y)
Base.push!(plt::Plot, x::Real, y::Real) = push!(plt, 1, x, y)
Base.push!(plt::Plot, x::Real, y::Real, z::Real) = push!(plt, 1, x, y, z)

# y only
function Base.push!(plt::Plot, i::Integer, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesByOne(xdata), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, y::AVec)
    xdata, ydata = plt[i]
    if !isa(xdata, UnitRange{Int})
        error("Expected x is a UnitRange since you're trying to push a y value only")
    end
    plt[i] = (extendSeriesByOne(xdata, length(y)), extendSeriesData(ydata, y))
    plt
end

# x and y
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec)
    @assert length(x) == length(y)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end

# x, y, and z
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real, z::Real)
    # @show i, x, y, z
    xdata, ydata, zdata = getxyz(plt, i)
    # @show xdata, ydata, zdata
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec, z::AVec)
    @assert length(x) == length(y) == length(z)
    xdata, ydata, zdata = getxyz(plt, i)
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end

# tuples
Base.push!{X,Y}(plt::Plot, xy::Tuple{X,Y})                  = push!(plt, 1, xy...)
Base.push!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z})             = push!(plt, 1, xyz...)
Base.push!{X,Y}(plt::Plot, i::Integer, xy::Tuple{X,Y})      = push!(plt, i, xy...)
Base.push!{X,Y,Z}(plt::Plot, i::Integer, xyz::Tuple{X,Y,Z}) = push!(plt, i, xyz...)

# -------------------------------------------------------
# push/append for all series

# push y[i] to the ith series
function Base.push!(plt::Plot, y::AVec)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, y[mod1(i,ny)])
    end
    plt
end

# push y[i] to the ith series
# same x for each series
function Base.push!(plt::Plot, x::Real, y::AVec)
    push!(plt, [x], y)
end

# push (x[i], y[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec)
    nx = length(x)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)])
    end
    plt
end

# push (x[i], y[i], z[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec, z::AVec)
    nx = length(x)
    ny = length(y)
    nz = length(z)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)], z[mod1(i,nz)])
    end
    plt
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------
# graphs detailing the features that each backend supports

function supportGraph(allvals, func)
    vals = reverse(sort(allvals))
    bs = sort(backends())
    x, y = map(string, bs), map(string, vals)
    nx, ny = map(length, (x,y))
    z = zeros(nx, ny)
    for i=1:nx, j=1:ny
        supported = func(Plots._backend_instance(bs[i]))
        z[i,j] = float(vals[j] in supported) * (0.4i/nx+0.6)
    end
    heatmap(x, y, z,
            color = ColorGradient([:white, :darkblue]),
            line = (1, :black),
            leg = false,
            size = (50nx+50, 35ny+100),
            xlim = (0.5, nx+0.5),
            ylim = (0.5, ny+0.5),
            xrotation = 60,
            aspect_ratio = :equal)
end

supportGraphArgs()    = supportGraph(_allArgs, supportedArgs)
supportGraphTypes()   = supportGraph(_allTypes, supportedTypes)
supportGraphStyles()  = supportGraph(_allStyles, supportedStyles)
supportGraphMarkers() = supportGraph(_allMarkers, supportedMarkers)
supportGraphScales()  = supportGraph(_allScales, supportedScales)
supportGraphAxes()    = supportGraph(_allAxes, supportedAxes)

function dumpSupportGraphs()
    for func in (supportGraphArgs, supportGraphTypes, supportGraphStyles,
               supportGraphMarkers, supportGraphScales, supportGraphAxes)
        plt = func()
        png(Pkg.dir("ExamplePlots", "docs", "examples", "img", "supported", "$(string(func))"))
    end
end

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

# Some conversion functions
# note: I borrowed these conversion constants from Compose.jl's Measure

const PX_PER_INCH   = 100
const DPI           = PX_PER_INCH
const MM_PER_INCH   = 25.4
const MM_PER_PX     = MM_PER_INCH / PX_PER_INCH

inch2px(inches::Real)   = float(inches * PX_PER_INCH)
px2inch(px::Real)       = float(px / PX_PER_INCH)
inch2mm(inches::Real)   = float(inches * MM_PER_INCH)
mm2inch(mm::Real)       = float(mm / MM_PER_INCH)
px2mm(px::Real)         = float(px * MM_PER_PX)
mm2px(mm::Real)         = float(px / MM_PER_PX)

"Smallest x in plot"
xmin(plt::Plot) = minimum([minimum(d[:x]) for d in plt.seriesargs])
"Largest x in plot"
xmax(plt::Plot) = maximum([maximum(d[:x]) for d in plt.seriesargs])

"Extrema of x-values in plot"
Base.extrema(plt::Plot) = (xmin(plt), xmax(plt))


text:
import NLPModelsKnitro: knitro
#Check https://github.com/JuliaSmoothOptimizers/NLPModelsKnitro.jl/blob/master/src/NLPModelsKnitro.jl
#=
*General options*
algorithm: Indicates which algorithm to use to solve the problem
blasoption: Specifies the BLAS/LAPACK function library to use for basic vector
            and matrix computations
cg_maxit: Determines the maximum allowable number of inner conjugate gradient
          (CG) iterations
cg_pmem: Specifies number of nonzero elements per hessian column when computing
         preconditioner
cg_precond: Specifies whether or not to apply preconditioning during
            CG iterations in barrier algorithms
cg_stoptol: Relative stopping tolerance for CG subproblems
convex: Identify convex models and apply specializations often beneficial for convex models
delta: Specifies the initial trust region radius scaling factor
eval_fcga: Specifies that gradients are provided together with functions in one callback
honorbnds: Indicates whether or not to enforce satisfaction of simple variable bounds
initpenalty: Initial penalty value used in Knitro merit function
linesearch_maxtrials: Indicates the maximum allowable number of trial points during the linesearch
linesearch: Indicates which linesearch strategy to use for the Interior/Direct or SQP algorithm
linsolver_ooc: Indicates whether to use Intel MKL PARDISO out-of-core solve of linear systems
linsolver: Indicates which linear solver to use to solve linear systems arising in Knitro algorithms
linsolver_maxitref
linsolver_pivottol
objrange: Specifies the extreme limits of the objective function for purposes of determining unboundedness
          Default value: 1.0e20
presolve: Determine whether or not to use the Knitro presolver
presolve_initpt: Controls whether Knitro presolver can shift user-supplied initial point
restarts: Specifies whether to enable automatic restarts
restarts_maxit: Maximum number of iterations before restarting when restarts are enabled
scale: Specifies whether to perform problem scaling
soc: Specifies whether or not to try second order corrections (SOC)
strat_warm_start: Specifies whether or not to invoke a warm-start strategy

*Derivatives options*
bfgs_scaling: Specifies the initial scaling for the BFGS or L-BFGS Hessian approximation
derivcheck: Determine whether or not to perform a derivative check on the model
gradopt: Specifies how to compute the gradients of the objective and constraint functions
hessian_no_f: Determines whether or not to allow Knitro to request Hessian
              evaluations without the objective component included.
hessopt: Specifies how to compute the (approximate) Hessian of the Lagrangian
         1 (exact) User provides a routine for computing the exact Hessian. (Default)
         4 (product_findiff) Knitro computes Hessian-vector products using finite-differences.
         5 (product) User provides a routine to compute the Hessian-vector products.
         6 (lbfgs) Knitro computes a limited-memory quasi-Newton BFGS Hessian
                   (its size is determined by the option lmsize).
lmsize: Specifies the number of limited memory pairs stored when approximating the Hessian

*Termination options*
feastol: Specifies the final relative stopping tolerance for the feasibility error.
         Default value: 1.0e-6
feastol_abs: Specifies the final absolute stopping tolerance for the feasibility error.
             Default value: 1.0e-3
fstopval: Used to implement a custom stopping condition based on the objective function value
ftol: The optimization process will terminate if feasible and the relative change
      in the objective function is less than ftol
      Default value: 1.0e-15
ftol_iters: The optimization process will terminate if the relative change in
            the objective function is less than ftol for ftol_iters consecutive
            feasible iterations. Default value: 5

maxfevals: Specifies the maximum number of function evaluations before termination.
           Default value: -1 (unlimited)
maxit: Specifies the maximum number of iterations before termination
       0 is default value, let Knitro set it (10000 for LP/NLP)
maxtime_cpu: Specifies, in seconds, the maximum allowable CPU time before termination.
             Default value: 1.0e8
maxtime_real: Specifies, in seconds, the maximum allowable real time before termination.
              Default value: 1.0e8
opttol: Specifies the final relative stopping tolerance for the KKT (optimality) error
        Default value: 1.0e-6
opttol_abs: Specifies the final absolute stopping tolerance for the KKT (optimality) error
            Default value: 1.0e-3
xtol: The optimization process will terminate if the relative change of the
      solution point estimate is less than xtol. Default value: 1.0e-12
xtol_iters: Number of consecutive iterations where change of the solution point
            estimate is less than xtol before Knitro stops.
            Default is 1

*Output options*
out_hints: Print diagnostic hints (e.g. on user option settings) after solving
           0 no prints, and 1 prints (default in Knitro)
outlev: Controls the level of output produced by Knitro
        0 (none) Printing of all output is suppressed.
        1 (summary) Print only summary information.
        2 (iter_10) Print basic information every 10 iterations. (default in Knitro)
=#
"""
knitro(nlp) DOESN'T CHECK THE WRONG KWARGS, AND RETURN AN ERROR.

knitro(::NLPStopping)

Selection of possible [options](https://www.artelys.com/docs/knitro/3_referenceManual/userOptions.html):
"""
function knitro(stp          :: NLPStopping;
                convex       :: Int  = -1, #let Knitro deal with it :)
                objrange     :: Real = stp.meta.unbounded_threshold,
                hessopt      :: Int  = 1,
                feastol      :: Real = stp.meta.rtol,
                feastol_abs  :: Real = stp.meta.atol,
                opttol       :: Real = stp.meta.rtol,
                opttol_abs   :: Real = stp.meta.atol,
                maxfevals    :: Int  = stp.meta.max_cntrs[:neval_sum],
                maxit        :: Int  = 0, #stp.meta.max_iter
                maxtime_real :: Real = stp.meta.max_time,
                out_hints    :: Int  = 0,
                outlev       :: Int  = 2,
                kwargs...)

  @assert -1 ≤ convex ≤ 1
  @assert 1  ≤ hessopt ≤ 7
  @assert 0  ≤ out_hints ≤ 1
  @assert 0  ≤ outlev ≤ 6
  @assert 0  ≤ maxit

  nlp = stp.pb
  #y0 = stp.current_state.lambda #si défini
  #z0 = stp.current_state.mu #si défini
  stats = knitro(nlp, x0           = stp.current_state.x,
                      objrange     = objrange,
                      feastol      = feastol,
                      feastol_abs  = feastol_abs,
                      opttol       = opttol,
                      opttol_abs   = opttol_abs,
                      maxfevals    = maxfevals,
                      maxit        = maxit,
                      maxtime_real = maxtime_real,
                      out_hints    = out_hints,
                      outlev       = outlev;
                      kwargs...)

  if stats.status ∉ (:exception, :unknwon) #∈ (:first_order, :acceptable)
    stp.current_state.x  = stats.solution
    stp.current_state.fx = stats.objective
    stp.current_state.gx = grad(nlp, stats.solution)#stats.dual_feas
    stp.current_state.current_score  = norm(stp.current_state.gx, Inf)#stats.dual_feas #TODO: this is for unconstrained problem!!
  end
  #Update the meta boolean with the output message
  stp = stats_status_to_meta!(stp, stats)

  if status(stp) == :Unknown
    @warn "Error in StoppingInterface statuses: return status is $(stats.status)"
    print(stats)
  end

  return stp #would be better to return the stats somewhere
end


text:
"""
Julia ain't quite perfect, so here are some tools to help make the ride a bit
smoother.
"""

using Logging

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

"""
As far as I can tell, the Julia logging library doesn't let us set the
logging level with a string, so let's fix that up.
"""
logLevels = Dict{String, Int32}(
    "debug" => Logging.Debug.level,
    "info"  => Logging.Info.level,
    "warn"  => Logging.Warn.level,
    "error" => Logging.Error.level,
    )

function initialise_logger(level::AbstractString)
    initialise_logger(logLevels[lowercase(level)])
end

function initialise_logger(level::Int32)
    global_logger(SimpleLogger(stdout, level))
    @debug "Set log level to $level"
end


text:
function transit_poly_integrate(r, b, u, N)
    n = length(u)
    if r ≤ b - 1
        return 1.0
    end
    s_1 = max(0, b - r)
    s_2 = min(1, b + r)
    ds = (s_2 - s_1) / N
    s = range(s_1 + 0.5 * ds, s_2 - 0.5 * ds, length=N)

    # normalization
    norm = mapreduce(i -> 2 * u[i] / (2 + 3 * i + i^2), -, 1:n, init=1.0)

    flux = sum(s) do s_j
        if s_j < r - b
            dϕ = 2π
        else
            f = (s_j^2 + b^2 - r^2) / (2 * b * s_j)
            if abs(f) ≈ 1 && abs(f) > 1
                f = clamp(f, -1, 1)
            end
            dϕ = 2 * acos(f)
        end
        z = sqrt(1 - s_j^2)
        iμ = mapreduce(i -> u[i] * (1 - z)^i, -, 1:n, init=1.0)

        return s_j * ds * dϕ * iμ
    end

    flux /= π * norm
    flux = 1 - flux
    return flux
end

@testset "polynomial limb darkening" begin
    r0 = [
        fill(0.01, 5)
        fill(0.1, 5)
        fill(1, 3)
        fill(10, 3)
        fill(100, 3)
    ]
    b0 = [
        0, 0.1, 0.99, 1, 1.01,
        0, 0.1, 0.9, 1, 1.1,
        0, 1, 2,
        9, 10, 11,
        99, 100, 101
    ]
    n = 2 + ceil(Int, rand(rng) * 20)
    u = rand(rng, n)
    u *= rand(rng) / sum(u)

    N = 5000

    ld = PolynomialLimbDark(u)
    flux = @. ld(b0, r0)
    f_num = @. transit_poly_integrate(r0, b0, (u,), N)

    @info "Maximum difference of lightcurve for N=$N: $(maximum(abs,flux-f_num)))"
    @test flux ≈ f_num atol=1e-5

    # test vs. quad limb dark
    ld_ = PolynomialLimbDark(u[1:2])
    ld_q = QuadLimbDark(u[1:2])
    @test ld_.(b0, r0) == ld_q.(b0, r0)

    npts = 1000
    iu = 30
    u_n = ones(iu) / iu
    b0 = @. 3 * ((1:npts) / npts) - 0.5
    ld = PolynomialLimbDark(u_n)
    lc_ana = @. ld(abs(b0), 0.1)
    lc_num = @. transit_poly_integrate(0.1, abs(b0), (u_n,), N)

    @info "Maximum difference of lightcurve for N=$iu: $(maximum(abs,lc_ana-lc_num)))"
    @test lc_ana ≈ lc_num atol=1e-5

end

@testset "QuadLimbDark interface" begin
    # test truncation
    ld = (@test_logs (:warn, "Higher-order terms will be ignored") QuadLimbDark(ones(10)))
    @test ld.n_max == 2
    @test ld.g_n == QuadLimbDark(ones(2)).g_n
end

@testset "IntegratedLimbDark interface" for basis in [:legendre, :radau, :lobatto]
    u = [0.4, 0.26]
    ld = PolynomialLimbDark(u)
    ldt = IntegratedLimbDark(ld, basis=basis)

    @test IntegratedLimbDark(u).driver isa PolynomialLimbDark

    orbit = SimpleOrbit(period=3, duration=1)
    # sharp edges
    ts = [-0.5, -0.25, 0, 0.25, 0.5]
    above_or_below = [-1, 1, 1, 1, -1]

    # test no exposure time is exactly same
    @test ld.(orbit, ts, 0.01) == ldt.(orbit, ts, 0.01)

    # test exposure "smears" signal causing it to be
    # lower or higher then the exact signal
    f = @. ldt(orbit, ts, 0.01, [0.3 0.5])
    base = @. ld(orbit, ts, 0.01)
    @test sign.(f[:, 1] .- base) == above_or_below
    @test sign.(f[:, 2] .- base) == above_or_below
    # test higher exposure time is more smearing
    @test sign.(f[:, 2] .- f[:, 1]) == above_or_below

end

@testset "secondary" begin
    u = [0.4, 0.26]
    d1 = PolynomialLimbDark(u)
    d2 = PolynomialLimbDark(ones(2))
    ld = SecondaryLimbDark(d1, d2; brightness_ratio=0.1)
    ld2 = SecondaryLimbDark(u, ones(2))
    @test ld.primary_driver.g_n == ld2.primary_driver.g_n == d1.g_n
    @test ld.secondary_driver.g_n == ld2.secondary_driver.g_n == d2.g_n

    orbit = SimpleOrbit(period=4, duration=1)
    r = 0.01
    f1 = ld(orbit, 0, r)
    f2 = ld(orbit, 2, r)
    @test f1 ≈ f2 rtol=1e-3
end

@testset "reciprocity" begin
    u = [0.4, 0.26]
    d1 = PolynomialLimbDark(u)
    d2 = PolynomialLimbDark(ones(2))
    ld1 = IntegratedLimbDark(SecondaryLimbDark(d1, d2))
    ld2 = SecondaryLimbDark(IntegratedLimbDark(d1), IntegratedLimbDark(d2))

    orbit = SimpleOrbit(period=4, duration=1)
    r = 0.01
    t = range(-1.5, 1.5, length=1000)
    for texp in [nothing, 0.1, 0.3]
        @test ld1.(orbit, t, r; texp) ≈ ld2.(orbit, t, r; texp)
    end
end


text:
############################################################################
### models w.o. meanstructure
############################################################################

model_ml = Sem(
    specification = spec,
    data = dat,
    diff = semdiff
)

model_ls_sym = Sem(
    specification = spec,
    data = dat,
    imply = RAMSymbolic,
    loss = SemWLS,
    start_val = start_simple,
    diff = semdiff
)

model_ml_sym = Sem(
    specification = spec,
    data = dat,
    imply = RAMSymbolic,
    diff = semdiff
)

model_ridge = Sem(
    specification = spec,
    data = dat,
    loss = (SemML, SemRidge),
    α_ridge = .001,
    which_ridge = 16:20,
    diff = semdiff
)

model_constant = Sem(
    specification = spec,
    data = dat,
    loss = (SemML, SemConstant),
    constant_loss = 3.465,
    diff = semdiff
)

model_ml_weighted = Sem(
    specification = partable,
    data = dat,
    loss_weights = (n_obs(model_ml),),
    diff = semdiff
)

############################################################################
### test gradients
############################################################################

models = [model_ml, model_ls_sym, model_ridge, model_constant, model_ml_sym, model_ml_weighted]
names = ["ml", "ls_sym", "ridge", "constant", "ml_sym", "ml_weighted"]

for (model, name) in zip(models, names)
    try
        @testset "$(name)_gradient" begin
            @test test_gradient(model, start_test; rtol = 1e-9)
        end
    catch
    end
end

############################################################################
### test solution
############################################################################

models = [model_ml, model_ls_sym, model_ml_sym, model_constant]
names = ["ml", "ls_sym", "ml_sym", "constant"]
solution_names = Symbol.("parameter_estimates_".*["ml", "ls", "ml", "ml"])

for (model, name, solution_name) in zip(models, names, solution_names)
    try
        @testset "$(name)_solution" begin
            solution = sem_fit(model)
            update_estimate!(partable, solution)
            @test compare_estimates(partable, solution_lav[solution_name]; atol = 1e-2)
        end
    catch
    end
end

@testset "ridge_solution" begin
    solution_ridge = sem_fit(model_ridge)
    solution_ml = sem_fit(model_ml)
    # solution_ridge_id = sem_fit(model_ridge_id)
    @test abs(solution_ridge.minimum - solution_ml.minimum) < 1
end

# test constant objective value
@testset "constant_objective_and_gradient" begin
    @test (objective!(model_constant, start_test) - 3.465) ≈ objective!(model_ml, start_test)
    @test gradient!(model_constant, start_test) ≈ gradient!(model_ml, start_test)
end

@testset "ml_solution_weighted" begin
    solution_ml = sem_fit(model_ml)
    solution_ml_weighted = sem_fit(model_ml_weighted)
    @test isapprox(solution(solution_ml), solution(solution_ml_weighted), rtol = 1e-3)
    @test isapprox(n_obs(model_ml)*StructuralEquationModels.minimum(solution_ml), StructuralEquationModels.minimum(solution_ml_weighted), rtol = 1e-6)
end

############################################################################
### test fit assessment
############################################################################

@testset "fitmeasures/se_ml" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml]; atol = 1e-3))

    update_partable!(partable, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3, col = :se, lav_col = :se)
end

@testset "fitmeasures/se_ls" begin
    solution_ls = sem_fit(model_ls_sym)
    fm = fit_measures(solution_ls)
    @test all(test_fitmeasures(fm, solution_lav[:fitmeasures_ls]; atol = 1e-3, fitmeasure_names = fitmeasure_names_ls))
    @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing)

    update_partable!(partable, identifier(model_ls_sym), se_hessian(solution_ls), :se)
    @test compare_estimates(partable, solution_lav[:parameter_estimates_ls]; atol = 1e-2, col = :se, lav_col = :se)
end

############################################################################
### test hessians
############################################################################

if semdiff == SemDiffOptim
    using Optim, LineSearches

    model_ls = Sem(
        specification = spec,
        data = dat,
        imply = RAMSymbolic,
        loss = SemWLS,
        hessian = true,
        algorithm = Newton(
            ;linesearch = BackTracking(order=3),
            alphaguess = InitialHagerZhang())
    )

    model_ml = Sem(
        specification = spec,
        data = dat,
        imply = RAMSymbolic,
        hessian = true,
        algorithm = Newton()
    )

    @testset "ml_hessians" begin
        @test test_hessian(model_ml, start_test; atol = 1e-4)
    end

    @testset "ls_hessians" begin
        @test test_hessian(model_ls, start_test; atol = 1e-4)
    end

    @testset "ml_solution_hessian" begin
        solution = sem_fit(model_ml)
        update_estimate!(partable, solution)
        @test compare_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-3)
    end

    @testset "ls_solution_hessian" begin
        solution = sem_fit(model_ls)
        update_estimate!(partable, solution)
        @test compare_estimates(partable, solution_lav[:parameter_estimates_ls]; atol = 1e-3)
    end

end

############################################################################
### meanstructure
############################################################################

# models
model_ls = Sem(
    specification = spec_mean,
    data = dat,
    imply = RAMSymbolic,
    loss = SemWLS,
    meanstructure = true,
    diff = semdiff
)

model_ml = Sem(
    specification = spec_mean,
    data = dat,
    meanstructure = true,
    diff = semdiff
)

model_ml_sym = Sem(
    specification = spec_mean,
    data = dat,
    imply = RAMSymbolic,
    meanstructure = true,
    start_val = start_test_mean,
    diff = semdiff
)

############################################################################
### test gradients
############################################################################

models = [model_ml, model_ls, model_ml_sym]
names = ["ml", "ls_sym", "ml_sym"]

for (model, name) in zip(models, names)
    try
        @testset "$(name)_gradient_mean" begin
            @test test_gradient(model, start_test_mean; rtol = 1e-9)
        end
    catch
    end
end

############################################################################
### test solution
############################################################################

solution_names = Symbol.("parameter_estimates_".*["ml", "ls", "ml"].*"_mean")

for (model, name, solution_name) in zip(models, names, solution_names)
    try
        @testset "$(name)_solution_mean" begin
            solution = sem_fit(model)
            update_estimate!(partable_mean, solution)
            @test compare_estimates(partable_mean, solution_lav[solution_name]; atol = 1e-2)
        end
    catch
    end
end

############################################################################
### test fit assessment
############################################################################

@testset "fitmeasures/se_ml_mean" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_ml_mean]; atol = 1e-3))

    update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_ml_mean]; atol = 1e-3, col = :se, lav_col = :se)
end

@testset "fitmeasures/se_ls_mean" begin
    solution_ls = sem_fit(model_ls)
    fm = fit_measures(solution_ls)
    @test all(test_fitmeasures(fm, solution_lav[:fitmeasures_ls_mean]; atol = 1e-3, fitmeasure_names = fitmeasure_names_ls))
    @test (fm[:AIC] === missing) & (fm[:BIC] === missing) & (fm[:minus2ll] === missing)

    update_partable!(partable_mean, identifier(model_ls), se_hessian(solution_ls), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_ls_mean]; atol = 1e-2, col = :se, lav_col = :se)
end

############################################################################
### fiml
############################################################################

# models
model_ml = Sem(
    specification = spec_mean,
    data = dat_missing,
    observed = SemObsMissing,
    loss = SemFIML,
    diff = semdiff
)

model_ml_sym = Sem(
    specification = spec_mean,
    data = dat_missing,
    observed = SemObsMissing,
    imply = RAMSymbolic,
    loss = SemFIML,
    start_val = start_test_mean,
    diff = semdiff
)

############################################################################
### test gradients
############################################################################

@testset "fiml_gradient" begin
    @test test_gradient(model_ml, start_test_mean; atol = 1e-6)
end

@testset "fiml_gradient_symbolic" begin
    @test test_gradient(model_ml_sym, start_test_mean; atol = 1e-6)
end

############################################################################
### test solution
############################################################################

@testset "fiml_solution" begin
    solution = sem_fit(model_ml)
    update_estimate!(partable_mean, solution)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
end

@testset "fiml_solution_symbolic" begin
    solution = sem_fit(model_ml_sym)
    update_estimate!(partable_mean, solution)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-2)
end

############################################################################
### test fit measures
############################################################################

@testset "fitmeasures/se_fiml" begin
    solution_ml = sem_fit(model_ml)
    @test all(test_fitmeasures(fit_measures(solution_ml), solution_lav[:fitmeasures_fiml]; atol = 1e-3))

    update_partable!(partable_mean, identifier(model_ml), se_hessian(solution_ml), :se)
    @test compare_estimates(partable_mean, solution_lav[:parameter_estimates_fiml]; atol = 1e-3, col = :se, lav_col = :se)
end


text:
using AMQPClient
const VIRTUALHOST ="/"
const HOST = "127.0.0.1"

function send(message)
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare a queue
            success, queue_name, message_count, consumer_count = queue_declare(chan, "hello")
            # 4.1 Create a message, AMQPCleint only accepts message in UInt8
            data = convert(Vector{UInt8}, codeunits(message))
            msg = Message(data, content_type="text/plain", delivery_mode=PERSISTENT)
            # 4.2 Send a message
            basic_publish(chan, msg; exchange="", routing_key="hello")
            println("Message sent: $message")
            # 5. Auto-closes the channel and connection
        end
    end
end

send("Hello World!")

max_forks_repo_name: tomashq/rabbitmq-tutorials


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

include("../src/Git-Work-Tracker.jl")

using Documenter, .GitWorkTracker

makedocs(sitename="Git Work Tracker Documentation", modules = [GitWorkTracker])

max_stars_repo_name: LuRenJiasWorld/Git-Work-Tracker

max_issues_repo_name: LuRenJiasWorld/Git-Work-Tracker

max_forks_repo_name: LuRenJiasWorld/Git-Work-Tracker


text:
module LibFEM

using StaticArrays

include("mesh.jl")

end # module


text:
using DataStructures
using Plots; pyplot(fmt="png")
PyPlot.matplotlib.rc("mathtext", fontset="cm");
PyPlot.matplotlib.rc("text", usetex=true);
PyPlot.matplotlib.rc("font", family="serif", size=12);
using LaTeXStrings
using LsqFit
using Statistics
using JLD

data_dir = "/home/user/thesis/data/scaling";
savepath = "/home/user/thesis/latex/images";
save_bool = true
dpi = 144
if save_bool
	dpi = 300
end
save_bool

# %%

Ns = 2 .^ collect(15:24);
save_file = joinpath(data_dir, "cluster_size_distribution_avgs.jld")
distribution_avgs = load(save_file)["distribution_avgs"]

# %% Plot t_0 data

distribution_avg = distribution_avgs["t_0"]
# distribution_avg = Dict(2^k => distributions["t_0"][2^k][1] for k in 15:24)
plot_0 = plot(dpi=dpi);
Ns = sort(collect(keys(distribution_avg)));
colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))];
reverse!(colors);
markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys);

for i in 1:length(Ns)
	N = Ns[i]
	distribution = distribution_avg[N]
	s = collect(keys(distribution))
	n_s = collect(values(distribution))
	x = log10.(s)
	y = log10.(n_s)
	scatter!(x, y,
		legend=:topright,
		xaxis=(L"\log_{10} (s)", (-0.2, 5)),
		yaxis=(L"\log_{10} (n_s(t_0))", (-12, 0)),
		marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])),
		label=latexstring("N = 2^{$(Int(log2(N)))}"),
	)
end

@. function f_linear(x, α)
	return α[1] * x + α[2]
end

αs    = [];
σs    = [];
ratio = 0.01;
for N in Ns[end:end]
	d = distribution_avg[N];
	x = log10.(collect(keys(d)));
	y = log10.(collect(values(d)));
	order = sortperm(x);
	x = x[order][1:Int(floor(ratio * length(x)))];
	y = y[order][1:Int(floor(ratio * length(y)))];
	fit = curve_fit(f_linear, x, y, [-0.1, 0])
	push!(αs, fit.param)
	push!(σs, stderror(fit))
end

α = sum(αs) / length(αs)
σ = sum(σs) / length(σs)
τ₀ = -α[1]

d = distribution_avg[2^24];
x = log10.(collect(keys(d)));
order = sortperm(x);
x = x[order][1:Int(floor(ratio * length(x)))];
y = f_linear(x, α);
plot!(x, y,
	label=latexstring("\\mathrm{Fit} (N = 2^{24}) \\sim $(round(α[1], digits=3)) ($(Int(round(σ[1] * 1e3, digits=0)))) \\cdot \\log_{10}(s)"),
	line=(1, :black),
)

if save_bool
	savefig(plot_0, joinpath(savepath, "n_s_t_0.png"));
end

plot_0

# %% Plot t_1 data

distribution_avg = distribution_avgs["t_1"];
plot_1 = plot(dpi=dpi);
Ns = sort(collect(keys(distribution_avg)));
colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))];
reverse!(colors);
markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys);

for i in 1:length(Ns)
	N = Ns[i]
	distribution = distribution_avg[N]
	s = collect(keys(distribution))
	n_s = collect(values(distribution))
	x = log10.(s ./ N)
	y = log10.(n_s)
	scatter!(x, y,
		legend=:topright,
		xaxis=(L"\log_{10} (s / N)", (-8, 0)),
		yaxis=(L"\log_{10} (n_s(t_1))", (-12, 0)),
		marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])),
		label=latexstring("N = 2^{$(Int(log2(N)))}"),
	)
end

@. function f_linear(x, α)
	return α[1] * x + α[2]
end

αs    = [];
σs    = [];
ratio = 0.10;
for N in Ns[end:end]
	d = distribution_avg[N];
	x = log10.(collect(keys(d)) ./ N);
	y = log10.(collect(values(d)));
	order = sortperm(x);
	x = x[order][1:Int(floor(ratio * length(x)))];
	y = y[order][1:Int(floor(ratio * length(y)))];
	fit = curve_fit(f_linear, x, y, [-0.1, 0])
	push!(αs, fit.param)
	push!(σs, stderror(fit))
end

α = sum(αs) / length(αs)
σ = sum(σs) / length(σs)
τ₁ = -α[1]

d = distribution_avg[2^24];
x = log10.(collect(keys(d)) ./ 2^24);
order = sortperm(x);
x = x[order][1:Int(floor(ratio * length(x)))];
y = f_linear(x, α)
plot!(x, y,
	label=latexstring("\\mathrm{Fit} (N = 2^{24}) \\sim $(round(α[1], digits=3)) ($(Int(round(σ[1] * 1e3, digits=0)))) \\cdot \\log_{10}(s/N)"),
	line=(1, :black),
)

if save_bool
	savefig(plot_1, joinpath(savepath, "n_s_t_1.png"));
end

plot_1

# %% Plot FSS collapse for τ₀

distribution_avg = distribution_avgs["t_0"]
plot_a = plot(dpi=dpi);
Ns = sort(collect(keys(distribution_avg)));
colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))]
reverse!(colors)
markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

for i in 1:length(Ns)
	N = Ns[i]
	distribution = distribution_avg[N]
	s = collect(keys(distribution))
	n_s = collect(values(distribution))
	x = log10.(s ./ N^(1/τ₁))
	y = log10.(s.^τ₁ .* n_s)
	scatter!(x, y,
		legend=:topleft,
		xaxis=(L"\log_{10} (s / N^{1/\tau})", (-4, 2), -4:1:2),
		yaxis=(L"\log_{10} (s^\tau n_s(t_0))", (-1.5, 0.5), -1.5:0.5:0.5),
		marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])),
		label=latexstring("N = 2^{$(Int(log2(N)))}"),
	)
	annotate!(-2, 0, text(latexstring("\\tau_1 \\approx $(round(τ₁, digits=3))")))
end

if save_bool
	savefig(plot_a, joinpath(savepath, "fss_collapse.png"));
end

plot_a

# %% Plot FSS collapse for τ₀

distribution_avg = distribution_avgs["t_0"]
plot_b = plot(dpi=dpi);
Ns = sort(collect(keys(distribution_avg)));
colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))]
reverse!(colors)
markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

for i in 1:length(Ns)
	N = Ns[i]
	distribution = distribution_avg[N]
	s = collect(keys(distribution))
	n_s = collect(values(distribution))
	x = log10.(s ./ N^(1/τ₀))
	y = log10.(s.^τ₀ .* n_s)
	scatter!(x, y,
		legend=false,
		xaxis=(L"\log_{10} (s / N^{1/\tau})", (-4, 2), -4:1:2),
		yaxis=(L"\log_{10} (s^\tau n_s(t_0))", (-1.5, 0.5), -1.5:0.5:0.5),
		marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])),
		label=latexstring("N = 2^{$(Int(log2(N)))}"),
	)
	annotate!(-2, 0, text(latexstring("\\tau_0 \\approx $(round(τ₀, digits=3))")))
end

plot_b

# %% Plot FSS collapse for τ₀

τ_ref = τ₁ + 0.05
distribution_avg = distribution_avgs["t_0"]
plot_c = plot(dpi=dpi);
Ns = sort(collect(keys(distribution_avg)));
colors = RGB[cgrad(:plasma)[z] for z=range(0, stop=1, length=length(Ns))]
reverse!(colors)
markers = filter((m->begin m in Plots.supported_markers() end), Plots._shape_keys)

for i in 1:length(Ns)
	N = Ns[i]
	distribution = distribution_avg[N]
	s = collect(keys(distribution))
	n_s = collect(values(distribution))
	x = log10.(s ./ N^(1/τ_ref))
	y = log10.(s.^τ_ref .* n_s)
	scatter!(x, y,
		legend=false,
		xaxis=(L"\log_{10} (s / N^{1/\tau})", (-4, 2), -4:1:2),
		yaxis=(L"\log_{10} (s^\tau n_s(t_0))", (-1.5, 0.5), -1.5:0.5:0.5),
		marker=(2, colors[i], markers[i], 0.8, stroke(colors[i])),
		label=latexstring("N = 2^{$(Int(log2(N)))}"),
	)
	annotate!(-2, 0, text(latexstring("\\tau_1 + 0.05")))
end

plot_c

# %% Plot FSS comparison

l = @layout [a{0.6h}; [b{0.5w} c]]
plot_triple = plot(plot_a, plot_b, plot_c, layout=l, size=(800, 600))

if save_bool
	savefig(plot_triple, joinpath(savepath, "fss_collapse_triple.png"));
end

plot_triple

# %% Maximual decay exponent calculation

δ = 0.302;
exponent_of_decay = -δ * (τ₁ - 2) / (τ₁ - 1)

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation