text:
struct AStarMCSPHEntry{D <: Number}
    v_idx::Int64
    parent_idx::Int64
    gvalue::D
    fvalue::D
    weights::Vector{D}
    parent_handle::Int64    # refers to the heap handle in CLOSED LIST that tracks this path
end

Base.isless(e1::AStarMCSPHEntry, e2::AStarMCSPHEntry) = (e1.fvalue, e1.gvalue) < (e2.fvalue, e2.gvalue) # Lexicographic

@with_kw mutable struct AStarMCSPStates{D <: Number}
    open_list::MutableBinaryMinHeap{AStarMCSPHEntry{D}}    = MutableBinaryMinHeap{AStarMCSPHEntry{D}}()
    closed_list::MutableBinaryMinHeap{AStarMCSPHEntry{D}}  = MutableBinaryMinHeap{AStarMCSPHEntry{D}}()
    open_list_hmap::Dict{Int64,Set{Int64}}                 = Dict{Int64,Set{Int64}}()
    closed_list_hmap::Dict{Int64,Set{Int64}}               = Dict{Int64,Set{Int64}}()
end

# A path is dominated if cost and weight are both worse than an existing entry
function is_dominated(state::AStarMCSPStates{D}, s::AStarMCSPHEntry{D}) where {D <: Number}

    if haskey(state.open_list_hmap, s.v_idx)
        for ol_idx in state.open_list_hmap[s.v_idx]
            ols_entry = state.open_list.nodes[state.open_list.node_map[ol_idx]].value
            if ols_entry.gvalue <= s.gvalue && ols_entry.weights <= s.weights
                return true
            end
        end
    end

    if haskey(state.closed_list_hmap, s.v_idx)
        for cl_idx in state.closed_list_hmap[s.v_idx]
            cls_entry = state.closed_list.nodes[state.closed_list.node_map[cl_idx]].value
            if cls_entry.gvalue <= s.gvalue && cls_entry.weights <= s.weights
                return true
            end
        end
    end

    return false
end

# A new entry is eligble if the weight-to-go based on this transition
# plus the heuristic weight-to-go is within the weight constraint
function is_eligible(graph::AbstractGraph{V},
                     entry::AStarMCSPHEntry{D},
                     nbr::Int64,
                     weight_functions::Vector{F1},
                     weight_heuristics::Vector{F2},
                     weight_constraints::Vector{D}) where {V, D <: Number, F1 <: Function, F2 <: Function}

    n_weights = length(weight_functions)
    new_weight_vector = zeros(D, n_weights)

    for i = 1:n_weights
        new_constr_wt = entry.weights[i] +  weight_functions[i](graph.vertices[entry.v_idx], graph.vertices[nbr])

        if new_constr_wt + weight_heuristics[i](graph.vertices[nbr])  > weight_constraints[i]
            return false, new_weight_vector
        else
            new_weight_vector[i] = new_constr_wt
        end
    end

    return true, new_weight_vector
end

function process_neighbors_implicit!(
    state::AStarMCSPStates{D},
    graph::AbstractGraph{V},
    edge_wt_fn::F1,
    neighbors::Vector{Int64},
    parent_entry::AStarMCSPHEntry{D},
    entry_cl_handle::Int64,    # For closed list
    visitor::AbstractDijkstraVisitor,
    heuristic::F2,
    weight_functions::Vector{F3},
    weight_heuristics::Vector{F4},
    weight_constraints::Vector{D}
    ) where {V, D <: Number, F1 <: Function, F2 <: Function, F3 <: Function, F4 <: Function}

    for nbr_idx in neighbors

        is_elig, new_weight_vector =  is_eligible(graph, parent_entry, nbr_idx, weight_functions, weight_heuristics, weight_constraints)

        if is_elig
            # Create new entry
            new_gvalue = parent_entry.gvalue + edge_wt_fn(graph.vertices[parent_entry.v_idx], graph.vertices[nbr_idx])
            new_fvalue = new_gvalue + heuristic(graph.vertices[nbr_idx])
            new_entry = AStarMCSPHEntry(nbr_idx, parent_entry.v_idx, new_gvalue, new_fvalue,
                                        new_weight_vector, entry_cl_handle)

            if ~(is_dominated(state, new_entry))
                new_entry_handle = push!(state.open_list, new_entry)
                if ~(haskey(state.open_list_hmap, nbr_idx))
                    Graphs.discover_vertex!(visitor, graph.vertices[parent_entry.v_idx], graph.vertices[nbr_idx], new_gvalue)
                    state.open_list_hmap[nbr_idx] = Set{Int64}(new_entry_handle)
                else
                    push!(state.open_list_hmap[nbr_idx], new_entry_handle)
                end
            end

        end
    end
end

"""
Runs the A*-MCSP (A* with Multiple Constraints Shortest Path) algorithm on the given list-of-indices graph.
As with all A* versions (except a_star_spath), the graph is implicit and the graph visitor must generate the neighbors of the expanded
vertex just-in-time.

Arguments:
    - `graph::AbstractGraph{V} where {V}`
    - `edge_wt_fn::F1 where {F1 <: Function}` Maps (u,v) to the edge weight of the u -> v edge
    - `source::Int64`
    - `visitor::AbstractDijkstraVisitor` The graph visitor that implements include_vertex!, which is called when
                                         a vertex is expanded
    - `heuristic::F2 where {F2 <: Function}` Maps u to h(u); an admissible heuristic for the cost-to-go
    - `weight_functions::Vector{F3} where {F3 <: Function}` A list of constraint weight functions for edges
    - `weight_heuristics::Vector{F4} where {F4 <: Function}` A list of admissible heuristics for the weight constraint for edges
    - `weight_constraints::Vector{D} where {D <: Number}` A list of correspoding weight constraint values
    - `::Type{D} = Float64` The type of the edge weight value

Returns:
    - `::AStarMCSPStates` The result of the A* search
    - `::AStarMCSPHEntry` The heap entry for the target node (or the source if no path found)
"""
function a_star_implicit_constrained_shortest_path!(
    graph::AbstractGraph{V},
    edge_wt_fn::F1,
    source::Int64,
    visitor::AbstractDijkstraVisitor,
    heuristic::F2,
    weight_functions::Vector{F3},
    weight_heuristics::Vector{F4},
    weight_constraints::Vector{D},
    ::Type{D} = Float64) where {V, D <: Number, F1 <: Function, F2 <: Function, F3 <: Function, F4 <: Function}

    state = AStarMCSPStates{D}()

    n_weights = length(weight_constraints)

    source_entry = AStarMCSPHEntry(source, source, zero(D), heuristic(graph.vertices[source]), zeros(D, n_weights), 0)
    source_handle = push!(state.open_list, source_entry)
    state.open_list_hmap[source] = Set{Int64}(source_handle)

    while ~(isempty(state.open_list))

        # Remove from open_list and open_list_hmap
        entry, entry_handle = top_with_handle(state.open_list)
        pop!(state.open_list)
        delete!(state.open_list_hmap[entry.v_idx], entry_handle)

        nbrs = Vector{Int64}(undef, 0)

        if Graphs.include_vertex!(visitor, graph.vertices[entry.parent_idx],
                                  graph.vertices[entry.v_idx], entry.gvalue, nbrs) == false
            return state, entry
        end

        # Insert into closed list to get new handle
        cl_entry_handle = push!(state.closed_list, entry)
        if ~(haskey(state.closed_list_hmap, entry.v_idx))
            state.closed_list_hmap[entry.v_idx] = Set{Int64}(cl_entry_handle)
        else
            push!(state.closed_list_hmap[entry.v_idx], cl_entry_handle)
        end

        # Process nbrs
        process_neighbors_implicit!(state, graph, edge_wt_fn, nbrs,
                                    entry, cl_entry_handle, visitor,
                                    heuristic, weight_functions, weight_heuristics, weight_constraints)
        Graphs.close_vertex!(visitor, graph.vertices[entry.v_idx])
    end

    return state, source_entry
end

"""
Given the solution state for the A*MCSP Heap, walk back using the closed list to actually get the shortest feasible path.

Argument:
    - `state::AStarMCSPStates{D} where D <: Number` The solution state returned by the shortest path call.
    - `graph::AbstractGraph{V} where V` The graph (as a vector of indices)
    - `source_idx::Int64`
    - `target_entry::AStarMCSPHEntry{D} where D <: Number` The heap entry of the target node returned by the shortest path call

Returns:
    - `::Vector{Int64}` The list of indices on the shortest feasible path in order.
    - `::Vector{D} where D <: Number` The g-values (cost-to-go) of the shortest path indices.
    - `::Vector{Vector{D}} where D <: Number` The constraint weight vectors for the shortest path indices.
"""
function shortest_path_cost_weights(state::AStarMCSPStates{D}, graph::AbstractGraph{V},
                                    source::Int64, target_entry::AStarMCSPHEntry{D}) where {V, D <: Number}

    target_idx = target_entry.v_idx

    sp = [target_idx]
    weights = [target_entry.weights]
    costs = [target_entry.gvalue]

    # Walk back USING closed list heap
    curr_idx = target_idx
    curr_par_handle = target_entry.parent_handle

    while curr_idx != source

        cl_list_entry = state.closed_list.nodes[state.closed_list.node_map[curr_par_handle]].value

        curr_idx = cl_list_entry.v_idx
        curr_par_handle = cl_list_entry.parent_handle
        pushfirst!(sp, curr_idx)
        pushfirst!(weights, cl_list_entry.weights)
        pushfirst!(costs, cl_list_entry.gvalue)
    end

    return sp, costs, weights
end


text:
"""
    Apply Laplace method for a given decision matrix

"""
function laplace(decisionMatrix::DataFrame)::LaplaceResult

    n, p = size(decisionMatrix)

    probs = [1.0 / p for i in 1:p]

    expecteds = zeros(Float64, n)

    # m = convert(Array{Float64,2}, decisionMatrix)
    m = Matrix{Float64}(decisionMatrix)

    for i in 1:n
        expecteds[i] = (m[i,:] .* probs)  |> sum
    end

    bestIndex = sortperm(expecteds) |> last

    result = LaplaceResult(
        expecteds,
        bestIndex
    )

    return result
end

"""
    Apply maximin method for a given decision matrix

"""
function maximin(decisionMatrix::DataFrame)::MaximinResult

    n, p = size(decisionMatrix)

    # m = convert(Array{Float64,2}, decisionMatrix)
    m = Matrix{Float64}(decisionMatrix)

    rmins = rowmins(decisionMatrix)

    bestIndex = sortperm(rmins) |> last

    result = MaximinResult(
        rmins,
        bestIndex
    )

    return result
end

"""
    Apply maximax method for a given decision matrix

"""
function maximax(decisionMatrix::DataFrame)::MaximaxResult

    n, p = size(decisionMatrix)

    # m = convert(Array{Float64,2}, decisionMatrix)
    m = Matrix{Float64}(decisionMatrix)

    rmaxs = rowmaxs(decisionMatrix)

    bestIndex = sortperm(rmaxs) |> last

    result = MaximaxResult(
        rmaxs,
        bestIndex
    )

    return result
end

"""
    Apply minimax method for a given decision matrix

"""
function minimax(decisionMatrix::DataFrame)::MinimaxResult

    n, p = size(decisionMatrix)

    # m = convert(Array{Float64,2}, decisionMatrix)
    m = Matrix{Float64}(decisionMatrix)

    rmaxs = rowmaxs(decisionMatrix)

    bestIndex = sortperm(rmaxs) |> first

    result = MinimaxResult(
        rmaxs,
        bestIndex
    )

    return result
end

"""
    Apply minimin method for a given decision matrix

"""
function minimin(decisionMatrix::DataFrame)::MiniminResult

    n, p = size(decisionMatrix)

    # m = convert(Array{Float64,2}, decisionMatrix)
    m = Matrix{Float64}(decisionMatrix)

    rmins = rowmins(decisionMatrix)

    bestIndex = sortperm(rmins) |> first

    result = MiniminResult(
        rmins,
        bestIndex
    )

    return result
end

"""
    Apply Savage method for a given decision matrix

"""
function savage(decisionMatrix::DataFrame)::SavageResult

    n, p = size(decisionMatrix)

    cmaxs = colmaxs(decisionMatrix)

    newDecisionMatrix = similar(decisionMatrix)

    @inbounds for i in 1:p
        newDecisionMatrix[:, i] = cmaxs[i] .- decisionMatrix[:, i]
    end

    rmaxs = scores = rowmaxs(newDecisionMatrix)

    bestIndex = rmaxs |> sortperm |> first

    result = SavageResult(
        newDecisionMatrix,
        scores,
        bestIndex
    )

    return result
end

"""
    Apply Hurwicz method for a given decision matrix

"""
function hurwicz(decisionMatrix::DataFrame; alpha::Float64=0.5)::HurwiczResult

    n, p = size(decisionMatrix)

    rmaxs = rowmaxs(decisionMatrix)
    rmins = rowmins(decisionMatrix)

    scores = alpha .* rmaxs .+ (1.0 - alpha) .* rmins

    bestIndex = scores |> sortperm |>  last

    result = HurwiczResult(
        scores,
        bestIndex
    )

    return result

end

"""
    Apply MLE (Maximum Likelihood) method for a given decision matrix

"""
function mle(decisionMatrix::DataFrame, weights::Array{Float64,1})::MLEResult

    w = unitize(weights)

    weightedMatrix = w * decisionMatrix

    scores = rowsums(weightedMatrix)

    bestIndex = scores |> sortperm |> last

    result = MLEResult(
        scores,
        bestIndex
    )

    return result

end

"""
    Apply Expected Regret method for a given decision matrix

"""
function expectedregret(decisionMatrix::DataFrame, weights::Array{Float64,1})

    _, p = size(decisionMatrix)

    w = unitize(weights)

    cmaxs = colmaxs(decisionMatrix)

    regretmat = similar(decisionMatrix)

    for i in 1:p
        regretmat[:, i] = cmaxs[i] .- decisionMatrix[:, i]
    end

    weightedregret = w * regretmat

    scores = rowsums(weightedregret)

    bestIndex = scores |> sortperm |> first

    result = ExpectedRegretResult(
        scores,
        bestIndex
    )

    return result
