text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), May 23, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Handling missing values
#
# A singelton type `Missings.Missing` allows us to deal with missing values.
# ------------------------------------------------------------------------------------------

missing, typeof(missing)

# ------------------------------------------------------------------------------------------
# Arrays automatically create an appropriate union type.
# ------------------------------------------------------------------------------------------

x = [1, 2, missing, 3]

# ------------------------------------------------------------------------------------------
# `ismissing` checks if passed value is missing.
# ------------------------------------------------------------------------------------------

ismissing(1), ismissing(missing), ismissing(x), ismissing.(x)

# ------------------------------------------------------------------------------------------
# We can extract the type combined with Missing from a `Union` via
#
# (This is useful for arrays!)
# ------------------------------------------------------------------------------------------

eltype(x), Missings.T(eltype(x))

# ------------------------------------------------------------------------------------------
# `missing` comparisons produce `missing`.
# ------------------------------------------------------------------------------------------

missing == missing, missing != missing, missing < missing

# ------------------------------------------------------------------------------------------
# This is also true when `missing`s are compared with values of other types.
# ------------------------------------------------------------------------------------------

1 == missing, 1 != missing, 1 < missing

# ------------------------------------------------------------------------------------------
# `isequal`, `isless`, and `===` produce results of type `Bool`.
# ------------------------------------------------------------------------------------------

isequal(missing, missing), missing === missing, isequal(1, missing), isless(1, missing)

# ------------------------------------------------------------------------------------------
# In the next few examples, we see that many (not all) functions handle `missing`.
# ------------------------------------------------------------------------------------------

map(x -> x(missing), [sin, cos, zero, sqrt]) # part 1

map(x -> x(missing, 1), [+, - , *, /, div]) # part 2

map(x -> x([1,2,missing]), [minimum, maximum, extrema, mean, any, float]) # part 3

# ------------------------------------------------------------------------------------------
# `skipmissing` returns iterator skipping missing values. We can use `collect` and
# `skipmissing` to create an array that excludes these missing values.
# ------------------------------------------------------------------------------------------

collect(skipmissing([1, missing, 2, missing]))

# ------------------------------------------------------------------------------------------
# Similarly, here we combine `collect` and `Missings.replace` to create an array that
# replaces all missing values with some value (`NaN` in this case).
# ------------------------------------------------------------------------------------------

collect(Missings.replace([1.0, missing, 2.0, missing], NaN))

# ------------------------------------------------------------------------------------------
# Another way to do this:
# ------------------------------------------------------------------------------------------

coalesce.([1.0, missing, 2.0, missing], NaN)

# ------------------------------------------------------------------------------------------
# Caution: `nothing` would also be replaced here (for Julia 0.7 a more sophisticated
# behavior of `coalesce` that allows to avoid this problem is planned).
# ------------------------------------------------------------------------------------------

coalesce.([1.0, missing, nothing, missing], NaN)

# ------------------------------------------------------------------------------------------
# You can use `recode` if you have homogenous output types.
# ------------------------------------------------------------------------------------------

recode([1.0, missing, 2.0, missing], missing=>NaN)

# ------------------------------------------------------------------------------------------
# You can use `unique` or `levels` to get unique values with or without missings,
# respectively.
# ------------------------------------------------------------------------------------------

unique([1, missing, 2, missing]), levels([1, missing, 2, missing])

# ------------------------------------------------------------------------------------------
# In this next example, we convert `x` to `y` with `allowmissing`, where `y` has a type that
# accepts missings.
# ------------------------------------------------------------------------------------------

x = [1,2,3]
y = allowmissing(x)

# ------------------------------------------------------------------------------------------
# Then, we convert back with `disallowmissing`. This would fail if `y` contained missing
# values!
# ------------------------------------------------------------------------------------------

z = disallowmissing(y)
x,y,z

# ------------------------------------------------------------------------------------------
# In this next example, we show that the type of each column in `x` is initially `Int64`.
# After using `allowmissing!` to accept missing values in columns 1 and 3, the types of
# those columns become `Union`s of `Int64` and `Missings.Missing`.
# ------------------------------------------------------------------------------------------

x = DataFrame(Int, 2, 3)
println("Before: ", eltypes(x))
allowmissing!(x, 1) # make first column accept missings
allowmissing!(x, :x3) # make :x3 column accept missings
println("After: ", eltypes(x))

# ------------------------------------------------------------------------------------------
# In this next example, we'll use `completecases` to find all the rows of a `DataFrame` that
# have complete data.
# ------------------------------------------------------------------------------------------

x = DataFrame(A=[1, missing, 3, 4], B=["A", "B", missing, "C"])
println(x)
println("Complete cases:\n", completecases(x))

# ------------------------------------------------------------------------------------------
# We can use `dropmissing` or `dropmissing!` to remove the rows with incomplete data from a
# `DataFrame` and either create a new `DataFrame` or mutate the original in-place.
# ------------------------------------------------------------------------------------------

y = dropmissing(x)
dropmissing!(x)
[x, y]

# ------------------------------------------------------------------------------------------
# When we call `showcols` on a `DataFrame` with dropped missing values, the columns still
# allow missing values.
# ------------------------------------------------------------------------------------------

showcols(x)

# ------------------------------------------------------------------------------------------
# Since we've excluded missing values, we can safely use `disallowmissing!` so that the
# columns will no longer accept missing values.
# ------------------------------------------------------------------------------------------

disallowmissing!(x)
showcols(x)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
Pkg.add.(["SampledSignals","DSP","PortAudio","JSON"])
Pkg.add.(["Interpolations","Colors","FixedPointNumbers"])
Pkg.clone("https://github.com/gaika/ArrayFire.jl")

ext: jl


text:
using Documenter
using ChargedParticleDynamics

makedocs(
    sitename = "ChargedParticleDynamics.jl",
    format = Documenter.HTML(
                prettyurls = get(ENV, "CI", nothing) == "true",
                assets = [asset("assets/style.css", class=:css, islocal=true)]),
    pages = ["Overview"                      => "index.md",
             "Normalization"                 => "normalization.md",
             "Initialization"                => "initialization.md",
             "Charged Particles in 3D"       => "charged_particle_3d.md",
             "Pauli Particles in 3D"         => "pauli_particle_3d.md",
             "Guiding Center Dynamics in 4D" => "guiding_center_4d.md",
            ]
)

deploydocs(
    repo   = "github.com/JuliaPlasma/ChargedParticleDynamics.jl",
)


text:
""" Creates and tests the axial flux motor model. """
function afpm_model()
    # Data preprocessing
    X = log.(CSV.read("data/afpm/afpm_inputs.csv", DataFrame,
                 copycols=true, delim=","));
    X_infeas = log.(CSV.read("data/afpm/afpm_infeas_inputs.csv", DataFrame,
                        copycols=true, delim=","));
    X = select!(X, Not(:Column1));
    X_infeas = select!(X_infeas, Not(:Column1));

    Y = CSV.read("data/afpm/afpm_outputs.csv", DataFrame, copycols=true, delim=","); # log taken later
    foms = ["P_shaft", "Torque", "Rotational Speed", "Efficiency", "Mass", "Mass Specific Power"]
    Y = select!(Y, foms);

    feas_idxs = findall(x -> !ismissing(x) && x .>= 0.5, Y.Efficiency);
    X_feas = X[feas_idxs, :]
    Y_feas = log.(Y[feas_idxs, :])

    # Creating model with relevant variables, and geometry constraints
    m = Model(with_optimizer(OCTHaGOn.SOLVER_SILENT))
    add_variables_from_data!(m, X)
    bound_to_data!(m, X_feas)
    add_variables_from_data!(m, Y)
    bound_to_data!(m, Y_feas)
    @assert length(all_variables(m)) == 15
    N_coils, TPC, p = m[:N_coils], m[:TPC], m[:p]

    # Geometry constraints (in logspace)
    D_out, D_in, N_coils, wire_w, wire_A, wire_h = m[:D_out], m[:D_in], m[:N_coils],
                                                   m[:wire_w], m[:wire_A], m[:wire_h]
    @constraint(m, D_out >= D_in)
    @constraint(m, log(pi) + D_in >= log(0.2) + wire_w + N_coils) #pi*D_in >= 2*0.1*wire_w*N_coils
    @constraint(m, wire_A == wire_h + wire_w)
    @constraint(m, N_coils >= p + 1e-3) # motor type 2

    # Objectives and FOMs
    P_shaft, Torque, Rotational_Speed, Efficiency, Mass, Mass_Specific_Power = [m[Symbol(fom)] for fom in foms]
    gm = GlobalModel(model = m)

    # Integer constraints (in logspace)
    restrict_to_set(N_coils, unique(X_feas[!, "N_coils"]));
    restrict_to_set(p, unique(X_feas[!, "p"]))
    restrict_to_set(TPC, unique(X_feas[!, "TPC"]))

    # Feasibility constraint
    add_datadriven_constraint(gm, X, Int64.([!ismissing(y) && y .>= 0.5 for y in Array(Y.Efficiency)]) .* 2 .- 1)
    # FOM regressions
    for fom in foms
        add_datadriven_constraint(gm, X_feas, Array(Y_feas[!, Symbol(fom)]), dependent_var = m[Symbol(fom)], equality=true)
    end

    # Requirements and objective
    @constraint(m, log(7900) <= Rotational_Speed <= log(8100))
    @constraint(m, log(9.8) <= P_shaft <= log(10.2))
    @objective(m, Min, Mass)

    return gm
end


text:
function get_vehicle_category(v_set::HomogeneousVehicleSet,
                              data::RvrpInstance,
                              computed_data::RvrpComputedData)
    return data.vehicle_categories[
        computed_data.vehicle_category_id_2_index[
            v_set.vehicle_category_id
        ]
    ]
end

function get_locations_from_group_id(group_id::String,
                                     data::RvrpInstance,
                                     computed_data::RvrpComputedData)
    return [
        data.locations[computed_data.location_id_2_index[id]]
        for id in data.location_groups[
            computed_data.location_group_id_2_index[group_id]
        ].location_ids
    ]
end

function get_departure_locations(v_set::HomogeneousVehicleSet,
                                 data::RvrpInstance,
                                 computed_data::RvrpComputedData)
    return get_locations_from_group_id(v_set.departure_location_group_id,
                                       data, computed_data)
end

function get_arrival_locations(v_set::HomogeneousVehicleSet,
                                 data::RvrpInstance,
                                 computed_data::RvrpComputedData)
    return get_locations_from_group_id(v_set.arrival_location_group_id,
                                       data, computed_data)
end

function get_pickup_locations(req::Request,
                              data::RvrpInstance,
                              computed_data::RvrpComputedData)
    return get_locations_from_group_id(req.pickup_location_group_id,
                                       data, computed_data)
end

function get_delivery_locations(req::Request,
                              data::RvrpInstance,
                              computed_data::RvrpComputedData)
    return get_locations_from_group_id(req.delivery_location_group_id,
                                       data, computed_data)
end

function get_closest_tw(tws::Vector{Range}, target_time::Float64)
    for tw in tws
        if target_time <= tw.ub
            return tw
        end
    end
    @show tws
    @show target_time
    error("No feasible time windows for target ", target_time, ". Time windows are : ", tws)
end


text:
using Test

const tests = [
    "Correlation",
    "PearsonMatching",
    "GeneralizedSDistribution",
    "RandomVector",
    "Utilities"
]

printstyled("Running tests:\n", color=:blue)

for t in tests
    @testset "Test $t" begin
        include("$t.jl")
    end
end


text:
module Day10

using Combinatorics

export day10

function day10(::Val{:a}, filename)
    input = Vector{Int}()
    for line in readlines(filename)
        push!(input, parse(Int, line))
    end
    sort!(input)
    differences = Dict{Int, Int}([1 => 0, 2 => 0, 3 => 0])
    differences[input[1]] = 1
    @inbounds for i in 1:length(input)-1
        d = input[i+1] - input[i]
        differences[d] += 1
    end
    differences[3] += 1
    differences[1] * differences[3]
end

function choices(slice, before, after)
    if isempty(slice)
        return 1
    end
    count = after - before <= 3 ? 1 : 0
    for t in Base.Iterators.drop(powerset(slice),1)
        if (t[1] - before <= 3) && (after - t[end] <= 3) && isvalid(t)
            count += 1
        end
    end
    count
end

