text:
using JLD2, Distributions, Plots, StatPlots, AverageShiftedHistograms, DynNets
 plotly()

N_est = 500
Tvals =  [50 500 5000]
N_ind = length(Tvals)[1]

N = 10 ;
GBAest =1; GWest = 10; useStartVal = false

useStartVal? startValFlag = "T" : startValFlag = "F";# false #false #
# save
save_fold = "./data/estimatesTest/"
file_name = "EstTest_$(N)_$(Tvals)_$(N_est)_$(useStartVal)_Sin_Miss_$(GWest)_$(GBAest).jld"
save_path = save_fold*file_name#
estDataT = @load(save_path,groupsIndsEst,ObsAllT,ParTvDgpAllT ,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)

##

RavgMSE(A,B) = sqrt.(mean((A .- B).^2,1))
AvgBias(A,B) = mean((A .- B),1)
tmpModelEst = DynNets.SdErgm1_UmBAgr(zeros(10,10),[ones(GWest),ones(GBAest),ones(GBAest)],groupsIndsEst)
rmseTotGas = zeros(N_est,N,N_ind)
biasTotGas = zeros(N_est,N,N_ind)
rmseTotSnap = zeros(N_est,N,N_ind)
biasTotSnap = zeros(N_est,N,N_ind)
filTvParT = Array{Array{Float64,3},1}(N_ind)
ParSnapT =  Array{Array{Float64,3},1}(N_ind)
# Compute MSE Bias and Snap estimates
@time for indT = 1:2#N_ind
    #indT=1
            T =Tvals[indT]
            tmpParDgp = ParTvDgpAllT[indT]
            tmpParFil = zeros(T,N,N_est)
            tmpParSnap = zeros(T,N,N_est)
        for n=1:N_est
    #n=100
    obsT = ObsAllT[indT][:,:,n]
    vecParEst =  [estPar_1[indT,n,:]; estPar_2[indT,n];  estPar_3[indT,n]]
    GroupIndsEst = Array{Array{<:Int,1},1}(2)
    GroupIndsEst[1] = groupsIndsEst[1]; GroupIndsEst[2] = groupsIndsEst[2]
    tmpParFil[:,:,n] =   DynNets.score_driven_filter_or_dgpAndLikeliood(tmpModelEst,vecParEst,
                                                degsT = obsT,
                                                groupsInds = GroupIndsEst )[1]

    fitSnap = DynNets.estSnapSeq(DynNets.SdErgm1(obsT)) # sequence of snapshots estimate
    # remove the infinites
    for i=1:5
        infInd = find(.!isfinite.(fitSnap))
        replInd = infInd[infInd .> 2]
        fitSnap[replInd] = fitSnap[replInd .- 1]
        infInd = find(.!isfinite.(fitSnap))
        replInd = infInd[infInd .< (Int(N*T) - 2) ]
        fitSnap[replInd] = fitSnap[replInd .+ 1]
    end
    tmpParSnap[:,:,n] = fitSnap

            rmseTotGas[n,:,indT] = RavgMSE(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            biasTotGas[n,:,indT] = AvgBias(tmpParFil[:,:,n],tmpParDgp[:,:,n] )
            rmseTotSnap[n,:,indT] = RavgMSE(fitSnap,tmpParDgp[:,:,n] )
            biasTotSnap[n,:,indT] = AvgBias(fitSnap,tmpParDgp[:,:,n] )
        end

    filTvParT[indT]  = tmpParFil
    ParSnapT[indT] = tmpParSnap
end
dgpTvParT = ParTvDgpAllT
## Example single realization plots
Nlines=10
indPl = round(Int,linspace(1,N,Nlines))
C(g::ColorGradient) = RGB[g[z] for z=linspace(0,1,Nlines)]
lineColors = C(cgrad(:rainbow))

indT = 2
n=4
T = Tvals[indT]

fit_T =   dgpTvParT[indT][:,:,n]
fil_fit_T =  filTvParT[indT][:,:,n]
fitSnap =  ParSnapT[indT][:,:,n]
 p=plot()
for i = 1:Nlines p = plot!(fit_T[:,indPl[i]],  color = lineColors[i] ) end #lineColors
for i = 1:Nlines p = plot!(fil_fit_T[:,indPl[i]],linestyle=:dash,  color = lineColors[i]  ) end #lineColors
for i = 1:Nlines p= plot!(fitSnap[:,indPl[i]],linestyle=:dot,markershape=:circle,
    markersize=1,linewidth = 0, markerstrokecolor=  lineColors[i], color = lineColors[i]  ) end

p =plot(p,titlefont = font(16),xlims = (0,T) ,titleloc=:center,
    title = "Gas DGP,N $(N) T = $(T) ",size=(1350,600))

##
p = plot()
 for i = 1:Nlines p = plot!(ObsAllT[indT][:,indPl[i]],  color = lineColors[i] ) end #lineColors
 plot(p)

## Box Plots
yLim = (0,1)
    pRT1 = boxplot( rmseTotGas[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N),  GWest = $(GWest)")
    pRT2 = boxplot( rmseTotGas[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
     pRT3 = boxplot( rmseTotGas[:,:,3]  , ylims = yLim , legend = false, title = "T = $(Tvals[3])")
     yLim = (-0.8,1.1)
    pBT1 = boxplot( biasTotGas[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
    pBT2 = boxplot( biasTotGas[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
    pBT3 = boxplot( biasTotGas[:,:,3] , ylims = yLim ,legend = false, title = "T = $(Tvals[3])")
    ptot = plot(pRT1,pRT2,pRT3,pBT1,pBT2,pBT3,titlefont = font(12) ,titleloc=:center,size=(1350,600) )
##
yLim = (0,1)
    pRT1 = boxplot( rmseTotSnap[:,:,1]  , ylims = yLim  , legend = false, title = "T = $(Tvals[1]), N = $(N)  , Single Snapshot ")
    pRT2 = boxplot( rmseTotSnap[:,:,2]  , ylims = yLim  , legend = false, title = "T = $(Tvals[2])")
     pRT3 = boxplot( rmseTotSnap[:,:,3]  , ylims = yLim , legend = false, title = "T = $(Tvals[3])")
     yLim = (-0.8,1.1)
    pBT1 = boxplot( biasTotSnap[:,:,1] , ylims = yLim , legend = false, title = "T = $(Tvals[1])")
    pBT2 = boxplot( biasTotSnap[:,:,2] , ylims = yLim , legend = false, title = "T = $(Tvals[2])")
    pBT3 = boxplot( biasTotSnap[:,:,3] , ylims = yLim ,legend = false, title = "T = $(Tvals[3])")
    ptot = plot(pRT1,pRT2,pRT3,pBT1,pBT2,pBT3,titlefont = font(12) ,titleloc=:center,size=(1350,600) )

#    title = ["($i)" for j = 1:1, i=1:11], titleloc = :right, titlefont = font(8)

## histograms of static parameters
#
# function Bias(A::Array{<:Real,2},B::Array{<:Real,2})
#     return  Bias =    (A.-B)./std(A .-B)
# end
# Bias(A::SharedArray{<:Real,2},B::SharedArray{<:Real,2}) = BiasRel( convert(Array{Real},A),convert(Array{Real},B))
#
# for ind = 1:length(Tvals)
#
# #ind=3
# errA = Bias(simPar_3[ind,:,:],estPar_3[ind,:,:])
# errB = Bias(simPar_2[ind,:,:],estPar_2[ind,:,:])
# errW =  (simPar_1[ind,:,:].- estPar_1[ind,:,:])./(repmat(std(estPar_1[ind,:,:],1),N_est,1))
# plot_tot =Vector{Plots.Plot{Plots.PlotlyBackend}}
# rng_space = -10:.1:10
#
# err= errW
# tmpAsh = ash(err;rng =rng_space )
# pW1 = plot(tmpAsh,
#     title ="  T= $(Tvals[ind]),
#     W   mu =  $(round(mean(err[1,:]),2)) ")
# plot!(Normal(0,1),linewidth = 5)
#
# err= errB
# tmpAsh = ash(err;rng = rng_space)
# pB1 = plot(tmpAsh,
# title =    " $(N_est) Est   B   mu =  $(round(mean(err[1,:]),2)) ")
# plot!(Normal(0,1),linewidth = 5)
#
# err= errA
# rng_space = -3:.1:3
# tmpAsh = ash(err;rng = rng_space)
# pA1 = plot(tmpAsh,
# title = startValFlag * "  $(GWest) / $(GBAest) / $(N) A  mu =  $(round(mean(err[1,:]),2)) ")
# plot!(Normal(0,1),linewidth = 5)
#
# exp = "pl_T_$(ind)=  plot(pW1,pB1,pA1,layout=(3,1) )"
# eval(parse(exp))
#
#
# end
#
#
# exp2 = " " ; for i = 1:ind exp2 *=  "pl_T_$(i) , " end
#
# exp3 = "plot( " * exp2 * " layout = (1,ind),size=(1350,600),legend =:none)"
# eval(parse(exp3))
# #plot(pl_T_1 , pl_T_2,pl_T_3 , pl_T_3, layout = (1,4),size=(1350,600),legend =:none)
#
# ##
#
#
#
#
#
#
# ##


text:
abstract type TargetStyle end
struct UnknownTargetStyle <: TargetStyle end
struct MixedTargetStyle <: TargetStyle end
struct Regression <: TargetStyle end
struct Classification{N} <: TargetStyle end
const BinaryClassification = Classification{2}

ext: jl


text:
function mdl_VScontrol_ode!(device_states,
        output_ode,
        device::DynInverter{C,O,CombinedVIwithVZ,DC,P,F}) where {C <: Converter,
                                                                O <: OuterControl,
                                                                DC<: DCSource,
                                                                P <: FrequencyEstimator,
                                                                F <: Filter}

    #Obtain external states inputs for component
    external_ix = device.input_port_mapping[device.vscontrol]
    iod = device_states[external_ix[1]] #TODO: Should be var referemce?
    ioq = device_states[external_ix[2]] #TODO: Should be state reference?
    icvd = device_states[external_ix[3]]
    icvq = device_states[external_ix[4]]
    vod = device_states[external_ix[5]]
    voq = device_states[external_ix[6]]

    #Obtain inner variables for component
       #vod = device.inner_vars[Vdo_var] #TODO: Should be state reference?
       #voq = device.inner_vars[Vqo_var] #TODO: Should be state reference?
     ω_vsm = device.inner_vars[ω_control_var]
    v_refr = device.inner_vars[v_control_var]
    vdc = device.inner_vars[Vdc_var]

    #Get Voltage Controller parameters
    kpv = device.vscontrol.kpv
    kiv = device.vscontrol.kiv
   kffi = device.vscontrol.kffi
     cf = device.filter.cf #TODO: Is this OK? Should be a getter?
     rv = device.vscontrol.rv
     lv = device.vscontrol.lv

    #Get Current Controller parameters
    kpc = device.vscontrol.kpc
    kic = device.vscontrol.kic
   kffv = device.vscontrol.kffv
     lf = device.filter.lf #TODO: Is this OK? Should be a getter?
    ωad = device.vscontrol.ωad
    kad = device.vscontrol.kad

    #Obtain indices for component w/r to device
    local_ix = device.local_state_ix[device.vscontrol]

    #Define internal states for frequency estimator
    internal_states = @view device_states[local_ix]
    ξ_d = internal_states[1]
    ξ_q = internal_states[2]
    γ_d = internal_states[3]
    γ_q = internal_states[4]
    ϕ_d = internal_states[5]
    ϕ_q = internal_states[6]

    #Inputs (control signals)

    ### Compute 6 states ODEs (D'Arco EPSR122 Model) ###
    ## SRF Voltage Control w/ Virtual Impedance ##
    #Virtual Impedance - Computation but not DAE
    #v_refr = V_ref + kq*(q_ref - qm)
    vod_ref = ( v_refr - rv*iod + ω_vsm*lv*ioq )
    voq_ref = (        - rv*ioq - ω_vsm*lv*iod )
    #Output Control Signal - Links to SRF Current Controller
    icvd_ref = ( kpv*(vod_ref-vod)
               + kiv*ξ_d
               - cf*ω_vsm*voq #j product - change axis
               + kffi*iod )

    icvq_ref = ( kpv*(voq_ref-voq)
               + kiv*ξ_q
               + cf*ω_vsm*vod #j product - change axis
               + kffi*ioq )
    #Voltage Control ODEs
    #PI Integrator (internal state)
    output_ode[local_ix[1]] = (vod_ref-vod)
    output_ode[local_ix[2]] = (voq_ref-voq)

    ## SRF Current Control ##
    #Active Damping
    #vad_d = kad*(vod-ϕ_d)
    #vad_q = kad*(voq-ϕ_q)
    #References for Converter Output Voltage
    vcvd_ref = ( kpc*(icvd_ref-icvd)
               + kic*γ_d
               - ω_vsm*lf*icvq #j product - change axis
               + kffv*vod
               - kad*(vod-ϕ_d))
    vcvq_ref = ( kpc*(icvq_ref-icvq)
               + kic*γ_q
               + ω_vsm*lf*icvd #j product - change axis
               + kffv*voq
               - kad*(voq-ϕ_q) )
    #Modulation Commands to Converter
    #md = vcvd_ref/vdc
    #mq = vcvq_ref/vdc
    #Current Control ODEs
    #PI Integrator (internal state)
    output_ode[local_ix[3]] = icvd_ref - icvd
    output_ode[local_ix[4]] = icvq_ref - icvq
    #Active Damping LPF (internal state)
    output_ode[local_ix[5]] = ωad*vod - ωad*ϕ_d
    output_ode[local_ix[6]] = ωad*voq - ωad*ϕ_q

    #Update inner_vars
    device.inner_vars[md_var] = vcvd_ref/vdc
    device.inner_vars[mq_var] = vcvq_ref/vdc

end


text:
function add_vars!(model::Model, obj::Vector, l_in::Bounds, u_in::Bounds)
    nvars = length(obj)
    obj = fvec(obj)
    l = fvecx(copy(l_in), nvars)
    u = fvecx(copy(u_in), nvars)
    for i = 1:nvars
        if l[i] == -Inf
            l[i] = -CPX_INFBOUND
        end
        if u[i] == Inf
            u[i] = CPX_INFBOUND
        end
    end
    if nvars > 0
        stat = @cpx_ccall(newcols, Cint, (
                          Ptr{Void},
                          Ptr{Void},
                          Cint,
                          Ptr{Cdouble},
                          Ptr{Cdouble},
                          Ptr{Cdouble},
                          Ptr{Cchar},
                          Ptr{Ptr{Cchar}}
                          ),
                          model.env.ptr, model.lp, nvars, float(obj), float(l), float(u), C_NULL, C_NULL)
        if stat != 0
            throw(CplexError(model.env, stat))
        end
    end
end

add_var!(model::Model, obj::Vector, l::Vector, u::Vector) = add_vars!(model, obj, l, u)

function add_var!(model::Model, constridx::IVec, constrcoef::FVec, l::FVec, u::FVec, objcoef::FVec)
    nvars = length(objcoef)
    (nvars == length(l) == length(u)) || error("Inconsistent dimensions when adding variables.")
    for i = 1:nvars
        if l[i] == -Inf
            l[i] = -CPX_INFBOUND
        end
        if u[i] == Inf
            u[i] = CPX_INFBOUND
        end
    end
    if nvars > 0
        stat = @cpx_ccall(addcols, Cint, (
                          Ptr{Void},
                          Ptr{Void},
                          Cint,
                          Cint,
                          Ptr{Cdouble},
                          Ptr{Cint},
                          Ptr{Cint},
                          Ptr{Cdouble},
                          Ptr{Cdouble},
                          Ptr{Cdouble},
                          Ptr{Ptr{Cchar}}
                          ),
                          model.env.ptr, model.lp, nvars, length(constridx), objcoef, Cint[0], constridx.-1, constrcoef, l, u, C_NULL)
        if stat != 0
            throw(CplexError(model.env, stat))
        end
    end
end

add_var!(model::Model, obj, l, u) = add_vars!(model, Cdouble[obj], Cdouble[l], Cdouble[u])

function add_var!(model::Model, constridx, constrcoef, l, u, objcoef)
    return add_var!(model,
                    convert(Vector{Cint},   vec(collect(constridx))),
                    convert(Vector{Cdouble},vec(collect(constrcoef))),
                    convert(Vector{Cdouble},vec(collect(l))),
                    convert(Vector{Cdouble},vec(collect(u))),
                    convert(Vector{Cdouble},vec(collect(objcoef))))
end

function add_var!(model::Model, constridx::Vector, constrcoef::Vector, l::Vector, u::Vector, objcoef::Vector)
    return add_var!(model, ivec(constridx), fvec(constrcoef), fvec(l), fvec(u), fvec(objcoef))
end

function get_varLB(model::Model)
    nvars = num_var(model)
    lb = Array(Cdouble, nvars)
    stat = @cpx_ccall(getlb, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, lb, 0, nvars-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return lb
end

function set_varLB!(model::Model, l::FVec)
    nvars = num_var(model)
    for i = 1:nvars
        if l[i] == -Inf
            l[i] = -CPX_INFBOUND
        end
    end
    stat = @cpx_ccall(chgbds, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cchar},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, nvars, Cint[0:nvars-1;], fill(convert(Cchar, 'L'), nvars), l)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function get_varUB(model::Model)
    nvars = num_var(model)
    ub = Array(Cdouble, nvars)
    stat = @cpx_ccall(getub, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, ub, 0, nvars-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return ub
end

function set_varUB!(model::Model, u::FVec)
    nvars = num_var(model)
        for i = 1:nvars
        if u[i] == Inf
            u[i] = CPX_INFBOUND
        end
    end
    stat = @cpx_ccall(chgbds, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cchar},
                      Ptr{Cdouble}
                      ),
                      model.env.ptr, model.lp, nvars, Cint[0:nvars-1;], fill(convert(Cchar, 'U'), nvars), u)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function set_vartype!(model::Model, vtype::Vector{Char})
    nvars = num_var(model)
    stat = @cpx_ccall(chgctype, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Cchar}
                      ),
                      model.env.ptr, model.lp, length(vtype), Cint[0:length(vtype)-1;], convert(Vector{Cchar},vtype))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    if !isempty(find(!(vtype.=='C')))
        model.has_int = true
    end
end

function get_vartype(model::Model)
    nvars = num_var(model)
    vartypes = Array(Cchar, nvars)
    stat = @cpx_ccall(getctype, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cchar},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, vartypes, 0, nvars-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return convert(Vector{Char},vartypes)
end

function num_var(model::Model)
    nvar = @cpx_ccall(getnumcols, Cint, (
                      Ptr{Void},
                      Ptr{Void}
                      ),
                      model.env.ptr, model.lp)
    return(nvar)
end

function set_varname!(model::Model, idx::Integer, name::AbstractString)
    s = bytestring(name)

    stat = @cpx_ccall(chgcolname, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Ptr{Cint},
                      Ptr{Ptr{UInt8}}
                      ),
                      model.env.ptr, model.lp, 1, Cint[idx-1], [pointer(s)])
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end


text:
export          OrthoPoly,
                MultiOrthoPoly,
                LegendreOrthoPoly,
                JacobiOrthoPoly,
                LaguerreOrthoPoly,
                genLaguerreOrthoPoly,
                HermiteOrthoPoly,
                genHermiteOrthoPoly,
                GaussOrthoPoly,
                Beta01OrthoPoly,
                GammaOrthoPoly,
                LogisticOrthoPoly,
                Uniform01OrthoPoly,
                MeixnerPollaczekOrthoPoly,
                InconsistencyError,
                Quad

struct InconsistencyError <: Exception
    var::String
end

Base.showerror(io::IO, err::InconsistencyError) = print(io, err.var)

function _checkConsistency(deg::Int, Nrec::Int)
    deg < 0 && throw(DomainError(deg, "degree has to be non-negative"))
    Nrec < deg + 1 && throw(DomainError(Nrec, "not enough recurrence coefficients specified (need >= $(deg + 1))"))
end

struct LegendreOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int            # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::LegendreMeasure
    quad::AbstractQuad

    # inner constructor
    function LegendreOrthoPoly(deg::Int; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_legendre(Nrec)
        new(deg, α, β, LegendreMeasure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct JacobiOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int            # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::JacobiMeasure
    quad::AbstractQuad

    # inner constructor
    function JacobiOrthoPoly(deg::Int, shape_a::Real, shape_b::Real; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_jacobi(Nrec, shape_a, shape_b)
        new(deg, α, β, JacobiMeasure(shape_a, shape_b), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct LaguerreOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::LaguerreMeasure
    quad::AbstractQuad

    # inner constructor
    function LaguerreOrthoPoly(deg::Int;Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_laguerre(Nrec)
        new(deg, α, β, LaguerreMeasure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct genLaguerreOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::genLaguerreMeasure
    quad::AbstractQuad

    # inner constructor
    function genLaguerreOrthoPoly(deg::Int, shape::Real ;Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_laguerre(Nrec, shape)
        new(deg, α, β, genLaguerreMeasure(shape), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct HermiteOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::HermiteMeasure
    quad::AbstractQuad

    # inner constructor
    function HermiteOrthoPoly(deg::Int; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_hermite(Nrec)
        new(deg, α, β, HermiteMeasure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct genHermiteOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::genHermiteMeasure
    quad::AbstractQuad

    # inner constructor
    function genHermiteOrthoPoly(deg::Int, mu::Real;Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_hermite(Nrec, mu)
        new(deg, α, β, genHermiteMeasure(mu), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct MeixnerPollaczekOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::MeixnerPollaczekMeasure
    quad::AbstractQuad

    # inner constructor
    function MeixnerPollaczekOrthoPoly(deg::Int, λ::Real, ϕ::Real; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = rm_meixner_pollaczek(Nrec, λ, ϕ)
        new(deg, α, β, MeixnerPollaczekMeasure(λ,ϕ), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct GaussOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::GaussMeasure
    quad::AbstractQuad

    # inner constructor
    function GaussOrthoPoly(deg::Int;Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = r_scale(1/sqrt(2pi), rm_hermite_prob(Nrec)...)
        new(deg, α, β, GaussMeasure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct Uniform01OrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::Uniform01Measure
    quad::AbstractQuad

    # inner constructor
    function Uniform01OrthoPoly(deg::Int;Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = r_scale(1., rm_legendre01(Nrec)...)
        new(deg, α, β, Uniform01Measure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct Beta01OrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::Beta01Measure
    quad::AbstractQuad

    # inner constructor
    function Beta01OrthoPoly(deg::Int, shape_a::Real, shape_b::Real; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = r_scale(1/beta(shape_a, shape_b), rm_jacobi01(Nrec, shape_b-1., shape_a-1.)...)
        new(deg, α, β, Beta01Measure(shape_a, shape_b), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct GammaOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::GammaMeasure
    quad::AbstractQuad

    # inner constructor
    function GammaOrthoPoly(deg::Int, shape::Real, rate::Real; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = r_scale((rate^shape)/gamma(shape), rm_laguerre(Nrec, shape-1.)...)
        new(deg, α, β, GammaMeasure(shape, rate), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct LogisticOrthoPoly <: AbstractCanonicalOrthoPoly
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::LogisticMeasure
    quad::AbstractQuad

    # inner constructor
    function LogisticOrthoPoly(deg::Int; Nrec::Int=deg+1, addQuadrature::Bool = true)
        _checkConsistency(deg, Nrec)
        α, β = r_scale(1., rm_logistic(Nrec)...)
        new(deg, α, β, LogisticMeasure(), addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

struct OrthoPoly <: AbstractOrthoPoly
    name::String
    deg::Int          # maximum degree
    α::Vector{<:Real}  # recurrence coefficients
    β::Vector{<:Real}  # recurrence coefficients
    measure::AbstractMeasure
    quad::AbstractQuad
    # inner constructor
    function OrthoPoly(name::String, deg::Int, α::Vector{<:Real}, β::Vector{<:Real}, measure::AbstractMeasure; addQuadrature::Bool = true)
        deg < 0 && throw(DomainError(deg, "degree has to be non-negative"))
        !(length(α) == length(β)) && throw(InconsistencyError("Inconsistent lengths"))
        new(lowercase(name), deg, α, β, measure, addQuadrature ?  Quad(length(α)-1,α,β) : EmptyQuad())
    end
end

# constructor for known Measure
function OrthoPoly(name::String, deg::Int, measure::AbstractMeasure; Nrec=deg+1, Nquad=10*Nrec, quadrature::Function=clenshaw_curtis, discretization::Function=stieltjes, addQuadrature::Bool=true)
    _checkConsistency(deg, Nrec)
    name = lowercase(name)
    α, β = rm_compute(measure, Nrec, Nquad, quadrature=quadrature, discretization=discretization)
    return OrthoPoly(name, deg, α, β, measure, addQuadrature=addQuadrature)
end

# general constructor
function OrthoPoly(name::String, deg::Int, w::Function, s::Tuple{<:Real,<:Real}, symm::Bool, d::Dict=Dict(); Nrec=deg+1, Nquad=10*Nrec, quadrature::Function=clenshaw_curtis, discretization::Function=stieltjes, addQuadrature::Bool = true)
  name = lowercase(name)
  measure = Measure(name, w, s, symm, d)
  OrthoPoly(name, deg, measure; Nrec=Nrec, Nquad=Nquad, quadrature=quadrature, discretization=discretization, addQuadrature=addQuadrature)
end

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

struct MultiOrthoPoly <: AbstractOrthoPoly
name::Vector{String}
deg::Int
dim::Int
ind::Matrix{<:Int} # multi-index
measure::ProductMeasure
uni::Vector{<:AbstractOrthoPoly}
    function MultiOrthoPoly(uniOrthoPolys::Vector{<:AbstractOrthoPoly}, deg::Int)
      degs = [ op.deg for op in uniOrthoPolys ]
      deg > minimum(degs) && throw(DomainError(deg, "Requested degree $deg is greater than smallest univariate degree $(minimum(degs))."))

      w(t) = prod([op.measure.w(t) for op in uniOrthoPolys])
      measures = [ op.measure for op in uniOrthoPolys ]
      measure = ProductMeasure(w, measures)

      names = [ supertype(typeof(op)) == AbstractCanonicalOrthoPoly ? string(typeof(op)) : op.name for op in uniOrthoPolys ]
      ind = calculateMultiIndices(length(uniOrthoPolys), deg)
      dim = size(ind,1)

      new(names, deg, dim, ind, measure, uniOrthoPolys)
    end
end


text:
using GeoPhyInv
#md using Plots; gr();

# ### Medium
medium_true=Medium(:pizza)

# ### SeisInvExpt
pa=SeisInvExpt(Fdtd(), LS(), :pizza);

# ### Initial plotting
#md p1=heatmap(pa.mediumm, :vp, title="Initial vp")
#md scatter!(pa.ageom, SSrcs())
#md scatter!(pa.ageom, Recs())
#md p2=heatmap(medium_true, :vp, title="True vp")
#md scatter!(pa.ageom, SSrcs())
#md scatter!(pa.ageom, Recs())
#md plot(p2, p1, size=(700,300))

# ### Final plotting
update!(pa,solver=:ipopt);

# ### Plotting #1
#md p1=heatmap(pa.mediumm, :vp, title="Inverted vp")
#md scatter!(pa.ageom, SSrcs())
#md scatter!(pa.ageom, Recs())
#md p2=heatmap(medium_true, :vp, title="True vp")
#md scatter!(pa.ageom, SSrcs())
#md scatter!(pa.ageom, Recs())
#md plot(p2, p1, size=(700,300))


text:
checkbandwidths(N, M, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    M == 1 || (length(u) == M && length(l) == M) ||
    throw(DimensionMismatch("For a matrix of $(N)×$(M) blocks, $(M) lower and upper column bandwidths are required"))

#### Routines for BandedSizes
function bb_blockstarts(ax, l::AbstractVector{Int}, u::AbstractVector{Int})
    N,M = blocksize.(ax,1)
    L,U = maximum(l), maximum(u)
    b_start = BandedMatrix{Int}(undef, (N, M), (L, U))
    -L > U && return b_start

    checkbandwidths(N, M, l, u)

    ind_shift = 0
    for J = 1:M
        KR = Block.(max(1,J-u[J]):min(J+l[J],N))
        if !isempty(KR)
            b_start[Int.(KR),J] .= ind_shift .+ first.(getindex.(Ref(ax[1]),KR)) .- first(ax[1][KR[1]]) .+ 1

            num_rows = length(ax[1][KR])
            num_cols = length(ax[2][Block(J)])
            ind_shift += num_rows*num_cols
        end
    end

    b_start
end

function bb_blockstrides(b_axes, l::AbstractVector{Int}, u::AbstractVector{Int})
    N, M = blocksize.(b_axes,1)
    L,U = maximum(l), maximum(u)
    checkbandwidths(N, M, l, u)
    b_strides = Vector{Int}(undef, M)
    for J=1:M
        KR = Block.(max(1,J-u[J]):min(J+l[J],N))
        if !isempty(KR)
            b_strides[J] = length(b_axes[1][KR])
        else
            b_strides[J] = 0
        end
    end
    b_strides
end

bb_blockstrides(b_axes, l::Integer, u::Integer) = bb_blockstrides(b_axes, Fill(l, blocksize(b_axes,1)), Fill(u, blocksize(b_axes,2)))

struct BlockSkylineSizes{BS<:NTuple{2,AbstractUnitRange{Int}}, LL<:AbstractVector{Int}, UU<:AbstractVector{Int}, BStarts, BStrides}
    axes::BS
    block_starts::BStarts # gives index where the blocks start, usually a BandedMatrix{Int}
    block_strides::BStrides # gives stride to next block for J-th column, usually a Vector{Int}
    l::LL
    u::UU
end

const BlockBandedSizes = BlockSkylineSizes{NTuple{2,BlockedUnitRange{Vector{Int}}}, Fill{Int,1,Tuple{OneTo{Int}}}, Fill{Int,1,Tuple{OneTo{Int}}},
                                            BandedMatrix{Int,Matrix{Int},OneTo{Int}}, Vector{Int}}

BlockSkylineSizes(b_axes::NTuple{2,AbstractUnitRange{Int}}, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    BlockSkylineSizes(b_axes, bb_blockstarts(b_axes, l, u), bb_blockstrides(b_axes, l, u), l, u)

BlockSkylineSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::AbstractVector{Int}, u::AbstractVector{Int}) =
    BlockSkylineSizes((blockedrange(rows),blockedrange(cols)), l, u)

BlockBandedSizes(b_axes::NTuple{2,AbstractUnitRange{Int}}, l::Int, u::Int) =
    convert(BlockBandedSizes, BlockSkylineSizes(b_axes, Fill(l, blocklength(b_axes[2])), Fill(u, blocklength(b_axes[2]))))

BlockBandedSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Int, u::Int) =
    convert(BlockBandedSizes, BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols))))

convert(::Type{BlockSkylineSizes{BS,LL,UU,BStarts,BStrides}}, b::BlockSkylineSizes) where {BS,LL,UU,BStarts,BStrides} =
    BlockSkylineSizes(convert(BS, b.axes), convert(BStarts, b.block_starts), convert(BStrides,b.block_strides), convert(LL,b.l), convert(UU,b.u))

colblockbandwidths(bs::BlockSkylineSizes) = (bs.l, bs.u)

blockstart(block_sizes::BlockSkylineSizes, K, J) = block_sizes.block_starts[K,J]
blockstride(block_sizes::BlockSkylineSizes, J) = block_sizes.block_strides[J]

blockstart(A::AbstractMatrix, K, J) = blockstart(A.block_sizes,K,J)
blockstride(A::AbstractMatrix, J) = blockstride(A.block_sizes,J)

blockbandwidths(bs::BlockSkylineSizes) = maximum.(colblockbandwidths(bs))
blockbandwidth(bs::BlockSkylineSizes, i::Int) = blockbandwidths(bs)[i]

size(bs::BlockSkylineSizes) = map(length,bs.axes)

==(A::BlockSkylineSizes, B::BlockSkylineSizes) = blockisequal(A.axes,B.axes) && A.block_starts == B.block_starts &&
    A.l == B.l && A.u == B.u

colrange(B::BlockSkylineSizes, J::Integer) = max(1, J-B.u[J]):min(blocklength(B.axes[1]), J+B.l[J])

function bb_numentries(B::BlockSkylineSizes)
    axes = B.axes
    numentries = 0
    N = blocklength(axes[1])
    for J = blockaxes(axes[2],1)
        KR = colrange(B, Int(J))
        num_rows = length(axes[1][Block.(KR)])
        num_cols = length(axes[2][J])
        numentries += num_rows*num_cols
    end
    numentries
end

function _BandedBlockMatrix end

#  A block matrix where only the bands are nonzero
#   isomorphic to BandedMatrix{Matrix{T}}
struct BlockSkylineMatrix{T, DATA<:AbstractVector{T}, BS<:BlockSkylineSizes} <: AbstractBlockBandedMatrix{T}
    data::DATA
    block_sizes::BS

    global function _BlockSkylineMatrix(data::DATA, block_sizes::BS) where {T,DATA<:AbstractVector{T}, BS<:BlockSkylineSizes}
        new{T,DATA,BS}(data, block_sizes)
    end
end

const BlockBandedMatrix{T} = BlockSkylineMatrix{T, Vector{T}, BlockBandedSizes}

# Auxiliary outer constructors
@inline _BlockBandedMatrix(data::AbstractVector, bs::BlockBandedSizes) =
    _BlockSkylineMatrix(data, bs)

@inline _BlockBandedMatrix(data::AbstractVector, kr::AbstractVector{Int}, jr::AbstractVector{Int}, (l,u)::NTuple{2, Int}) =
    _BlockBandedMatrix(data, BlockBandedSizes(kr,jr, l,u))

@inline BlockSkylineMatrix{T}(::UndefInitializer, block_sizes::BlockSkylineSizes) where T =
    _BlockSkylineMatrix(Vector{T}(undef, bb_numentries(block_sizes)), block_sizes)

@inline BlockBandedMatrix{T}(::UndefInitializer, block_sizes::BlockBandedSizes) where T =
    _BlockSkylineMatrix(Vector{T}(undef, bb_numentries(block_sizes)), block_sizes)

"""
    BlockSkylineMatrix{T,LL,UU}(M::Union{UndefInitializer,UniformScaling,AbstractMatrix},
                                rows, cols, (l::LL, u::UU))

returns a `sum(rows)`×`sum(cols)` block-banded matrix `A` of type `T`
with block-bandwidths `(l,u)` and where `A[Block(K,J)]` is a
`Matrix{T}` of size `rows[K]`×`cols[J]`.

`(l,u)` may be integers for constant bandwidths or integer vectors of
lengths `rows` and `cols`, respectively, for ragged bands.

# Examples

```jldoctest
julia> BlockSkylineMatrix(I, [2,3,4], [1,2,3], ([2,0,1],[1,1,1]))
3×3-blocked 9×6 BlockSkylineMatrix{Bool,Array{Bool,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 1  │  0  0  │  ⋅  ⋅  ⋅
 0  │  1  0  │  ⋅  ⋅  ⋅
 ───┼────────┼─────────
 0  │  0  1  │  0  0  0
 0  │  0  0  │  1  0  0
 0  │  0  0  │  0  1  0
 ───┼────────┼─────────
 0  │  ⋅  ⋅  │  0  0  1
 0  │  ⋅  ⋅  │  0  0  0
 0  │  ⋅  ⋅  │  0  0  0
 0  │  ⋅  ⋅  │  0  0  0

julia> BlockSkylineMatrix(Ones(9,6), [2,3,4], [1,2,3], ([2,0,1],[1,1,1]))
3×3-blocked 9×6 BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 1.0  │  1.0  1.0  │   ⋅    ⋅    ⋅
 1.0  │  1.0  1.0  │   ⋅    ⋅    ⋅
 ─────┼────────────┼───────────────
 1.0  │  1.0  1.0  │  1.0  1.0  1.0
 1.0  │  1.0  1.0  │  1.0  1.0  1.0
 1.0  │  1.0  1.0  │  1.0  1.0  1.0
 ─────┼────────────┼───────────────
 1.0  │   ⋅    ⋅   │  1.0  1.0  1.0
 1.0  │   ⋅    ⋅   │  1.0  1.0  1.0
 1.0  │   ⋅    ⋅   │  1.0  1.0  1.0
 1.0  │   ⋅    ⋅   │  1.0  1.0  1.0
```
"""
BlockSkylineMatrix

@inline BlockBandedMatrix{T}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2, Int}) where T =
    BlockSkylineMatrix{T}(undef, BlockBandedSizes(axes, lu...))

@inline BlockSkylineMatrix{T}(::UndefInitializer, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2, AbstractVector{Int}}) where T =
    BlockSkylineMatrix{T}(undef, BlockSkylineSizes(axes, lu...))

@inline BlockBandedMatrix{T}(::UndefInitializer, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2, Int}) where T =
    BlockSkylineMatrix{T}(undef, BlockBandedSizes(rdims, cdims, lu...))

@inline BlockSkylineMatrix{T}(::UndefInitializer, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2, AbstractVector{Int}}) where T =
    BlockSkylineMatrix{T}(undef, BlockSkylineSizes(rdims, cdims, lu...))

function BlockSkylineMatrix{T}(A::AbstractMatrix, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
    for J = blockaxes(ret,2), K = blockcolsupport(ret, Int(J))
        kr, jr = getindex.(block_sizes.axes, (K, J))
        view(ret, K, J) .= view(A, kr, jr)
    end
    ret
end

function BlockSkylineMatrix{T}(A::AbstractBlockBandedMatrix, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
    blockisequal(axes(A), block_sizes.axes) || throw(ArgumentError())
    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
        view(ret, K, J) .= view(A, K, J)
    end
    ret
end

BlockBandedMatrix{T}(A::AbstractMatrix, block_sizes::BlockBandedSizes) where T =
    BlockSkylineMatrix{T}(A, block_sizes)

##
# Special cases
##

function BlockSkylineMatrix{T}(Z::Zeros, block_sizes::BlockSkylineSizes) where T
    if size(Z) ≠ map(length,block_sizes.axes)
        throw(DimensionMismatch("Size of input $(size(Z)) must be consistent with $(size(block_sizes))"))
    end
    _BlockSkylineMatrix(zeros(T, bb_numentries(block_sizes)), block_sizes)
end

function BlockSkylineMatrix{T}(E::Eye, block_sizes::BlockSkylineSizes) where T
    if size(E) ≠ size(block_sizes)
        throw(DimensionMismatch("Size of input $(size(E)) must be consistent with $(sum.(dims))"))
    end
    ret = BlockSkylineMatrix(Zeros{T}(size(E)), block_sizes)
    ret[diagind(ret)] .= one(T)
    ret
end

function BlockSkylineMatrix{T}(A::UniformScaling, block_sizes::BlockSkylineSizes) where T
    ret = BlockSkylineMatrix(Zeros{T}(size(block_sizes)), block_sizes)
    ret[diagind(ret)] .= convert(T, A.λ)
    ret
end

BlockSkylineMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::BlockSkylineSizes) = BlockSkylineMatrix{eltype(A)}(A, block_sizes)
BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::BlockBandedSizes) = BlockBandedMatrix{eltype(A)}(A, block_sizes)

BlockSkylineMatrix{T}(A::Union{AbstractMatrix,UniformScaling},
                           rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,AbstractVector{Int}}) where T =
                               BlockSkylineMatrix{T}(A, BlockSkylineSizes(rdims, cdims, lu...))

BlockBandedMatrix{T}(A::Union{AbstractMatrix,UniformScaling},
                           rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,Int}) where T =
                               BlockSkylineMatrix{T}(A, BlockBandedSizes(rdims, cdims, lu...))

BlockSkylineMatrix(A::Union{AbstractMatrix,UniformScaling},
                        rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                        lu::NTuple{2,AbstractVector{Int}}) = BlockSkylineMatrix{eltype(A)}(A, rdims, cdims, lu)
BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
                        lu::NTuple{2,Int}) = BlockBandedMatrix{eltype(A)}(A, rdims, cdims, lu)

BlockBandedMatrix(A::AbstractMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))

function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
    block_sizes = BlockSkylineSizes(axes(A), colblockbandwidths(A)...)

    copyto!(BlockSkylineMatrix{eltype(A)}(undef, block_sizes), A)
end

function convert(::Type{BlockBandedMatrix}, A::AbstractMatrix)
    convert(BlockSkylineMatrix, A)
end

BlockSkylineMatrix(A::AbstractMatrix) = convert(BlockSkylineMatrix, A)
BlockBandedMatrix(A::AbstractMatrix) = convert(BlockBandedMatrix, A)

similar(A::BlockSkylineMatrix, T::Type=eltype(A), bs::BlockSkylineSizes=A.block_sizes) =
    BlockSkylineMatrix{T}(undef, bs)

axes(A::BlockSkylineMatrix) = A.block_sizes.axes

################################
# BlockSkylineMatrix Interface #
################################

MemoryLayout(::Type{<:BlockSkylineMatrix}) = BlockBandedColumnMajor()
colblockbandwidths(A::BlockSkylineMatrix) = (A.block_sizes.l, A.block_sizes.u)
blockbandwidths(A::BlockSkylineMatrix) = maximum.(colblockbandwidths(A))
BroadcastStyle(::Type{<:BlockSkylineMatrix}) = BlockSkylineStyle()
BroadcastStyle(::Type{<:BlockBandedMatrix}) = BlockBandedStyle()

################################
# AbstractBlockArray Interface #
################################

zeroblock(A::BlockSkylineMatrix, K::Int, J::Int) =
    Matrix(Zeros{eltype(A)}(length.(getindex.(axes(A),(Block(K),Block(J))))))

###########################
# AbstractArray Interface #
###########################

# @inline function Base.similar(block_array::BlockSkylineMatrix{T}, ::Type{T2}) where {T,N,T2}
#     BlockArray(similar(block_array.blocks, Array{T2, N}), copy(block_array.block_sizes))
# end

Base.size(arr::BlockSkylineMatrix) =
    @inbounds return map(length,axes(arr))

@inline function getindex(A::BlockSkylineMatrix, i::Int, j::Int)
    @boundscheck checkbounds(A, i, j)
    bi = findblockindex.(axes(A), (i,j))
    @inbounds v = view(A, block.(bi)...)[blockindex.(bi)...]
    return v
end

@inline function setindex!(A::BlockSkylineMatrix{T}, v, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = findblockindex.(axes(A), (i,j))
    V = view(A, block.(bi)...)
    @inbounds V[blockindex.(bi)...] = convert(T, v)::T
    return v
end

############
# Indexing #
############

# function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int}) where {T,N}
#     for i in 1:N
#         if size(v, i) != blocksize(block_arr.block_sizes, (i, block[i]))
#             throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
#         end
#     end
# end
#
#
@inline function setblock!(A::BlockSkylineMatrix, v, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)

    V = view(A, Block(K), Block(J))
    V .= v
    return A
end

##################
# BlockBandedBlock
#
#   views of the blocks satisfy the Matrix interface, and are memory-compatible
#   with StridedMatrix.
##################

const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{BlockSlice1,BlockSlice1},false}

# gives the columns of parent(V).data that encode the block
_parent_blocks(V::BlockBandedBlock)::Tuple{Int,Int} =
    first(first(parentindices(V)).block.n),first(last(parentindices(V)).block.n)

######################################
# Matrix interface  for Blocks #
######################################

MemoryLayout(::Type{<:BlockBandedBlock}) = ColumnMajor()
Base.elsize(::Type{<:BlockSkylineMatrix{T,R}}) where {T,R} = Base.elsize(R)

function Base.unsafe_convert(::Type{Ptr{T}}, V::BlockBandedBlock{T}) where T
    A = parent(V)
    K,J = _parent_blocks(V)
    Base.unsafe_convert(Ptr{T}, A.data) + sizeof(T)*(blockstart(A,K,J)-1)
end

strides(V::BlockBandedBlock) = (1,parent(V).block_sizes.block_strides[_parent_blocks(V)[2]])

@propagate_inbounds function getindex(V::BlockBandedBlock, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = _parent_blocks(V)
    if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
        b_start = blockstart(A,K,J)
        b_start == 0 && return zero(eltype(V))
        b_stride = blockstride(A,J)
        A.data[b_start + k-1 + (j-1)*b_stride ]
    else
        zero(eltype(V))
    end
end

@propagate_inbounds function setindex!(V::BlockBandedBlock, v, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = _parent_blocks(V)
    if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
        b_start = blockstart(A,K,J)
        # TODO: What to do if b_start == 0 ?
        b_stride = A.block_sizes.block_strides[J]
        A.data[b_start + k-1 + (j-1)*b_stride ] = v
    elseif iszero(v) # allow setindex for 0 data
        v
    else
        throw(BandError(A, J-K))
    end
end

"""
    copy_accommodating_diagonals(A::BlockSkylineMatrix, diagonals::UnitRange{<:Integer})

Return copy of `A`, ensuring blocks are present covering the
`diagonals` as well.
"""
function copy_accommodating_diagonals(A::BlockSkylineMatrix, diagonals::UnitRange{<:Integer}, ::Type{T}=eltype(A)) where T
    checksquareblocks(A)
    bs = A.block_sizes
    l,u = bs.l,bs.u
    ax = first(bs.axes)

    if 0 ∈ diagonals
        l = max.(l, 0)
        u = max.(u, 0)
    end

    for d = extrema(diagonals)
        d == 0 && continue # Already taken care of above
        v = d > 0 ? u : l
        # The j:th element of rows is the row index which the diagonal
        # d covers in the j:th column.
        rows = clamp.(ax .- d, 1, last(ax))
        for (j,i) in enumerate(rows)
            # First we find which block the j:th element of the main
            # diagonal would occupy.
            md_block = findfirst(≥(j), ax.lasts)

            # Next we find which block covers row i
            d_block = findfirst(≥(i), ax.lasts)

            # Finally, we increase the block-bandwidth as necessary
            v[md_block] = max(v[md_block], abs(d_block-md_block))
        end
    end

    BlockSkylineMatrix{T}(A, BlockSkylineSizes((ax,ax), l, u))
end

for op in (:-, :+)
    @eval begin
        function $op(A::BlockSkylineMatrix, I::UniformScaling)
            B = copy_accommodating_diagonals(A, 0:0, Base._return_type(+, Tuple{eltype(A), eltype(I)}))
            @inbounds for i in axes(A, 1)
                B[i,i] = $op(B[i,i], I.λ)
            end
            B
        end
        function $op(I::UniformScaling, A::BlockSkylineMatrix)
            B = copy_accommodating_diagonals($op(A), 0:0, Base._return_type(+, Tuple{eltype(A), eltype(I)}))
            @inbounds for i in axes(A, 1)
                B[i,i] += I.λ
            end
            B
        end

        function $op(A::BlockSkylineMatrix, D::Diagonal)
            B = copy_accommodating_diagonals(A, 0:0, Base._return_type(+, Tuple{eltype(A), eltype(D)}))
            @inbounds for i in axes(A, 1)
                B[i,i] = $op(B[i,i], D.diag[i])
            end
            B
        end
        function $op(D::Diagonal, A::BlockSkylineMatrix)
            B = copy_accommodating_diagonals($op(A), 0:0, Base._return_type(+, Tuple{eltype(A), eltype(D)}))
            @inbounds for i in axes(A, 1)
                B[i,i] += D.diag[i]
            end
            B
        end

        function $op(A::BlockSkylineMatrix, Bd::Bidiagonal)
            B = copy_accommodating_diagonals(A, Bd.uplo == 'U' ? (0:1) : (-1:0),
                                             Base._return_type(+, Tuple{eltype(A), eltype(Bd)}))
            @inbounds for i in axes(A, 1)
                B[i,i] = $op(B[i,i], Bd.dv[i])
            end
            @inbounds for i in 1:size(A, 1)-1
                Bd.uplo == 'U' && (B[i,i+1] = $op(B[i,i+1], Bd.ev[i]))
                Bd.uplo == 'L' && (B[i+1,i] = $op(B[i+1,i], Bd.ev[i]))
            end
            B
        end
        function $op(Bd::Bidiagonal, A::BlockSkylineMatrix)
            B = copy_accommodating_diagonals($op(A), Bd.uplo == 'U' ? (0:1) : (-1:0),
                                             Base._return_type(+, Tuple{eltype(A), eltype(Bd)}))
            @inbounds for i in axes(A, 1)
                B[i,i] += Bd.dv[i]
            end
            @inbounds for i in 1:size(A, 1)-1
                Bd.uplo == 'U' && (B[i,i+1] += Bd.ev[i])
                Bd.uplo == 'L' && (B[i+1,i] += Bd.ev[i])
            end
            B
        end

        function $op(A::BlockSkylineMatrix, T::Tridiagonal)
            B = copy_accommodating_diagonals(A, -1:1, Base._return_type(+, Tuple{eltype(A), eltype(T)}))
            @inbounds for i in axes(A, 1)
                B[i,i] = $op(B[i,i], T.d[i])
            end
            @inbounds for i in 1:size(A, 1)-1
                B[i,i+1] = $op(B[i,i+1], T.du[i])
                B[i+1,i] = $op(B[i+1,i], T.dl[i])
            end
            B
        end
        function $op(T::Tridiagonal, A::BlockSkylineMatrix)
            B = copy_accommodating_diagonals($op(A), -1:1, Base._return_type(+, Tuple{eltype(A), eltype(T)}))
            @inbounds for i in axes(A, 1)
                B[i,i] += T.d[i]
            end
            @inbounds for i in 1:size(A, 1)-1
                B[i,i+1] += T.du[i]
                B[i+1,i] += T.dl[i]
            end
            B
        end

        function $op(A::BlockSkylineMatrix, T::SymTridiagonal)
            B = copy_accommodating_diagonals(A, -1:1, Base._return_type(+, Tuple{eltype(A), eltype(T)}))
            @inbounds for i in axes(A, 1)
                B[i,i] = $op(B[i,i], T.dv[i])
            end
            @inbounds for i in 1:size(A, 1)-1
                B[i,i+1] = $op(B[i,i+1], T.ev[i])
                B[i+1,i] = $op(B[i+1,i], T.ev[i])
            end
            B
        end
        function $op(T::SymTridiagonal, A::BlockSkylineMatrix)
            B = copy_accommodating_diagonals($op(A), -1:1, Base._return_type(+, Tuple{eltype(A), eltype(T)}))
            @inbounds for i in axes(A, 1)
                B[i,i] += T.dv[i]
            end
            @inbounds for i in 1:size(A, 1)-1
                B[i,i+1] += T.ev[i]
                B[i+1,i] += T.ev[i]
            end
            B
        end
    end
end


text:
# runtests.jl

using Test

include("../src/infrastructure.jl")

struct AppliPerson
   id::String
   name::String
   email::String
    #constructors
   AppliPerson(id, name) = new(id, name, "")
   AppliPerson(id, name, email) = new(id, name, email)
end

# TEST MODEL
@testset "AppliPerson id, and name test" begin
    scrooge = AppliPerson("smcd", "Scrooge McDuck", "")
    @test scrooge.id == "smcd"
    @test scrooge.name == "Scrooge McDuck"
    @test scrooge.email == ""
end

@testset "AppliPerson id, name, and email test" begin
    scrooge = AppliPerson("smcd", "Scrooge McDuck", "scrooge@duckcity.com")
    @test scrooge.id == "smcd"
    @test scrooge.name == "Scrooge McDuck"
    @test scrooge.email == "scrooge@duckcity.com"
end

@testset "In-memory database test" begin
    db = connect()
    @test db.file == ":memory:"
end

@testset "Archive AppliPerson in database" begin
    db = DBInterface.connect(SQLite.DB, "test.sqlite")
    daisy = AppliPerson("dd", "Daisy", "")
    id = daisy.id
    name = daisy.name
    email = daisy.email
    archive(db, "persons", [daisy])
    r = retrieve(db, "persons", "key = '" * daisy.id * "'").item
    #r = SQLite.DBInterface.execute( db, "select * from persons where key = '$(daisy.id)'") |> DataFrame
    @test r[1].id == id
    @test r[1].name == name
    @test r[1].email == email

    cmd = `rm test.sqlite`
    run(cmd)
end

@testset "Retrieve AppliPerson from database" begin
    db = DBInterface.connect(SQLite.DB)
    daisy = AppliPerson("dd", "Daisy", "")
    id = daisy.id
    name = daisy.name
    email = daisy.email
    store(db, "persons", [daisy])
    r = retrieve(db, "persons", "id = '" * daisy.id * "'")
    #r = SQLite.DBInterface.execute(db, "select * from persons where id = '$(daisy.id)'") |> DataFrame
    @test r.id[1] == id
    @test r.name[1] == name
    @test r.email[1] == email
end


text:
#### Precipitation component in atmosphere model
abstract type PrecipitationModel end

export NoPrecipitation, RainModel, RainSnowModel

eq_tends(pv::PV, m::PrecipitationModel, ::Flux{O}) where {PV, O} = ()

using ..Microphysics

vars_state(::PrecipitationModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::PrecipitationModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_first_order!(
    ::PrecipitationModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    direction,
) end
function compute_gradient_flux!(
    ::PrecipitationModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end
function flux_second_order!(
    ::PrecipitationModel,
    flux::Grad,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    diffusive::Vars,
    hyperdiffusive::Vars,
) end
function flux_first_order!(::PrecipitationModel, _...) end
function compute_gradient_argument!(
    ::PrecipitationModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

source!(::PrecipitationModel, args...) = nothing

"""
    PrecipitationFlux{PV <: Union{Rain, Snow}} <: TendencyDef{Flux{FirstOrder}, PV}

Computes the precipitation flux as a sum of air velocity and terminal velocity
multiplied by the advected variable.
"""
struct PrecipitationFlux{PV <: Union{Rain, Snow}} <:
       TendencyDef{Flux{FirstOrder}, PV} end

PrecipitationFlux() = (PrecipitationFlux{Rain()}, PrecipitationFlux{Snow()})

function flux(::PrecipitationFlux{Rain}, m, state, aux, t, ts, direction)
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_rai = state.precipitation.ρq_rai / state.ρ

    v_term_rai::FT = FT(0)
    if q_rai > FT(0)
        v_term_rai = terminal_velocity(
            m.param_set,
            m.param_set.microphys.rai,
            state.ρ,
            q_rai,
        )
    end

    k̂ = vertical_unit_vector(m, aux)
    return state.precipitation.ρq_rai * (u - k̂ * v_term_rai)
end
function flux(::PrecipitationFlux{Snow}, m, state, aux, t, ts, direction)
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_sno = state.precipitation.ρq_sno / state.ρ

    v_term_sno::FT = FT(0)
    if q_sno > FT(0)
        v_term_sno = terminal_velocity(
            m.param_set,
            m.param_set.microphys.sno,
            state.ρ,
            q_sno,
        )
    end

    k̂ = vertical_unit_vector(m, aux)
    return state.precipitation.ρq_sno * (u - k̂ * v_term_sno)
end

"""
    NoPrecipitation <: PrecipitationModel

No precipitation.
"""
struct NoPrecipitation <: PrecipitationModel end

"""
    RainModel <: PrecipitationModel

Precipitation model with rain.
"""
struct RainModel <: PrecipitationModel end

vars_state(::RainModel, ::Prognostic, FT) = @vars(ρq_rai::FT)
vars_state(::RainModel, ::Gradient, FT) = @vars(q_rai::FT)
vars_state(::RainModel, ::GradientFlux, FT) = @vars(∇q_rai::SVector{3, FT})

function atmos_nodal_update_auxiliary_state!(
    precip::RainModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
end

function compute_gradient_flux!(
    precip::RainModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
end

function flux_first_order!(
    precip::RainModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    direction,
)
    tend = Flux{FirstOrder}()
    args = (atmos, state, aux, t, ts, direction)
    flux.precipitation.ρq_rai = Σfluxes(eq_tends(Rain(), atmos, tend), args...)
end

function flux_second_order!(
    precip::RainModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
)
    d_q_rai = (-D_t) .* diffusive.precipitation.∇q_rai

    flux_second_order!(precip, flux, state, d_q_rai)
end
function flux_second_order!(precip::RainModel, flux::Grad, state::Vars, d_q_rai)
    flux.precipitation.ρq_rai += d_q_rai * state.ρ
end

function source!(m::RainModel, source::Vars, atmos::AtmosModel, args)
    tend = Source()
    source.precipitation.ρq_rai =
        Σsources(eq_tends(Rain(), atmos, tend), atmos, args)
end

"""
    RainSnowModel <: PrecipitationModel

Precipitation model with rain and snow.
"""
struct RainSnowModel <: PrecipitationModel end

vars_state(::RainSnowModel, ::Prognostic, FT) = @vars(ρq_rai::FT, ρq_sno::FT)
vars_state(::RainSnowModel, ::Gradient, FT) = @vars(q_rai::FT, q_sno::FT)
vars_state(::RainSnowModel, ::GradientFlux, FT) =
    @vars(∇q_rai::SVector{3, FT}, ∇q_sno::SVector{3, FT})

function atmos_nodal_update_auxiliary_state!(
    precip::RainSnowModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function compute_gradient_argument!(
    precip::RainSnowModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.precipitation.q_rai = state.precipitation.ρq_rai * ρinv
    transform.precipitation.q_sno = state.precipitation.ρq_sno * ρinv
end

function compute_gradient_flux!(
    precip::RainSnowModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.precipitation.∇q_rai = ∇transform.precipitation.q_rai
    diffusive.precipitation.∇q_sno = ∇transform.precipitation.q_sno
end

function flux_first_order!(
    precip::RainSnowModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    direction,
)
    tend = Flux{FirstOrder}()
    args = (atmos, state, aux, t, ts, direction)
    flux.precipitation.ρq_rai = Σfluxes(eq_tends(Rain(), atmos, tend), args...)
    flux.precipitation.ρq_sno = Σfluxes(eq_tends(Snow(), atmos, tend), args...)
end

function flux_second_order!(
    precip::RainSnowModel,
    flux::Grad,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
    ts,
    diffusive::Vars,
    hyperdiffusive::Vars,
)
    tend = Flux{SecondOrder}()
    args = (atmos, state, aux, t, ts, diffusive, hyperdiffusive)
    flux.precipitation.ρq_rai = Σfluxes(eq_tends(Rain(), atmos, tend), args...)
    flux.precipitation.ρq_sno = Σfluxes(eq_tends(Snow(), atmos, tend), args...)
end

function source!(m::RainSnowModel, source::Vars, atmos::AtmosModel, args)
    tend = Source()
    source.precipitation.ρq_rai =
        Σsources(eq_tends(Rain(), atmos, tend), atmos, args)
    source.precipitation.ρq_sno =
        Σsources(eq_tends(Snow(), atmos, tend), atmos, args)
end

#####
##### Tendency specifications
#####

eq_tends(pv::PV, ::RainModel, ::Flux{FirstOrder}) where {PV <: Rain} =
    (PrecipitationFlux{PV}(),)
eq_tends(pv::PV, ::RainModel, ::Flux{SecondOrder}) where {PV <: Rain} =
    (Diffusion{PV}(),)

eq_tends(
    pv::PV,
    ::RainSnowModel,
    ::Flux{FirstOrder},
) where {PV <: Precipitation} = (PrecipitationFlux{PV}(),)
eq_tends(
    pv::PV,
    ::RainSnowModel,
    ::Flux{SecondOrder},
) where {PV <: Precipitation} = (Diffusion{PV}(),)


text:
using StructuralEquationModels, Test, Statistics
import StructuralEquationModels: obs_cov, obs_mean, get_data

### model specification --------------------------------------------------------------------

spec = ParameterTable(nothing)
spec.variables[:observed_vars] = [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8]
spec.variables[:latent_vars] = [:ind60, :dem60, :dem65]

### data -----------------------------------------------------------------------------------

dat = example_data("political_democracy")
dat_missing = example_data("political_democracy_missing")[:, names(dat)]

dat_matrix = Matrix(dat)
dat_missing_matrix = Matrix(dat_missing)

dat_cov = Statistics.cov(dat_matrix)
dat_mean = vcat(Statistics.mean(dat_matrix, dims = 1)...)

############################################################################################
### tests - SemObservedData
############################################################################################

# w.o. means -------------------------------------------------------------------------------

# errors
@test_throws ArgumentError(
    "You passed your data as a `DataFrame`, but also specified `obs_colnames`. "*
    "Please make shure the column names of your data frame indicate the correct variables "*
    "or pass your data in a different format.") begin
    SemObservedData(specification = spec, data = dat, obs_colnames = Symbol.(names(dat)))
end

@test_throws ArgumentError(
    "Your `data` can not be indexed by symbols. "*
    "Maybe you forgot to provide column names via the `obs_colnames = ...` argument.") begin
    SemObservedData(specification = spec, data = dat_matrix)
end

@test_throws ArgumentError("please specify `obs_colnames` as a vector of Symbols") begin
    SemObservedData(specification = spec, data = dat_matrix, obs_colnames = names(dat))
end

@test_throws UndefKeywordError(:data) SemObservedData(specification = spec)

@test_throws UndefKeywordError(:specification) SemObservedData(data = dat_matrix)

# should work
observed = SemObservedData(
    specification = spec,
    data = dat
)

observed_nospec = SemObservedData(
    specification = nothing,
    data = dat_matrix
)

observed_matrix = SemObservedData(
    specification = spec,
    data = dat_matrix,
    obs_colnames = Symbol.(names(dat))
)

all_equal_cov =
    (obs_cov(observed) == obs_cov(observed_nospec)) &
    (obs_cov(observed) == obs_cov(observed_matrix))

all_equal_data =
    (get_data(observed) == get_data(observed_nospec)) &
    (get_data(observed) == get_data(observed_matrix))

@testset "unit tests | SemObservedData | input formats" begin
    @test all_equal_cov
    @test all_equal_data
end

# shuffle variables
new_order = [3,2,7,8,5,6,9,11,1,10,4]

shuffle_names = Symbol.(names(dat))[new_order]

shuffle_dat = dat[:, new_order]

shuffle_dat_matrix = dat_matrix[:, new_order]

observed_shuffle = SemObservedData(
    specification = spec,
    data = shuffle_dat
)

observed_matrix_shuffle = SemObservedData(
    specification = spec,
    data = shuffle_dat_matrix,
    obs_colnames = shuffle_names
)

all_equal_cov_suffled =
    (obs_cov(observed) == obs_cov(observed_shuffle)) &
    (obs_cov(observed) == obs_cov(observed_matrix_shuffle))

all_equal_data_suffled =
    (get_data(observed) == get_data(observed_shuffle)) &
    (get_data(observed) == get_data(observed_matrix_shuffle))

@testset "unit tests | SemObservedData | input formats shuffled " begin
    @test all_equal_cov_suffled
    @test all_equal_data_suffled
end

# with means -------------------------------------------------------------------------------

# errors
@test_throws ArgumentError(
    "You passed your data as a `DataFrame`, but also specified `obs_colnames`. "*
    "Please make shure the column names of your data frame indicate the correct variables "*
    "or pass your data in a different format.") begin
    SemObservedData(specification = spec, data = dat, obs_colnames = Symbol.(names(dat)), meanstructure = true)
end

@test_throws ArgumentError(
    "Your `data` can not be indexed by symbols. "*
    "Maybe you forgot to provide column names via the `obs_colnames = ...` argument.") begin
    SemObservedData(specification = spec, data = dat_matrix, meanstructure = true)
end

@test_throws ArgumentError("please specify `obs_colnames` as a vector of Symbols") begin
    SemObservedData(specification = spec, data = dat_matrix, obs_colnames = names(dat), meanstructure = true)
end

@test_throws UndefKeywordError(:data) SemObservedData(specification = spec, meanstructure = true)

@test_throws UndefKeywordError(:specification) SemObservedData(data = dat_matrix, meanstructure = true)

# should work
observed = SemObservedData(
    specification = spec,
    data = dat,
    meanstructure = true
)

observed_nospec = SemObservedData(
    specification = nothing,
    data = dat_matrix,
    meanstructure = true
)

observed_matrix = SemObservedData(
    specification = spec,
    data = dat_matrix,
    obs_colnames = Symbol.(names(dat)),
    meanstructure = true
)

all_equal_mean =
    (obs_mean(observed) == obs_mean(observed_nospec)) &
    (obs_mean(observed) == obs_mean(observed_matrix))

@testset "unit tests | SemObservedData | input formats - means" begin
    @test all_equal_mean
end

# shuffle variables
new_order = [3,2,7,8,5,6,9,11,1,10,4]

shuffle_names = Symbol.(names(dat))[new_order]

shuffle_dat = dat[:, new_order]

shuffle_dat_matrix = dat_matrix[:, new_order]

observed_shuffle = SemObservedData(
    specification = spec,
    data = shuffle_dat,
    meanstructure = true
)

observed_matrix_shuffle = SemObservedData(
    specification = spec,
    data = shuffle_dat_matrix,
    obs_colnames = shuffle_names,
    meanstructure = true
)

all_equal_mean_suffled =
    (obs_mean(observed) == obs_mean(observed_shuffle)) &
    (obs_mean(observed) == obs_mean(observed_matrix_shuffle))

@testset "unit tests | SemObservedData | input formats shuffled - mean" begin
    @test all_equal_mean_suffled
end

############################################################################################
### tests - SemObservedCovariance
############################################################################################

# w.o. means -------------------------------------------------------------------------------

# errors

@test_throws ArgumentError("observed means were passed, but `meanstructure = false`") begin
    SemObservedCovariance(specification = nothing, obs_cov = dat_cov, obs_mean = dat_mean)
end

@test_throws UndefKeywordError(:specification) SemObservedCovariance(obs_cov = dat_cov)

@test_throws ArgumentError("no `obs_colnames` were specified") begin
    SemObservedCovariance(specification = spec, obs_cov = dat_cov)
end

@test_throws ArgumentError("please specify `obs_colnames` as a vector of Symbols") begin
    SemObservedCovariance(specification = spec, obs_cov = dat_cov, obs_colnames = names(dat))
end

# should work
observed = SemObservedCovariance(
    specification = spec,
    obs_cov = dat_cov,
    obs_colnames = obs_colnames = Symbol.(names(dat))
)

observed_nospec = SemObservedCovariance(
    specification = nothing,
    obs_cov = dat_cov,
    n_obs = 75.0
)

all_equal_cov = (obs_cov(observed) == obs_cov(observed_nospec))

@testset "unit tests | SemObservedCovariance | input formats" begin
    @test all_equal_cov
    @test isnothing(n_obs(observed))
    @test n_obs(observed_nospec) == 75.0
end

# shuffle variables
new_order = [3,2,7,8,5,6,9,11,1,10,4]

shuffle_names = Symbol.(names(dat))[new_order]

shuffle_dat_matrix = dat_matrix[:, new_order]

shuffle_dat_cov = Statistics.cov(shuffle_dat_matrix)

observed_shuffle = SemObservedCovariance(
    specification = spec,
    obs_cov = shuffle_dat_cov,
    obs_colnames = shuffle_names
)

all_equal_cov_suffled = (obs_cov(observed) ≈ obs_cov(observed_shuffle))

@testset "unit tests | SemObservedCovariance | input formats shuffled " begin
    @test all_equal_cov_suffled
end

# with means -------------------------------------------------------------------------------

# errors
@test_throws ArgumentError("`meanstructure = true`, but no observed means were passed") begin
    SemObservedCovariance(specification = spec, obs_cov = dat_cov, meanstructure = true)
end

@test_throws UndefKeywordError SemObservedCovariance(data = dat_matrix, meanstructure = true)

@test_throws UndefKeywordError SemObservedCovariance(obs_cov = dat_cov, meanstructure = true)

@test_throws ArgumentError("`meanstructure = true`, but no observed means were passed") begin
    SemObservedCovariance(
        specification = spec,
        obs_cov = dat_cov,
        obs_colnames = Symbol.(names(dat)),
        meanstructure = true
    )
end

# should work
observed = SemObservedCovariance(
    specification = spec,
    obs_cov = dat_cov,
    obs_mean = dat_mean,
    obs_colnames = Symbol.(names(dat)),
    n_obs = 75.0,
    meanstructure = true
)

observed_nospec = SemObservedCovariance(
    specification = nothing,
    obs_cov = dat_cov,
    obs_mean = dat_mean,
    meanstructure = true,
    n_obs = 75.0
)

all_equal_mean = (obs_mean(observed) == obs_mean(observed_nospec))

@testset "unit tests | SemObservedCovariance | input formats - means" begin
    @test all_equal_mean
end

# shuffle variables
new_order = [3,2,7,8,5,6,9,11,1,10,4]

shuffle_names = Symbol.(names(dat))[new_order]

shuffle_dat = dat[:, new_order]

shuffle_dat_matrix = dat_matrix[:, new_order]

shuffle_dat_cov = Statistics.cov(shuffle_dat_matrix)
shuffle_dat_mean = vcat(Statistics.mean(shuffle_dat_matrix, dims = 1)...)

observed_shuffle = SemObservedCovariance(
    specification = spec,
    obs_cov = shuffle_dat_cov,
    obs_mean = shuffle_dat_mean,
    obs_colnames = shuffle_names,
    n_obs = 75.0,
    meanstructure = true
)

all_equal_mean_suffled = (obs_mean(observed) == obs_mean(observed_shuffle))

@testset "unit tests | SemObservedCovariance | input formats shuffled - mean" begin
    @test all_equal_mean_suffled
end

############################################################################################
### tests - SemObservedMissing
############################################################################################

# errors
@test_throws ArgumentError(
    "You passed your data as a `DataFrame`, but also specified `obs_colnames`. "*
    "Please make shure the column names of your data frame indicate the correct variables "*
    "or pass your data in a different format.") begin
    SemObservedMissing(specification = spec, data = dat_missing, obs_colnames = Symbol.(names(dat)))
end

@test_throws ArgumentError(
    "Your `data` can not be indexed by symbols. "*
    "Maybe you forgot to provide column names via the `obs_colnames = ...` argument.") begin
    SemObservedMissing(specification = spec, data = dat_missing_matrix)
end

@test_throws ArgumentError("please specify `obs_colnames` as a vector of Symbols") begin
    SemObservedMissing(specification = spec, data = dat_missing_matrix, obs_colnames = names(dat))
end

@test_throws UndefKeywordError(:data) SemObservedMissing(specification = spec)

@test_throws UndefKeywordError(:specification) SemObservedMissing(data = dat_missing_matrix)

# should work
observed = SemObservedMissing(
    specification = spec,
    data = dat_missing
)

observed_nospec = SemObservedMissing(
    specification = nothing,
    data = dat_missing_matrix
)

observed_matrix = SemObservedMissing(
    specification = spec,
    data = dat_missing_matrix,
    obs_colnames = Symbol.(names(dat))
)

all_equal_data =
    isequal(get_data(observed), get_data(observed_nospec)) &
    isequal(get_data(observed), get_data(observed_matrix))

@testset "unit tests | SemObservedMissing | input formats" begin
    @test all_equal_data
end

# shuffle variables
new_order = [3,2,7,8,5,6,9,11,1,10,4]

shuffle_names = Symbol.(names(dat))[new_order]

shuffle_dat_missing = dat_missing[:, new_order]

shuffle_dat_missing_matrix = dat_missing_matrix[:, new_order]

observed_shuffle = SemObservedMissing(
    specification = spec,
    data = shuffle_dat_missing
)

observed_matrix_shuffle = SemObservedMissing(
    specification = spec,
    data = shuffle_dat_missing_matrix,
    obs_colnames = shuffle_names
)

all_equal_data_shuffled =
    isequal(get_data(observed), get_data(observed_shuffle)) &
    isequal(get_data(observed), get_data(observed_matrix_shuffle))

@testset "unit tests | SemObservedMissing | input formats shuffled " begin
    @test all_equal_data_suffled
end


text:
"""
    Separable{Tl<:Kernel, Tr<:Kernel} <: Kernel

The kernel `k` given by
```julia
k((xl, xr), (yl, yr)) = k.l(xl, yl) * k.r(xr, yr)
```
"""
struct Separable{Tl<:Kernel, Tr<:Kernel} <: Kernel
    l::Tl
    r::Tr
end

# Unary methods.
function KernelFunctions.kernelmatrix_diag(
    k::Separable, x::AbstractVector{<:Tuple{Any, Any}},
)
    return kernelmatrix_diag(k.l, first.(x)) .* kernelmatrix_diag(k.r, last.(x))
end
function KernelFunctions.kernelmatrix(
    k::Separable, x::AbstractVector{<:Tuple{Any, Any}},
)
    return kernelmatrix(k.l, first.(x)) .* kernelmatrix(k.r, last.(x))
end

# Binary methods.
function KernelFunctions.kernelmatrix_diag(
    k::Separable,
    x::AbstractVector{<:Tuple{Any, Any}},
    y::AbstractVector{<:Tuple{Any, Any}},
)
    return kernelmatrix_diag(k.l, first.(x), first.(y)) .*
        kernelmatrix_diag(k.r, last.(x), last.(y))
end
function KernelFunctions.kernelmatrix(
    k::Separable,
    x::AbstractVector{<:Tuple{Any, Any}},
    y::AbstractVector{<:Tuple{Any, Any}},
)
    return kernelmatrix(k.l, first.(x), first.(y)) .* kernelmatrix(k.r, last.(x), last.(y))
end


text:
using Documenter, DocumenterMarkdown
using MonteCarloObservable

makedocs(modules = [MonteCarloObservable], doctest=false, format = :markdown)

deploydocs(
    deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4", "python-markdown-math",
        "pygments", "pymdown-extensions"),
    repo   = "github.com/crstnbr/MonteCarloObservable.jl.git",
    target = "site",
    make   = () -> run(`mkdocs build`)
)


text:
struct MockField{T<:Number} <: Field
  v::T
end

function evaluate!(c,f::MockField,x::Point)
  f.v
end

function evaluate!(cache,f::FieldGradient{1,<:MockField},x::Point)
  zero(outer(x,f.object.v))
end

function evaluate!(cache,f::FieldGradient{2,<:MockField},x::Point)
  zero(outer(x,outer(x,f.object.v)))
end

testvalue(::Type{MockField{T}}) where T = MockField(zero(T))

# This is to emulate MonomialBasis that only implements at the level of array of fields.
struct MockFieldArray{N,A} <: AbstractArray{GenericField{Nothing},N}
  values::A
  function MockFieldArray(values::AbstractArray{<:Number})
    N = ndims(values)
    A = typeof(values)
    new{N,A}(values)
  end
end

Base.size(a::MockFieldArray) = size(a.values)
Base.IndexStyle(::Type{<:MockFieldArray}) = IndexLinear()
@inline Base.getindex(a::MockFieldArray,i::Integer) = GenericField(nothing)

function return_cache(f::MockFieldArray,x::Point)
  nothing
end

function evaluate!(c,f::MockFieldArray,x::Point)
  f.values
end

function return_cache(f::MockFieldArray,x::AbstractVector{<:Point})
  T = eltype(f.values)
  CachedArray(zeros(T,length(x),length(f.values)))
end

function evaluate!(c,f::MockFieldArray,x::AbstractVector{<:Point})
  s = (length(x),length(f.values))
  setsize!(c,s)
  r = c.array
  for p in 1:length(x)
    for i in eachindex(f.values)
      r[p,i] = f.values[i]
    end
  end
  r
end

function return_cache(f::FieldGradientArray{1,<:MockFieldArray},x::Point)
  T = return_type(outer,x,testitem(f.fa.values))
  r = zeros(T,size(f.fa.values))
  CachedArray(r)
end

function evaluate!(cache,f::FieldGradientArray{1,<:MockFieldArray},x::Point)
  if size(cache) != size(f.fa.values)
    setsize!(cache,size(f.fa.values))
    r = cache.array
    fill!(r,zero(eltype(r)))
  end
  cache.array
end

function return_cache(f::FieldGradientArray{1,<:MockFieldArray},x::AbstractVector{<:Point})
  T = return_type(outer,testitem(x),testitem(f.fa.values))
  r = zeros(T,length(x),length(f.fa.values))
  CachedArray(r)
end

function evaluate!(cache,f::FieldGradientArray{1,<:MockFieldArray},x::AbstractVector{<:Point})
  setsize!(cache,(length(x),length(f.fa.values)))
  r = cache.array
  fill!(r,zero(eltype(r)))
  r
end


text:
#=
test_plots:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-12
=#

module TestPlots
    # packages
    using Test

    # external modules
    include("plots_2d.jl")
    include("plots_3d.jl")
    include("subplot_2d.jl")
    include("animation_2d.jl")

    # methods
    function test()
        @testset "Plots" begin
            @testset "Plots2D" begin
                Plots2D.set_show_plot(false)
                @test_nowarn Plots2D.main()
            end
            @testset "Plots3D" begin
                Plots3D.set_show_plot(false)
                @test_nowarn Plots3D.main()
            end
            @testset "SubPlot2D" begin
                SubPlot2D.set_show_plot(false)
                @test_nowarn SubPlot2D.main()
            end
            @testset "Anime2D" begin
                Anime2D.set_show_plot(false)
                @test_nowarn Anime2D.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestPlots
    TestPlots.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
import Base: show
import GPUArrays: device

struct CLContext
    context::cl.Context
    device::cl.Device
    queue::cl.CmdQueue
end

function show(io::IO, ctx::CLContext)
    println(io, "OpenCL context with:")
    println(io, "CL version: ", cl.info(ctx.device, :version))
    GPUArrays.device_summary(io, ctx.device)
end

function CLContext(device::cl.Device)
    context = cl.Context(device)
    queue = cl.CmdQueue(context)
    CLContext(context, device, queue)
end

new_context(dev::cl.Device) = CLContext(dev)
function default_device()
    devs = sort(devices(), by = x-> !is_gpu(x))
    if isempty(devs)
        error("No OpenCL devices found")
    end
    first(devs)
end

global all_contexts, global_context, current_device, getcontext!

global_context(A::CLArray) = context(A)

global_queue(A::CLArray) = global_queue(context(A))

function global_queue(cl_ctx::cl.Context = global_context())
    getcontext!(cl_ctx).queue
end

function device(x::CLArray)
    getcontext!(context(x)).device
end
function device(x::cl.Context)
    getcontext!(x).device
end

global init
let contexts = Dict{cl.Device, CLContext}(), active_device = cl.Device[], clcontext2context = Dict{cl.Context, CLContext}()
    function getcontext!(ctx::cl.Context)
        if haskey(clcontext2context, ctx)
            return clcontext2context[ctx]
        else
            device = first(cl.devices(ctx))
            getcontext!(device)
        end
    end
    function getcontext!(device::cl.Device)
        get!(contexts, device) do
            ctx = new_context(device)
            clcontext2context[ctx.context] = ctx
            ctx
        end
    end
    all_contexts() = values(contexts)
    function current_device()
        device = default_device()
        if isempty(active_device)
            push!(active_device, device)
        end
        active_device[]
    end
    function global_context()
        device = current_device()
        getcontext!(device).context
    end
    function init(device::cl.Device)
        if isempty(active_device)
            push!(active_device, device)
        else
            active_device[] = device
        end
        ctx = get!(()-> new_context(device), contexts, device)
        ctx
    end

    function destroy!(context::CLContext)
        # don't destroy primary device context
        dev = context.device
        if haskey(contexts, dev) && contexts[dev] == context
            error("Trying to destroy primary device context which is prohibited. Please use reset!(context)")
        end
        finalize(context.ctx)
        return
    end
end

function reset!(context::CLContext)
    device = context.device
    finalize(context.context)
    context.context = cl.Context(device)
    context.queue = cl.CmdQueue(context.context)
    return
end


text:
using ProbabilityDistributions
using Test

using ProbabilityDistributions: Normal, ∂Normal, Normalc, ∂Normalc
using DistributionParameters, PaddedMatrices, LinearAlgebra
using SIMDPirates: vsum

# quarter of a GiB
# up to 1<<(28-3) == 33_554_432 doubles
const STACK_POINTER = PaddedMatrices.StackPointer(Libc.malloc(1<<28));

function clean_ret(f::F, args::Vararg{<:Any,K}) where {F,K}
    track = Val{(ntuple(Val(K)) do k true end)}()
    sp2, vt = f(STACK_POINTER, args..., track)
    vsum(vt)
end
function clean_ret(f::F, args::Tuple, ::Val{track}) where {F,track}
    sp2, vt = f(STACK_POINTER, args..., Val{track}())
    vsum(vt)
end
@generated function ∂clean_ret(f::F, args::Vararg{<:Any,K}) where {F,K}
    ∂args = Expr(:tuple, (gensym() for k ∈ 1:K)...)
    ∂argsc = copy(∂args)
    pushfirst!(∂argsc.args, :vt)
    pushfirst!(∂args.args, :(vsum(vt)))
    quote
        track = Val{$(Expr(:tuple,(true for k ∈ 1:K)...))}()
        sp2, $∂argsc = f(STACK_POINTER, args..., track)
        $∂args
    end
end
@generated function ∂clean_ret(f::F, args::Tuple, ::Val{track}) where {F,track}
    K = length(args.parameters)
    num_track = sum(track)
    ∂args = Expr(:tuple, (gensym() for k ∈ 1:num_track)...)
    ∂argsc = copy(∂args)
    pushfirst!(∂argsc.args, :vt)
    pushfirst!(∂args.args, :(vsum(vt)))
    quote
        sp2, $∂argsc = f(STACK_POINTER, args..., Val{track}())
        $∂args
    end
end

rel_error(x, y) = (x == y ? zero(x) : (x - y) / y)
@generated function test_grad(f, ∂f, t1, a, t2, tol = 1e-5)
    K1 = length(t1.parameters)
    K2 = length(t2.parameters)
    track = Expr(:tuple,append!(push!(fill(false, K1), true), fill(false, K2))...)
    argstup = Expr(:tuple, [:(t1[$k]) for k ∈ 1:K1]..., :a, [:(t2[$k]) for k ∈ 1:K2]...)
    call = Expr(:call, :clean_ret, :f, :argstuple, :(Val{$track}()))
    ∂call = Expr(:call, :∂clean_ret, :∂f, :argstuple, :(Val{$track}()))
    quote
        @inbounds argstuple = $argstup
        dbase, gptr = $∂call
        g = copy(gptr)
        for i ∈ eachindex(a)
            aᵢ = a[i]
            step = max(cbrt(eps(aᵢ)), 1e-20)
            a[i] = aᵢ + step
            dforward = $call
            a[i] = aᵢ - step
            dback = $call

            gad = g[i]
            gforward = (dforward - dbase) / step
            gback = (dbase - dback) / step
            gcenter = (dforward - dback) / (2step)

            if abs(gcenter) < 1e-20 && abs(gad) < 1e-20
                continue
            end
            r1, r2, r3 = rel_error.(gad, (gforward, gback, gcenter))
            if abs(r3) < tol
            else
                if sign(r1) != sign(r2)
                    print("; i = $i: $r3")
                else
                    @show i, (gad, (gforward,gback,gcenter)),(r1,r2,r3)
                end
            end
            a[i] = aᵢ
        end
    end
end

@testset "ProbabilityDistributions.jl" begin
    # Write your own tests here.

    M = 100
    Σ = DistributionParameters.MutableFixedSizeCovarianceMatrix{M,Float64}(undef);
    fill!(Σ, 0.0);
    Σcopy = DistributionParameters.MutableFixedSizeCovarianceMatrix{M,Float64}(undef);
    A = randn(M, 2M);
    BLAS.syrk!('L', 'N', 1.0, A, 0.0, Σ);

    L = MutableFixedSizeMatrix{M,M,Float64,M}(undef);
    copyto!(L, Σ);
    LAPACK.potrf!('L', L);

    μ = @Mutable randn(M);
    Y = LowerTriangular(L) * randn(M, 50) .+ μ;

    function mvnormdirect(Y, μ, Σ0)
        Σ1 = Symmetric(Array(Σ), :L)
        Lc = cholesky!(Σ1)
        L = LowerTriangular(Lc.L)
        LldivY = L \ (Y .- μ)
        -size(Y,2)*logdet(L) - 0.5dot(LldivY, LldivY)
    end
    t0 = mvnormdirect(Y, μ, Σ)

    t1 = clean_ret(Normal, Y, μ, copyto!(Σcopy, Σ))
    t2, ∂Y, ∂μ, ∂Σ = ∂clean_ret(∂Normal, Y, μ, copyto!(Σcopy, Σ))

    t0, t1, t2

    test_grad(Normalc, ∂Normalc, (), Y, (μ, Σ, Σcopy))
    test_grad(Normalc, ∂Normalc, (Y,), μ, (Σ, Σcopy))
    test_grad(Normalc, ∂Normalc, (Y, μ), Σ, (Σcopy,), 1e-2)

end

using LoopVectorization, SIMDPirates, VectorizationBase
function vdot(δ::AbstractArray{T}) where {T}
    starget = vbroadcast(SVec{(VectorizationBase.pick_vector_width(T)),T}, zero(T))
    @vectorize for i ∈ 1:length(δ)
        vδ = δ[i]
        starget = vmuladd( vδ, vδ, starget )
    end
    target = extract_data(starget)
end


text:
"""
    rampedHalfHalfGenerator(gpExperimentInfo::GPExperimentInfo, popSize::Integer,
                            rng::Random.AbstractRNG, RHHFactor::AbstractFloat = 0.5)

Initializes a population for a GP problem by the \"ramped half and half\" method,
which creates part of it using the full method, and the other part using the grow
method.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `gpExperimentInfo::GPExperimentInfo`: information about the GP experiment.
- `popSize::Integer`: size of the population that is going to be created, that is
    to say, the number of individuals of the population.
- `rng::Random.AbstractRNG`: random number generator for random number consistency
    along an experiment.

# User arguments
- `RHHFactor::AbstractFloat = 0.5`: percentage of the population that is going to
    be created by the \"full\" method, the remaining percentage will be created by
    the \"grow\" method.

# Returns
A population of GP individuals generated by the \"ramped half and half\" method.
"""
function rampedHalfHalfGenerator(gpExperimentInfo::GPExperimentInfo, popSize::Integer,
                                 rng::Random.AbstractRNG, RHHFactor::AbstractFloat = 0.5)

    indType = typeof(growGenerator(gpExperimentInfo, rng))
    population = Array{indType}(undef, popSize)
    N = popSize - round(Integer, popSize * RHHFactor)

    for i = 1:(popSize-N)
        population[i] = fullGenerator(gpExperimentInfo, rng)
    end

    for i = (popSize-N+1):popSize
        population[i] = growGenerator(gpExperimentInfo, rng)
    end

    return Random.shuffle(rng, population)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
# # Tutorial 🐧

# This is a gentle and lighthearted tutorial on how to use tools from AlgebraOfGraphics,
# using as example dataset a collection of measurements on penguins[^1]. See
# the [Palmer penguins website](https://allisonhorst.github.io/palmerpenguins/index.html)
# for more information.
#
# [^1]: Gorman KB, Williams TD, Fraser WR (2014) Ecological Sexual Dimorphism and Environmental Variability within a Community of Antarctic Penguins (Genus Pygoscelis). PLoS ONE 9(3): e90081. [DOI](https://doi.org/10.1371/journal.pone.0090081)

using PalmerPenguins, DataFrames

penguins = dropmissing(DataFrame(PalmerPenguins.load()))
first(penguins, 6)

# ## Frequency plots
#
# Let us start by getting a rough idea of how the data is distributed

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

axis = (width = 225, height = 225)
penguin_frequency = data(penguins) * frequency() * mapping(:species)
draw(penguin_frequency; axis)

# Next, let us see whether the distribution is the same across islands.

plt = penguin_frequency * mapping(color = :island)
draw(plt; axis)

# Ups! The bars are in the same spot and are hiding each other. We need to specify
# how we want to fix this. Bars can either `dodge` each other, or be `stack`ed on top
# of each other.

plt = penguin_frequency * mapping(color = :island, dodge = :island)
draw(plt; axis)

# This is our first finding. `Adelie` is the only species of penguins that can be
# found on all three islands. To be able to see both which species is more numerous
# and how different species are distributed across islands in a unique plot,
# we could have used `stack`.

plt = penguin_frequency * mapping(color = :island, stack = :island)
draw(plt; axis)

# ## Correlating two variables
#
# Now that we have understood the distribution of these three penguin species, we can
# start analyzing their features.

penguin_bill = data(penguins) * mapping(:bill_length_mm, :bill_depth_mm)
draw(penguin_bill; axis)

# We would actually prefer to visualize these measures in centimeters, and to have
# cleaner axes labels. As we want this setting to be preserved in all of our `bill`
# visualizations, let us save it in the variable `penguin_bill`, to be reused in
# subsequent plots.

penguin_bill = data(penguins) * mapping(
    :bill_length_mm => (t -> t / 10) => "bill length (cm)",
    :bill_depth_mm => (t -> t / 10) => "bill depth (cm)",
)
draw(penguin_bill; axis)

# Much better! Note the parentheses around the function `t -> t / 10`. They are
# necessary to specify that the function maps `t` to `t / 10`, and not to
# `t / 10 => "bill length (cm)"`.

# There does not seem to be a strong correlation between the two dimensions, which
# is odd. Maybe dividing the data by species will help.

plt = penguin_bill * mapping(color = :species)
draw(plt; axis)

# Ha! Within each species, penguins with a longer bill also have a deeper bill.
# We can confirm that with a linear regression

plt = penguin_bill * linear() * mapping(color = :species)
draw(plt; axis)

# This unfortunately no longer shows our data!
# We can use `+` to plot both things on top of each other:

plt = penguin_bill * linear() * mapping(color = :species) + penguin_bill * mapping(color = :species)
draw(plt; axis)

# Note that the above expression seems a bit redundant, as we wrote the same thing twice.
# We can "factor it out" as follows

plt = penguin_bill * (linear() + mapping()) * mapping(color = :species)
draw(plt; axis)

# where `mapping()` is a neutral multiplicative element.
# Of course, the above could be refactored as

layers = linear() + mapping()
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# We could actually take advantage of the spare `mapping()` and use it to pass some
# extra info to the scatter, while still using all the species members to compute
# the linear fit.

layers = linear() + mapping(marker = :sex)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# This plot is getting a little bit crowded. We could instead analyze female and
# male penguins in separate subplots.

layers = linear() + mapping(col = :sex)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# ## Smooth density plots
#
# An alternative approach to understanding how two variables interact is to consider
# their joint probability density distribution (pdf).

using AlgebraOfGraphics: density
plt = penguin_bill * density(npoints=50) * mapping(col = :species)
draw(plt; axis)

# The default colormap is multi-hue, but it is possible to pass single-hue colormaps as well:

plt *= visual(colormap = :grayC)
draw(plt; axis)

# A `Heatmap` (the default visualization for a 2D density) is a bit unfortunate if
# we want to mark species by color. In that case, one can use `visual` to change
# the default visualization and, optionally, fine tune some arguments.
# In this case, a `Wireframe` with thin lines looks quite nice. (Note that, for the
# time being, we must specify explicitly that we require a 3D axis.)

axis = (type = Axis3, width = 300, height = 300)
layer = density() * visual(Wireframe, linewidth=0.05)
plt = penguin_bill * layer * mapping(color = :species)
draw(plt; axis)

# Of course, a more traditional approach would be to use a `Contour` plot instead:

axis = (width = 225, height = 225)
layer = density() * visual(Contour)
plt = penguin_bill * layer * mapping(color = :species)
draw(plt; axis)

# The data and the linear fit can also be added back to the plot:

layers = density() * visual(Contour) + linear() + mapping()
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# In the case of many layers (contour, density and scatter) it is important to think
# about balance. In the above plot, the markers are quite heavy and can obscure the linear
# fit and the contour lines.
# We can lighten the markers using alpha transparency.

layers = density() * visual(Contour) + linear() + visual(alpha = 0.5)
plt = penguin_bill * layers * mapping(color = :species)
draw(plt; axis)

# ## Correlating three variables
#
# We are now mostly up to speed with `bill` size, but we have not consider how
# it relates to other penguin features, such as their weight.
# For that, a possible approach is to use a continuous color
# on a gradient to denote weight and different marker shapes to denote species.
# Here we use `group` to split the data for the linear regression without adding
# any additional style.

body_mass = :body_mass_g => (t -> t / 1000) => "body mass (kg)"
layers = linear() * mapping(group = :species) + mapping(color = body_mass, marker = :species)
plt = penguin_bill * layers
draw(plt; axis)

# Naturally, within each species, heavier penguins have bigger bills, but perhaps
# counter-intuitively the species with the shallowest bills features the heaviest penguins.
# We could also try and see the interplay of these three variables in a 3D plot.

axis = (type = Axis3, width = 300, height = 300)
plt = penguin_bill * mapping(body_mass, color = :species)
draw(plt; axis)

#

plt = penguin_bill * mapping(body_mass, color = :species, layout = :sex)
draw(plt; axis)

# Note that static 3D plot can be misleading, as they only show one projection
# of 3D data. They are mostly useful when shown interactively.


text:
# Tables.jl interface

Tables.isrowtable(::Type{<:TableRowIterator}) = true
Tables.columnnames(tr::TableRow) = Tuple(tr.index.column_labels)
Tables.getcolumn(tr::TableRow, nm::Symbol) = getdata(tr, nm)

function writetable(filename::AbstractString, x; kw...)

	_as_vector(y::AbstractVector) = y
	_as_vector(y) = collect(y)

    writetable(filename, Any[_as_vector(c) for c in Tables.Columns(x)], collect(Symbol, Tables.columnnames(x)); kw...)
end

function writetable(filename::AbstractString, data_vector::Vector{Tuple{S, T}}; kw...) where {S<:AbstractString, T}
	_as_vector(y::AbstractVector) = y
	_as_vector(y) = collect(y)

	tables = Tuple{AbstractString, Vector{Any}, Vector{Any}}[]
	for sheet in data_vector
		push!(tables, (sheet[1], Any[_as_vector(c) for c in Tables.Columns(sheet[2])], collect(Symbol, Tables.columnnames(sheet[2]))))
	end
	writetable(filename, tables; kw...)
end

function writetable(filename::AbstractString, data_dict::Dict{S, T}; kw...) where {S<:AbstractString, T}
	writetable(filename, [(sheetname, data) for (sheetname, data) in data_dict]; kw...)
end


text:
@time a = rand([randstring(rand(0:16)) for i in 1:100_000], 10_000_000);

@time as = StrF{16}.(a);

@time sort(as, alg=RadixSort) |> issorted

@time asb = radixsort(as, rev=true)
@test issorted(asb, rev=true)

@time asb = radixsort(as)
@test issorted(asb)


text:
# This file is a part of HistTools.jl, licensed under the MIT License (MIT).

using HistTools
using Test