end

ext: jl

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
module ACTransforms
export
    acwt,
    iacwt,
    acwpt,
    iacwpt
using ..ACWT
using LinearAlgebra, Wavelets

"""
	acwt(x, wt[, L=maxtransformlevels(x)])

Perform a forward ac wavelet transform of the array `x`.
This method works for the 2D case as well.
The wavelet type `wt` determines the transform type.
Refer to Wavelet.jl for a list of available methods.

# Examples
```julia
acwt(x, wavelet(WT.db4))
```

**See also:** `iacwt`
"""
function acwt end

"""
    acwpt(x, wt[, L=maxtransformlevels(x)])

Perform a ac wavelet packet transform of the array `x`.
The wavelet type `wt` determines the transform type.

# Examples
```julia
acwpt(x, wavelet(WT.db4))
```

**See also:** `iacwpt`
"""
function acwpt end

"""
	iacwt

The inverse of `acwt`

**See also:** `acwt`
"""
function iacwt end

"""
    iacwpt

The inverse of `acwpt`

**See also:** `acwpt`
"""

function iacwt!(x::AbstractArray{<:Number,2})
    n,m = size(x)
    for i = 2:m
        @inbounds x[:,1] = (x[:,1] + x[:,i])/sqrt(2)
    end
end

function iacwt(x::AbstractArray{<:Number,2})
    y = deepcopy(x)
    iacwt!(y)
    return y[:,1]
end

function acwt_step(v::AbstractVector{T}, j::Integer, h::Array{T,1}, g::Array{T,1}) where T <: Number
    N = length(v)
    L = length(h)
    v1 = zeros(T, N)
    w1 = zeros(T, N)

    @inbounds begin
        for i in 1:N
            t = i
            i = mod1(i + (L÷2+1) * 2^(j-1),N) # Need to shift by half the filter size because of periodicity assumption
            for n in 1:L
                t += 2^(j-1)
                t = mod1(t, N)
                w1[i] += h[n] * v[t]
                v1[i] += g[n] * v[t]
            end
        end
    end
    return v1, w1
end

function acwt(x::AbstractVector{<:Number}, wt::OrthoFilter, L::Integer=maxtransformlevels(x))

    @assert L <= maxtransformlevels(x) || throw(ArgumentError("Too many transform levels (length(x) < 2^L"))
    @assert L >= 1 || throw(ArgumentError("L must be >= 1"))

    # Setup
	n = length(x)
    Pmf, Qmf = ACWT.makereverseqmfpair(wt)
	wp = zeros(n,L+1)
	wp[:,1] = x

    for j in 1:L
        @inbounds wp[:,1], wp[:,L+2-j] = acwt_step(wp[:,1],j,Qmf,Pmf)
    end

	return wp
end

"""Computes the column-wise forward acwat coeficients for 2D signals."""
function hacwt(x::AbstractArray{<:Number,2},
                         wt::OrthoFilter,
                         L::Integer=maxtransformlevels(x[:,1]))
    nrow, ncol = size(x)
    W = Array{Float64,3}(undef,nrow,L+1,ncol)
    for i in 1:ncol
        @inbounds W[:,:,i] = acwt(x[:,i],wt,L)
    end
    return W
end

"""Computes the row-wise 1D acwt coeficients for 2D signals."""
function vacwt(x::AbstractArray{<:Number,2},
                      wt::OrthoFilter,
                      L::Integer=maxtransformlevels(x[1,:]))
    nrow, ncol = size(x)
    W = Array{Number,3}(undef,ncol,L+1,nrow)
    for i in 1:nrow
        W[:,:,i] = acwt(x[i,:],wt,L)
    end
    return W
end

function acwt(x::AbstractArray{<:Number,2}, wt::OrthoFilter,
              Lrow::Integer=maxtransformlevels(x[1,:]),
              Lcol::Integer=maxtransformlevels(x[:,1]))
    nrow, ncol = size(x)
    W3d = hacwt(x,wt,Lcol)
    W4d = Array{Number,4}(undef,Lcol+1,ncol,Lrow+1,nrow)
    for i in 1:Lcol+1
        @inbounds W4d[i,:,:,:] = vacwt(W3d[:,i,:],wt,Lrow)
    end
    W4d = permutedims(W4d, [4,2,3,1])
    return W4d
end

# Inverse Transforms
function iacwt(x::AbstractArray{<:Number,4})
    nrow, ncol, Lrow, Lcol = size(x)
    W4d = permutedims(x,[4,2,3,1])
    W3d = Array{Number,3}(undef, nrow, Lcol, ncol)
    for i in 1:Lcol
        for j in 1:nrow
            @inbounds W3d[j,i,:] = iacwt(W4d[i,:,:,j])
        end
    end
    y = Array{Number,2}(undef, nrow, ncol)
    for i in 1:ncol
        @inbounds y[:,i] = iacwt(W3d[:,:,i])
    end
    return y
end

function acwpt_step(W::AbstractArray{T,2}, i::Integer, d::Integer, Qmf::Vector{T}, Pmf::Vector{T}) where T <: Number
    n,m = size(W)
    if i<<1+1 <= m
      W[:,i<<1], W[:,i<<1+1] = acwt_step(W[:,i],d,Qmf,Pmf)
      acwpt_step(W, leftchild(i), d+1, Qmf, Pmf) # left
      acwpt_step(W, rightchild(i), d+1, Qmf, Pmf) # right
    end
  end

  function acwpt(x::AbstractVector{T}, wt::OrthoFilter, L::Integer=maxtransformlevels(x)) where T<:Number
    W = Array{Float64,2}(undef,length(x),2<<L-1)
    W[:,1] = x
    Pmf, Qmf = ACWT.makereverseqmfpair(wt)
    acwpt_step(W,1,1,Qmf,Pmf)
    return W
  end

function iacwpt(xw::AbstractArray{<:Number,2}, tree::BitVector, i::Integer=1)

    @assert i <= size(xw, 2)
    @assert isvalidtree(xw[:,1], tree)
    n₀ = length(tree)
    if i > n₀ || tree[i] == false       # leaf node
        return xw[:,i]
    end

    v₀ = iacwpt(xw,tree,leftchild(i))
    v₁ = iacwpt(xw,tree,rightchild(i))

    return (v₀ + v₁) / √2
end

end # End module

max_stars_repo_name: ShozenD/AutocorrelationShell

max_issues_repo_name: ShozenD/AutocorrelationShell

max_forks_repo_name: ShozenD/AutocorrelationShell


text:
#=
n2,i2 are light, n3,i3 are regular, n4,i4 are medium, n6,i6 are semi-bold
<script src="/use.edgefonts.net/source-sans-pro:n2,i2,n3,i3,n4,i4,n6,i6:all.js"></script>
=#


text:
using Polynomials: Poly, polyder, polyint

struct LegendreBasis{T, N, fT, pT, dT, iT} <: PolynomialBasis{T,N}
    factors::fT
    polys::pT
    derivs::dT
    ints::iT
end

function LegendreBasis(T,N)

    factors = Tuple{}()
    polys   = Tuple{}()
    derivs  = Tuple{}()
    ints    = Tuple{}()

    for p in 0:N-1
        tfac, tpol = legendre_polynomial(p)
        tder = polyder(tpol)
        tint = polyint(tpol)

        factors = (factors..., tfac)
        polys   = (polys..., tpol)
        derivs  = (derivs..., tder)
        ints    = (ints..., tint)
    end

    LegendreBasis{T, N, typeof(factors), typeof(polys), typeof(derivs), typeof(ints)}(factors, polys, derivs, ints)
end

nodes(b::LegendreBasis{T,N})  where {T,N} = Vector{T}([])

Base.hash(b::LegendreBasis{T,N}, h::UInt) where {T,N} = hash(T, hash(N, h))

Base.:(==)(b1::LegendreBasis{T1,N1}, b2::LegendreBasis{T2,N2}) where {T1,N1,T2,N2} = (T1 == T2 && N1 == N2)

Base.isequal(b1::LegendreBasis{T1,N1}, b2::LegendreBasis{T2,N2}) where {T1,N1,T2,N2} = (T1 == T2 && N1 == N2)

function legendre_polynomial(p, T=Float64)
    if p == 0
        (T(sqrt(1)), Poly([1]))
    elseif p == 1
        (T(sqrt(3)), Poly([-1, +2]))
    elseif p == 2
        (T(sqrt(5)), Poly([+1, -6, +6]))
    elseif p == 3
        (T(sqrt(7)), Poly([-1, +12, -30, +20]))
    elseif p == 4
        (T(sqrt(9)), Poly([+1, -20, +90, -140, +70]))
    elseif p == 5
        (T(sqrt(11)), Poly([-1, +30, -210, +560, -630, +252]))
    elseif p == 6
        (T(sqrt(13)), Poly([+1, -42, +420, -1680, +3150, -2772, +924]))
    elseif p == 7
        (T(sqrt(15)), Poly([-1, +56, -756, +4200, -11550, +16632, -12012, +3432]))
    elseif p == 8
        (T(sqrt(17)), Poly([+1, -72, +1260, -9240, +34650, -72072, +84084, -51480, +12870]))
    elseif p == 9
        (T(sqrt(19)), Poly([-1, +90, -1980, +18480, -90090, +252252, -420420, +411840, -218790, +48620]))
    elseif p == 10
        (T(sqrt(21)), Poly([+1, -110, +2970, -34320, +210210, -756756, +1681680, -2333760, +1969110, -923780, +184756]))
    else
        error("Legendre polynomials only implemented for 0 ≤ p ≤ 10.")
    end
end