function isvalid(t)
    for i in 1:length(t)-1
        if t[i+1] - t[i] > 3
            return false
        end
    end
    true
end

function day10(::Val{:b}, filename)
    input = Vector{Int}()
    for line in readlines(filename)
        push!(input, parse(Int, line))
    end
    sort!(input)
    indices = Vector{Int}()
    if input[1] == 3
        push!(indices, 1)
    end
    for i in 2:length(input)
        if input[i] - input[i-1] == 3
            push!(indices, i)
            push!(indices, i-1)
        end
        if input[i] - input[i-1] == 2 && input[i+1] - input[1] == 2
            push!(indices, i)
        end
    end
    push!(indices, length(input))
    indices = sort(unique(indices))
    count = choices(input[1:indices[1]-1], 0, input[indices[1]])
    for i in 1:length(indices)-1
        count *= choices(input[indices[i]+1:indices[i+1]-1], input[indices[i]], input[indices[i+1]])
    end
    count
end

end


text:
#=
dynamic_factor_model.jl

    Dynamic factor model specification in state space formulation

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/07
=#

"""
    DynamicFactorModel <: StateSpaceModel

Constructor for a dynamic factor model instance of the state space model type
with hyperparameters `Λ` and `ϕ`, error specification, and mean
specification.
"""
struct DynamicFactorModel{Ty, Tr, Tm, TΛ, Tϕ, Te} <: StateSpaceModel
    y::Ty       # data
    r::Tr       # number of factors
    mean::Tm    # mean specification
    Λ::TΛ       # Loading matrix
    ϕ::Tϕ       # autoregressive parameters
    error::Te   # error specification
end
# Constructors
function DynamicFactorModel(y::AbstractMatrix, r::Integer, mean::AbstractMeanModel, error::AbstractErrorModel)
    # Get dims
    n= size(y,1)

    # hyper paremeters
    Λ= similar(y, n, r)
    ϕ= Diagonal(similar(y, r))

    return DynamicFactorModel(y, r, mean, Λ, ϕ, error)
end

# Methods
mean(model::DynamicFactorModel)= mean(model.mean)
cov(model::DynamicFactorModel)= cov(model.error)
prec(model::DynamicFactorModel)= prec(model.error)
resid(model::DynamicFactorModel)= resid(model.error)

# State space system and hyperparameters
function nparams(model::DynamicFactorModel)
    # Get dims
    (n,r)= size(model.Λ)

    n_params= n*r + r + nparams(model.error) + nparams(model.mean)

    return n_params
end

function get_params!(ψ::AbstractVector, model::DynamicFactorModel)
    # Get dims
    (n,r)= size(model.Λ)
    idx= 0  # index counter

    # Store values
    # Λ
    ψ[1:n*r].= vec(model.Λ)
    idx+= n*r

    # ϕ
    ψ[idx+1:idx+r].= model.ϕ.diag
    idx+= r

    # mean model
    get_params!(view(ψ,idx+1:idx+nparams(model.mean)), model.mean)
    idx+= nparams(model.mean)

    # error model
    get_params!(view(ψ,idx+1:idx+nparams(model.error)), model.error)

    return nothing
end

function get_system!(sys::LinearTimeInvariant, model::DynamicFactorModel, method::Symbol)
    # Store values
    sys.T.= model.ϕ
    if method === :univariate
        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Store values
        ldiv!(sys.y, C.L, model.y)
        ldiv!(sys.Z, C.L, model.Λ)
        if !isa(model.mean, NoConstant)
            ldiv!(sys.d, C.L, mean(model))
        end
    elseif method === :collapsed
        # Collapsing transformation
        # Z'×H⁻¹×Z
        tmp_pn= transpose(model.Λ) * prec(model)
        tmp= tmp_pn * model.Λ
        # (Z'×H⁻¹×Z)⁻¹
        # perform pseudo inverse as Λ can contain zero columns
        pseudo= pinv(tmp)
        # pivoted Cholesky decomposition
        C= cholesky!(Hermitian(pseudo), Val(true), check=false)
        # transformation
        ip= invperm(C.p)
        U= C.U[:,ip]
        A= U * tmp_pn

        # Store values
        mul!(sys.y, A, model.y)
        sys.Z.= pinv(transpose(U))
        if !isa(model.mean, NoConstant)
            mul!(sys.d, A, mean(model))
        end
    else
        # Store values
        sys.Z.= model.Λ
        sys.H.= cov(model)
        if !isa(model.mean, NoConstant)
            sys.d.= mean(model)
        end
    end

    return nothing
end

function get_system!(sys::LinearTimeVariant, model::DynamicFactorModel, method::Symbol)
    # get dims
    T_len= size(model.y,2)

    # Store values
    @inbounds for t in 1:T_len
        sys.T[t].= model.ϕ
    end
    if method === :univariate
        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Store values
        ldiv!(sys.y, C.L, model.y)
        ldiv!(sys.Z[1], C.L, model.Λ)
        @inbounds for t in 2:T_len
            sys.Z[t].= sys.Z[1]
        end
        if !isa(model.mean, NoConstant)
            ldiv!(sys.d, C.L, mean(model))
        end
    elseif method === :collapsed
        # Collapsing transformation
        # Z'×H⁻¹×Z
        tmp_pn= transpose(model.Λ) * prec(model)
        tmp= tmp_pn * model.Λ
        # (Z'×H⁻¹×Z)⁻¹
        # perform pseudo inverse as Λ can contain zero columns
        pseudo= pinv(tmp)
        # pivoted Cholesky decomposition
        C= cholesky!(Hermitian(pseudo), Val(true), check=false)
        # transformation
        ip= invperm(C.p)
        U= C.U[:,ip]
        A= U * tmp_pn

        # Store values
        mul!(sys.y, A, model.y)
        sys.Z[1].= pinv(transpose(U))
        @inbounds for t in 2:T_len
            sys.Z[t].= sys.Z[1]
        end
        if !isa(model.mean, NoConstant)
            mul!(sys.d, A, mean(model))
        end
    else
        # Store values
        @inbounds for t in 1:T_len
            sys.Z[t].= model.Λ
            sys.H[t].= cov(model)
        end
        if !isa(model.mean, NoConstant)
            sys.d.= mean(model)
        end
    end

    return nothing
end

function store_params!(model::DynamicFactorModel, ψ::AbstractVector)
    # Get dims
    (n,r)= size(model.Λ)
    idx= 0  # index counter

    # Store values
    # Λ
    vec(model.Λ).= view(ψ,1:n*r)
    idx+= n*r

    # ϕ
    model.ϕ.diag.= view(ψ,idx+1:idx+r)
    idx+= r

    # mean model
    store_params!(model.mean, view(ψ,idx+1:idx+nparams(model.mean)))
    idx+= nparams(model.mean)

    # error model
    store_params!(model.error, view(ψ,idx+1:idx+nparams(model.error)))

    return nothing
end

# Initialization
function init!(model::DynamicFactorModel, init::NamedTuple, method::Symbol)
    # Model
    init_model!(model, init)

    # State space system
    sys= create_system(model, method)
    fix_system!(sys, model, method)
    init_system!(sys, model)

    return sys
end

"""
    init_ϕ!(ϕ, pc)

Initialize autoregressive parameters `ϕ` using OLS with principal components,
storing the result in `ϕ`.

#### Arguments
  - `pc::AbstractMatrix`: principal components

#### Returns
  - `ϕ::Diagonal`   : autoregressive parameters
"""
function init_ϕ!(ϕ::Diagonal, pc::AbstractMatrix)
    # Get dims
    (r,T)= size(pc)

    # OLS
    @inbounds @fastmath for i in 1:r
        y= view(pc,i,2:T)
        x= view(pc,i,1:T-1)
        ϕ.diag[i]= dot(x, y) * inv(sum(abs2, x))
    end

    return nothing
end

function init_model!(model::DynamicFactorModel, init::NamedTuple)
    # Principal component analysis
    M= fit(PCA, model.y, maxoutdim=model.r, pratio=1.)
    loadings= projection(M)
    pc= transform(M, model.y)

    # residuals
    resid(model).= model.y
    mul!(resid(model), loadings, pc, -1., 1.)

    # Λ
    model.Λ.= haskey(init, :Λ) ? init.Λ : loadings
    # ϕ
    haskey(init, :ϕ) ? model.ϕ.= init.ϕ : init_ϕ!(model.ϕ, pc)

    # mean model
    init_mean!(model.mean, model.y, init)

    # update residuals
    if !isa(model.mean, NoConstant)
        resid(model).-= mean(model)
    end

    # error model
    init_error!(model.error, init)

    return nothing
end

function create_system(model::DynamicFactorModel, method::Symbol)
    # get dims
    (n,T_len)= method === :collapsed ? (model.r, size(model.y,2)) : size(model.y)

    # types
    Te= eltype(model.y)
    Tv= Vector{Te}
    Tm= Matrix{Te}
    Td= Diagonal{Te}
    if method === :univariate || method === :collapsed
        TH= Td
    else
        TH= cov(model) isa Symmetric ? Symmetric{Te} : Td
    end

    # State space system
    if model.mean isa Exogeneous
        sys= LinearTimeVariant{Tm, Tm, Td, Tm, Tm, TH, Td, Tv, Tm}(n, model.r, T_len)
    else
        sys= LinearTimeInvariant{Tm, Tm, Td, Tv, Tv, TH, Td, Tv, Tm}(n, model.r, T_len)
    end

    return sys
end

function fix_system!(sys::LinearTimeInvariant, model::DynamicFactorModel, method::Symbol)
    # infer type
    T= eltype(model.y)

    if method === :multivariate || method === :woodbury
        sys.y.= model.y
    end
    if method === :univariate || method === :collapsed
        sys.H.diag.= one(T)
    end
    sys.Q.diag.= one(T)
    if model.mean isa NoConstant
        sys.d.= zero(T)
    end
    sys.c.= zero(T)
end

function fix_system!(sys::LinearTimeVariant, model::DynamicFactorModel, method::Symbol)
    #  get dims
    T_len= size(model.y,2)

    # infer type
    T= eltype(model.y)

    if method === :multivariate || method === :woodbury
        sys.y.= model.y
    end
    if method === :univariate || method === :collapsed
        @inbounds for t in 1:T_len
            sys.H[t].diag.= one(T)
        end
    end
    @inbounds for t in 1:T_len
        sys.Q[t].diag.= one(T)
    end
    if model.mean isa NoConstant
        sys.d.= zero(T)
    end
    sys.c.= zero(T)
end

function init_system!(sys::StateSpaceSystem, model::DynamicFactorModel)
    T= eltype(model.y)
    # a
    sys.a1.= zero(T)
    # P
    sys.P1.= zero(T)
    @inbounds @fastmath for i in 1:model.r
        sys.P1[i,i]= one(T)
    end

    return nothing
end

# Estimation
# Log-likelihood
function loglik(filter::KalmanFilter, model::DynamicFactorModel, method::Symbol)
    # Compute log-likelihood
    ll= loglik(filter)

    # Account for transformation
    if method === :univariate
        # get dims
        T_len= size(model.y,2)

        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Add Jacobian determinant term
        d,s= logabsdet(C.L)
        ll-= T_len * (d + log(s))
    elseif method === :collapsed
        # get dims
        (n,T_len)= size(model.y)

        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Add Jacobian determinant term
        ll-= .5 * T_len * logdet(C)

        # Z'×H⁻¹×Z
        tmp_pn= transpose(model.Λ) * prec(model)
        tmp= tmp_pn * model.Λ
        # (Z'×H⁻¹×Z)⁻¹
        # perform pseudo inverse as Λ can contain zero columns
        pseudo= pinv(tmp)
        # Projection
        M= I - model.Λ * pseudo * tmp_pn

        # Add projected out term
        ll-= .5 * (n - model.r) * T_len * log(2*π)
        e= similar(model.y, n)
        @inbounds @fastmath for t in 1:T_len
            mul!(e, M, view(model.y,:,t))
            ll-= .5 * dot(e, prec(model), e)
        end
    end

    return ll
end

# EM
"""
    update_ϕ!(model, state)

Update autoregressive parameters `ϕ`, storing the result in `model`.

#### Arguments
  - `model::DynamicFactorModel` : state space model
  - `state:::EMOptimizerState`  : state variables
"""
function update_ϕ!(model::DynamicFactorModel, state::EMOptimizerState)
    @inbounds @fastmath for i in axes(model.ϕ,1)
        model.ϕ.diag[i]= state.V_01[i,i] * inv(state.V_1[i,i])
    end

    return nothing
end

"""
    update_Λ!(model, state, smoother, pen)

Update loading matrix parameters `Λ`, storing the result in `model`.

#### Arguments
  - `model::DynamicFactorModel` : state space model
  - `state:::EMOptimizerState`  : state variables
  - `smoother::Smoother`        : Kalman smoother output
  - `pen::Penalization`         : penalization parameters
"""
function update_Λ!(model::DynamicFactorModel, state::EMOptimizerState,
                    smoother::Smoother, pen::Penalization)
    # residuals
    resid(model).= model.y .- mean(model)
    # precision matrix
    Ω= prec(model)

    # Gram matrix
    α_1= view(smoother.α,:,1)
    state.V_0.+= view(smoother.V,:,:,1,1) .+ α_1 .* transpose(α_1)

    # linear coefficient b
    b= -inv(prod(size(model.y))) * vec(Ω * resid(model) * transpose(smoother.α))
    # I + A, with A quadratic coefficient
    tmp= inv(prod(size(model.y))) * kron(state.V_0, Ω)
    @inbounds @fastmath for i in axes(tmp,1)
        tmp[i,i]+= one(eltype(tmp))
    end
    # Cholesky decomposition
    C= cholesky!(Hermitian(tmp))

    # Proximal operators
    prox_g!(x::AbstractVector, λ::Real)= prox!(x, λ, pen)
    prox_f!(x::AbstractVector, λ::Real)= shrinkage!(x, λ, C, b)

    # Penalized estimation via admm
    vec(model.Λ).= admm!(vec(model.Λ), prox_f!, prox_g!)

    return nothing
end

function update_Λ!(model::DynamicFactorModel, state::EMOptimizerState,
                    smoother::Smoother, pen::NoPen)
    # residuals
    resid(model).= model.y .- mean(model)

    # Gram matrix
    α_1= view(smoother.α,:,1)
    state.V_0.+= view(smoother.V,:,:,1,1) .+ α_1 .* transpose(α_1)

    # Generalized Ridge regression
    mul!(model.Λ, resid(model), transpose(smoother.α))
    C= cholesky!(Hermitian(state.V_0))
    rdiv!(model.Λ, C)

    return nothing
end

"""
    resid_factors!(ε, model, smoother)

Calculate residuals based on factors only.

#### Arguments
  - `model::DynamicFactorModel` : state space model
  - `smoother::Smoother`        : Kalman smoother output

#### Returns
  - `ε::AbstractMatrix` : residuals
"""
function resid_factors!(ε::AbstractMatrix, model::DynamicFactorModel, smoother::Smoother)
    # Update residuals
    ε.= model.y
    mul!(ε, model.Λ, smoother.α, -1., 1.)

    return nothing
end

function update_model!(model::DynamicFactorModel, state::EMOptimizerState,
                        smoother::Smoother, pen::NamedTuple)
    # Closure of function
    init_resid!(ε::AbstractMatrix)= resid_factors!(ε, model, smoother)

    # Quadratic ridge component
    quad= model.Λ * state.V_sum * transpose(model.Λ)

    # Update hyper parameters
    update_ϕ!(model, state)
    update_error!(model.error, quad, pen.error)
    update_Λ!(model, state, smoother, pen.Λ)
    update_mean!(model.mean, init_resid!, model.error, pen.mean)

    # Update residuals
    ε= resid(model) # retrieve residual container
    ε.= model.y .- mean(model)
    mul!(ε, model.Λ, smoother.α, -1., 1.)

    return nothing
end

# Forecast
function reinstantiate(model::DynamicFactorModel, h::Integer)
    # number of time series
    n= size(model.y,1)

    # expand
    forecast_mean= forecast(model.mean, h)
    y_f= hcat(model.y, fill(NaN, n, h))

    return DynamicFactorModel(y_f, model.r, forecast_mean, model.Λ, model.ϕ, model.error)
end

function forecast(model::DynamicFactorModel, h::Integer)
    # reinstantiate
    forecast_model= reinstantiate(model, h)

    # State space system
    sys= create_system(forecast_model, :multivariate)
    fix_system!(sys, forecast_model, :multivariate)
    init_system!(sys, forecast_model)
    get_system!(sys, forecast_model, :multivariate)

    # forecasts
    f= forecast(sys, h)

    return f
end


text:
#
# Correctness Tests
#

using Base.Test
using MachineLearningMetrics

my_tests = [
    "classification.jl",
    "regression.jl"
  ]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\tPASSED: $(my_test)")
    catch e
        println("\tFAILED: $(my_test)")
    end
end


text:
import Bio.Align.BAM

# Count the number of mapped records.
function countmapped(reader)
    ret = 0
    record = BAM.Record()
    while !eof(reader)
        # in-place reading
        read!(reader, record)
        if BAM.ismapped(record)
            ret += 1
        end
    end
    return ret
end

println(open(countmapped, BAM.Reader, ARGS[1]))

## USAGE : C:\Program\Julia-1.2.0\bin\julia.exe .\test_prog.jl E:\bam_test\X001\1870157\g_1870157_chr8_48711871.bam


text:
"""
    HyperRectangle{N, T}

A `HyperRectangle` is a generalization of a rectangle into N-dimensions.
Formally it is the cartesian product of intervals, which is represented by the
`origin` and `widths` fields, whose indices correspond to each of the `N` axes.
"""
struct HyperRectangle{N,T} <: GeometryPrimitive{N,T}
    origin::Vec{N,T}
    widths::Vec{N,T}
end

##
# Constructors & typealiases

const Rect{N,T} = HyperRectangle{N,T}
const Rect2D{T} = Rect{2,T}
const Rect3D{T} = Rect{3,T}
const TRect{T} = Rect{N,T} where {N}

const FRect{N} = Rect{N,Float32}
const FRect2D = Rect2D{Float32}
const FRect3D = Rect3D{Float32}

const IRect{N} = HyperRectangle{N,Int}
const IRect2D = Rect2D{Int}
const IRect3D = Rect3D{Int}

Rect() = Rect{2,Float32}()

TRect{T}() where {T} = Rect{2,T}()

Rect{N}() where {N} = Rect{N,Float32}()

function Rect{N,T}() where {T,N}
    # empty constructor such that update will always include the first point
    return Rect{N,T}(Vec{N,T}(typemax(T)), Vec{N,T}(typemin(T)))
end

# conversion from other Rects
function Rect{N,T1}(a::Rect{N,T2}) where {N,T1,T2}
    return Rect(Vec{N,T1}(minimum(a)), Vec{N,T1}(widths(a)))
end

function Rect(v1::Vec{N,T1}, v2::Vec{N,T2}) where {N,T1,T2}
    T = promote_type(T1, T2)
    return Rect{N,T}(Vec{N,T}(v1), Vec{N,T}(v2))
end

function TRect{T}(v1::VecTypes{N,T1}, v2::VecTypes{N,T2}) where {N,T,T1,T2}
    return if T <: Integer
        Rect{N,T}(round.(T, v1), round.(T, v2))
    else
        return Rect{N,T}(Vec{N,T}(v1), Vec{N,T}(v2))
    end
end

function Rect{N}(v1::VecTypes{N,T1}, v2::VecTypes{N,T2}) where {N,T1,T2}
    T = promote_type(T1, T2)
    return Rect{N,T}(Vec{N,T}(v1), Vec{N,T}(v2))
end

"""
    Rect(vals::Number...)

```
Rect(vals::Number...)
```
Rect constructor for indidually specified intervals.
e.g. Rect(0,0,1,2) has origin == Vec(0,0) and
width == Vec(1,2)
"""
@generated function Rect(vals::Number...)
    # Generated so we get goodish codegen on each signature
    n = length(vals)
    @assert iseven(n)
    mid = div(n, 2)
    v1 = Expr(:call, :Vec)
    v2 = Expr(:call, :Vec)
    # TODO this can be inbounds
    append!(v1.args, [:(vals[$i]) for i in 1:mid])
    append!(v2.args, [:(vals[$i]) for i in (mid + 1):length(vals)])
    return Expr(:call, :Rect, v1, v2)
end

Rect3D(a::Vararg{Number,6}) = Rect3D(Vec{3}(a[1], a[2], a[3]), Vec{3}(a[4], a[5], a[6]))
Rect3D(args::Vararg{Number,4}) = Rect3D(Rect{2}(args...))
#=
From different args
=#
function (Rect)(args::Vararg{Number,4})
    args_prom = promote(args...)
    return Rect2D{typeof(args_prom[1])}(args_prom...)
end

function (Rect2D)(args::Vararg{Number,4})
    args_prom = promote(args...)
    return Rect2D{typeof(args_prom[1])}(args_prom...)
end

function (Rect{2,T})(args::Vararg{Number,4}) where {T}
    x, y, w, h = T <: Integer ? round.(T, args) : args
    return Rect2D{T}(Vec{2,T}(x, y), Vec{2,T}(w, h))
end

function TRect{T}(args::Vararg{Number,4}) where {T}
    x, y, w, h = T <: Integer ? round.(T, args) : args
    return Rect2D{T}(Vec{2,T}(x, y), Vec{2,T}(w, h))
end

function FRect3D(x::Rect2D{T}) where {T}
    return Rect{3,T}(Vec{3,T}(minimum(x)..., 0), Vec{3,T}(widths(x)..., 0.0))
end

function Rect2D{T}(a::Rect2D) where {T}
    return Rect2D{T}(minimum(a), widths(a))
end

function TRect{T}(a::Rect2D) where {T}
    return Rect2D{T}(minimum(a), widths(a))
end

function Rect{N,T}(a::GeometryPrimitive) where {N,T}
    return Rect{N,T}(Vec{N,T}(minimum(a)), Vec{N,T}(widths(a)))
end

function Rect2D(xy::VecTypes{2}, w::Number, h::Number)
    return Rect2D(xy..., w, h)
end

function Rect2D(x::Number, y::Number, wh::VecTypes{2})
    return Rect2D(x, y, wh...)
end

function TRect{T}(xy::VecTypes{2}, w::Number, h::Number) where {T}
    return Rect2D{T}(xy..., w, h)
end

function TRect{T}(x::Number, y::Number, wh::VecTypes{2}) where {T}
    return Rect2D{T}(x, y, wh...)
end

# TODO These are kinda silly
function Rect2D(xy::NamedTuple{(:x, :y)}, wh::NamedTuple{(:width, :height)})
    return Rect2D(xy.x, xy.y, wh.width, wh.height)
end

function FRect3D(x::Tuple{Tuple{<:Number,<:Number},Tuple{<:Number,<:Number}})
    return FRect3D(Vec3f0(x[1]..., 0), Vec3f0(x[2]..., 0))
end

function FRect3D(x::Tuple{Tuple{<:Number,<:Number,<:Number},
                          Tuple{<:Number,<:Number,<:Number}})
    return FRect3D(Vec3f0(x[1]...), Vec3f0(x[2]...))
end

origin(prim::Rect) = prim.origin
Base.maximum(prim::Rect) = origin(prim) + widths(prim)
Base.minimum(prim::Rect) = origin(prim)
Base.length(prim::Rect{N,T}) where {T,N} = N
widths(prim::Rect) = prim.widths

width(prim::Rect) = prim.widths[1]
height(prim::Rect) = prim.widths[2]

"""
    split(rectangle, axis, value)

Splits an Rect into two along an axis at a given location.
"""
split(b::Rect, axis, value::Integer) = _split(b, axis, value)
split(b::Rect, axis, value::Number) = _split(b, axis, value)
function _split(b::H, axis, value) where {H<:Rect}
    bmin = minimum(b)
    bmax = maximum(b)
    b1max = setindex(bmax, value, axis)
    b2min = setindex(bmin, value, axis)

    return H(bmin, b1max - bmin), H(b2min, bmax - b2min)
end

###
# Algebraic operations