function eval_basis(b::LegendreBasis{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.polys[j](x)
end

function deriv_basis(b::LegendreBasis{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.derivs[j](x)
end

function int_basis(b::LegendreBasis{T,N}, j::Int, x::T) where {T,N}
    return b.factors[j] * b.ints[j](x)
end


text:
module BigCombinatorics

using Combinatorics

export Fibonacci
export Factorial, DoubleFactorial, FallingFactorial, RisingFactorial
export Binomial, Catalan
export Derangements, MultiChoose, Multinomial
export Bell, Stirling1, Stirling2
export IntPartitions, IntPartitionsDistinct
export Euler, PowerSum

_master_table = Dict{Function,Dict}()

function _save(f::Function, x, val::BigInt)
    d = _master_table[f]
    d[x] = val
    nothing
end

function _has(f::Function, x)::Bool
    d = _master_table[f]
    return haskey(d,x)
end

function _get(f::Function, x)::BigInt
    d = _master_table[f]
    return d[x]
end

function _make(f::Function, T::Type)
    _master_table[f] = Dict{T,BigInt}()
    nothing
end

"""
`BigCombinatorics.cache_report()` reports how many
entries are saved for each function in the `BigCombinatorics`
module.
"""
function cache_report()
    total = 0
    for func in keys(_master_table)
        parse_func = split(string(func),".")
        func_name = last(parse_func)
        sz = length(_master_table[func])
        total += sz
        println("$sz\t$func_name")
    end
    println()
    println("$total\tTotal entries")
    nothing
end

"""
`BigCombinatorics.cache_clear(f)` clears the cached values for
the function `f`. Returns `false` if `f` values are not cached.

`BigCombinatorics.cache_clear()` clears all cached values.
"""
function cache_clear(f::Function)::Bool
    if !haskey(_master_table,f)
        return false
    end
    tab = _master_table[f]
    for k in keys(tab)
        delete!(tab,k)
    end
    return true
end

function cache_clear()
    for f in keys(_master_table)
        cache_clear(f)
    end
    true
end

"""
`Fibonacci(n)` returns the `n`-th Fibonacci number.
We begin with `Fibonacci(0)==0` and `Fibonacci(1)==1`.
"""
function Fibonacci(n::Integer)::BigInt
  if n<0
    throw(DomainError(n,"argument must be nonngative"))
  end
  if n==0
    return big(0)
  end
  if n==1
    return big(1)
  end
  if _has(Fibonacci,n)
      return _get(Fibonacci,n)
  end
  val = Fibonacci(n-1)+Fibonacci(n-2)
  _save(Fibonacci,n,val)
  return val
end

_make(Fibonacci,Integer)

"""
`Factorial(n)` returns `n!` for nonnegative integers `n`.
`Factorial(n,k)` returns `n!/k!` (to be consistent with Julia's
`factorial`.) Requires `0 <= k <= n`.

See also `FallingFactorial` and `RisingFactorial`.
"""
function Factorial(n::Integer)::BigInt
    if n < 0
        throw(DomainError(n,"argument must be nonngative"))
    end
    if n <= 1
        return big(1)
    end
    if _has(Factorial,n)
        return _get(Factorial,n)
    end
    val = big(n) * Factorial(n-1)
    _save(Factorial,n,val)
    return val
end

function Factorial(n::Integer,k::Integer)::BigInt
    if k>n
        throw(DomainError((n,k),"$k cannot exceed $n"))
    end
    return div(Factorial(n),Factorial(k))
end
_make(Factorial,Integer)

"""
`FallingFactorial(n,k)` returns `n*(n-1)*(n-2)*...*(n-k+1)`
(with a total of `k` factors). Requires `n,k >= 0`.
If `k>n` then `0` is returned.
"""
function FallingFactorial(n::Integer, k::Integer)::BigInt
  if n<0 || k<0
    throw(DomainError((n,k)," arguments must be nonnegative"))
  end
  if k>n
    return big(0)
  end
  return Factorial(n,n-k)
end

"""
`RisingFactorial(n,k)` returns `n*(n+1)*(n+2)*...*(n+k-1)`
(with a total of `k` factors). Requires `n,k >= 0`.
"""
function RisingFactorial(n::Integer,k::Integer)::BigInt
  if n<0 || k<0
    throw(DomainError((n,k),"arguments must be nonnegative"))
  end
  if k==0
    return big(1)
  end
  if n==0
    return big(0)
  end

  return FallingFactorial(n+k-1,k)
end

"""
`DoubleFactorial(n)` returns `n!!`, i.e.,
`n*(n-2)*...` with `(-1)!! == 0!! == 1!! == 1`.
"""
function DoubleFactorial(n::Integer)::BigInt
  if n<-1
    throw(DomainError(n,"argument must be at least -1"))
  end
  if n<2
    return big(1)
  end
  if _has(DoubleFactorial,n)
      return _get(DoubleFactorial,n)
  end
  val = n * DoubleFactorial(n-2)
  _save(DoubleFactorial,n,val)
  return val
end
_make(DoubleFactorial,Integer)

"""
`Binomial(n,k)` returns the binomial coefficient `n`-choose-`k`.
This is the number of `k`-element subsets of an `n`-element set.
"""
Binomial(n::Integer, k::Integer) = binomial(big(n),big(k))::BigInt

"""
`Multinomial(vec)` returns the multinomial coefficient whose
top index is the sum of `vec` (an array of `Int`s) and whose
bottom indices are given by `vec`.

This may also be called with a common-separated list of arguments,
that is, either of `Multinomial([1,2,3])` or `Multinomial(1,2,3)`.
The result is `60` in both cases as these equal `6!/(1! 2! 3!)`.

**Warning**: This is not the same as `MultiChoose`.
"""
function Multinomial(v...)::BigInt
  nv = length(v)
  for i=1:nv
    typeof(v[i])<:Integer || throw(DomainError(v,"arguments must be integers"))
    v[i]>=0 || throw(DomainError(v,"arguments must be nonngative"))
  end
  vals = [t for t in v]
  return Multinomial(vals)
end

Multinomial() = big(1)::BigInt

function Multinomial(vals::Vector{T})::BigInt where T<:Integer
  if any([t<0 for t in vals])
    throw(DomainError(vals,"arguments must be nonnegative"))
  end

  nv = length(vals)
  n  = sum(vals)
  # base cases
  if nv<=1 || n==0
    return big(1)
  end
  # reduce
  return Binomial(n,vals[end]) * Multinomial(vals[1:nv-1])
end

"""
`MultiChoose(n,k)` returns the number of `k`-element
*multisets* that can be formed using the elements of an
`n`-element set.

**Warning**: This is not the same as `Multinomial`.
"""
function MultiChoose(n::Integer,k::Integer)::BigInt
  return Binomial(n+k-1,k)
end

"""
`Catalan(n)` returns the `n`-th Catalan number.
"""
function Catalan(n::Integer)::BigInt
  n >= 0 || throw(DomainError(n,"argument must be nonnegative"))
  return div(Binomial(2n,n),n+1)
end

"""
`Derangements(n)` returns the number of permutations of
an `n`-set that have no fixed point.
"""
function Derangements(n::Integer)::BigInt
  if n<0
    throw(DomainError(n,"argument must be nonnegative"))
  end
  if n==0
    return big(1)
  end
  if n==1
    return big(0)
  end
  if _has(Derangements,n)
      return _get(Derangements,n)
  end

  val = (n-1)*(Derangements(n-1)+Derangements(n-2))
  _save(Derangements,n,val)
  return val
end

_make(Derangements,Integer)

"""
`Bell(n)` gives the `n`-th Bell number, that is,
the number of partitions of an `n`-element set.
"""
function Bell(n::Integer)::BigInt
  if n<0
    throw(DomainError(n,"argument must be nonnegative"))
  end
  if n==1 || n==0
    return big(1)
  end
  if _has(Bell,n)
      return _get(Bell,n)
  end
  N1 = n-1
  result = big(0)
  for k=0:n-1
    result += Binomial(n-1,k) * Bell(k)
  end
  _save(Bell,n,result)
  return result
end
_make(Bell,Integer)

"""
`Stirling2(n,k)` gives the Stirling number of the second kind,
that is, the number of paritions of an `n`-set into `k`-parts."
"""
function Stirling2(n::Integer,k::Integer)::BigInt
  # special cases
  if k<0 || n<0
    throw(DomainError((n,k),"arguments must be nonnegative"))
  end

  if k>n
    return big(0)
  end

  if n==0  # and by logic, k==0
    return big(1)
  end

  if k==0
    return big(0)
  end

  if n==k
    return big(1)
  end
  # END OF SPECIAL CASES, invoke recursion
  if _has(Stirling2,(n,k))
      return _get(Stirling2,(n,k))
  end

  val = Stirling2(n-1,k-1) + Stirling2(n-1,k)*k
  _save(Stirling2,(n,k),val)
  return val
end
_make(Stirling2,Tuple{Integer,Integer})

"""
`Stirling1(n,k)` gives the (signed) Stirling number
of the first kind, that is, the coefficient of `x^k`
in the poynomial `x(x-1)(x-2)...(x-n+1)`.
"""
function Stirling1(n::Integer,k::Integer)::BigInt
  # special cases
  if k<0 || n<0
    throw(DomainError((n,k),"arguments must be nonnegative"))
  end

  if k>n
    return big(0)
  end

  if n==0  # and, by logic, k==0
    return big(1)
  end

  if k==0  # and, by logic, n>0
    return big(0)
  end

  if _has(Stirling1,(n,k))
      return _get(Stirling1,(n,k))
  end

  # end of special cases, invoke recursion

  val = Stirling1(n-1,k-1) - (n-1)*Stirling1(n-1,k)
  _save(Stirling1,(n,k),val)
  return val
end

_make(Stirling1,Tuple{Integer,Integer})

"""
`IntPartitions(n)` is the number of partitions of the integer `n`.

`IntPartitions(n,k)` is the number of partitions of the integer
`n` with exactly `k` (nonzero) parts.
"""
function IntPartitions(n::Integer,k::Integer)::BigInt
  if n<0 || k<0
    throw(DomainError((n,k),"arguments must be nonnegative"))
  end
  # lots of special cases
  if k>n
    return big(0)
  end
  if n==0
    return big(1)
  end
  if k==0
    return big(0)
  end
  if k==n || k==1
    return big(1)
  end
  if _has(IntPartitions,(n,k))
      return _get(IntPartitions,(n,k))
  end

  val =  sum([IntPartitions(n-k,i) for i=0:k])
  _save(IntPartitions,(n,k),val)
  return val
end

function IntPartitions(n::Integer)::BigInt
    if n<0
        throw(DomainError(n,"argument must be nonnegative"))
    end
    if _has(IntPartitions,n)
        return _get(IntPartitions,n)
    end
    val = sum([IntPartitions(n,k) for k=0:n])
    _save(IntPartitions,n,val)
    return val
end

_make(IntPartitions,Union{Tuple{Integer,Integer},Integer})

"""
`IntPartitionsDistinct(n,k)` is the number of partitions of
the integer `n` into exactly `k` *distinct* parts.

`IntPartitionsDistinct(n)` is the number of partitions of `n`
into *distinct* parts.
"""
function IntPartitionsDistinct(n::Integer,k::Integer)::BigInt
  if n<0 || k<0
    throw(DomainError((n,k),"arguments must be nonnegative"))
  end
  Ck2 = div(k*(k-1),2)
  if n < Ck2
    return big(0)
  end
  return IntPartitions(n-Ck2,k)
end

function IntPartitionsDistinct(n::Integer)::BigInt
  if n<0
    throw(DomainError(n,"argument must be nonngative"))
  end
  result = big(0)
  for k=1:n
    s = IntPartitionsDistinct(n,k)
    if s==0
      break
    end
    result += s
  end
  return result
end

"""
`Euler(n)` returns the `n`-th Euler number. Starting with `n=0`
this is the sequence
1, 0, -1, 0, 5, 0, -61, 0, 1385 and so on.
"""
function Euler(n::Integer)::BigInt
  n>=0 || throw(DomainError(n,"argument must be nonngative"))
  if n%2 == 1
    return big(0)
  end
  if n==0
    return big(1)
  end
  if _has(Euler,n)
      return _get(Euler,n)
  end

  last = div(n,2)-1
  val =  -sum([ Binomial(n,2k)*Euler(2k) for k=0:last])
  _save(Euler,n,val)
  return val
end
_make(Euler,Integer)

"""
`PowerSum(n,k)` returns the sum of the `k`-th powers of the
integers `1` through `n`, i.e.,
`1^k + 2^k + 3^k + ... + n^k`.
"""
function PowerSum(n::Integer, k::Integer)::BigInt
  (n>=0 && k>=0) || throw(DomainError((n,k),"arguments must be nonngative"))
  # Base and special cases
  if n==0
    return big(0)
  end
  if k==0
    return big(n)
  end
  if k==1
    return Binomial(n,2)
  end
  if _has(PowerSum,(n,k))
      return _get(PowerSum,(n,k))
  end

  val =  big(n)^k + PowerSum(n-1,k)
  _save(PowerSum,(n,k),val)
  return val
end
_make(PowerSum,Tuple{Integer,Integer})

include("eulerian.jl")

end  #end of module


text:
# Automatically generated using Clang.jl wrap_c

const leveldb_t = Cvoid
const leveldb_cache_t = Cvoid
const leveldb_comparator_t = Cvoid
const leveldb_env_t = Cvoid
const leveldb_filelock_t = Cvoid
const leveldb_filterpolicy_t = Cvoid
const leveldb_iterator_t = Cvoid
const leveldb_logger_t = Cvoid
const leveldb_options_t = Cvoid
const leveldb_randomfile_t = Cvoid
const leveldb_readoptions_t = Cvoid
const leveldb_seqfile_t = Cvoid
const leveldb_snapshot_t = Cvoid
const leveldb_writablefile_t = Cvoid
const leveldb_writebatch_t = Cvoid
const leveldb_writeoptions_t = Cvoid


text:
#Matrices related to stochastic processes

import Base: start, next, done
export AiryProcess, BrownianProcess, WhiteNoiseProcess, next!

abstract StochasticProcess{T<:Real}

immutable WhiteNoiseProcess{T<:Real} <: StochasticProcess{T}
    dt::T
end

immutable BrownianProcess{T<:Real} <: StochasticProcess{T}
    dt::T
end

immutable AiryProcess{S<:Real, T<:Real} <: StochasticProcess{T}
    dt::T
    beta::S
end

done{T}(p::StochasticProcess{T}, x...) = false #Processes can always go on forever

#######################
# White noise process #
#######################

start{T}(p::WhiteNoiseProcess{T}) = nothing
next{T}(p::WhiteNoiseProcess{T}, x::Void=nothing) = (randn()*sqrt(p.dt), nothing)

####################
# Brownian process #
####################

start{T}(p::BrownianProcess{T}) = zero(T)
function next{T}(p::BrownianProcess{T}, x::T)
    newx = x + randn()*sqrt(p.dt)
    newx, newx
end

################
# Airy process #
################
start{T}(p::AiryProcess{T}) = SymTridiagonal(T[-(2/p.dt^2)], T[])

"""
Like next, but update only the state of the AiryProcess

Skip the eigenvalue computation, which gets expensive
"""
function next!{T}(p::AiryProcess{T}, S::SymTridiagonal{T})
    t = (size(S, 1)-1)*p.dt

    #Discretized Airy operator plus diagonal noise
    x = inv(p.dt^2)
    push!(S.dv, -2x - t + 2/sqrt(p.dt*p.beta)*randn())
    push!(S.ev, x)

    S
end
function next{T}(p::AiryProcess{T}, S::SymTridiagonal{T})
    t = (size(S, 1)-1)*p.dt

    #Discretized Airy operator plus diagonal noise
    x = inv(p.dt^2)
    push!(S.dv, -2x - t + 2/sqrt(p.dt*p.beta)*randn())
    push!(S.ev, x)

    (eigmax(S), S)
end


text:
@testset "766.toeplitz-matrix" begin
    @test is_toeplitz_matrix([[1, 2, 3, 4], [5, 1, 2, 3], [9, 5, 1, 2]]) == true
    @test is_toeplitz_matrix([[1, 2], [2, 2]]) == false
end


text:
module ShellModule
export AbstractShell, Shell, expandShell
using ..BaseModule
using ..BasisFunctionsModule
using ..BasisSetModule
using ..GeometryModule
import Base.display

abstract type AbstractShell end

"""
    Shell(center::Position, lqn::LQuantumNumber, exponents::[Float], coefficients::[Float64]; renorm=true)

A shell of basis functions is the collection of basis functions that share a center and radial definition (e.g.
the px, py and pz orbital that are defined together by their center, exponents and coefficients. Contrary to other codes we also consider a shell to have a unique LQuantumNumber not allowing for shared-sp tricks. Working with shells
instead of the individual basis functions allows for some performance gains in integral evaluation.
If renorm is set to false, then the coefficients are read in "as is". With renorm set to true the coefficients are expected to be given as with a basis set file and normalized accordingly. This is the default in accordance with LibInt2Shell, where this is necessary due to performance reasons. (As such the normalization scaling is taken from the renorm() function in libint2/shell.h.

In a shell, there is no discrimination between MQuantumNumbers. This allows for efficient recursion schemes, but poses a challenge to normalization, because the normalization factor depends explicitly on the MQuantumNumber. Following the choice of libint2 (see: https://github.com/evaleev/libint/blob/80d829b0e9932cd118a85e8aab03a41fe0699a7b/doc/progman/progman.tex#L498) to scale all coefficients according to the axial MQuantumNumbers (e.g. x, yyy, zz, ...). After evaluation of the integrals, those including non-axial (xy, xz, ...) MQuantumNumbers need to be rescaled by the relative factors (`renormNonAxial!`).

[Fermann,Valeev: Fundamentals of Molecular Integrals Evaluation, Eq. 2.25]
"""
type Shell <: AbstractShell
  center::Position
  lqn::LQuantumNumber
  exponents::Vector{Float64}
  coefficients::Vector{Float64}
  function Shell(center,lqn,exponents,coefficients;renorm=true)
    L = lqn.exponent
    if renorm
      nprim = length(exponents)
      # renorm primitive-by-primitive
      for p in 1:nprim
	coefficients[p] *= sqrt( (2^L*(2*exponents[p])^(L+3/2)) / (π^(3/2)*doublefactorial(2*L-1)) )
      end
      # normalize total shell
      norm = 0.
      for p in 1:nprim, q in 1:nprim
        norm += (doublefactorial(2*L-1)*π^(3/2)*coefficients[p]*coefficients[q])/
	              ((2^L)*(exponents[p]+exponents[q])^(L+3/2))
      end
      coefficients /= sqrt(norm)
    end
    new(center,lqn,exponents,coefficients)
  end
end

function computeFactorsNonAxial(lqn::LQuantumNumber)
  L = lqn.exponent
  result = Vector{Float64}()
  for mqn in MQuantumNumbers(lqn)
    na = doublefactorial(2*mqn.x-1)*doublefactorial(2*mqn.y-1)*doublefactorial(2*mqn.z-1)/doublefactorial(2*L-1)
    push!(result,1/√na)
  end
  return result
end

function computeFactorsNonAxial(shells::Vector{Shell})
  naxfacs = Vector{Float64}()
  for sh in shells
    push!(naxfacs,computeFactorsNonAxial(sh.lqn)...)
  end
  return naxfacs
end

function renormNonAxial!(vec::Vector{Float64},shells::Vector{<:AbstractShell})
  naxfacs = computeFactorsNonAxial(shells)
  vec .*= naxfacs
end

function renormNonAxial!(cgbfs::Vector{ContractedGaussianBasisFunction},shells::Vector{<:AbstractShell})
  naxfacs = computeFactorsNonAxial(shells)
  for (cgbf,fac) in zip(cgbfs,naxfacs)
    scale!(cgbf.coefficients,fac)
  end
end

function renormNonAxial!(mat::Matrix{Float64},shells::Vector{<:AbstractShell})
  naxfacs = computeFactorsNonAxial(shells)
  scale!(naxfacs,mat)
  scale!(mat,naxfacs)
end

function renormNonAxial!(mat::Matrix{Float64},shells1::Vector{<:AbstractShell},shells2::Vector{<:AbstractShell})
  naxfacs1 = computeFactorsNonAxial(shells1)
  scale!(naxfacs1,mat)
  naxfacs2 = computeFactorsNonAxial(shells2)
  scale!(mat,naxfacs2)
end

function renormNonAxial!(tensor::Array{Float64,N},shells::Vararg{Vector{<:AbstractShell},N}) where {N}
  naxfacs = computeFactorsNonAxial.(shells)
  for idx in CartesianRange(size(tensor))
    for dim in 1:length(shells)
      tensor[idx] *= naxfacs[dim][idx.I[dim]]
    end
  end
end

function display(sh::Shell)
  println("$(sh.lqn.symbol)-Shell   @    $(sh.center)")
  print("  Exponents:    ")
  for exp in sh.exponents
    @printf("  %f",exp)
  end
  println("")
  print("  Coefficients: ")
  for coeff in sh.coefficients
    @printf("  %f",coeff)
  end
  println("")
end

"""
Compute a list of ContractedGaussianBasisFunctions which are described within the input shell.
"""
function expandShell(sh::Shell)
  result = ContractedGaussianBasisFunction[]
  for mqn in MQuantumNumbers(sh.lqn)
    pgbfs = [PrimitiveGaussianBasisFunction(sh.center,exp,mqn) for exp in sh.exponents]
    cgbf = ContractedGaussianBasisFunction(copy(sh.coefficients),pgbfs)
    push!(result,cgbf)
  end
  renormNonAxial!(result,[sh])
  return result
end

"""
    nbf(shell)
computes the number of basisfunctions (cartesian) described by the shell.
"""
nbf(shell::Shell) = BaseModule.numberMQNsCartesian(shell.lqn)
end # module


text:
"""
    plot_confidence_timeseries!(p::Plot, data::AbstractMatrix{<:Real}; label="", kwargs...)

Plots confidence intervals for the time-series represented by `data`.
Assumes each row corresponds to the samples for a single time-step.
"""
function plot_confidence_timeseries!(p::Plots.AbstractPlot, data::AbstractMatrix; no_label = false, label="", kwargs...)
    intervals = [0.025, 0.25, 0.5, 0.75, 0.975]

    qs = [quantile(v, intervals) for v in eachrow(data)]
    llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )
    plot!(mq, ribbon=(mq - llq, uuq - mq), linewidth=0, label=(no_label ? "" : "$(label) (95% quantiles)"), kwargs...)
    plot!(mq, ribbon=(mq - lq, uq - mq), linewidth=2, label=(no_label ? "" : "$(label) (50% quantiles)"), kwargs...)

    return p
end

"""
    plot_confidence_timeseries!(p::Plot, time::AbstractVector, data::AbstractMatrix{<:Real}; label="", kwargs...)

Plots confidence intervals for the time-series represented by `data` and a vector of Dates.
Assumes each row corresponds to the samples for a single time-step.
"""
function plot_confidence_timeseries!(p::Plots.AbstractPlot, time::AbstractVector{<:Dates.TimeType}, data::AbstractMatrix; no_label = false, label="", c=:peru, kwargs...)
    intervals = [0.025, 0.25, 0.5, 0.75, 0.975]

    qs = [quantile(v, intervals) for v in eachrow(data)]
    llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )
    strtime = time .|> Symbol .|> String
    #plot!(time, mq, ribbon=(mq - llq, uuq - mq), c=c, α=0.5, linewidth=0, label=(no_label ? "" : "$(label) (95% quantiles)"), hover=strtime, kwargs...)
    plot!(time, mq, ribbon=(mq - lq, uq - mq), c=c, α=0.5, linewidth=2, label=(no_label ? "" : "$(label) (50% quantiles)"), hover=mq, kwargs...)

    return p
end

"""
    plot_confidence_timeseries(data::AbstractMatrix{<:Real}; label="", kwargs...)

See `plot_confidence_timeseries!`.
"""
plot_confidence_timeseries(data::AbstractVector; kwargs...) = plot_confidence_timeseries(hcat(data...) ; kwargs...)
plot_confidence_timeseries(data::AbstractMatrix; kwargs...) = plot_confidence_timeseries!(plot(), data; kwargs...)

"""
    plot_confidence_timeseries(data::AbstractMatrix{<:Real}; label="", kwargs...)

See `plot_confidence_timeseries!`.
"""
plot_confidence_timeseries(
    time::AbstractVector{<:Dates.TimeType},
    data::AbstractVector;
    kwargs...
) = plot_confidence_timeseries!(plot(), time, hcat(data...); kwargs...)
plot_confidence_timeseries!(
    p   ::Plots.AbstractPlot,
    time::AbstractVector{<:Dates.TimeType},
    data::AbstractVector;
    kwargs...
) = plot_confidence_timeseries!(p, time, hcat(data...); kwargs...)
plot_confidence_timeseries(
      time::AbstractVector{<:Dates.TimeType}
    , data::AbstractMatrix; kwargs...
) = plot_confidence_timeseries!(plot(), time, data; kwargs...)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# Hacky until the package is registered, then there should be a Project.toml in docs with dependencies listed
push!(LOAD_PATH, "../")

using Documenter, MPCCModelManager

makedocs(
    sitename = "MPCCModelManager.jl",
    authors = "Peter Maxwell",
    modules = [MPCCModelManager],
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "manual/function_list.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/PeterMaxwell/MPCCModelManager.jl.git",
)


text:
## Test π₀ estimator methods ##
module Test_pi0

using MultipleTesting
using Test
using StatsBase

@testset "π₀ estimators" begin

    ## test case: deterministic
    p0 = collect(0.01:0.01:1)
    p1 = p0 .^ 10
    p = [p0; p1] ## unsorted
    pi0 = length(p0) / length(p)

    lambdas = collect(0.05:0.05:0.95)

    unsort = MultipleTesting.unsort

    @testset "Storey π₀" begin

        @test typeof(Storey()) <: Pi0Estimator
        @test typeof(Storey(0.5)) <: Pi0Estimator
        @test estimate_pi0(p, Storey(0.2)) ≈ 0.6
        @test estimate_pi0(p, Storey(0.0)) ≈ 1.0
        @test estimate_pi0(p, Storey(1.0)) ≈ 1.0

        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test estimate_pi0(p_unsort, Storey(0.2)) ≈ 0.6
        @test !issorted(p_unsort)

        @test_throws DomainError Storey(-0.1)
        @test_throws DomainError Storey(1.1)
        #@test_throws DomainError estimate_pi0(p, Storey(1.0)) ## CHCK

    end

    @testset "StoreyBootstrap π₀" begin

        q = 0.1

        @test typeof(StoreyBootstrap()) <: Pi0Estimator

        @test estimate_pi0(p, StoreyBootstrap()) ≈ 0.6
        @test estimate_pi0(p0, StoreyBootstrap()) ≈ 1.0
        @test estimate_pi0(p1, StoreyBootstrap()) ≈ 0.15

        @test estimate_pi0(p, StoreyBootstrap(lambdas, q)) ≈ 0.6
        @test estimate_pi0(p0, StoreyBootstrap(lambdas, q)) ≈ 1.0
        @test estimate_pi0(p1, StoreyBootstrap(lambdas, q)) ≈ 0.15

        ## unsorted lambdas
        lambdas_unsort = unsort(lambdas)
        @test !issorted(lambdas_unsort)
        @test estimate_pi0(p, StoreyBootstrap(lambdas_unsort, q)) ≈ 0.6
        @test !issorted(lambdas_unsort)

        @test_throws DomainError StoreyBootstrap(lambdas, -0.1)
        @test_throws DomainError StoreyBootstrap(lambdas, 1.1)
        @test_throws MethodError StoreyBootstrap(0.5)
        @test_throws MethodError StoreyBootstrap(lambdas)
        @test_throws MethodError StoreyBootstrap(0.5, lambdas)

    end

    @testset "LeastSlope π₀" begin

        @test typeof(LeastSlope()) <: Pi0Estimator

        ## checked against structSSI::pi0.lsl
        @test isapprox( estimate_pi0(p, LeastSlope()), 0.62, atol = 1e-2 )
        @test isapprox( estimate_pi0(p0, LeastSlope()), 1.0, atol = 1e-2 )
        @test isapprox( estimate_pi0(p1, LeastSlope()), 0.16, atol = 1e-2 )

        ## check against internal reference implementation
        @test estimate_pi0(p, LeastSlope()) ≈ MultipleTesting.lsl_pi0_vec(p)
        @test estimate_pi0(p0, LeastSlope()) ≈ MultipleTesting.lsl_pi0_vec(p0)
        @test estimate_pi0(p1, LeastSlope()) ≈ MultipleTesting.lsl_pi0_vec(p1)

        @test_throws MethodError LeastSlope(0.1)

        ## unsorted p-values
        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test isapprox( estimate_pi0(p_unsort, LeastSlope()), 0.62, atol = 1e-2 )
        @test !issorted(p_unsort)

        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test isapprox( MultipleTesting.lsl_pi0_vec(p_unsort), 0.62, atol = 1e-2 )
        @test !issorted(p_unsort)

    end

    @testset "Oracle π₀" begin

        @test estimate_pi0(p, Oracle(0.5)) == 0.5
        @test estimate_pi0(p0, Oracle(0.6)) == 0.6
        @test estimate_pi0(p1, Oracle()) == 1.0

    end

    @testset "TwoStep π₀" begin

        alpha = 0.05

        @test typeof(TwoStep(alpha)) <: Pi0Estimator

        ## checked against mutoss::TSBKY_pi0_est
        @test estimate_pi0(p, TwoStep()) ≈ 0.665
        @test estimate_pi0(p, TwoStep(alpha)) ≈ 0.665
        @test estimate_pi0(p0, TwoStep(alpha)) ≈ 1.0
        @test estimate_pi0(p1, TwoStep(alpha)) ≈ 0.29
        @test estimate_pi0(p, TwoStep(alpha, BenjaminiHochberg())) ≈ 0.665

        @test estimate_pi0(p, TwoStep(0.1)) ≈ 0.63
        @test estimate_pi0(p, TwoStep(0.0)) ≈ 1.0
        @test estimate_pi0(p, TwoStep(1.0)) ≈ 0.415
        @test estimate_pi0(p, TwoStep(0.1, BenjaminiHochberg())) ≈ 0.63

        ## unsorted p-values
        p_unsort = unsort(p)
        @test !issorted(p_unsort)
        @test estimate_pi0(p_unsort, TwoStep(alpha)) ≈ 0.665
        @test !issorted(p_unsort)

    end

    @testset "RightBoundary π₀" begin

        # R code tested against (note this uses equidistant grid so we have to use
        # λseq as in Storey for comparison
        # ```R
        # library(pi0)
        # p0 <- seq(0.01,1, by=0.01)
        # histf1(p0, rightBoundary=TRUE)
        # [1] 1
        # p1 <- p0^10
        # histf1(p1, rightBoundary=TRUE)
        # [1] 0.1428571
        # p <- c(p0,p1)# histf1(p, rightBoundary=TRUE)
        # [1] 0.5714286
        # ```

        @test typeof(RightBoundary()) <: Pi0Estimator
        @test typeof(RightBoundary(lambdas)) <: Pi0Estimator

        # only eps because pi0 package uses right closed histograms
        @test isapprox( estimate_pi0(p, RightBoundary(lambdas)), 0.5714286, atol = 0.02 )
        @test estimate_pi0(p0, RightBoundary(lambdas)) ≈ 1.0
        @test isapprox( estimate_pi0(p1, RightBoundary(lambdas)), 0.1428571, atol = 1e-7 )

        @test estimate_pi0(p, RightBoundary(lambdas)) ≈ estimate_pi0(p, RightBoundary(unsort(lambdas)))

        # not checked against R implementation but should hold (check for default lambda grid)
        @test estimate_pi0(p0, RightBoundary()) ≈ 1.0

        p_unsort = unsort(p1)
        @test !issorted(p_unsort)
        @test isapprox( estimate_pi0(p_unsort, RightBoundary(lambdas)), 0.1428571, atol = 1e-7 )
        @test !issorted(p_unsort)

    end

    @testset "CensoredBUM π₀" begin

        @test typeof(CensoredBUM()) <: Pi0Estimator
        @test typeof(CensoredBUM(0.2, 0.1)) <: Pi0Estimator

        @test isapprox( estimate_pi0(p, CensoredBUM()), 0.55797, atol = 1e-5 )
        @test isapprox( estimate_pi0(p0, CensoredBUM()), 1.0, atol = 1e-5 )
        @test isapprox( estimate_pi0(p1, CensoredBUM()), 0.11608, atol = 2e-5 )

        # test against internal reference implementation
        # p0 case is not handled well by naive implementation
        @test estimate_pi0(p, CensoredBUM()) ≈ MultipleTesting.cbum_pi0_naive(p)[1]
        @test estimate_pi0(p1, CensoredBUM()) ≈ MultipleTesting.cbum_pi0_naive(p1)[1]

        ## test case that does not converge
        @test isnan(estimate_pi0(p, CensoredBUM(0.5, 0.05, 1e-6, 2)))

        @test_throws DomainError CensoredBUM(-0.5, 0.05)
        @test_throws DomainError CensoredBUM(1.5, 0.05)
        @test_throws DomainError CensoredBUM(0.5, -0.05)
        @test_throws DomainError CensoredBUM(0.5, 1.05)

        @test_throws DomainError CensoredBUM(0.5, 0.05, -1e-6, 100)
        @test_throws DomainError CensoredBUM(0.5, 0.05, 1.1, 100)
        @test_throws DomainError CensoredBUM(0.5, 0.05, 1e-6, -10)

        f = fit(CensoredBUM(), p)
        @test typeof(f) <: CensoredBUMFit
        @test isapprox( f.π0, 0.55797, atol = 1e-5 )

        # denominator becomes 0 if all p-values are 1
        @test isnan( estimate_pi0(ones(20), CensoredBUM()) )

        pi0_est, pars, is_converged = MultipleTesting.cbum_pi0(ones(20))
        @test isnan( pi0_est )
        @test !is_converged

        pi0_est, pars, is_converged = MultipleTesting.cbum_pi0_naive(ones(20))
        @test isnan( pi0_est )
        @test !is_converged

    end

    @testset "BUM π₀" begin

        @test isapprox( estimate_pi0(p, BUM(0.5)), 0.55528, atol = 1e-5 )
        @test isapprox( estimate_pi0(p, BUM()), 0.55528, atol = 1e-5 )
        @test isapprox( estimate_pi0(p0, BUM()), 1.0, atol = 1e-5 )
        @test isapprox( estimate_pi0(p1, BUM()), 0.10874, atol = 1e-5 )

        ## test case that does not converge
        @test isnan(estimate_pi0(p, BUM(0.5, 1e-6, 2)))

        @test typeof(BUM()) <: Pi0Estimator

        @test_throws DomainError BUM(-0.5)
        @test_throws DomainError BUM(1.5)

    end

    @testset "FlatGrenander π₀" begin

        @test typeof(FlatGrenander()) <: Pi0Estimator

        pu = collect(0.1:0.05:0.9)
        @test estimate_pi0(pu, FlatGrenander()) ≈ 1.0
        @test estimate_pi0(pu.^0.5, FlatGrenander()) ≈ 1.0

        # no reference values from literature/software available
        @test estimate_pi0(p0, FlatGrenander()) ≈ 1.0
        @test isapprox( estimate_pi0(p1, FlatGrenander()), 0.10458, atol = 1e-4 )
        @test isapprox( estimate_pi0(p, FlatGrenander()), pi0, atol = 0.1)

        ## longest constant interval: low level
        lci = MultipleTesting.longest_constant_interval

        px = [0:0.1:1;];
        f = [1.5, 1.5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.5

        f = [1.5, 1.5, 1.5, 1.5, 1.5, 1.2, 1.2, 1.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.2

        f = [0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.2

        f = [0.5, 0.5, 0.5, 0.5, 0.5, 0.2, 0.2, 0.2, 0.2, 0.2, 0.1];
        @test lci(px, f) ≈ 0.5

        px = [0.1, 0.3, 0.5, 0.9];
        f = [0.5, 0.5, 0.2, 0.2];
        @test lci(px, f) ≈ 0.2

        px = [0.1, 0.5, 0.7, 0.9];
        f = [0.5, 0.5, 0.2, 0.2];
        @test lci(px, f) ≈ 0.5

    end

    @testset "ConvexDecreasing π₀" begin

        @test typeof(ConvexDecreasing()) <: Pi0Estimator
        @test typeof(ConvexDecreasing(100, 1e-6, 1000)) <: Pi0Estimator

        @test isapprox( estimate_pi0(p, ConvexDecreasing()), 0.5739054, atol = 1e-4)
        @test isapprox( estimate_pi0(p0, ConvexDecreasing()), 1.0, atol = 1e-6 )
        @test isapprox( estimate_pi0(p1, ConvexDecreasing()), 0.1323284, atol = 1e-3)

        @test isapprox( estimate_pi0(ones(20), ConvexDecreasing()), 1.0, atol = 1e-6 )

        ## test case that does not converge
        @test isnan(estimate_pi0(p, ConvexDecreasing(100, 1e-6, 2)))

        @test_throws DomainError ConvexDecreasing(100, 2.0, 1000)
        @test_throws DomainError ConvexDecreasing(100, -1.0, 1000)
        @test_throws DomainError ConvexDecreasing(-100, 0.01, 1000)
        @test_throws DomainError ConvexDecreasing(100, 0.01, -1000)

        f = fit(ConvexDecreasing(), p)
        @test typeof(f) <: ConvexDecreasingFit
        @test isapprox( f.π0, 0.5739054, atol = 1e-4 )

    end

end

end


text:
using ErrorLineNumber
using ErrorLineNumber.@errorlinenumber
using Base.Test

@errorlinenumber function erroneous()
    x = [1,2]
    return x[3]
end

function runtest3()
    (outRead, outWrite) = redirect_stderr()
    msg = ""
    try
        erroneous()
    catch
        msg = readavailable(outRead)
    finally
        close(outWrite)
        close(outRead)
    end
    matchstring = "    !!!!!!!!!!!!!!!!!"
    @test ASCIIString(msg[1:length(matchstring)]) == matchstring
end

runtest3()


text:
using Quaternions
using DualNumbers
using ForwardDiff
using LinearAlgebra
using Random
using Test

function fdquat_to_dualquat(q::Quaternion{<:ForwardDiff.Dual})
    q0 = quat(
        ForwardDiff.value(q.s),
        ForwardDiff.value(q.v1),
        ForwardDiff.value(q.v2),
        ForwardDiff.value(q.v3),
    )
    qe = quat(
        ForwardDiff.partials(q.s, 1),
        ForwardDiff.partials(q.v1, 1),
        ForwardDiff.partials(q.v2, 1),
        ForwardDiff.partials(q.v3, 1),
    )
    return dualquat(q0, qe)
end

function dualquat_to_fdquat(dq::DualQuaternion)
    return quat(
        ForwardDiff.Dual(dq.q0.s, dq.qe.s),
        ForwardDiff.Dual(dq.q0.v1, dq.qe.v1),
        ForwardDiff.Dual(dq.q0.v2, dq.qe.v2),
        ForwardDiff.Dual(dq.q0.v3, dq.qe.v3),
    )
end

@testset "DualQuaternion" begin
    @testset "type aliases" begin
        @test DualQuaternionF16 === DualQuaternion{Float16}
        @test DualQuaternionF32 === DualQuaternion{Float32}
        @test DualQuaternionF64 === DualQuaternion{Float64}
    end

    @testset "Constructors" begin
        @testset "from quaternions" begin
            q0 = Quaternion{Int}(1, 2, 3, 4, false)
            qe = QuaternionF32(5, 6, 7, 8, false)
            @test @inferred(DualQuaternion(q0, qe)) isa DualQuaternionF32
            @test DualQuaternion(q0, qe) === DualQuaternionF32(
                QuaternionF32(1, 2, 3, 4, false), QuaternionF32(5, 6, 7, 8, false), false
            )
            @test @inferred(DualQuaternionF64(q0, qe, false)) === DualQuaternionF64(
                QuaternionF64(1, 2, 3, 4, false), QuaternionF64(5, 6, 7, 8, false), false
            )
            @test DualQuaternionF64(q0, qe, true) === DualQuaternionF64(
                QuaternionF64(1, 2, 3, 4, false), QuaternionF64(5, 6, 7, 8, false), true
            )

            qnorm = Quaternion(0//1, 1//1, 0//1, 0//1, true)
            @test @inferred(DualQuaternion(q0)) isa DualQuaternion{Int}
            @test DualQuaternion(qnorm) === DualQuaternion(qnorm, zero(qnorm), true)
            @test @inferred(DualQuaternionF64(q0)) ===
                DualQuaternionF64(QuaternionF64(q0), zero(QuaternionF64), false)
        end
        @testset "from dual" begin
            coef = (dual(1, 2), dual(3.0, 4.0), dual(5.0f0, 6.0f0), dual(7//1, 8//1))
            @testset for T in (Float32, Float64, Int), norm in (true, false)
                dq = @inferred DualQuaternion{T}(coef..., norm)
                @test dq isa DualQuaternion{T}
                @test dq.norm === norm
                @test dq === DualQuaternion{T}(convert.(Dual{T}, coef)..., norm)
                @test dq == DualQuaternion(
                    Quaternion(DualNumbers.value.(coef)...),
                    Quaternion(DualNumbers.epsilon.(coef)...),
                )
                dq2 = @inferred DualQuaternion(convert.(Dual{T}, coef)..., norm)
                @test DualQuaternion(convert.(Dual{T}, coef)..., norm) === dq
                @test DualQuaternion(coef[1]) ==
                    DualQuaternion(coef[1], fill(zero(coef[1]), 3)...)
                @test DualQuaternion{T}(coef[1]) ==
                    DualQuaternion(convert(Dual{T}, coef[1]), fill(zero(Dual{T}), 3)...)
                if !norm
                    @test DualQuaternion(convert.(Dual{T}, coef)...) === dq
                end
            end
        end
        @testset "from real" begin
            @testset for x in (-1//1, 1.0, 2.0), T in (Float32, Float64, Int, Rational{Int})
                coef = (Quaternion{T}(x, 0, 0, 0, isone(abs(x))), zero(Quaternion{T}))
                @test @inferred(DualQuaternion{T}(x)) ===
                    DualQuaternion{T}(coef..., isone(abs(x)))
                @test @inferred(DualQuaternion(T(x))) ===
                    DualQuaternion{T}(coef..., isone(abs(x)))
            end
        end
        @testset "from dual quaternion" begin
            dq = DualQuaternion(
                QuaternionF32(1, 2, 3, 4, false), QuaternionF32(4, 5, 6, 7, false), false
            )
            @test @inferred(DualQuaternion(dq)) === dq
            @test @inferred(DualQuaternionF64(dq)) ===
                DualQuaternionF64(dq.q0, dq.qe, false)
        end
        @testset "from vector" begin
            v = randn(3)
            @test @inferred(DualQuaternion(v)) ===
                DualQuaternion(Quaternion(zero(v)), Quaternion(v))
        end
    end

    @testset "==" begin
        @test DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) ==
            DualQuaternion(Quaternion(1.0, 2, 3, 4), Quaternion(5, 6, 7, 8))
        @test DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) !=
            DualQuaternion(Quaternion(1.0, 2, 3, 4), Quaternion(1, 2, 3, 4))
    end

    @testset "convert" begin
        @test convert(DualQuaternion{Float64}, 1) === DualQuaternion(1.0)
        @test convert(DualQuaternion{Float64}, DualNumbers.Dual(1, 2)) ===
            DualQuaternion(Quaternion(1.0), Quaternion(2.0))
        @test convert(DualQuaternion{Float64}, Quaternion(1, 2, 3, 4)) ===
            DualQuaternion(Quaternion(1.0, 2.0, 3.0, 4.0))
        @test convert(
            DualQuaternion{Float64},
            DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)),
        ) === DualQuaternion(Quaternion(1.0, 2.0, 3.0, 4.0), Quaternion(5.0, 6.0, 7.0, 8.0))
    end

    @testset "promote" begin
        @test promote(DualQuaternion(1.0), 1.0) ===
            (DualQuaternion(1.0), DualQuaternion(1.0))
        @test promote(DualQuaternion(1.0f0), 2.0) ===
            (DualQuaternion(1.0), DualQuaternion(2.0))
        @test promote(DualQuaternion(1.0f0), dual(1, 2)) ===
            (DualQuaternion(1.0f0), DualQuaternion(dual(1.0f0, 2.0f0)))
        @test promote(DualQuaternion(1.0f0), Quaternion(3//1)) ===
            (DualQuaternion(1.0f0), DualQuaternion(3.0f0))
        @test promote(DualQuaternion(1.0f0), DualQuaternion(2.0)) ===
            (DualQuaternion(1.0), DualQuaternion(2.0))

        @test Quaternion(1, 2, 3, 4) == DualQuaternion(Quaternion(1, 2, 3, 4))
        @test Quaternion(1, 2, 3, 4) !=
            DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8))
        @test DualQuaternion(1) == 1.0
    end

    @testset "shorthands" begin
        @test dualquat(Quaternion(1, 0, 0, 0)) == Quaternion(1, 0, 0, 0)
        @test dualquat(Quaternion(1, 2, 3, 4)) == Quaternion(1, 2, 3, 4)
        @test dualquat(Quaternion(1, 0, 0, 0)) === DualQuaternion(Quaternion(1, 0, 0, 0)) # checking the .norm field in particular
        @test dualquat(Quaternion(1, 2, 3, 4)) === DualQuaternion(Quaternion(1, 2, 3, 4))
        @test dualquat(1) === DualQuaternion(1)
        @test dualquat(Dual(1, 2)) === DualQuaternion(Dual(1, 2))
        @test dualquat(Dual(1, 2), Dual(0), Dual(0), Dual(0)) ===
            DualQuaternion(Dual(1, 2), Dual(0), Dual(0), Dual(0))
        @test dualquat(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8)) ==
            DualQuaternion(Quaternion(1, 2, 3, 4), Quaternion(5, 6, 7, 8))
        @test dualquat(Quaternion(1, 0, 0, 0), Quaternion(0)).norm == false
        @test dualquat(Quaternion(1, 0, 0, 0), Quaternion(0), false).norm == false # respect the .norm input (even if wrong)
        @test dualquat(Quaternion(1, 2, 3, 4), Quaternion(0), true).norm == true # respect the .norm input (even if wrong)
        @test dualquat(Dual(2, 0), Dual(0), Dual(0), Dual(0), true).norm == true # respect the .norm input (even if wrong)
        @test dualquat(Dual(1, 0), Dual(0), Dual(0), Dual(0), false).norm == false # respect the .norm input (even if wrong)
    end

    @testset "random generation" begin
        @testset "dualquatrand" begin
            dq = dualquatrand()
            @test dq isa DualQuaternionF64
            @test !dq.norm
        end

        @testset "ndualquatrand" begin
            dq = ndualquatrand()
            @test dq isa DualQuaternionF64
            @test dq.norm
        end

        @testset "rand($H)" for H in (DualQuaternionF32, DualQuaternionF64)
            rng = Random.MersenneTwister(42)
            dq = rand(rng, H)
            @test dq isa H
            @test !dq.norm

            dqs = rand(rng, H, 1000)
            @test eltype(dqs) === H
            @test length(dqs) == 1000
            xs = map(dqs) do dq
                return [
                    real(dq.q0)
                    imag_part(dq.q0)...
                    real(dq.qe)
                    imag_part(dq.qe)...
                ]
            end
            xs_mean = sum(xs) / length(xs)
            xs_var = sum(x -> abs2.(x .- xs_mean), xs) / (length(xs) - 1)
            @test all(isapprox.(xs_mean, 0.5; atol=0.1))
            @test all(isapprox.(xs_var, 1 / 12; atol=0.01))
        end
    end

    @testset "basic" begin
        q = rand(DualQuaternionF64)
        qnorm = normalize(q)
        @test_throws MethodError imag(q)
        @test_throws MethodError Quaternions.imag(q)
        @test_throws MethodError imag_part(q)
        @test conj(q) === dualquat(conj(q.q0), conj(q.qe), q.norm)
        @test conj(qnorm) === dualquat(conj(qnorm.q0), conj(qnorm.qe), qnorm.norm)
        @test conj(conj(q)) === q
        @test conj(conj(qnorm)) === qnorm
        @test float(dualquat(dual.(1:4, 5:8)...)) === dualquat(dual.(1.0:4.0, 5.0:8.0)...)
        @test dconj(q) === dualquat(q.q0, -q.qe, q.norm)
        @test dconj(qnorm) === dualquat(qnorm.q0, -qnorm.qe, qnorm.norm)
        @test Quaternions.Q0(q) === q.q0
        @test Quaternions.Qe(q) === q.qe
    end

    @testset "algebraic properties" begin
        @testset "addition/subtraction" begin
            for _ in 1:100
                dq1, dq2 = rand(DualQuaternionF64, 2)
                @test (dq1 + dq2).q0 ≈ dq1.q0 + dq2.q0
                @test (dq1 + dq2).qe ≈ dq1.qe + dq2.qe
                @test (dq1 - dq2).q0 ≈ dq1.q0 - dq2.q0
                @test (dq1 - dq2).qe ≈ dq1.qe - dq2.qe
                @test (+dq1).q0 === dq1.q0
                @test (+dq1).qe === dq1.qe
                @test (-dq1).q0 === -dq1.q0
                @test (-dq1).qe === -dq1.qe
            end
        end
        @testset "division" begin
            for _ in 1:100
                dq = rand(DualQuaternionF64)
                dq2 = dq / dq
                @test dq2.q0 ≈ 1
                @test dq2.qe ≈ zero(dq2.qe) atol = 1e-6
                dq2 = dq \ dq
                @test dq2.q0 ≈ 1
                @test dq2.qe ≈ zero(dq2.qe) atol = 1e-6
                s = randn()
                @test (dq / s).q0 ≈ dq.q0 / s
                @test (dq / s).qe ≈ dq.qe / s
                d = dual(randn(2)...)
                @test (dq / d).q0 ≈ (dq / dualquat(d)).q0
                @test (dq / d).qe ≈ (dq / dualquat(d)).qe
            end
        end
        @testset "multiplication is associative" begin
            for _ in 1:100
                dq1, dq2, dq3 = rand(DualQuaternionF64, 3)
                dq4 = (dq1 * dq2) * dq3
                dq5 = dq1 * (dq2 * dq3)
                @test dq4.q0 ≈ dq5.q0
                @test dq4.qe ≈ dq5.qe
            end
        end
        @testset "multiplication by dual" begin
            for _ in 1:100
                dq = rand(DualQuaternionF64)
                d = dual(randn(2)...)
                @test (dq * d).q0 ≈ (dq * dualquat(d)).q0
                @test (dq * d).qe ≈ (dq * dualquat(d)).qe
                @test (d * dq).q0 ≈ (dualquat(d) * dq).q0
                @test (d * dq).qe ≈ (dualquat(d) * dq).qe
            end
        end
    end

    @testset "iszero" begin
        @test iszero(dualquat(0))
        @test !iszero(dualquat(1))
        @test !iszero(dualquat(quat(0, 1, 0, 0)))
        @test !iszero(dualquat(quat(0, 0, 1, 0)))
        @test !iszero(dualquat(quat(0, 0, 0, 1)))
        @test !iszero(dualquat(quat(0), quat(1)))
        @test !iszero(dualquat(quat(0), quat(0, 1, 0, 0)))
        @test !iszero(dualquat(quat(0), quat(0, 0, 1, 0)))
        @test !iszero(dualquat(quat(0), quat(0, 0, 0, 1)))
    end

    @testset "isone" begin
        @test isone(dualquat(1))
        @test !isone(dualquat(-1))
        @test !isone(dualquat(quat(0, 1, 0, 0)))
        @test !isone(dualquat(quat(1, 1, 0, 0)))
        @test !isone(dualquat(quat(1, 0, 1, 0)))
        @test !isone(dualquat(quat(1, 0, 0, 1)))
        @test !isone(dualquat(quat(0), quat(1)))
        @test !isone(dualquat(quat(0), quat(0, 1, 0, 0)))
        @test !isone(dualquat(quat(0), quat(0, 0, 1, 0)))
        @test !isone(dualquat(quat(0), quat(0, 0, 0, 1)))
    end

    @testset "^" begin
        @testset "^(::DualQuaternion, ::Real)" begin
            for _ in 1:100
                dq = rand(DualQuaternionF64)
                @test dq^2 === dq * dq
                @test dq^1 === dq
                @test_broken (dq^2.0).q0 ≈ (dq * dq).q0
                @test_broken (dq^2.0).qe ≈ (dq * dq).qe
                @test_broken (dq^1.0).q0 ≈ dq.q0
                @test_broken (dq^1.0).qe ≈ dq.qe
                @test_broken (dq^-1.0).q0 ≈ inv(dq).q0
                @test_broken (dq^-1.0).qe ≈ inv(dq).qe
                for p in (1.3, 7.8, 1.3f0, 7.8f0)
                    @test (dq^p).q0 ≈ exp(p * log(dq)).q0
                    @test (dq^p).qe ≈ exp(p * log(dq)).qe
                end
            end
        end
        @testset "^(::DualQuaternion, ::DualQuaternion)" begin
            for _ in 1:100
                q, p = randn(QuaternionF64, 2)
                @test q^p ≈ exp(p * log(q))
            end
        end
    end

    @testset "non-analytic functions" begin
        unary_funs = [conj, abs, abs2, norm, sign]
        # since every dual quaternion is conjugate to a dual complex number,
        # one can establish correctness as follows:
        @testset for fun in unary_funs
            for _ in 1:100
                dq1, dq2 = rand(DualQuaternionF64, 2)
                c = dual(randn(ComplexF64), randn(ComplexF64))
                dq = dualquat(reim(c)..., reim(zero(c))...)
                p = dualquat(@inferred(fun(dq)))
                @test p.q0 ≈ DualNumbers.value(fun(c))
                @test p.qe ≈ DualNumbers.epsilon(fun(c))
                p2 = dq2 * fun(dq1) * inv(dq2)
                p3 = dualquat(fun(dq2 * dq1 * inv(dq2)))
                @test p2.q0 ≈ p3.q0
                @test p2.qe ≈ p3.qe
            end
        end
    end

    @testset "analytic functions" begin
        unary_funs = [sqrt, inv, exp, log]
        # since every dual quaternion is conjugate to a dual complex number,
        # one can establish correctness as follows:
        @testset for fun in unary_funs
            (fun === log || fun === sqrt) && continue  # log is currently broken
            for _ in 1:100
                dq1, dq2 = rand(DualQuaternionF64, 2)
                c = dual(randn(ComplexF64), randn(ComplexF64))
                dq = dualquat(reim(c)..., reim(zero(c))...)
                p = dualquat(@inferred(fun(dq)))
                @test p.q0 ≈ DualNumbers.value(fun(c))
                @test p.qe ≈ DualNumbers.epsilon(fun(c))
                p2 = dq2 * fun(dq1) * inv(dq2)
                p3 = dualquat(fun(dq2 * dq1 * inv(dq2)))
                @test p2.q0 ≈ p3.q0
                @test p2.qe ≈ p3.qe
            end
        end

        @testset "identities" begin
            for _ in 1:100
                dq = rand(DualQuaternionF64)
                @test (inv(dq) * dq).q0 ≈ one(dq.q0)
                @test (inv(dq) * dq).qe ≈ zero(dq.qe) atol = 1e-6
                @test (dq * inv(dq)).q0 ≈ one(dq.q0)
                @test (dq * inv(dq)).qe ≈ zero(dq.qe) atol = 1e-6
                @test_broken (sqrt(dq) * sqrt(dq)).q0 ≈ dq.q0
                @test_broken (sqrt(dq) * sqrt(dq)).qe ≈ dq.qe
                @test_broken exp(log(dq)).q0 ≈ dq.q0
                @test_broken exp(log(dq)).qe ≈ dq.qe
                @test exp(zero(dq)).q0 ≈ one(dq.q0)
                @test exp(zero(dq)).qe ≈ zero(dq.qe) atol = 1e-6
                @test log(one(dq)).q0 ≈ zero(dq.q0) atol = 1e-6
                @test log(one(dq)).qe ≈ zero(dq.qe) atol = 1e-6
            end
            @test_broken log(zero(DualQuaternionF64)) === dualquat(-Inf)
        end
    end

    @testset "normalize" begin
        for _ in 1:100
            dq = rand(DualQuaternionF64)
            q = dualquat_to_fdquat(dq)
            dqnorm = normalize(dq)
            qnorm = normalize(q)
            dqnorm2 = fdquat_to_dualquat(qnorm)
            @test dqnorm.q0 ≈ dqnorm2.q0
            @test dqnorm.qe ≈ dqnorm2.qe
            @test normalize(dqnorm) === dqnorm
            @test iszero(normalize(zero(dq)))
        end
    end

    @testset "normalizea" begin
        for _ in 1:100
            dq = rand(DualQuaternionF64)
            q = dualquat_to_fdquat(dq)
            dqnorm, dqa = normalizea(dq)
            qnorm, qa = normalizea(q)
            dqnorm2 = fdquat_to_dualquat(qnorm)
            @test dqnorm.q0 ≈ dqnorm2.q0
            @test dqnorm.qe ≈ dqnorm2.qe
            @test DualNumbers.value(dqa) ≈ ForwardDiff.value(qa)
            @test DualNumbers.epsilon(dqa) ≈ ForwardDiff.partials(qa, 1)
        end
    end

    @testset "angle/axis/angleaxis" begin
        dq1 = rand(DualQuaternionF64)
        dq2 = dq1 + (1 - dq1.q0.s)
        for dq in (dq1, dq2)
            t1, p1 = angleaxis(dq)
            t2, v = angleaxis(dualquat_to_fdquat(dq))
            p2 = fdquat_to_dualquat(quat(v))
            @test DualNumbers.value(t1) ≈ ForwardDiff.value(t2)
            @test_broken DualNumbers.epsilon(t1) ≈ ForwardDiff.partials(t2, 1)
            @test p2.q0 ≈ p1.q0
            @test_broken p2.qe ≈ p1.qe
            t3 = angle(dq)
            @test DualNumbers.value(t3) ≈ DualNumbers.value(t1)
            @test DualNumbers.epsilon(t3) ≈ DualNumbers.epsilon(t1)
            p3 = axis(dq)
            @test p3.q0 ≈ p1.q0
            @test p3.qe ≈ p1.qe
        end
    end
end


text:
# ## Histogram Plot
# ### Histogram Plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Random
Random.seed!(2018)

x = randn(1000)
y = randn(1000)
z = randn(1000)

histogram(x, bins=20, alpha=0.4, label="A")
histogram!(y, bins=20, alpha=0.6, label="B")
histogram!(z, bins=20, alpha=0.8, label="C")


text:
using Plots
plotlyjs();

include("../julia/intensity_gamma.jl");
include("../julia/utils.jl");

x_true = readcsv("../data/spikes.csv")[:];
rs_true = readcsv("../data/firing_rate_true.csv")[:];

rs_mle = intensity_gamma(x_true);

plot(ts, rs_true, label="True")
plot!(ts, rs_mle, label="Estimated")

xlabel!("Time (s)")
ylabel!("Firing rate (Hz)")

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
struct HessianDFP{T,VT,MT,OBJ} <: Hessian{T}
    objective::OBJ

    x̄::VT    # previous solution
    x::VT    # current solution
    δ::VT    # difference of current and previous solution

    ḡ::VT    # previous gradient
    g::VT    # current gradient
    γ::VT    # difference of current and previous gradient

    Q::MT

    T1::MT
    T2::MT
    γγ::MT
    δδ::MT

    function HessianDFP(objective::MultivariateObjective, x::AbstractVector{T}) where {T}
        Q = alloc_h(x)

        T1 = zero(Q)
        T2 = zero(Q)
        γγ = zero(Q)
        δδ = zero(Q)

        new{T,typeof(x),typeof(Q),typeof(objective)}(objective, zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), Q, T1, T2, γγ, δδ)
    end
end

HessianDFP(F, x) = HessianDFP(MultivariateObjective(F, x), x)

function initialize!(H::HessianDFP, x::AbstractVector)
    H.Q .= Matrix(1.0I, size(H.Q)...)

    H.x̄ .= eltype(x)(NaN)
    H.δ .= eltype(x)(NaN)
    H.ḡ .= eltype(x)(NaN)
    H.γ .= eltype(x)(NaN)

    H.x .= x
    H.g .= gradient!(H.objective, x)

    return H
end

function update!(H::HessianDFP, x::AbstractVector)
    # copy previous data and compute new gradient
    H.ḡ .= H.g
    H.x̄ .= H.x
    H.x .= x
    H.g .= gradient!(H.objective, x)

    # δ = x - x̄
    H.δ .= H.x .- H.x̄

    # γ = g - ḡ
    H.γ .= H.g .- H.ḡ

    # γQγ = γᵀQγ
    γQγ = dot(H.γ, H.Q, H.γ)

    # δγ = δᵀγ
    δγ = H.δ ⋅ H.γ

    # DFP
    # Q = Q - ... + ...
    # H.Q .-= H.Q * H.γ * H.γ' * H.Q / (H.γ' * H.Q * H.γ) .-
    #         H.δ * H.δ' ./ δγ

    if δγ ≠ 0 && γQγ ≠ 0
        outer!(H.γγ, H.γ, H.γ)
        outer!(H.δδ, H.δ, H.δ)

        mul!(H.T1, H.γγ, H.Q)
        mul!(H.T2, H.Q, H.T1)

        H.Q .-= H.T2 ./ γQγ
        H.Q .+= H.δδ ./ δγ
    end
end

Base.inv(H::HessianDFP) = H.Q

Base.:\(H::HessianDFP, b) = inv(H) * b

LinearAlgebra.ldiv!(x, H::HessianDFP, b) = mul!(x, inv(H), b)


text:
# Tests for the functions library in deps/examples

using CxxWrap
using Base.Test
using Compat
cxx_available = true
try
  using Cxx
catch
  cxx_available = false
end

const functions_lib_path = CxxWrap._l_functions

# Wrap the functions defined in C++
wrap_modules(functions_lib_path)

# Test functions from the CppHalfFunctions module
@test CppHalfFunctions.half_d(3) == 1.5
@show methods(CppHalfFunctions.half_d)
@test CppHalfFunctions.half_i(-2) == -1
@test CppHalfFunctions.half_u(3) == 1
@test CppHalfFunctions.half_lambda(2.) == 1.
if get(ENV, "MSYSTEM", "") == "" || get(ENV, "APPVEYOR", "") == "" # Disabled on MINGW CI due to https://discourse.julialang.org/t/c-and-bits-types-confusion/1168/8
  @test CppHalfFunctions.strict_half(3.) == 1.5
end
@test_throws MethodError CppHalfFunctions.strict_half(3)

# Test functions from the CppTestFunctions module
@test CppTestFunctions.concatenate_numbers(4, 2.) == "42"
@test method_exists(CppTestFunctions.concatenate_numbers, (Union{Cint,CxxWrap.argument_overloads(Cint)...},Union{Cdouble,CxxWrap.argument_overloads(Cdouble)...}))
@test CppTestFunctions.concatenate_strings(2, "ho", "la") == "holahola"
@test CppTestFunctions.test_int32_array(Int32[1,2])
@test CppTestFunctions.test_int64_array(Int64[1,2])
@test CppTestFunctions.test_float_array(Float32[1.,2.])
@test CppTestFunctions.test_double_array([1.,2.])
if !(is_windows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException CppTestFunctions.test_exception()
end
ta = [1.,2.]
@test CppTestFunctions.test_array_len(ta) == 2
@test CppTestFunctions.test_array_get(ta, Int64(0)) == 1.
@test CppTestFunctions.test_array_get(ta, Int64(1)) == 2.
CppTestFunctions.test_array_set(ta, Int64(0), 3.)
CppTestFunctions.test_array_set(ta, Int64(1), 4.)
@test ta[1] == 3.
@test ta[2] == 4.
@test CppTestFunctions.test_type_name("IO") == "IO"

@test CppTestFunctions.test_long_long() == 42
@test CppTestFunctions.test_short() == 43

# Test GC protection array
a = "str1"
b = "str2"
c = "str3"
protect_arr = CxxWrap._gc_protected
start_len = length(protect_arr)
gcprotect(a)
gcprotect(b)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == a
@test protect_arr[end] == b
gcunprotect(a)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == nothing
@test protect_arr[end] == b
gcprotect(c)
@test length(protect_arr) == start_len + 2
@test protect_arr[end-1] == c
@test protect_arr[end] == b
@test CppTestFunctions.test_julia_call(1.,2.) == 2
@test CppTestFunctions.test_string_array(["first", "second"])
darr = [1.,2.]
CppTestFunctions.test_append_array!(darr)
@test darr == [1.,2.,3.]

testf(x,y) = x+y
@show c_func = safe_cfunction(testf, Float64, (Float64,Float64))
CppTestFunctions.test_safe_cfunction(c_func)
CppTestFunctions.test_safe_cfunction2(c_func)

function testf_arf(v::Vector{Float64}, s::String)
  r = sum(v)
  print_with_color(:green, "callback in Julia: $(s) = $(r)\n")
  return r
end

c_func_arf = safe_cfunction(testf_arf, Float64, (Any,Any))

CppTestFunctions.fn_clb(c_func_arf)
CppTestFunctions.fn_clb2(testf_arf)

function testf2(p::ConstPtr{Float64}, n_elems::Int)
  arr = ConstArray(p, n_elems)
  @test arr[1] == 1.0
  @test arr[2] == 2.0
  return
end
c_func2 = safe_cfunction(testf2, Void, (ConstPtr{Float64},Int))
CppTestFunctions.test_safe_cfunction3(c_func2)

dref = Ref(0.0)
CppTestFunctions.test_double_ref(dref)
@test dref[] == 1.0

@test CppTestFunctions.get_test_double() == 0.0
cppdref = CppTestFunctions.get_test_double_ref()
@test unsafe_load(cppdref) == 0.0
unsafe_store!(cppdref, 1.0)
@test CppTestFunctions.get_test_double() == 1.0

@test CppTestFunctions.test_const_string_return() == "test"
@test CppTestFunctions.test_datatype_conversion(Float64) == Float64

@test typeof(CppTestFunctions.test_double_pointer()) == Ptr{Float64}
@test CppTestFunctions.test_double_pointer() == C_NULL
@test typeof(CppTestFunctions.test_double2_pointer()) == Ptr{Ptr{Float64}}
@test CppTestFunctions.test_double2_pointer() == C_NULL
@test typeof(CppTestFunctions.test_double3_pointer()) == Ptr{Ptr{Ptr{Float64}}}
@test CppTestFunctions.test_double3_pointer() == C_NULL

@test CppTestFunctions.test_wstring_to_julia() == "šČô_φ_привет_일보"
@test CppTestFunctions.test_wstring_to_cpp("šČô_φ_привет_일보")

@test CppTestFunctions.real_part(2.0 + 1.0*im) == 2.0
@test CppTestFunctions.imag_part(2.0 + 1.0*im) == 1.0
@test CppTestFunctions.make_complex(Float32(3.0), Float32(4.0)) == 3.0 + 4.0*im
@test typeof(CppTestFunctions.make_complex(Float32(3.0), Float32(4.0))) == Complex{Float32}

@test CppTestFunctions.process_irrational(φ, 2) == 2*φ

# Performance tests
const test_size = Sys.ARCH == :armv7l ? 1000000 : 50000000
const numbers = rand(test_size)
output = zeros(test_size)

# Build a function to loop over the test array
function make_loop_function(name)
    fname = Symbol(:half_loop_,name,:!)
    inner_name = Symbol(:half_,name)
    @eval begin
        function $(fname)(n::Array{Float64,1}, out_arr::Array{Float64,1})
            test_length = length(n)
          for i in 1:test_length
                out_arr[i] = $(inner_name)(n[i])
          end
        end
    end
end

# Julia version
half_julia(d::Float64) = d*0.5

# C version
half_c(d::Float64) = ccall((:half_c, functions_lib_path), Cdouble, (Cdouble,), d)

# Bring C++ versions into scope
using CppHalfFunctions.half_d, CppHalfFunctions.half_lambda, CppHalfFunctions.half_loop_cpp!, CppHalfFunctions.half_loop_jlcall!, CppHalfFunctions.half_loop_cfunc!

@static if cxx_available
  # Cxx.jl version
  cxx"""
  double half_cxx(const double d)
  {
    return 0.5*d;
  }
  """
  half_cxxjl(d::Float64) = @cxx half_cxx(d)
end

# Make the looping functions
make_loop_function(:julia)
make_loop_function(:c)
make_loop_function(:d) # C++ with regular C++ function pointer
make_loop_function(:lambda) # C++ lambda, so using std::function
if cxx_available
  make_loop_function(:cxxjl) # Cxx.jl version
end

# test that a "half" function does what it should
function test_half_function(f)
  input = [2.]
  output = [0.]
  f(input, output)
  @test output[1] == 1.
end
test_half_function(half_loop_julia!)
test_half_function(half_loop_c!)
test_half_function(half_loop_d!)
test_half_function(half_loop_lambda!)
test_half_function(half_loop_cpp!)
if cxx_available
  test_half_function(half_loop_cxxjl!)
end

# Run timing tests
println("---- Half test timings ----")
println("Julia test:")
@time half_loop_julia!(numbers, output)
@time half_loop_julia!(numbers, output)
@time half_loop_julia!(numbers, output)

println("C test:")
@time half_loop_c!(numbers, output)
@time half_loop_c!(numbers, output)
@time half_loop_c!(numbers, output)

println("C++ test:")
@time half_loop_d!(numbers, output)
@time half_loop_d!(numbers, output)
@time half_loop_d!(numbers, output)

if cxx_available
  println("Cxx.jl test:")
  @time half_loop_cxxjl!(numbers, output)
  @time half_loop_cxxjl!(numbers, output)
  @time half_loop_cxxjl!(numbers, output)
end

println("C++ lambda test:")
@time half_loop_lambda!(numbers, output)
@time half_loop_lambda!(numbers, output)
@time half_loop_lambda!(numbers, output)

println("C++ test, loop in the C++ code:")
@time half_loop_cpp!(numbers, output)
@time half_loop_cpp!(numbers, output)
@time half_loop_cpp!(numbers, output)

println("cfunction in C++ loop")
half_cfunc = safe_cfunction(half_julia, Float64, (Float64,))
@time half_loop_cfunc!(numbers, output, half_cfunc)
@time half_loop_cfunc!(numbers, output, half_cfunc)
@time half_loop_cfunc!(numbers, output, half_cfunc)

const small_in = rand(test_size÷100)
small_out = zeros(test_size÷100)

println("jl_call inside C++ loop (array is 100 times smaller than other tests):")
@time half_loop_jlcall!(small_in, small_out)
@time half_loop_jlcall!(small_in, small_out)
@time half_loop_jlcall!(small_in, small_out)


text:
module TutorialsForDevelopers

files = [
 "arrays.jl"=>"Working with arrays",
]

end

max_stars_repo_name: gridap/TutorialsForDevelopers

max_issues_repo_name: gridap/TutorialsForDevelopers

max_forks_repo_name: gridap/TutorialsForDevelopers


text:
using Test
using Suppressor: @suppress_out

@testset "Example Directory Tests" begin

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/mapreduce.jl")

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_divide.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_mod.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_discrete.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_arith.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_condint.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/gcd_mod.jl"
    )

    @test_nowarn @suppress_out include(
        "$(@__DIR__)/../examples/test_ifelse.jl"
    )

end


text:
function no_whitening(
    dataset::DataSet{T, I})::WhiteningResult{T} where {T<:AbstractFloat, I<:Integer}

    datamean = zeros(T, dataset.P)

    transform_data!(dataset, Array{T, 2}(LinearAlgebra.I, dataset.P, dataset.P), datamean)
end

function cholesky_whitening!(
    dataset::DataSet{T, I})::WhiteningResult{T} where {T<:AbstractFloat, I<:Integer}

    datamean = zeros(T, dataset.P)

    for p in eachindex(datamean)
        datamean[p] = mean(view(dataset.data, p, :))
    end

    for n in 1:dataset.N
        buffer = view(dataset.data, :, n) - datamean
        setindex!(dataset.data, buffer, :, n)
    end

    #cov of a Transpose object yields the wrong (0.0) result in 1D, if all weights are equal (e.g. direct sampling)!
    covmatrix = cov(convert(typeof(dataset.data), transpose(dataset.data)), FrequencyWeights(dataset.weights), corrected=true)
    symcovmatrix = Symmetric(covmatrix)
    covmatrix_inv = inv(symcovmatrix)
    w = cholesky(covmatrix_inv).U
    wres = convert(Matrix{T}, w)

    transform_data!(dataset, wres, datamean)
end

function cholesky_partial_whitening!(
    dataset::DataSet{T, I})::WhiteningResult{T} where {T<:AbstractFloat, I<:Integer}

    datamean = zeros(T, dataset.P)
	avoid_whitening_dims = find_steep_edges(dataset)

    for p in eachindex(datamean)
        datamean[p] = mean(view(dataset.data, p, :))
    end

    for n in 1:dataset.N
        buffer = view(dataset.data, :, n) - datamean
        setindex!(dataset.data, buffer, :, n)
    end

    covmatrix = cov(convert(typeof(dataset.data), transpose(dataset.data)), FrequencyWeights(dataset.weights), corrected=true)

	for param in avoid_whitening_dims # replace non-diagonal elements with zeros to preserve rectangular bounds.
		covmatrix[1:param-1, param] .= 0.0
		covmatrix[param+1:end, param] .= 0.0
		covmatrix[param, param+1:end] .= 0.0
		covmatrix[param, 1:param-1] .= 0.0
	end

    symcovmatrix = Symmetric(covmatrix)
    covmatrix_inv = inv(symcovmatrix)
    w = cholesky(covmatrix_inv).U
    wres = convert(Matrix{T}, w)

    transform_data!(dataset, wres, datamean)
end

function find_steep_edges(dataset::DataSet{T, I})::Array{Int64,1} where {T<:AbstractFloat, I<:Integer}

    N_bins::Int64 = 100
    treshold::Float64 = 0.00
    avoid_whitening_dims::Array{Int64,1} = []

    for param in 1:dataset.P
        hist = fit(Histogram, dataset.data[param, :], FrequencyWeights(dataset.weights), nbins = N_bins, closed = :left)
        treshold = 0.01 * maximum(hist.weights)
		if ((hist.weights[1] < treshold) && (hist.weights[end] < treshold)) == false
            push!(avoid_whitening_dims, param)
        end
    end

    return avoid_whitening_dims
end

function statistical_whitening!(
    dataset::DataSet{T, I})::WhiteningResult{T} where {T<:AbstractFloat, I<:Integer}

    datamean = zeros(T, dataset.P)

    for p in eachindex(datamean)
        datamean[p] = mean(view(dataset.data, p, :))
    end

    dataset.data .-= datamean

    covmatrix = cov(transpose(dataset.data), FrequencyWeights(dataset.weights), corrected=true)

    E = eigfact(covmatrix).vectors
    w_d = transpose(E)
    D = Diagonal(inv(E) * covmatrix * E)
    w = inv(full(sqrt.(D))) * w_d
    wres = convert(Matrix{T}, w)

    transform_data!(dataset, wres, datamean)
end

function transform_data!(
    dataset::DataSet{T, I},
    W::Matrix{T},
    datamean::Vector{T},
    substractmean::Bool = false)::WhiteningResult{T} where {T<:AbstractFloat, I<:Integer}

    local determinant::T

    if substractmean
        dataset.data .-= datamean
    end

    if diag(W) == ones(T, dataset.P)
        determinant = 1.0
    else
        dataset.data = W * dataset.data

        determinant = abs(det(W))
    end

    maxP::T = maximum(dataset.logprob)
    suggTargetProb::T = 500 #exp(maxP - partialsort(dataset.logprob, floor(Int64, dataset.N * 0.2)))

    dataset.iswhitened = true

    @debug "Determinant:\t" * string(determinant)
    @debug "Suggested Target Probability Factor:\t" * string(suggTargetProb)

    WhiteningResult(determinant, suggTargetProb, W, datamean)
end


text:
import FinanceLib.Derivatives.Forwards as Fd
using Dates: Date
using FinanceLib: yearFrac

@testset "FinanceLib.Derivatives.Forwards                             " begin

  fw = Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50, benf = -1)
  @test fw.FT == 51.37827066066438
  @test Fd.getS0(fw) == 50.0
  @test Fd.valuation(fw, 52, 1/12) == 1.8766865113848255

  fw = Fd.Forward_S0(r = 0.03, T = 0.25, S0 = 50)
  @test fw.FT == 50.370853588886646
  @test Fd.getS0(fw) == 50.0

  T0 = Date(2015,10,5)
  T1 = Date(2016,1,8)
  TX = yearFrac(T0, T1)

  @test TX == 0.2600958247775496

  fw = Fd.Forward_S0(r = 0.03, T = TX, S0 = 50, benf = -1)
  # println(fw)

  dw = Fd.XForward_S0(r = 0.03, T0 = T0, T = T1, S0 = 50, benf = -1.0)
  # println(dw)

  @test fw.FT == dw.FT
  @test Fd.getS0(fw) == Fd.getS0(dw)

  TM = Date(2015, 12, 18)
  TZ = yearFrac(T0, TM)
  # println(TZ)

  @test Fd.valuation(fw, 55, TZ) == Fd.valuation(dw, 55, TM)

end

import FinanceLib.Derivatives.Options as Op

@testset "FinanceLib.Derivatives.Options                              " begin

  ct = Op.Call(2_000, 81.75)

  @test Op.valuation(ct, ST = 1_900) == 0.0
  @test Op.profit(ct, ST = 1_900) == -81.75

  @test Op.valuation(ct, ST = 2_100) == 100.0
  @test Op.profit(ct, ST = 2_100) == 18.25

  pt = Op.Put(2_000, 79.25)

  @test Op.valuation(pt, ST = 1_900) == 100.0
  @test Op.profit(pt, ST = 1_900) == 20.75

  @test Op.valuation(pt, ST = 2_100) == 0.0
  @test Op.profit(pt, ST = 2_100) == -79.25

  cv = Op.CoveredCall(98, Op.Call(105, 8))

  @test Op.valuation(cv, ST = 110) == 105
  @test Op.profit(cv, ST = 110) == 15

  @test Op.valuation(cv, ST = 88) == 88
  @test Op.profit(cv, ST = 88) == -2

  pt = Op.ProtectivePut(0.875, Op.Put(0.90, 0.075))

  @test Op.valuation(pt, ST = 0.96) == 0.96
  @test Op.profit(pt, ST = 0.96) ≈ 0.01

  @test Op.valuation(pt, ST = 0.75) == 0.90
  @test Op.profit(pt, ST = 0.75) ≈ -0.05

end


text:
# This script shows how to compute risk-adjusted linearizations using sparse Jacobian methods
# for calls to nlsolve
using RiskAdjustedLinearizations, LinearAlgebra, SparseArrays, BenchmarkTools

# Settings
define_functions = true
time_methods     = true
algorithm        = :relaxation
N_approx         = 10          # Number of periods ahead used for forward-difference equations

if define_functions
    include(joinpath(dirname(@__FILE__), "..", "nk_with_capital", "nk_with_capital.jl"))
end

# Set up
m_nk = NKCapital(; N_approx = N_approx) # create parameters
m = nk_capital(m_nk) # instantiate risk-adjusted linearization
zinit = deepcopy(m.z)
yinit = deepcopy(m.y)
Ψinit = deepcopy(m.Ψ)

# Compute steady state
solve!(m; algorithm = :deterministic, verbose = :none)
zdet = deepcopy(m.z)
ydet = deepcopy(m.y)
Ψdet = deepcopy(m.Ψ)

# Using sparsity pattern and matrix coloring to find deterministic steady state
sparsity_det, colorvec_det = compute_sparsity_pattern(m, :deterministic) # Create sparsity pattern/matrix coloring vector
update!(m, zinit, yinit, Ψinit) # need to re-initialize or else the system of equations will be solved at the initial guess
solve!(m; algorithm = :deterministic, ftol = 1e-10, # this call to solve! tries to infer the sparsity pattern
       sparse_jacobian = true, verbose = :none)     # by computing the Jacobian one time w/finite differences
update!(m, zinit, yinit, Ψinit)
solve!(m; algorithm = :deterministic, ftol = 1e-10,                              # need to add tighter tolerance to avoid
       sparse_jacobian = true, sparsity = sparsity_det, colorvec = colorvec_det, # a LAPACK exception when computing the
       verbose = :none)                                                          # Schur decomposition for Ψ

# Using sparsity pattern and matrix coloring to find stochastic steady state via relaxation
# Note that the syntax is essentially the same as the deterministic case, except that
# compute_sparsity_pattern's second argument is now :relaxation
sparsity_rel, colorvec_rel = compute_sparsity_pattern(m, :relaxation)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :relaxation, autodiff = autodiff_method,
       sparse_jacobian = true, verbose = :none)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :relaxation, autodiff = autodiff_method,
       sparse_jacobian = true, sparsity = sparsity_rel, colorvec = colorvec_rel,
       verbose = :none)