"""
    *(m::Mat, h::Rect)

Transform a `Rect` using a matrix. Maintains axis-align properties
so a significantly larger Rect may be generated.
"""
function Base.:(*)(m::Mat{N1,N1,T1}, h::Rect{N2,T2}) where {N1,N2,T1,T2}

    # TypeVar constants
    T = promote_type(T1, T2)
    D = N1 - N2

    # get all points on the Rect
    d = decompose(Point, h)
    # make sure our points are sized for the tranform
    pts = (Vec{N1,T}[vcat(pt, ones(Vec{D,T})) for pt in d]...,)::NTuple{2^N2,Vec{N1,T}}

    vmin = Vec{N1,T}(typemax(T))
    vmax = Vec{N1,T}(typemin(T))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * pt
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    return Rect{N2,T}(vmin, vmax - vmin)
end

# equal dimension case
function Base.:(*)(m::Mat{N,N,T1}, h::Rect{N,T2}) where {N,T1,T2}

    # TypeVar constants
    T = promote_type(T1, T2)

    # get all points on the Rect
    pts = decompose(Point, h)

    # make sure our points are sized for the tranform
    vmin = Vec{N,T}(typemax(T))
    vmax = Vec{N,T}(typemin(T))

    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * Vec(pt)
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    return Rect{N,T}(vmin, vmax - vmin)
end

# fast path. TODO make other versions fast without code duplications like now
function Base.:(*)(m::Mat{4,4,T}, h::Rect{3,T}) where {T}
    # equal dimension case

    # get all points on the Rect
    pts = (Vec{4,T}(0.0, 0.0, 0.0, 1.0), Vec{4,T}(1.0, 0.0, 0.0, 1.0),
           Vec{4,T}(0.0, 1.0, 0.0, 1.0), Vec{4,T}(1.0, 1.0, 0.0, 1.0),
           Vec{4,T}(0.0, 0.0, 1.0, 1.0), Vec{4,T}(1.0, 0.0, 1.0, 1.0),
           Vec{4,T}(0.0, 1.0, 1.0, 1.0), Vec{4,T}(1.0, 1.0, 1.0, 1.0))

    # make sure our points are sized for the tranform
    vmin = Vec{4,T}(typemax(T))
    vmax = Vec{4,T}(typemin(T))
    o, w = origin(h), widths(h)
    _o = Vec{4,T}(o[1], o[2], o[3], T(0))
    _w = Vec{4,T}(w[1], w[2], w[3], T(1))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * (_o + (pt .* _w))
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    _vmin = Vec{3,T}(vmin[1], vmin[2], vmin[3])
    _vmax = Vec{3,T}(vmax[1], vmax[2], vmax[3])
    return Rect{3,T}(_vmin, _vmax - _vmin)
end

Base.:(-)(h::Rect{N,T}, move::Number) where {N,T} = h - Vec{N,T}(move)
Base.:(+)(h::Rect{N,T}, move::Number) where {N,T} = h + Vec{N,T}(move)

function Base.:(-)(h::Rect{N,T}, move::StaticVector{N}) where {N,T}
    return Rect{N,T}(minimum(h) .- move, widths(h))
end

function Base.:(+)(h::Rect{N,T}, move::StaticVector{N}) where {N,T}
    return Rect{N,T}(minimum(h) .+ move, widths(h))
end

function Base.:(*)(rect::Rect, scaling::Union{Number,StaticVector})
    return Rect(minimum(rect) .* scaling, widths(rect) .* scaling)
end

# Enables rectangular indexing into a matrix
function Base.to_indices(A::AbstractMatrix{T}, I::Tuple{Rect2D{IT}}) where {T,IT<:Integer}
    rect = I[1]
    mini = minimum(rect)
    wh = widths(rect)
    return ((mini[1] + 1):(mini[1] + wh[1]), (mini[2] + 1):(mini[2] + wh[2]))
end

function minmax(p::StaticVector, vmin, vmax)
    any(isnan, p) && return (vmin, vmax)
    return min.(p, vmin), max.(p, vmax)
end

# Annoying special case for view(Vector{Point}, Vector{Face})
function minmax(tup::Tuple, vmin, vmax)
    for p in tup
        any(isnan, p) && continue
        vmin = min.(p, vmin)
        vmax = max.(p, vmax)
    end
    return vmin, vmax
end

function positive_widths(rect::Rect{N,T}) where {N,T}
    mini, maxi = minimum(rect), maximum(rect)
    realmin = min.(mini, maxi)
    realmax = max.(mini, maxi)
    return Rect{N,T}(realmin, realmax .- realmin)
end

###
# set operations

"""
Perform a union between two Rects.
"""
function Base.union(h1::Rect{N}, h2::Rect{N}) where {N}
    m = min.(minimum(h1), minimum(h2))
    mm = max.(maximum(h1), maximum(h2))
    return Rect{N}(m, mm - m)
end

"""
    diff(h1::Rect, h2::Rect)

Perform a difference between two Rects.
"""
diff(h1::Rect, h2::Rect) = h1

"""
    intersect(h1::Rect, h2::Rect)

Perform a intersection between two Rects.
"""
function intersect(h1::Rect{N}, h2::Rect{N}) where {N}
    m = max.(minimum(h1), minimum(h2))
    mm = min.(maximum(h1), maximum(h2))
    return Rect{N}(m, mm - m)
end

function update(b::Rect{N,T}, v::Vec{N,T2}) where {N,T,T2}
    return update(b, Vec{N,T}(v))
end

function update(b::Rect{N,T}, v::Vec{N,T}) where {N,T}
    m = min.(minimum(b), v)
    maxi = maximum(b)
    mm = if any(isnan, maxi)
        v - m
    else
        max.(v, maxi) - m
    end
    return Rect{N,T}(m, mm)
end

# Min maximum distance functions between hrectangle and point for a given dimension
function min_dist_dim(rect::Rect{N,T}, p::Vec{N,T}, dim::Int) where {N,T}
    return max(zero(T), max(minimum(rect)[dim] - p[dim], p[dim] - maximum(rect)[dim]))
end

function max_dist_dim(rect::Rect{N,T}, p::Vec{N,T}, dim::Int) where {N,T}
    return max(maximum(rect)[dim] - p[dim], p[dim] - minimum(rect)[dim])
end

function min_dist_dim(rect1::Rect{N,T}, rect2::Rect{N,T}, dim::Int) where {N,T}
    return max(zero(T),
               max(minimum(rect1)[dim] - maximum(rect2)[dim],
                   minimum(rect2)[dim] - maximum(rect1)[dim]))
end

function max_dist_dim(rect1::Rect{N,T}, rect2::Rect{N,T}, dim::Int) where {N,T}
    return max(maximum(rect1)[dim] - minimum(rect2)[dim],
               maximum(rect2)[dim] - minimum(rect1)[dim])
end

# Total minimum maximum distance functions
function min_euclideansq(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    minimum_dist = T(0.0)
    for dim in 1:length(p)
        d = min_dist_dim(rect, p, dim)
        minimum_dist += d * d
    end
    return minimum_dist
end

function max_euclideansq(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    maximum_dist = T(0.0)
    for dim in 1:length(p)
        d = max_dist_dim(rect, p, dim)
        maximum_dist += d * d
    end
    return maximum_dist
end

function min_euclidean(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    return sqrt(min_euclideansq(rect, p))
end

function max_euclidean(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    return sqrt(max_euclideansq(rect, p))
end

# Functions that return both minimum and maximum for convenience
function minmax_dist_dim(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}},
                         dim::Int) where {N,T}
    minimum_d = min_dist_dim(rect, p, dim)
    maximum_d = max_dist_dim(rect, p, dim)
    return minimum_d, maximum_d
end

function minmax_euclideansq(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    minimum_dist = min_euclideansq(rect, p)
    maximum_dist = max_euclideansq(rect, p)
    return minimum_dist, maximum_dist
end

function minmax_euclidean(rect::Rect{N,T}, p::Union{Vec{N,T},Rect{N,T}}) where {N,T}
    minimumsq, maximumsq = minmax_euclideansq(rect, p)
    return sqrt(minimumsq), sqrt(maximumsq)
end

# http://en.wikipedia.org/wiki/Allen%27s_interval_algebra
function before(b1::Rect{N}, b2::Rect{N}) where {N}
    for i in 1:N
        maximum(b1)[i] < minimum(b2)[i] || return false
    end
    return true
end

meets(b1::Rect{N}, b2::Rect{N}) where {N} = maximum(b1) == minimum(b2)

function overlaps(b1::Rect{N}, b2::Rect{N}) where {N}
    for i in 1:N
        maximum(b2)[i] > maximum(b1)[i] > minimum(b2)[i] &&
            minimum(b1)[i] < minimum(b2)[i] || return false
    end
    return true
end

function starts(b1::Rect{N}, b2::Rect{N}) where {N}
    return if minimum(b1) == minimum(b2)
        for i in 1:N
            maximum(b1)[i] < maximum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

function during(b1::Rect{N}, b2::Rect{N}) where {N}
    for i in 1:N
        maximum(b1)[i] < maximum(b2)[i] && minimum(b1)[i] > minimum(b2)[i] || return false
    end
    return true
end

function finishes(b1::Rect{N}, b2::Rect{N}) where {N}
    return if maximum(b1) == maximum(b2)
        for i in 1:N
            minimum(b1)[i] > minimum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

###
# Containment

"""
    in(b1::Rect, b2::Rect)

Check if Rect `b1` is contained in `b2`. This does not use
strict inequality, so Rects may share faces and this will still
return true.
"""
function Base.in(b1::Rect{N}, b2::Rect{N}) where {N}
    for i in 1:N
        maximum(b1)[i] <= maximum(b2)[i] && minimum(b1)[i] >= minimum(b2)[i] || return false
    end
    return true
end

"""
    in(pt::VecTypes, b1::Rect{N, T})

Check if a point is contained in a Rect. This will return true if
the point is on a face of the Rect.
"""
function Base.in(pt::VecTypes, b1::Rect{N,T}) where {T,N}
    for i in 1:N
        pt[i] <= maximum(b1)[i] && pt[i] >= minimum(b1)[i] || return false
    end
    return true
end

#
# Equality
#
Base.:(==)(b1::Rect, b2::Rect) = minimum(b1) == minimum(b2) && widths(b1) == widths(b2)

Base.isequal(b1::Rect, b2::Rect) = b1 == b2

centered(R::Type{Rect{N,T}}) where {N,T} = R(Vec{N,T}(-0.5), Vec{N,T}(1))
centered(R::Type{Rect{N}}) where {N} = R(Vec{N,Float32}(-0.5), Vec{N,Float32}(1))
centered(R::Type{Rect}) where {N} = R(Vec{2,Float32}(-0.5), Vec{2,Float32}(1))

##
# Rect2D decomposition

function faces(rect::Rect2D, nvertices=(2, 2))
    w, h = nvertices
    idx = LinearIndices(nvertices)
    quad(i, j) = QuadFace{Int}(idx[i, j], idx[i + 1, j], idx[i + 1, j + 1], idx[i, j + 1])
    return ivec((quad(i, j) for i in 1:(w - 1), j in 1:(h - 1)))
end

function coordinates(rect::Rect2D, nvertices=(2, 2))
    mini, maxi = extrema(rect)
    xrange, yrange = LinRange.(mini, maxi, nvertices)
    return ivec(((x, y) for x in xrange, y in yrange))
end

function texturecoordinates(rect::Rect2D, nvertices=(2, 2))
    xrange, yrange = LinRange.((0, 1), (1, 0), nvertices)
    return ivec(((x, y) for x in xrange, y in yrange))
end

function normals(rect::Rect2D, nvertices=(2, 2))
    return Iterators.repeated((0, 0, 1), prod(nvertices))
end

##
# Rect3D decomposition
function coordinates(rect::Rect3D)
    # TODO use n
    w = widths(rect)
    o = origin(rect)
    points = Point{3,Int}[(0, 0, 0), (0, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 0), (1, 0, 0),
                          (1, 0, 1), (0, 0, 1), (0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 0, 0),
                          (1, 1, 1), (0, 1, 1), (0, 0, 1), (1, 0, 1), (1, 1, 1), (1, 0, 1),
                          (1, 0, 0), (1, 1, 0), (1, 1, 1), (1, 1, 0), (0, 1, 0), (0, 1, 1)]
    return ((x .* w .+ o) for x in points)
end

function texturecoordinates(rect::Rect3D)
    return coordinates(Rect3D(0, 0, 0, 1, 1, 1))
end

function faces(rect::Rect3D)
    return QuadFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8), (9, 10, 11, 12), (13, 14, 15, 16),
                         (17, 18, 19, 20), (21, 22, 23, 24),]
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    KLIEPWeighter(tdata, vars=nothing;
                  bwidth=1.0, nbasis=100)

Kullback-Leibler importance estimation procedure based on empirical
target distribution of variables `vars` in spatial data `tdata`. The
procedure depends on the kernel bandwidth `bwidth` and the number of
basis functions `nbasis`.
"""
struct KLIEPWeighter{DΩ<:AbstractData} <: AbstractWeighter
  tdata::DΩ
  vars::Vector{Symbol}
  bwidth::Float64
  nbasis::Int
end

function KLIEPWeighter(tdata::DΩ, vars=nothing;
                       bwidth=1.0, nbasis=100) where {DΩ<:AbstractData}
  validvars = collect(keys(variables(tdata)))
  wvars = vars ≠ nothing ? vars : validvars
  @assert wvars ⊆ validvars "invalid variables ($wvars) for spatial data"
  @assert nbasis ≤ npoints(tdata) "invalid number of basis functions"
  KLIEPWeighter{DΩ}(tdata, wvars, bwidth, min(nbasis, 100))
end

function weight(sdata::AbstractData, weighter::KLIEPWeighter)
  # retrieve method parameters
  tdata = weighter.tdata
  vars  = weighter.vars
  b     = weighter.nbasis
  σ     = weighter.bwidth

  @assert vars ⊆ keys(variables(sdata)) "invalid variables for spatial data"

  # Gaussian kernel
  kern(x, y) = exp(-norm(x - y)^2 / (2σ^2))

  # number of points in source and target
  nₛ = npoints(sdata)
  nₜ = npoints(tdata)

  # source and target features
  xₛ = i -> sdata[i,vars]
  xₜ = j -> tdata[j,vars]

  # basis for kernel approximation
  basis = sample(1:nₜ, b, replace=false)

  # constants for objective
  Φ = Matrix{Float64}(undef, nₜ, b)
  for l in 1:b
    xₗ = xₜ(basis[l])
    for k in 1:nₜ
      xₖ = xₜ(k)
      Φ[k,l] = kern(xₖ, xₗ)
    end
  end

  # constants for equality constraints
  A = Matrix{Float64}(undef, 1, b)
  for l in 1:b
    xₗ = xₜ(basis[l])
    A[l] = sum(kern(xₛ(k), xₗ) for k in 1:nₛ)
  end
  lc = uc = [nₛ]

  # constants for inequality constraints
  lx = fill(0.0, b)
  ux = fill(Inf, b)

  # objective
  f(α) = -sum(log, Φ*α)
  function ∇f!(g, α)
    p = Φ*α
    for l in 1:b
      g[l] = -sum(Φ[j,l] / p[j] for j in 1:nₛ)
    end
  end
  function ∇²f!(h, α)
    p = Φ*α
    for k in 1:b, l in 1:b
      h[k,l] = sum(view(Φ,:,k) .* view(Φ,:,l) ./ p)
    end
  end

  # equality constraint
  c!(c, α)    = c  .= A*α
  J!(J, α)    = J  .= A
  H!(H, α, λ) = H .+= 0.0

  # initial guess
  αₒ = fill(nₛ/sum(A), b)

  # optimization problem
  objective   = TwiceDifferentiable(f, ∇f!, ∇²f!, αₒ)
  constraints = TwiceDifferentiableConstraints(c!, J!, H!, lx, ux, lc, uc)
  initguess   = αₒ

  # solve problem with interior-point primal-dual Newton
  solution = optimize(objective, constraints, initguess, IPNewton())

  # optimal weights
  α = solution.minimizer
  weights = map(1:nₛ) do i
    sum(α[l] * kern(xₛ(i), xₜ(c)) for (l, c) in enumerate(basis))
  end

  SpatialWeights(domain(sdata), weights)
end


text:
#=##############################################################################
# DESCRIPTION
    Methods for generation of conic sections
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Sep 2019
  * License   : MIT License
=###############################################################################

"""
    `line_intersect(x1, y1, x2, y2, x3, y3, x4, y4)`

    Returns the intersect between two lines, with the first line
passing through the points `(x1, y1), (x2, y2)`, and the second line passing
through the points `(x3, y3), (x4, y4)`.
"""
function line_intersect(x1, y1, x2, y2, x3, y3, x4, y4)
    Px = (
            (x1*y2-y1*x2)*(x3-x4) - (x1-x2)*(x3*y4 - y3*x4)
         ) / (
            (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
         )
    Py = (
            (x1*y2-y1*x2)*(y3-y4) - (y1-y2)*(x3*y4-y3*x4)
         ) / (
            (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
         )

    return Px, Py
end

"""
    `generate_conic_fun(Ax::Real, Ay::Real, Bx::Real, By::Real, Cx::Real, Cy::Real, Sx::Real, Sy::Real)`

     Returns a conic curve function with control point A, B, and C, and shoulder point S.
        See Raymer's Aircraft Design, lofting chapter (p. 128).
"""
function generate_conic_fun(Ax::Real, Ay::Real, Bx::Real, By::Real,
                            Cx::Real, Cy::Real, Sx::Real, Sy::Real)

    """
        Receives a value between 0 and 1, and returns the
    corresponding point along the conic curve, with 0 being A
    and 1 B.
    """
    function conic_fun(s)
        if s<0 || s>1
            error("Invalid conic parameter $s. Value between 0 and 1 expected.")
        end

        # D: Point in between A and B
        Dx = Ax + s*(Bx-Ax)
        Dy = Ay + s*(By-Ay)

        # E: Intersect between AS and CD
        Ex, Ey = line_intersect(Ax, Ay, Sx, Sy, Cx, Cy, Dx, Dy)

        # F: Intersect between BS and CD
        Fx, Fy = line_intersect(Bx, By, Sx, Sy, Cx, Cy, Dx, Dy)

        # P: Intersect between AF and BE
        Px, Py = line_intersect(Ax, Ay, Fx, Fy, Bx, By, Ex, Ey)

        return Px, Py
    end

    return conic_fun
end

"""
    `generate_conic_fun(Ax::Real, Ay::Real, Bx::Real, By::Real, Cx::Real, Cy::Real, rho::Real)`

    Returns a conic curve function with control point A, B, and C, and shape parameter rho.
        See Raymer's Aircraft Design, lofting chapter (p. 132).

    Hyperbola: rho > 0.5
    Parabola:  rho = 0.5
    Ellipse:   rho < 0.5
    Circle:    rho = 0.4142 and AC = BC
"""
function generate_conic_fun(Ax::Real, Ay::Real, Bx::Real, By::Real,
                            Cx::Real, Cy::Real, rho::Real)

    # D: Midpoint between A and B
    Dx = Ax + 0.5*(Bx-Ax)
    Dy = Ay + 0.5*(By-Ay)

    # S: Shoulder point in between D and C
    Sx = Dx + rho*(Cx-Dx)
    Sy = Dy + rho*(Cy-Dy)

    return generate_conic_fun(Ax, Ay, Bx, By, Cx, Cy, Sx, Sy)
end

"""
    `conic_cross_section(Ps::Array{Array{T1, 1},1}, CPs::Array{Array{T2, 1},1},
rhos::Array{T3, 1}, ss::Array{Array{T4, 1},1})`

    Receives a collection of points `Ps` along an open contour, and
stretching control points `CPs`, shape parameters `rhos`, and probing
parameters `ss` associated to every section of the contour, and it returns
a compound-conic cross section (it is a closed loop, so the number of
sections is equal to the number of points)
"""
function conic_cross_section(Ps::Array{Array{T1, 1},1},
                             CPs::Array{Array{T2, 1},1},
                             rhos::Array{T3, 1},
                             ss::Array{Array{T4, 1},1}
                            ) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real}
    nPs = size(Ps, 1)

    if size(CPs, 1)!=nPs
        error("Invalid CPs. Expected $nPs, got $(size(CPs, 1)).")
    elseif size(rhos, 1)!=nPs
        error("Invalid rhos. Expected $nPs, got $(size(rhos, 1)).")
    elseif size(ss, 1)!=nPs
        error("Invalid ss. Expected $nPs, got $(size(ss, 1)).")
    end

    points = []

    for i in 1:nPs
        A = Ps[i]
        B = Ps[i%nPs + 1]

        conic_fun = generate_conic_fun(Ps[i]..., Ps[i%nPs + 1]...,
                                                CPs[i]..., rhos[i])
        for s in ss[i]
            push!(points, [conic_fun(s)...])
        end
    end

    return points
end


text:
abstract type Endpoint end

path(e::Union{Type{T},T}) where T <: Endpoint =  Symbol(split(repr(T), '.')[end])
Base.hasproperty(::Type{T}, x::Symbol) where T <: Endpoint = hasproperty(T(), x)

function Base.show(io::IO, ::MIME"text/plain", e::T) where T <: Endpoint
    println(io, path(e))
    for (nm, prop) in applied_filters(e)
        println(io, "  - $nm: $prop")
    end
end

function applied_filters(ep::Endpoint)::Dict{Symbol,Any}
    out = Dict{Symbol,Any}()
    for nm in propertynames(ep)
        prop = getproperty(ep, nm)
        ismissing(prop) && continue
        out[nm] = prop
    end
    out
end

## Metaprogramming stuff
function get_endpoint_params(fn::Symbol, sw::Dict)
    thekey = "\$ref"
    prefix = "#/parameters/"
    params = sw["paths"]["/$(fn)"]["get"]["parameters"]
    specs = Dict[]

    for param in params
        length(param) != 1 && continue
        haskey(param, thekey) || continue

        pname = param[thekey]
        startswith(pname, prefix) || continue

        param_spec = sw["parameters"][replace(pname, prefix => "")]
        param_spec["in"] === "query" && push!(specs, param_spec)
    end
    specs
end

function make_endpoint(fn::Symbol, sw::Dict)
    specs = get_endpoint_params(fn, sw)
    fields = Expr[]
    for spec in specs
        param_sym = Symbol(spec["name"])
        push!(fields, Expr(:(=), param_sym, :missing))
        param_sym == :location && push!(fields, Expr(:(=), :state, :missing))
    end
    out = :(
        @with_kw struct $(fn) <: Endpoint
        end
    )
    out.args[3].args[3] = Expr(:block, fields...)
    out
end

function make_endpoint_docs(fn::Symbol, sw::Dict, )
    desc = get(sw["paths"][string("/", fn)]["get"], "description", missing)
    docs = IOBuffer()
    println(docs, "CovidCountyData endpoint $fn\n")
    if !ismissing(desc)
        println(docs, desc)
        println(docs, "")
    end

    filters = [x["name"] for x in get_endpoint_params(fn, sw)]
    print(docs, "Available filters are: ")
    join(docs, filters, ", ", " and ")
    println(docs, "\n\nSet any filter as keyword argument to $(fn)() function")
    String(take!(docs))
end


text:
# GENERATING A LIST OF PARTICLES

newParticle( ::Val{true} , y0, x0, z0, ys, xs, zs ) = ( x0 + rand()*xs, y0 + rand()*ys, z0 + rand()*zs );
newParticle( ::Val{false}, y0, x0, z0, ys, xs, zs ) = ( y0 + rand()*ys, x0 + rand()*xs, z0 + rand()*zs );

function generateParticles( p::SP; pivparams=setPIVParameters()::PIVParameters )
    if p.dens > 0
        return generateDensityParticles( p.dens, p.w, p.h, p.d, pivparams.interSize, p.mode );
    else
        return generateTotalNumParticles( p.n, p.w, p.h, p.d, p.mode )
    end
end

function generateTotalNumParticles( n::I, w::I, h::I, d::I, mode::Bool )

    particles = [ newParticle( Val(mode), 0, 0, 0, h, w, d ) for x in 1:n ];
    return sort!( particles, by = x -> x[3] )
end

function generateDensityParticles( dens::I, w::I, h::I, d::I, ISize::III, mode::Bool )

    particles  = Array{ Tuple{Float32,Float32,Float32}, 1}(undef,0);
    Y0, X0, Z0 = min.( ISize, ( h, w, d ) );
    thicknessZ = ( d == 1 ) ? 1 : ISize[3];

    for z_max in Z0:ISize[3]:d, x_max in X0:ISize[2]:w, y_max in Y0:ISize[1]:h

		for p in 1:dens
			push!( particles, newParticle( Val(mode), y_max-Y0, x_max-X0, z_max-Z0, ISize[1], ISize[2], thicknessZ ) );
		end
    end

    return sort!( particles, by = x -> x[3] );
end

# RENDERING THE PARTICLES TO IMAGES/VOLUMES

function computeIntensityGaussian( params::SP, x::I, y::I, z::I, x0::T, y0::T, z0::T ) where {T<:Real}

    return computeIntensityGaussian( params.i0, params.dt, params.th, x, y, z , x0, y0, z0 )
end

function computeIntensityGaussian( i0::F, dt::F, th::F, x::I, y::I, z::I, x0::T, y0::T, z0::T ) where {T<:Real}

    return i0*exp( -8*( (z - z0)^2/th^2 + ((x - x0)^2 + (y - y0)^2)/dt^2 ) )
end

#=
    Instead of calculating the Gaussian intensity of each particle over the whole image, it is
    more efficient to limit the area over which we compute the gaussian for each particle. This
	method calculates the radius ( in pixels ) at which the Gaussian reaches a certain intensity,
	named "error". The formula is obtained by considering a 1D gaussian, and the distance "d" at
	which the gaussian equals "error".

    I( d ) = i0 * e^( ( - d² ) / ( 1/8 * dt²) ) = error
                       ln(i0) - d² / ( dt²/8 )  = ln(error)
                                              d = sqrt( dt²( ln(i0) - ln(error) )/8 )
=#

function calculateGaussianRadius( params::SP )
	return calculateGaussianRadius( params.err, params.i0, params.dt, params.rad )
end

function calculateGaussianRadius( error::F, i0::F, dt::F, rad::I )
    return ( rad > 0 ) ? rad : round( Int32, sqrt( dt*dt*( log(i0) - log(error) )/8 ) )
end

function calculateGaussianZRadius( params::SP )
	return calculateGaussianZRadius( params.err, params.i0, params.th, params.d )
end

function calculateGaussianZRadius( error::F, i0::F, th::F, depth::I )
    return ( th == Inf ) ? depth : round( Int32, sqrt( th*th*( log(i0) - log(error) )/8 ) )
end

function renderParticles( particles::Array{Tuple{T,T,T},1}, params::SP ) where {T<:Real}
    if params.d == 1
        return renderSlice( particles, params );
    else
        return renderVolume( particles, params );
    end
end

renderSlice( particles::A{Tuple{T,T,T},1}, p::SP ) where {T<:Real} = renderSlice( particles, p.w, p.h, p.d, p.i0, p.dt, p.th, p.z, p.err, p.noise, p.rad, p.mode )

function renderSlice( particles::A{Tuple{T,T,T},1}, w::I, h::I, d::I, i0::F, dt::F, th::F, z::I,
                      err::F, noise::I, rad::I, mode::Bool ) where {T<:Real}

    img = zeros( Float64, ( h, w ) );
    rz  = calculateGaussianZRadius( err, i0, th,   d );
    r   =  calculateGaussianRadius( err, i0, dt, rad );

    for p in particles
        if abs( p[3] - z ) <= rz

            xmin, xmax = max( 1, floor( Int, p[1] - r ) ), min( w, ceil( Int, p[1] + r ) );
            ymin, ymax = max( 1, floor( Int, p[2] - r ) ), min( h, ceil( Int, p[2] + r ) );

            for x in xmin:xmax
                for y in ymin:ymax
                    @inbounds img[ y, x ] += computeIntensityGaussian( i0, dt, th, x, y, z, p[1], p[2], p[3] )
                end
            end
        end
    end

    # Allow for noise. Only Poisson noise is implemented by now.
    if ( noise > 0 )
        poi = Poisson( noise );
        for e in eachindex( img )
            img[e] += rand( poi, 1 )[1]
        end
    end

    return img
end

function renderVolume( particles::A{Tuple{T,T,T},1}, p::SP ) where {T<:Real}
	renderVolume( particles, p.w, p.h, p.d, p.i0, p.dt, p.th, p.err, p.rad )
end

function renderVolume( particles::A{Tuple{T,T,T},1}, w::I, h::I, d::I, i0::F, dt::F, th::F, err::F, rad::I ) where { T<:Real }

    out = zeros( Float64, h, w, d );
    r   = calculateGaussianRadius( err, i0, dt, rad );

    #println( "[h,w,d] = ", [h,w,d], " radius = ", r );

    for p in particles

        p_lo, p_hi = floor.( Int, p .- r ), floor.( Int, p .+ r );
        ymin, ymax = max( 1, p_lo[1] ), min( h, p_hi[1] );
        xmin, xmax = max( 1, p_lo[2] ), min( w, p_hi[2] );
        zmin, zmax = max( 1, p_lo[3] ), min( d, p_hi[3] );

        for z in zmin:zmax
            for x in xmin:xmax
                for y in ymin:ymax
                    @inbounds out[ y, x, z ] += computeIntensityGaussian( i0, dt, dt, y, x, z, p[1], p[2], p[3] )
                end
            end
        end
    end

    return out
end

"""
    normalize( img::Array{Float64,2}, range=(0,0)::Tuple{R, R} )
    ---
    Normalizes the intensities values between 0.0 and 255.0.

    Normalization is not automatically applied because it can change the intensity values of pixels when there is
    out-of-frame loss of the particles with the maximum intensities. Therefore, when normalizing the original and
    tranformed image, the minimum and maximum values of the original image should be used in the normalization
    of the transformed image. Example:

                     slice     = syntheticParticles.renderSlice( coords, synParams );
                     slice_t   = syntheticParticles.renderSlice( transformed_coords, synParams );

                     min, max  = extrema(slice)
                     slice      = syntheticParticles.normalize( slice );
                     slice_t    = syntheticParticles.normalize( slice_t , ( min, max ) );
"""
function normalize( img::AbstractArray{T,N} ) where {T<:Real,N}

    min, max = extrema( img );
    return normalize( img, ( min, max ) )
end

function normalize( img::AbstractArray{T,N}, range::Tuple{T,T} ) where {T<:Real,N}

    min, max = range;
    norm_img = zeros( eltype(img), size(img) );
    diff     = max - min;
    @simd for i in 1:length(norm_img)
       @inbounds norm_img[i] = (img[i] - min)/diff
    end
    return norm_img
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
#=
Created on 31/01/2021 17:57:19
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Illustration of performance types on a simple sum.
=#

using BenchmarkTools

function mysum(A)
    total = 0.0
    n, m = size(A)
    for i in 1:n
        for j in 1:m
            total += A[i,j]
        end
    end
    return total
end

function mysum_rowwise(A)
    total = 0.0
    n, m = size(A)
    for j in 1:m
        for i in 1:n
            total += A[i,j]
        end
    end
    return total
end

function mysum_inbounds(A)
    total = 0.0
    n, m = size(A)
    for i in 1:n
        for j in 1:m
            @inbounds total += A[i,j]
        end
    end
    return total
end

function mysum_simd(A)
    total = 0.0
    n, m = size(A)
    @simd for i in 1:n
        @simd for j in 1:m
            total += A[i,j]
        end
    end
    return total
end

A = randn(10000, 10000)

@btime sum($A);  # reference

@btime mysum($A);  # standard implementation

@btime mysum_rowwise($A);  # changing the order of crossing the matrix

@btime mysum_inbounds($A);  # not performing checks

@btime mysum_simd($A);  # making loops independent


text:
module SupportMod

export helper_func1

"""
  This is a helper function

  This is a new paragraph

  Inputs:

    x: some value
    y: an Array
"""
function helper_func1(x, y::Array)
  return x + y
end

end  # end module


text:
function simulateimpedance_noiseless(circuitfunc,parameters,frequencies)
    return [circuitfunc(parameters,fr) for fr in frequencies]
end

function simulateimpedance_noiseless(circuit,frequencies)
    circuitfunc = tree_to_function(get_circuit_tree(circuit))
    return [circuitfunc(circuit.parameters,fr) for fr in frequencies]
end

function simulateimpedance(circuitfunc,parameters,frequencies,noise_ratio=0.01)
    impedances = [circuitfunc(parameters,fr) for fr in frequencies]
    R = real(impedances)
    I = imag(impedances)
    Z = abs.(impedances)
    return [R[i]+rand(Normal(0,Z[i]*noise_ratio))+I[i]*im+rand(Normal(0,Z[i]*noise_ratio))im for i in 1:length(impedances)]
end

function simulateimpedance(circuit::Circuit,frequencies,noise_ratio=0.01)
    circuitfunc = tree_to_function(get_circuit_tree(circuit))
    impedances = [circuitfunc(circuit.parameters,fr) for fr in frequencies]
    R = real(impedances)
    I = imag(impedances)
    Z = abs.(impedances)
    return [R[i]+rand(Normal(0,Z[i]*noise_ratio))+I[i]*im+rand(Normal(0,Z[i]*noise_ratio))im for i in 1:length(impedances)]
end


text:
#eliminates missing values
function cleanpanel(panel::DataFrame, ξs::Vector{Characteristic})
  essentialfields =  [ξ.Fξ for ξ ∈ ξs]

  pushmaybe(s::Symbol) = s ∈ essentialfields ? essentialfields : push!(essentialfields, s)
  pushmaybe(::Nothing) = essentialfields
  for ξ ∈ ξs
    if ξ.ws.specname === PARAM[:caraname]
      pushmaybe(PARAM[:iivol])
    elseif ξ.ws.specname === PARAM[:sortedname]
      pushmaybe(PARAM[:sorteddefaultconditioning])
    elseif ξ.ws.specname === PARAM[:doublesortedname]
      pushmaybe(PARAM[:doublesorteddefaultconditioning])
      pushmaybe(PARAM[:doublesorteddefaultsecondary])
    elseif ξ.ws.specname === PARAM[:zname]
      continue
    elseif ξ.ws.specname === PARAM[:fixedname]
      continue
    else
      @assert false
    end
  end

  for f ∈ essentialfields
    panel[!, f] .= finiteormissing.(panel[!, f])
  end

  panel = panel[completecases(panel, essentialfields), :]

  return panel
end

function basiccharacteristics!(panel::DataFrame;
  strategies=error("strategies is required"),
  strategyprefix::Symbol = error("strategyprefix is required"),
  weightspec::Symbol=error("weightspec is required"),
  preconditionertype::Symbol=error("preconditioner is required"))

  #form the weight spec
  local preconditioner
  if preconditionertype === :quantiledemean
    preconditioner = demeanξ ∘ quantiles
  elseif preconditionertype === :quantile
    preconditioner = quantiles
  elseif preconditionertype === :identity
    preconditioner = (x)->x
  else
    error("unrecognized preconditionertype $preconditionertype")
  end

  local ws::WeightSpec
  if weightspec === :z
    ws = zweightspec(; preconditioner)
  elseif weightspec === :cara
    ws = caraweightspec(; preconditioner)
  elseif weightspec === :sorted
    ws = sortedweightspec(; preconditioner)
  elseif weightspec === :sorted2
    ws = doublesortedweightspec(; preconditioner)
  elseif weightspec === :fixed
    ws = fixedweightspec(; preconditioner)
  else
    error("unrecognized weightspec $weightspec")
  end

  Fstrategies = Symbol.(strategyprefix, strategies)

  ξs = (s->Characteristic(ws, s)).(Fstrategies)

  return ξs
end
basiccharacteristics!(panel::DataFrame, ::Val{:momentum};
  strategies=PARAM[:momentumstrategiesused],
  strategyprefix::Symbol = PARAM[:momentumstrategyprefix],
  preconditionertype::Symbol = PARAM[:momentumpreconditionertype],
  weightspec::Symbol=PARAM[:momentumweightspec]) = basiccharacteristics!(panel;
    strategies, strategyprefix, preconditionertype, weightspec)

basiccharacteristics!(panel::DataFrame ,::Val{:placebo};
  strategies=PARAM[:placebostrategiesused],
  strategyprefix::Symbol = Symbol(),
  preconditionertype::Symbol = PARAM[:placebopreconditionertype],
  weightspec::Symbol=PARAM[:placeboweightspec]) = basiccharacteristics!(panel;
    strategies, strategyprefix, preconditionertype, weightspec)

basiccharacteristics!(panel::DataFrame, ::Val{:fundamental};
  strategies=PARAM[:fundamentalstrategiesused],
  strategyprefix::Symbol = PARAM[:fundamentalstrategyprefix],
  preconditionertype::Symbol = PARAM[:fundamentalpreconditionertype],
  weightspec::Symbol=PARAM[:fundamentalweightspec]) = basiccharacteristics!(panel;
    strategies, strategyprefix, preconditionertype, weightspec)

basiccharacteristics!(panel::DataFrame, ::Val{:fixed};
  strategies=PARAM[:fixedstrategiesused],
  strategyprefix::Symbol = PARAM[:fixedstrategyprefix],
  preconditionertype::Symbol = :identity,
  weightspec::Symbol=:fixed) = basiccharacteristics!(panel;
    strategies, strategyprefix, preconditionertype, weightspec)

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

#forms a measure from a combination of different characteristic categories
#this is the main function for forming the weights and measure spec,
#though we can make others as needed with the measuretype parameter
function measure!(panel::DataFrame, ::Val{:combinedbasic};
  strategycategories::Vector{Symbol} = PARAM[:combinedstrategycategories])

  ξs = reduce(vcat, (C->basiccharacteristics!(panel, Val(C))).(strategycategories))

  panel = cleanpanel(panel, ξs)

  #assign weights
  (ξ->assignweights!(panel, ξ)).(ξs)
  #error("Got through it! make sure all the testing code is reverted")

  #form the measure
  ms::MeasureSpec = MeasureSpec(ξs...)

  return (panel, ms)
end

#a simple version for grabbing the MeasureSpec
function MeasureSpec(::Val{:combinedbasic};
  strategycategories::Vector{Symbol} = PARAM[:combinedstrategycategories])

  ξs = reduce(vcat, (C->basiccharacteristics!(panel, Val(C))).(strategycategories))
  return MeasureSpec(ξs...)
end

##########################################
#Simulated data from excel
#replaces actual data with simulated data
function measure!(::Any, ::Val{:testsimulated};
  Fxis = [:xi1,:xi2],
  testpath::String = PARAM[:testpath],
  testdataname::String = PARAM[:testdataname],
  testdateformat=PARAM[:wrdsdateformat])

  #read in the test dataset
  panel = CSV.File("$testpath\\$testdataname.csv") |> DataFrame
  panel.date = (s->Dates.Date("$s", testdateformat)).(panel.date)
  sort!(panel, [:permno,:date])

  panel.Wdvol = panel.dvol #do not winsorize since we are trying to match excel
  #differencewithin2!(panel, [:Wdvol], :permno, date=:date)

  wsz::WeightSpec = zweightspec(preconditioner=demeanξ)
  #ξs = (Fxi->Characteristic(wsz, Symbol(:x_,Fxi))).(Fxis)
  ξs = (Fxi->Characteristic(wsz, Fxi)).(Fxis)
  (ξ->assignweights!(panel, ξ)).(ξs)
  ms::MeasureSpec = MeasureSpec(ξs...)

  return (panel, ms)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    ThinAlgo

A thinning algorithm.
"""
abstract type ThinAlgo end

"""
    struct GuoAlgo <: ThinAlgo end

The Guo algorithm evaluates three conditions in order to determine which pixels of
the image should be removed.

The three conditions are explained in the page 361 of **Guo, Z., & Hall, R. W. (1989).
Parallel thinning with two-subiteration algorithms. Communications of the ACM, 32(3), 359-373.**
"""
struct GuoAlgo <: ThinAlgo end

"""
    thinning(img::AbstractArray{Bool}; algo::ThinAlgo=GuoAlgo())

Applies a binary blob thinning operation to achieve a skeletization of the input image.

See also: [`GuoAlgo`](@ref)
"""
function thinning(img::AbstractArray{Bool}; algo::ThinAlgo=GuoAlgo())
    # dispatch appropriate implementation
    thinning_impl(img, algo)
end

function thinning_impl(img::AbstractArray{Bool,2}, algo::GuoAlgo)
    # pad input image
    prev = falses(size(img).+2)
    prev[2:end-1,2:end-1] = img

    # preallocate memory for the mask
    mask = falses(size(prev))

    # perform single iteration
    it = 1; curr = copy(prev)
    guo_iteration!(mask, curr, isodd(it))
    curr[mask] .= false

    # loop if necessary
    while prev != curr
        it += 1; prev .= curr
        guo_iteration!(mask, curr, isodd(it))
        curr[mask] .= false
    end

    # unpad result
    curr[2:end-1,2:end-1]
end

# update mask with Guo iteration on padded image
function guo_iteration!(mask::AbstractArray{Bool,2}, img::AbstractArray{Bool,2}, odd_iteration::Bool)
    # start clean
    mask .= false

    # loop over pixels and update mask
    for j=2:size(img,2)-1, i=2:size(img,1)-1
        img[i,j] || continue
        p1 = img[i-1,j-1]
        p2 = img[i-1,j]
        p3 = img[i-1,j+1]
        p4 = img[i,j+1]
        p5 = img[i+1,j+1]
        p6 = img[i+1,j]
        p7 = img[i+1,j-1]
        p8 = img[i,j-1]
        C = (!p2 && (p3 || p4)) + (!p4 && (p5 || p6)) + (!p6 && (p7 || p8)) + (!p8 && (p1 || p2))
        N1 = (p1 || p2) + (p3 || p4) + (p5 || p6) + (p7 || p8)
        N2 = (p2 || p3) + (p4 || p5) + (p6 || p7) + (p8 || p1)
        N = min(N1, N2)
        if odd_iteration
            O = (p2 || p3 || !p5) && p4
        else
            O = (p6 || p7 || !p1) && p8
        end
        if C == 1 && (2 ≤ N ≤ 3) && !O
            mask[i,j] = true
        end
    end
end


text:
typealias FunctionOrNothing Union(Function, Nothing)
typealias F64OrVectorF64 Union(Float64, Vector{Float64})

abstract MCModel

### Sampler types hold the components that fully specify a Monte Carlo sampler
### The fields of sampler types are user-inputed

abstract MCSampler

abstract MHSampler <: MCSampler # Family of samplers based on Metropolis-Hastings
abstract HMCSampler <: MCSampler # Family of Hamiltonian Monte Carlo samplers
abstract LMCSampler <: MCSampler # Family of Langevin Monte Carlo samplers

### Monte Carlo sample types hold a single Monte Carlo sample
### A typical MCSample includes at least the vector of parameter values and the respective log-target value

abstract DerivativeOrder
type NullOrder <: DerivativeOrder end
type FirstOrder <: DerivativeOrder end
type SecondOrder <: DerivativeOrder end
type ThirdOrder <: DerivativeOrder end

abstract MCSample{O<:DerivativeOrder}

### Runner types indicate what type of simulation will be run (ex serial or sequential Monte Carlo)
### Their fields fully specify the simulation details (ex total number or number of burn-in iterations)

abstract ParallelismLevel
type Serial <: ParallelismLevel end
type Parallel <: ParallelismLevel end

abstract MCRunner{P<:ParallelismLevel}

typealias SerialMCRunner MCRunner{Serial}
typealias ParrallelMCRunner MCRunner{Parallel}

### Tuner types

abstract MCTuner

### Tune types hold the temporary output of the sampler that is used for tuning the sampler

abstract MCTune

### Heap types hold the temporary components used by a Monte Carlo sampler during its run
### This means that heap types represent the internal state ("local variables") of a Monte Carlo sampler

abstract MCHeap{S<:MCSample}

### Monte Carlo Jobs (ex plain jobs, jobs using tasks or MPI jobs)

abstract MCJob


text:
using LinearAlgebra
function sigmoid(a)
    return 1.0 / (1 + exp(-a))
end
function fit(x, t; alpha = 0.01, tau_max = 1000)
    function CEE(w, x, t)
        grad = zeros(size(w))
        for i in 1:length(t)
            ti = (t[i] > 0) ? 1 : 0
            h = sigmoid(dot(w, x[i, :]))
            grad += ((h - ti) * x[i, :])'
        end
        return grad / length(t)
    end
    c = size(t)[2]
    w = []
    x = hcat(ones(size(x)[1], 1), x)
    for i in 1:c
        if size(x)[1] != size(t[:, i])[1]#Processing when the matrix is organized by dependent variable
            x = x'
        end
        w0 = ones(1, length(x[1, :]))
        for j in 1:tau_max
            grad = CEE(w0, x, t[:, i])
            w0 -= alpha * grad
        end
        push!(w, w0)
    end
    return w
end

function forecast(x, w)
    x = hcat(ones(size(x)[1], 1), x)
    p = nothing
    for i in 1 : size(x)[1]
        tp = [sigmoid(dot(w[j, :][1], x[i, :])) for j in 1:size(w)[1]]
        tp = tp'
        (p == nothing) ? p = tp : p = vcat(p, tp)
    end
    return p
end

function predict(x, w)
    x = hcat(ones(size(x)[1], 1), x)
    p = []
    for i in 1 : size(x)[1]
        tp = [sigmoid(dot(w[j, :][1], x[i, :])) for j in 1:size(w)[1]]
        tp = argmax(tp)
        push!(p, tp)
    end
    return p
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ThingType <: AbstractThingType
    label::Label
    is_root::Bool
end

Instanciate_Types = Union{<:AbstractAttribute,
                            <:AbstractRelation,
                            <:AbstractEntity,
                            <:AbstractThingType,
                            <:AbstractEntityType,
                            <:AbstractAttributeType,
                            <:AbstractRelationType}

# Porting note: Java client calls into RequestBuilder but it really
# has nothing to to with requests... I think it's probably better
# migrating the function here.
function proto(t::AbstractThingType)
    Proto._Type(
        label=t.label.name,
        encoding=encoding(t)
    )
    # return ThingTypeRequestBuilder.proto_thing_type(t.label, encoding(t))
end

# a convinient function to prevent code to work around a argument nothing
proto(t::Nothing) = t

# Contract: all subtypes of AbstractThingType should have these two fields
is_root(t::AbstractThingType) = t.is_root
label(t::AbstractThingType) = t.label

# ------------------------------------------------------------------------
# Remote functions
# ------------------------------------------------------------------------

"""
    set_supertype(r::RemoteConcept{<: AbstractThingType,<: AbstractCoreTransaction},
        super_type::AbstractThingType)

Here we have the chance to set a given ThingType as the supertype of an other Thingtype.
Mixing of diffent root types eg. AttributeType with EntityType is not allowed.
"""
function set_supertype(r::RemoteConcept{<:AbstractThingType},
                        super_type::AbstractThingType)

    req = ThingTypeRequestBuilder.set_supertype_req(r.concept.label, proto(super_type))
    execute(r.transaction, req)
end

"""
    get_supertype(r::RemoteConcept{<:AbstractType})

Returns one step upwards in the chain of types and returns the direct supertype of a
given type.
"""
function get_supertype(r::RemoteConcept{<:AbstractType})
    req = TypeRequestBuilder.get_supertype_req(r.concept.label)
    res = execute(r.transaction, req)
    typ = res.type_res.type_get_supertype_res._type
    return instantiate(typ)
end

"""
    get_supertypes(r::RemoteConcept{<:AbstractType})

Here all supertypes in the chain upwards from the given type will be returned.
"""
function get_supertypes(r::RemoteConcept{<:AbstractType})

    req = TypeRequestBuilder.get_supertypes_req(r.concept.label)
    res = execute(r.transaction, req)
    if res !== nothing
        typs = res.type_res_part.type_get_supertypes_res_part.types
        return instantiate.(typs)
    else
        return Instanciate_Types[]
    end
end

"""
    get_subtypes(r::RemoteConcept{<:AbstractType})

get_subtypes returns all subtypes for a given type. A specialzation of this function
is made for AttributeType.
"""
function get_subtypes(r::RemoteConcept{<:AbstractType})
    req = TypeRequestBuilder.get_subtypes_req(r.concept.label)
    res = execute(r.transaction, req)
    if res !== nothing
        typs = res.type_res_part.type_get_subtypes_res_part.types
        return instantiate.(typs)
    else
        return []
    end
end

"""
    get_instances(r::RemoteConcept{<:AbstractThingType})

Here you can get all instances for a given Type. The has to be given as a RemoteConcept
"""
function get_instances(r::RemoteConcept{<:AbstractThingType})

    req = ThingTypeRequestBuilder.get_instances_req(r.concept.label)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_instances_res_part.things for r in res)))