# Using sparsity pattern and matrix coloring to find stochastic steady state via homotopy
sparsity_hom, colorvec_hom = compute_sparsity_pattern(m, :homotopy)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :homotopy, autodiff = autodiff_method,
       sparse_jacobian = true, verbose = :none)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :homotopy, autodiff = autodiff_method,
       sparse_jacobian = true, sparsity = sparsity_hom, colorvec = colorvec_hom,
       verbose = :none)

# Using Jacobian cache to find deterministic steady state
# Like compute_sparsity_pattern, the user only needs to select the
# algorithm for which the Jacobian cache will be used.
jac_cache_det = preallocate_jac_cache(m, :deterministic)
update!(m, zinit, yinit, Ψinit)
solve!(m; algorithm = :deterministic, autodiff = autodiff_method,
       ftol = 1e-10, sparse_jacobian = true, jac_cache = jac_cache_det,
       verbose = :none)

# Using Jacobian cache to find stochastic steady state via relaxation
jac_cache_rel = preallocate_jac_cache(m, :relaxation)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :relaxation, ftol = 1e-10,
       sparse_jacobian = true, jac_cache = jac_cache_rel,
       verbose = :none)

# Using Jacobian cache to find stochastic steady state via homotopy
jac_cache_hom = preallocate_jac_cache(m, :homotopy)
update!(m, zdet, ydet, Ψdet)
solve!(m; algorithm = :homotopy, autodiff = autodiff_method,
       sparse_jacobian = true, jac_cache = jac_cache_hom,
       verbose = :none)