end

"""
    set_abstract(r::RemoteConcept{<:AbstractThingType})

In the chain of fine grained functions in the concept section here we can determine
a given ThingType is abstract.
"""
function set_abstract(r::RemoteConcept{<:AbstractThingType})
    req = ThingTypeRequestBuilder.set_abstract_req(r.concept.label)
    execute(r.transaction, req)
end

"""
    unset_abstract(r::RemoteConcept{<:AbstractThingType})

With this function it is possible to change an ThingType from abstact to normal behavior
which means that the given ThingType can be instanceated.
"""
function unset_abstract(r::RemoteConcept{<:AbstractThingType})

    req = ThingTypeRequestBuilder.unset_abstract_req(r.concept.label)
    execute(r.transaction, req)
end

function is_abstract(r::RemoteConcept{<:AbstractType})

    req = TypeRequestBuilder.is_abstract_req(r.concept.label)
    res = execute(r.transaction, req)
    return res.type_res.type_is_abstract_res._abstract
end

"""
    set_plays(
        r::RemoteConcept{<:AbstractThingType},
        role_type::AbstractRoleType,
        overridden_role_type::Optional{AbstractRoleType}=nothing
    )

With set_play we are able to set a roletype to an ThingType and it is possible to
set a new RoleType instead of the old. Be cautious, this is only allowed for not instantiated
types.
"""
function set_plays(
    r::RemoteConcept{<:AbstractThingType},
    role_type::AbstractRoleType,
    overridden_role_type::Optional{AbstractRoleType}=nothing
)
    req = ThingTypeRequestBuilder.set_plays_req(
        r.concept.label,
        proto(role_type),
        proto(overridden_role_type)
    )
    execute(r.transaction, req)
end

function unset_plays(r::RemoteConcept{<:AbstractThingType}, role_type::AbstractRoleType)
    unset_req = ThingTypeRequestBuilder.unset_plays_req(r.concept.label, proto(role_type))
    execute(r.transaction, unset_req)
end

"""
    set_owns(
        r::RemoteConcept{<:AbstractType},
        attribute_type::AbstractType,
        is_key::Bool=false,
        overriden_type::Optional{AbstractType}=nothing
    )

With set_owns it is possible to assign an attribute to a type. If it is needed it is
possible to set the attribute as unique. This means for the type where the attribute is set
as key only one instance of the type can have a specific attribute:
E.g. person entity has a unique email as an attribute. So it is not possible to have two
entities with the same email address.
"""
function set_owns(
    r::RemoteConcept{<:AbstractType},
    attribute_type::AbstractType,
    is_key::Bool=false,
    overriden_type::Optional{AbstractType}=nothing
)
    req = ThingTypeRequestBuilder.set_owns_req(
        r.concept.label,
        is_key,
        proto(attribute_type),
        proto(overriden_type)
    )
    execute(r.transaction, req)
end

function unset_owns(
    r::RemoteConcept{<: AbstractType},
    attribute_type::AbstractType)
    req = ThingTypeRequestBuilder.unset_owns_req(
        r.concept.label,
        proto(attribute_type)
    )
    execute(r.transaction, req)