if time_methods
    println("Deterministic steady state with dense Jacobian")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = :deterministic, verbose = :none)
    end

    println("Deterministic steady state with sparsity pattern of Jacobian and matrix coloring vector")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = :deterministic, ftol = 1e-10,
               sparse_jacobian = true, sparsity = sparsity_det, colorvec = colorvec_det,
               verbose = :none)
    end

    println("Deterministic steady state with sparsity pattern of Jacobian, matrix coloring vector, and caching")
    @btime begin
        update!(m, zinit, yinit, Ψinit)
        solve!(m; algorithm = :deterministic, autodiff = autodiff_method,
               ftol = 1e-10, sparse_jacobian = true, jac_cache = jac_cache_det,
               verbose = :none)
    end

    println("Relaxation with dense Jacobian")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :relaxation, verbose = :none)
    end

    println("Relaxation with sparsity pattern of Jacobian and matrix coloring vector")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :relaxation, autodiff = autodiff_method,
               sparse_jacobian = true, sparsity = sparsity_rel, colorvec = colorvec_rel,
               verbose = :none)
    end

    println("Relaxation with sparsity pattern of Jacobian, matrix coloring vector, and caching")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :relaxation, autodiff = autodiff_method,
               sparse_jacobian = true, jac_cache = jac_cache_rel,
               verbose = :none)
    end

    println("Homotopy with dense Jacobian")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :homotopy, verbose = :none)
    end

    println("Homotopy with sparsity pattern of Jacobian and matrix coloring vector")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :homotopy,
               sparse_jacobian = true, sparsity = sparsity_hom, colorvec = colorvec_hom,
               verbose = :none)
    end

    println("Homotopy with sparsity pattern of Jacobian, matrix coloring vector, and caching")
    @btime begin
        update!(m, zdet, ydet, Ψdet)
        solve!(m; algorithm = :homotopy, autodiff = autodiff_method,
               sparse_jacobian = true, jac_cache = jac_cache_hom,
               verbose = :none)
    end
end


text:
const liblapack = Base.liblapack_name
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, chkstride1, Givens

function dorcsd2by1!(X11::AbstractMatrix{Float64}, X21::AbstractMatrix{Float64})
    # chkstride1(X11, X21)
    jobu1 = 'Y'
    jobu2 = 'Y'
    jobv1t= 'Y'
    p, q = size(X11)
    m = p + size(X21, 1)
    ldx11 = max(1, p)
    ldx21 = max(1, m-p)
    theta = Vector{Float64}(undef, min(p, m-p, q, m-q))
    U1 = similar(X11, Float64, ldx11, p)
    U2 = similar(X21, Float64, ldx21, m-p)
    ldu1 = ldx11
    ldu2 = ldx21
    V1T = Matrix{Float64}(undef, q, q)
    ldv1t = max(1, q)

    work = Vector{Float64}(undef, 1)
    lwork = BlasInt(-1)
    iwork = Vector{BlasInt}(undef, m-min(p, m-p, q, m-q))
    info  = Ref{BlasInt}()

    for i in 1:2  # first call returns lwork as work[1]
        ccall((@blasfunc(dorcsd2by1_), liblapack), Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
        Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
        Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
        Ptr{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
        Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}
        ),
        jobu1, jobu2, jobv1t,
        m, p, q,
        X11, ldx11, X21, ldx21,
        theta, U1, ldu1, U2, ldu2, V1T, ldv1t,
        work, lwork, iwork, info)
        # chklapackerror(info[])
        if i == 1
            lwork = BlasInt(real(work[1]))
            resize!(work, lwork)
        end
    end
    # return X11, X21, theta, U1, U2, V1T
    return theta, U1, U2, V1T
end