end

"""
    get_owns(
        r::RemoteConcept{<:AbstractThingType},
        value_type::Optional{EnumType}=nothing,
        keys_only::Bool=false
    )

get_owns will give back all thingtypes which owns a given AttributeType
"""
function get_owns(
    r::RemoteConcept{<:AbstractThingType},
    value_type::Optional{EnumType}=nothing,
    keys_only::Bool=false
)
    req = ThingTypeRequestBuilder.get_owns_req(r.concept.label, value_type, keys_only)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_owns_res_part.attribute_types for r in res)))
end

"""
    get_plays(r::RemoteConcept{<: AbstractThingType})

Returns which roles are played by the given ThingType
"""
function get_plays(r::RemoteConcept{<:AbstractThingType})
    req = ThingTypeRequestBuilder.get_plays_req(r.concept.label)
    res = stream(r.transaction, req)
    return instantiate.(collect(Iterators.flatten(
        r.type_res_part.thing_type_get_plays_res_part.roles for r in res)))
end

"""
    set_label(r::RemoteConcept{<:AbstractType},
        new_label_name::AbstractString)

With this function we are able to set the label for a given Type. This gives us the
chance to rename a given type. But be prepared, this is only allowed if the type is
not instantiated by inserted data.
"""
function set_label(r::RemoteConcept{<:AbstractType},
    new_label_name::AbstractString)

    set_label_req = TypeRequestBuilder.set_label_req(r.concept.label, new_label_name)
    execute(r.transaction, set_label_req)

    return nothing
end

"""
    delete(r::RemoteConcept{<:AbstractThingType})
To delete a type in the database pack the type with the transaction to a RemoteConcept.
Be aware that a type can only be deleted if no Entity, Attribute or Relation is in the
database which is based on this type.
"""
function delete(r::RemoteConcept{<:AbstractThingType})
    del_req = TypeRequestBuilder.delete_req(r.concept.label)
    execute(r.transaction, del_req)
end


text:
invalid_imgs = [
    ("too_few_dimensions", MethodError, rand(127)),
    ("too_many_channels", AssertionError, rand(127, 257, 5)),
    ("too_many_dimensions", MethodError, rand(127, 257, 3, 1)),
]

invalid_palette_arguments = [
    ("palette_too_large", ArgumentError, IndirectArray(rand(1:257, 127, 256), rand(RGB{N0f8}, 257))),
    ("palette_index_too_few_dimensions", MethodError, IndirectArray(rand(1:4, 127), rand(RGB{N0f8}, 4))),
    ("palette_index_too_many_dimensions", ArgumentError, IndirectArray(rand(1:4, 127, 256, 1), rand(RGB{N0f8}, 4))),
    ("palette_eltype_too_wide", ArgumentError, IndirectArray(rand(1:4, 127, 256), rand(RGB{N0f16}, 4))),
    ("palette_non_RGB_colorant", ArgumentError, IndirectArray(rand(1:4, 127, 256), rand(Gray{N0f8}, 4))),
]

incompatible_background_arguments = [
    ("graya_img_rgb_background", ArgumentError, rand(127, 257, 2), "temp/PngSuite/bgbn4a08.png", RGB{N0f8}(1.0, 0.5, 0.0)),
    ("graya_img_index_background", ArgumentError, rand(127, 257, 2), "temp/PngSuite/bgbn4a08.png", 0x00),
    ("rgba_img_index_background", ArgumentError, rand(127, 257, 4), "temp/PngSuite/bgwn6a08.png", 0x00),
]

warn_background_arguments = [
    ("rgb_no_alpha_background", rand(127, 257, 3), "temp/PngSuite/basn2c08.png", RGB{N0f8}(0.0)),
    ("gray_no_alpha_background", rand(127, 257, 1), "temp/PngSuite/basn0g08.png", Gray{N0f8}(0.0)),
]

@testset "Invalid inputs" begin
    @testset "Dimensionality" begin
        for (case, exception, image) in invalid_imgs
            @testset "$(case) throws" begin
                @test_throws exception PNGFiles.save(joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"),image)
                @test_throws exception open(io->PNGFiles.save(io,image), joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"), "w")
            end
        end
    end

    @testset "Paletted" begin
        for (case, exception, image) in invalid_palette_arguments
            @testset "$(case) throws" begin
                @test_throws exception PNGFiles.save(joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"),image)
                @test_throws exception open(io->PNGFiles.save(io,image), joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"), "w")
            end
        end
    end

    @testset "Background" begin
        for (case, exception, image, load_path, bg) in incompatible_background_arguments
            @testset "$(case) throws save" begin
                @test_throws exception PNGFiles.save(joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"), image, background=bg)
                @test_throws exception open(io->PNGFiles.save(io, image, background=bg), joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"), "w")
            end
            @testset "$(case) throws load" begin
                @test_throws exception PNGFiles.load(load_path, background=bg)
            end
        end
        @testset "palette_img_index_background_not_expandded" begin
            @test_throws ArgumentError PNGFiles.load("test_images/paletted_implicit_transparency.png", background=0x00)
        end
    end

    warn_msg = "Background color for non-transparent image: PNGFiles._png_color_16_struct(0x00, 0x0000, 0x0000, 0x0000, 0x0000)"
    for (case, image, load_path, bg) in warn_background_arguments
        @testset "$(case) throws save" begin
            @test PNGFiles.save(joinpath(PNG_TEST_PATH, "test_img_err_$(case).png"), image, background=bg) === nothing
        end
        @testset "$(case) throws load" begin
            @test_logs (:warn, warn_msg) PNGFiles.load(load_path, background=bg)
        end
    end
    @testset "non_transparent_image_with_background_warns" begin
        @test_throws ArgumentError PNGFiles.load("test_images/nontransparent_image_with_background.png", background=0x00)
    end
end


text:
export Bernoulli

"""
Description:

    Bernoulli factor node

    out ∈ {0, 1}
    p ∈ [0, 1]

    f(out, p) = Ber(out|p) = p^out (1 - p)^{1 - out}

Interfaces:

    1. out
    2. p

Construction:

    Bernoulli(id=:some_id)
"""
mutable struct Bernoulli <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Bernoulli(out, p; id=generateId(Bernoulli))
        @ensureVariables(out, p)
        self = new(id, Array{Interface}(undef, 2), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:p] = self.interfaces[2] = associate!(Interface(self), p)

        return self
    end
end

slug(::Type{Bernoulli}) = "Ber"

format(dist::ProbabilityDistribution{Univariate, Bernoulli}) = "$(slug(Bernoulli))(p=$(format(dist.params[:p])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{Bernoulli}; p=0.5) = ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>p))
ProbabilityDistribution(::Type{Bernoulli}; p=0.5) = ProbabilityDistribution{Univariate, Bernoulli}(Dict(:p=>p))

dims(dist::ProbabilityDistribution{Univariate, Bernoulli}) = 1

vague(::Type{Bernoulli}) = ProbabilityDistribution(Univariate, Bernoulli, p=0.5)

isProper(dist::ProbabilityDistribution{Univariate, Bernoulli}) = (0 <= dist.params[:p] <= 1)

unsafeMean(dist::ProbabilityDistribution{Univariate, Bernoulli}) = dist.params[:p]

unsafeMode(dist::ProbabilityDistribution{Univariate, Bernoulli}) = round(dist.params[:p])

unsafeMeanVector(dist::ProbabilityDistribution{Univariate, Bernoulli}) = [dist.params[:p], 1 - dist.params[:p]]

unsafeVar(dist::ProbabilityDistribution{Univariate, Bernoulli}) = dist.params[:p]*(1 - dist.params[:p])

logPdf(dist::ProbabilityDistribution{Univariate, Bernoulli}, x) = x*log(dist.params[:p]) + (1.0-x)*log(1.0-dist.params[:p])

sample(dist::ProbabilityDistribution{Univariate, Bernoulli}) = 1.0*(rand() < dist.params[:p])

function prod!( x::ProbabilityDistribution{Univariate, Bernoulli},
                y::ProbabilityDistribution{Univariate, Bernoulli},
                z::ProbabilityDistribution{Univariate, Bernoulli}=ProbabilityDistribution(Univariate, Bernoulli, p=0.5))

    norm = x.params[:p] * y.params[:p] + (1 - x.params[:p]) * (1 - y.params[:p])
    (norm > 0) || error("Product of $(x) and $(y) cannot be normalized")
    z.params[:p] = (x.params[:p] * y.params[:p]) / norm

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate, Bernoulli},
    y::ProbabilityDistribution{Univariate, PointMass},
    z::ProbabilityDistribution{Univariate, PointMass}=ProbabilityDistribution(Univariate, PointMass, m=NaN))

    z.params[:m] = y.params[:m]

    return z
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, Bernoulli})
    p = clamp(dist.params[:p], tiny, 1.0 - tiny)

    -(1.0 - p)*log(1.0 - p) -
    p*log(p)
end

# Average energy functional
function averageEnergy(::Type{Bernoulli}, marg_out::ProbabilityDistribution{Univariate}, marg_p::ProbabilityDistribution{Univariate})
    -unsafeMean(marg_out)*unsafeLogMean(marg_p) -
    (1.0 - unsafeMean(marg_out))*unsafeMirroredLogMean(marg_p)
end