text:
# This function finds a new node to start making a tree from
function findNewWork(successor, numberOfThreads)

    # Search from (threadID / Number of threads) * successors to the end of successors
    for i in (round(Int64, ((Threads.threadid() - 1) / numberOfThreads) * length(successor)) + 1) : length(successor)
        if successor[i][] == 0
            return i
        end
    end

    # Search from 1 to (threadID / Number of threads) * successor
    for i in 1 : round(Int64, ((Threads.threadid() - 1) / numberOfThreads) * length(successor)) + 1
        if successor[i][] == 0
            return i
        end
    end

    return 0
end

# Function intended to be called from different threads to make a subset of the MST of graph g
function partialPrims(g::AbstractGraph{U}, successorArray, numberOfThreads) where {U}

    graphWidth = LightGraphs.nv(g)

    # Grab the first root
    root = findNewWork(successorArray, numberOfThreads);

    # Create an array of all the edges in the msf
    edgeIterator = 1;
    arrayOfEdges = Array{Edge{Int64}, 1}(undef, graphWidth - 1)
    for edge in 1 : length(arrayOfEdges)
        arrayOfEdges[edge] = Edge(0, 0)
    end

    # While there is new work to do
    while (root != 0)

        # Add this root as its own parent, if that was already done find a new node
        ret = Threads.atomic_cas!(successorArray[root], 0, root);
        if (ret != 0)
            root = findNewWork(successorArray, numberOfThreads);
            continue
        end

        # Priority queue of edges
        edgeQueue = DataStructures.PriorityQueue{Tuple{Int64,Int64},Float64}()

        # Create a set which represents the current tree
        subtree = Set{Int64}()
        push!(subtree, root)

        # Add edges from root
        for neighbor in LightGraphs.outneighbors(g, root)
            if(!(in(neighbor, subtree)))
                DataStructures.enqueue!(edgeQueue, (root, neighbor), g.weights[root, neighbor])
            end
        end

        # Go through the queue of all edges adding small ones
        while (true)

            nextEdge = Tuple{Int64, Int64}

            somethingFound = false

            # Find the next edge which connects subtree to a node outside of it
            while((DataStructures.length(edgeQueue)) > 0)
                nextEdge = DataStructures.dequeue!(edgeQueue)
                if((in(nextEdge[1], subtree)) && (!(in(nextEdge[2], subtree))))
                    somethingFound = true
                    break;
                end
            end

            # If there aren't any edges left break loop
            if(((DataStructures.length(edgeQueue)) == 0) && !somethingFound)
                break;
            end

            # Add the edge to the array of MSF
            arrayOfEdges[edgeIterator] = Edge(nextEdge[1], nextEdge[2])
            edgeIterator += 1

            # If our destination doesn't have a successor set out roots successor to it successor
            parent = successorArray[root][]
            ret = Threads.atomic_cas!(successorArray[nextEdge[2]], 0, successorArray[parent][]);

            # If we do have a successor already
            if(ret != 0)
                parent = successorArray[nextEdge[2]][]
                Threads.atomic_xchg!(successorArray[root], successorArray[parent][])
                break;

            # Add edge to list and add neighbours to explorables
            else
                push!(subtree, nextEdge[2])
                for neighbor in LightGraphs.outneighbors(g, nextEdge[2])
                    if(!(in(neighbor, subtree)))
                        DataStructures.enqueue!(edgeQueue, (nextEdge[2], neighbor), g.weights[nextEdge[2], neighbor])
                    end
                end
            end
        end
        root = findNewWork(successorArray, numberOfThreads);
    end
    return arrayOfEdges
end

# Trace through the successor array recursively to find the true root
function trueRoot(successorArray::AbstractArray, root::Int64, prev::Int64)
    if(successorArray[root][] == prev)
        return prev
    elseif(successorArray[root][] != root)
        return trueRoot(successorArray, successorArray[root][], root)
    else
        return root
    end
end

# Step to connect all subsets created during the partial prims step
function unification(g::AbstractGraph, subtrees)
    subtree = subtrees[1]
    edgeArray = []

    for i in 1 : length(subtrees)

        # Find lowest edge that links to another subtree
        lowestEdge = Edge(0, 0)
        lowestWeight = 100000.0
        for vertex in subtree
            for neighbor in LightGraphs.outneighbors(g, vertex)
                if !(in(neighbor, subtree))
                    if (g.weights[vertex, neighbor] < lowestWeight)
                        lowestEdge = Edge(vertex, neighbor)
                        lowestWeight = g.weights[vertex, neighbor]
                    end
                end
            end
        end

        # Find the subtree the destination belongs too
        subTreeIndex = 1
        for subtreeCheck in 1 : length(subtrees)
            if (in(lowestEdge.dst, subtrees[subtreeCheck]))
                subTreeIndex = subtreeCheck
            end
        end

        # add that subtree to our one
        subtree = union(subtree, subtrees[subTreeIndex])
        push!(edgeArray, lowestEdge)
    end
    return edgeArray
end

# Function to run prims on multiple processors
function prims(g::AbstractGraph, numberOfThreads::Int64)

    # Create a successor array shared between threads, used to unify
    graphWidth = LightGraphs.nv(g)
    successorArray = Vector(undef, graphWidth)
    for i in 1 : (graphWidth)
        successorArray[i] = Threads.Atomic{Int}(0)
    end

    totalTree = Array{Array{Edge{Int64}, 1}, 1}(undef, numberOfThreads + 1)

    # Run partial prims on multiple threads
    Threads.@threads for i in 1 : numberOfThreads
        subTree = partialPrims(g, successorArray, numberOfThreads)
        totalTree[i] = subTree
    end

    theSets = Array{Array{Int64}, 1}()
    setIterator = 1
    rootDict = Dict()
    doneRoots = Set()

    # Find which unconnected subtree each vertex belongs to
    for i in 1 : length(successorArray)
        trueRootValue = trueRoot(successorArray, successorArray[i][], 0)
        if !(in(trueRootValue, doneRoots))
            push!(doneRoots, trueRootValue)

            push!(theSets, [])
            push!(theSets[setIterator], i)

            rootDict[trueRootValue] = setIterator

            setIterator += 1
        else
            push!(theSets[rootDict[trueRootValue]], i)
        end
    end

    # Run through the subtrees to find edges needed to connect islands
    unificationEdges = unification(g, theSets)
    totalTree[numberOfThreads + 1] = unificationEdges

    # Add together all valid edges to construct overall graph
    edges = Array{Edge{Int64}, 1}(undef, graphWidth - 1)
    for edge in 1 : length(edges)
        edges[edge] = Edge(0, 0)
    end
    edgeIterator = 0
    for i in 1 : length(totalTree)
        for edge in totalTree[i]
            if(!(in(edge, edges)) && !(in(Edge(edge.dst, edge.src), edges)) && edge.src != 0)
                edgeIterator += 1
                edges[edgeIterator] = edge
            end
        end
    end

    return edges
end


text:
module TestHelpers

type FakeTerminal <: Base.Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    FakeTerminal(stdin,stdout,stderr,hascolor=true) =
        new(stdin,stdout,stderr,hascolor,false)
end

Base.Terminals.hascolor(t::FakeTerminal) = t.hascolor
Base.Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
Base.Terminals.size(t::FakeTerminal) = (24, 80)

end

ext: jl

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
import Base.deepcopy_internal
deepcopy_internal(x::Module,stackdict::IdDict) = x

# Model Types.  The template args are required for inference for cache/perturbation solutions
struct PerturbationModel{MaxOrder,N_y,N_x,N_ϵ,N_z,N_p,HasΩ,T1,T2}
    mod::Module

    # Could extract from type, but here for simplicity
    max_order::Int64
    n_y::Int64
    n_x::Int64
    n_p::Int64
    n_ϵ::Int64
    n_z::Int64

    # Variations
    has_Ω::Bool
    η::T1
    Q::T2
end

# Construct from a module.  Inherently type unstable, so use function barrier from return type
function PerturbationModel(mod)
    return PerturbationModel{mod.max_order,mod.n_y,mod.n_x,mod.n_ϵ,mod.n_z,mod.n_p,
                             mod.has_Ω,typeof(mod.η),typeof(mod.Q)}(mod, mod.max_order,
                                                                    mod.n_y, mod.n_x,
                                                                    mod.n_p, mod.n_ϵ,
                                                                    mod.n_z, mod.has_Ω,
                                                                    mod.η, mod.Q)
end

# TODO: Add in latex stuff for the mod.H_latex,
function Base.show(io::IO, ::MIME"text/plain", m::PerturbationModel) where {T}
    return print(io,
                 "Perturbation Model: n_y = $(m.n_y), n_x = $(m.n_x), n_p = $(m.n_p), n_ϵ = $(m.n_ϵ), n_z = $(m.n_z)\n y = $(m.mod.y_symbols) \n x = $(m.mod.x_symbols) \n p = $(m.mod.p_symbols)")
end

# Buffers for the solvers to reduce allocations
# General rule for cache vs. buffers
# 1. If something should be used in multiple parts of the algorithm, put it in the cache
# 2. Otherwise, use the buffers, which you can "trash" with inplace operations as required
# 3. The cache should never be modified after it has been filled in a given sequence of events, buffers can be

# Note: When fully generic, this pattern has an issue when the exact number fields meets the number of arguments to the convenience constructor.  Forcing the first type bypasses the issue
Base.@kwdef struct FirstOrderSolverBuffers{ComplexMatrixType<:AbstractMatrix,RealMatrixType,
                                           UpperTriangularMatrixType}
    A::ComplexMatrixType
    B::ComplexMatrixType
    Z::RealMatrixType  # real version, transposed relative to the schur
    Z_ll::RealMatrixType
    S_bb::UpperTriangularMatrixType
    T_bb::UpperTriangularMatrixType
end
function FirstOrderSolverBuffers(n_y, n_x, n_p_d, n_ϵ, n_z)
    n = n_x + n_y
    return FirstOrderSolverBuffers(; A = zeros(Complex{Float64}, n, n),
                                   B = zeros(Complex{Float64}, n, n), Z = zeros(n, n),
                                   Z_ll = zeros(n_y, n_y),
                                   S_bb = UpperTriangular(zeros(n_x, n_x)),
                                   T_bb = UpperTriangular(zeros(n_x, n_x)))
end
Base.@kwdef struct FirstOrderDerivativeSolverBuffers{RealMatrixType}
    R::RealMatrixType
    A::RealMatrixType
    C::RealMatrixType
    D::RealMatrixType
    E::RealMatrixType
    dH::RealMatrixType
    bar::RealMatrixType
end
function FirstOrderDerivativeSolverBuffers(n_y, n_x, n_p_d, n_ϵ, n_z)
    n = n_x + n_y
    return FirstOrderDerivativeSolverBuffers(; R = zeros(2 * n, n_x), A = zeros(n, n),
                                             C = zeros(n, n), D = zeros(n_x, n_x),
                                             dH = zeros(n, 2n), bar = zeros(2n, 1),
                                             E = zeros(n, n_x))
end
Base.@kwdef struct SecondOrderSolverBuffers{RealMatrixType}
    A::RealMatrixType
    C::RealMatrixType
    B::RealMatrixType
    D::RealMatrixType
    E::RealMatrixType
    R::RealMatrixType
    A_σ::RealMatrixType
    R_σ::RealMatrixType
end
function SecondOrderSolverBuffers(n_y, n_x, n_p_d, n_ϵ, n_z)
    n = n_x + n_y
    return SecondOrderSolverBuffers(; A = zeros(n, n), B = zeros(n_x^2, n_x^2),
                                    C = zeros(n, n), D = zeros(n_x^2, n_x^2),
                                    E = zeros(n, n_x^2), R = zeros(2 * n, n_x),
                                    A_σ = zeros(n, n), R_σ = zeros(2 * n, n_x))
end
Base.@kwdef struct SecondOrderDerivativeSolverBuffers{RealMatrixType,VectorOfMatrixType,VectorOfVectorOfMatrixType}
    A::RealMatrixType
    B::RealMatrixType
    C::RealMatrixType
    D::RealMatrixType
    E::RealMatrixType
    R::RealMatrixType
    dH::RealMatrixType
    dΨ::VectorOfMatrixType
    gh_stack::RealMatrixType
    g_xx_flat::RealMatrixType
    Ψ_x_sum::VectorOfVectorOfMatrixType
    Ψ_y_sum::VectorOfVectorOfMatrixType
    bar::RealMatrixType
    kron_h_x::RealMatrixType
    R_p::RealMatrixType
    A_σ::RealMatrixType
    R_σ::RealMatrixType
end

function SecondOrderDerivativeSolverBuffers(n_y, n_x, n_p_d, n_ϵ, n_z)
    n = n_x + n_y
    return SecondOrderDerivativeSolverBuffers(; A = zeros(n, n), B = zeros(n_x^2, n_x^2),
                                              C = zeros(n, n), D = zeros(n_x^2, n_x^2),
                                              E = zeros(n, n_x^2), R = zeros(2n, n_x),
                                              dH = zeros(n, 2n), dΨ = [zeros(2n, 2n) for _ in 1:n],
                                              gh_stack = zeros(n, n_x^2), g_xx_flat = zeros(n_y, n_x^2),
                                              Ψ_x_sum =  [[zeros(2n, 2n) for _ in 1:n] for _ in 1:n_x],
                                              Ψ_y_sum =  [[zeros(2n, 2n) for _ in 1:n] for _ in 1:n_y],
                                              bar = zeros(2n, 1), kron_h_x = zeros(n_x^2, n_x^2),
                                              R_p = zeros(2n, n_x),
                                              A_σ = zeros(n, n), R_σ = zeros(2 * n, n_x),
                                              )
end

# The cache if for both 1st and 2nd order
# Constructors set values to nothing as appropriate
abstract type AbstractSolverCache{Order} end
Base.@kwdef struct SolverCache{Order,MatrixType,MatrixType2,MatrixType3,MatrixType4,
                               MatrixType5,VectorType,VectorType2,VectorOfVectorType,
                               VectorOfMatrixType,VectorOfMatrixType2,VectorOfMatrixType3,
                               VectorOrNothingType,VectorOfVectorOrNothingType,
                               MatrixScalingOrNothingType,SymmetricMatrixType,
                               SymmetricVectorOfMatrixType,VectorOfVectorOfMatrixType,
                               ThreeTensorType,CholeskyType,ChangeVarianceType,
                               VectorOfThreeTensorType,FirstOrderSolverBuffersType,
                               FirstOrderSolverDerivativesBuffersType,
                               SecondOrderSolverBuffersOrNothingType,
                               SecondOrderSolverDerivativesBuffersOrNothingType} <:
                   AbstractSolverCache{Order}
    order::Val{Order}  # allows inference in construction
    p_d_symbols::Vector{Symbol}
    H::VectorType
    H_yp::MatrixType
    H_y::MatrixType
    H_xp::MatrixType
    H_x::MatrixType
    H_yp_p::VectorOfMatrixType
    H_y_p::VectorOfMatrixType
    H_xp_p::VectorOfMatrixType
    H_x_p::VectorOfMatrixType
    H_p::VectorOfVectorType
    Γ::MatrixType2
    Γ_p::VectorOfMatrixType2
    Σ::SymmetricMatrixType
    Σ_p::SymmetricVectorOfMatrixType
    Ω::VectorOrNothingType
    Ω_p::VectorOfVectorOrNothingType
    Ψ::VectorOfMatrixType

    # Used in solution
    x::VectorType
    y::VectorType
    y_p::VectorOfVectorType
    x_p::VectorOfVectorType
    g_x::MatrixType4
    h_x::MatrixType4
    g_x_p::VectorOfMatrixType3
    h_x_p::VectorOfMatrixType3
    B::MatrixType2
    B_p::VectorOfMatrixType2
    Q::MatrixScalingOrNothingType
    η::MatrixType3
    A_1_p::VectorOfMatrixType3
    C_1::MatrixType4
    C_1_p::VectorOfMatrixType3
    V::CholeskyType
    V_p::ChangeVarianceType
    η_Σ_sq::SymmetricMatrixType

    # Additional for 2nd order
    Ψ_p::VectorOfVectorOfMatrixType
    Ψ_yp::VectorOfVectorOfMatrixType
    Ψ_y::VectorOfVectorOfMatrixType
    Ψ_xp::VectorOfVectorOfMatrixType
    Ψ_x::VectorOfVectorOfMatrixType
    g_xx::ThreeTensorType
    h_xx::ThreeTensorType
    g_σσ::VectorType2
    h_σσ::VectorType2
    g_xx_p::VectorOfThreeTensorType
    h_xx_p::VectorOfThreeTensorType
    g_σσ_p::MatrixType5
    h_σσ_p::MatrixType5

    # Additional for solution type 2nd order
    A_0_p::MatrixType5
    A_2_p::VectorOfThreeTensorType
    C_0::VectorType2
    C_2::ThreeTensorType
    C_0_p::MatrixType5
    C_2_p::VectorOfThreeTensorType

    # Buffers for additional calculations
    first_order_solver_buffer::FirstOrderSolverBuffersType
    first_order_solver_p_buffer::FirstOrderSolverDerivativesBuffersType
    second_order_solver_buffer::SecondOrderSolverBuffersOrNothingType
    second_order_solver_p_buffer::SecondOrderSolverDerivativesBuffersOrNothingType
    I_x::MatrixType  #dense identity matrices
    I_x_2::MatrixType
    zeros_x_x::MatrixType
    zeros_y_x::MatrixType
end

# The Val(2), etc. for the order required for inference to function
# Note that the n_p_d is the number of differentiated parameters to allocate for
# Extracts the field names from the order of the p_d argument (which could be dict, named tuple, etc.)
function SolverCache(m::PerturbationModel{MaxOrder,N_y,N_x,N_ϵ,N_z,N_p,HasΩ,T1,T2},
                     ::Val{Order},
                     p_d) where {Order,MaxOrder,N_y,N_x,N_ϵ,N_z,N_p,HasΩ,T1,T2}
    n_p_d = length(p_d)
    p_d_symbols = collect(Symbol.(keys(p_d)))
    return SolverCache(; order = Val(Order), p_d_symbols, H = zeros(N_x + N_y),
                       H_yp = zeros(N_x + N_y, N_y), H_y = zeros(N_x + N_y, N_y),
                       H_xp = zeros(N_x + N_y, N_x), H_x = zeros(N_x + N_y, N_x),
                       Γ = zeros(N_ϵ, N_ϵ), Ω = !HasΩ ? nothing : zeros(N_z),
                       Ψ = [zeros(2(N_x + N_y), 2(N_x + N_y)) for i in 1:(N_x + N_y)],
                       H_p = [zeros(N_x + N_y) for i in 1:n_p_d],
                       H_yp_p = [zeros(N_x + N_y, N_y) for i in 1:n_p_d],
                       H_y_p = [zeros(N_x + N_y, N_y) for i in 1:n_p_d],
                       H_xp_p = [zeros(N_x + N_y, N_x) for i in 1:n_p_d],
                       H_x_p = [zeros(N_x + N_y, N_x) for i in 1:n_p_d],
                       Γ_p = [zeros(N_ϵ, N_ϵ) for i in 1:n_p_d],
                       Ω_p = !HasΩ ? nothing : [zeros(N_z) for i in 1:n_p_d],
                       x = zeros(N_x), y = zeros(N_y), y_p = [zeros(N_y) for i in 1:n_p_d],
                       x_p = [zeros(N_x) for i in 1:n_p_d], g_x = zeros(N_y, N_x),
                       h_x = zeros(N_x, N_x), g_x_p = [zeros(N_y, N_x) for _ in 1:n_p_d],
                       h_x_p = [zeros(N_x, N_x) for _ in 1:n_p_d],
                       Σ = Symmetric(zeros(N_ϵ, N_ϵ)), η_Σ_sq = Symmetric(zeros(N_x, N_x)),
                       Σ_p = [Symmetric(zeros(N_ϵ, N_ϵ)) for _ in 1:n_p_d], m.Q, m.η,
                       B = zeros(N_x, N_ϵ), B_p = [zeros(N_x, N_ϵ) for _ in 1:n_p_d],
                       C_1 = zeros(N_z, N_x), C_1_p = [zeros(N_z, N_x) for _ in 1:n_p_d],
                       A_1_p = [zeros(N_x, N_x) for _ in 1:n_p_d],
                       V = cholesky(Array(I(N_x))),
                       V_p = [zeros(N_x, N_x) for _ in 1:n_p_d],

                       # Stuff for 2nd order
                       Ψ_p = (Order == 1) ? nothing :
                             [[zeros(2 * (N_x + N_y), 2 * (N_x + N_y))
                               for _ in 1:(N_x + N_y)] for _ in 1:n_p_d],
                       Ψ_yp = (Order == 1) ? nothing :
                              [[zeros(2 * (N_x + N_y), 2 * (N_x + N_y))
                                for _ in 1:(N_x + N_y)] for _ in 1:N_y],
                       Ψ_y = (Order == 1) ? nothing :
                             [[zeros(2 * (N_x + N_y), 2 * (N_x + N_y))
                               for _ in 1:(N_x + N_y)] for _ in 1:N_y],
                       Ψ_xp = (Order == 1) ? nothing :
                              [[zeros(2 * (N_x + N_y), 2 * (N_x + N_y))
                                for _ in 1:(N_x + N_y)] for _ in 1:N_x],
                       Ψ_x = (Order == 1) ? nothing :
                             [[zeros(2 * (N_x + N_y), 2 * (N_x + N_y))
                               for _ in 1:(N_x + N_y)] for _ in 1:N_x],
                       g_xx = (Order == 1) ? nothing : zeros(N_y, N_x, N_x),
                       h_xx = (Order == 1) ? nothing : zeros(N_x, N_x, N_x),
                       g_σσ = (Order == 1) ? nothing : zeros(N_y),
                       h_σσ = (Order == 1) ? nothing : zeros(N_x),
                       g_xx_p = (Order == 1) ? nothing :
                                [zeros(N_y, N_x, N_x) for _ in 1:n_p_d],
                       h_xx_p = (Order == 1) ? nothing :
                                [zeros(N_x, N_x, N_x) for _ in 1:n_p_d],
                       g_σσ_p = (Order == 1) ? nothing : zeros(N_y, n_p_d),
                       h_σσ_p = (Order == 1) ? nothing : zeros(N_x, n_p_d),
                       A_0_p = (Order == 1) ? nothing : zeros(N_x, n_p_d),
                       A_2_p = (Order == 1) ? nothing :
                               [zeros(N_x, N_x, N_x) for _ in 1:n_p_d],
                       C_0 = (Order == 1) ? nothing : zeros(N_z),
                       C_0_p = (Order == 1) ? nothing : zeros(N_z, n_p_d),
                       C_2 = (Order == 1) ? nothing : zeros(N_z, N_x, N_x),
                       C_2_p = (Order == 1) ? nothing :
                               [zeros(N_z, N_x, N_x) for _ in 1:n_p_d],

                       # buffers for algorithms
                       first_order_solver_buffer = FirstOrderSolverBuffers(N_y, N_x, n_p_d,
                                                                           N_ϵ, N_z),
                       first_order_solver_p_buffer = FirstOrderDerivativeSolverBuffers(N_y,
                                                                                N_x,
                                                                                n_p_d,
                                                                                N_ϵ,
                                                                                N_z),
                       second_order_solver_buffer = (Order == 1) ? nothing :
                                                    SecondOrderSolverBuffers(N_y, N_x,
                                                                             n_p_d, N_ϵ,
                                                                             N_z),
                       second_order_solver_p_buffer = (Order == 1) ? nothing :
                                                      SecondOrderDerivativeSolverBuffers(N_y,
                                                                                N_x,
                                                                                n_p_d,
                                                                                N_ϵ,
                                                                                N_z),
                       I_x = Matrix{Float64}(I(N_x)), I_x_2 = Matrix{Float64}(I(N_x^2)),
                       zeros_x_x = zeros(N_x, N_x), zeros_y_x = zeros(N_y, N_x))
end
Base.@kwdef struct PerturbationSolverSettings{T1,T2,T3,T4,T5,T6}
    print_level::Int64 = 1  # 0 is no output at all
    ϵ_BK::Float64 = 1e-6 # For checking Blanchard-Kahn condition
    nlsolve_method::Symbol = :trust_region
    nlsolve_iterations::Int64 = 1000
    nlsolve_show_trace::Bool = false
    nlsolve_ftol::Float64 = 1e-8
    use_solution_cache::Bool = true
    evaluate_functions_callback::T1 = nothing
    calculate_steady_state_callback::T2 = nothing
    solve_first_order_callback::T3 = nothing
    solve_first_order_p_callback::T4 = nothing
    solve_second_order_callback::T5 = nothing
    solve_second_order_p_callback::T6 = nothing
end

function nlsolve_options(s::PerturbationSolverSettings)
    return (method = s.nlsolve_method, iterations = s.nlsolve_iterations,
            show_trace = s.nlsolve_show_trace, ftol = s.nlsolve_ftol)
end

# State Space types
abstract type AbstractPerturbationSolution end
abstract type AbstractFirstOrderPerturbationSolution <: AbstractPerturbationSolution end

# For this, all are dense due to schur decomposition
# All are dense due to schur decomposition
Base.@kwdef struct FirstOrderPerturbationSolution{T1<:AbstractVector,T2<:AbstractVector,
                                                  T3<:AbstractMatrix,T4<:AbstractMatrix,
                                                  T5<:AbstractMatrix,
                                                  T6<:Union{Nothing,Distribution},
                                                  T7<:Union{Nothing,AbstractMatrix,
                                                            UniformScaling},
                                                  T8<:AbstractMatrix,T9<:AbstractMatrix,
                                                  T10<:Distribution,T11<:AbstractMatrix} <:
                   AbstractFirstOrderPerturbationSolution
    retcode::Symbol
    x_symbols::Vector{Symbol}
    y_symbols::Vector{Symbol}
    p_symbols::Vector{Symbol}
    p_d_symbols::Vector{Symbol}
    u_symbols::Vector{Symbol}
    # TODO: differentiated parameter ordering?
    n_y::Int64
    n_x::Int64
    n_p::Int64
    n_ϵ::Int64
    n_z::Int64
    y::T1
    x::T2
    g_x::T3
    A::T4
    B::T5
    D::T6  # current a matrix or nothing, later could make more general
    Q::T7  # can be nothing
    η::T8
    C::T9  # i.e. Q * g_x
    x_ergodic::T10
    Γ::T11
end

maybe_diagonal(x::AbstractVector) = TuringDiagMvNormal(zero(x), x)
maybe_diagonal(x) = x # otherwise, just return raw.  e.g. nothing

function FirstOrderPerturbationSolution(retcode, m::PerturbationModel, c::SolverCache)
    return FirstOrderPerturbationSolution(; retcode, m.mod.x_symbols, m.mod.y_symbols,
                                          m.mod.u_symbols, m.mod.p_symbols, c.p_d_symbols,
                                          m.n_x, m.n_y, m.n_p, m.n_ϵ, m.n_z, c.Q, c.η, c.y,
                                          c.x, c.B, D = maybe_diagonal(c.Ω), c.g_x,
                                          A = c.h_x, C = c.C_1,
                                          x_ergodic = TuringDenseMvNormal(zeros(m.n_x),
                                                                          c.V), c.Γ)
end

Base.@kwdef struct SecondOrderPerturbationSolution{T1<:AbstractVector,T2<:AbstractVector,
                                                   T3<:AbstractMatrix,T4<:AbstractMatrix,
                                                   T5<:AbstractMatrix,
                                                   T6<:Union{Nothing,Distribution},
                                                   T7<:Union{Nothing,AbstractMatrix,
                                                             UniformScaling},
                                                   T8<:AbstractMatrix,T9<:AbstractMatrix,
                                                   T10,T11<:AbstractArray,T12,
                                                   T13<:AbstractVector,T14<:AbstractMatrix,
                                                   T15<:AbstractVector,
                                                   T16<:AbstractArray} <:
                   AbstractPerturbationSolution
    retcode::Symbol
    x_symbols::Vector{Symbol}
    y_symbols::Vector{Symbol}
    p_symbols::Vector{Symbol}
    p_d_symbols::Vector{Symbol}
    u_symbols::Vector{Symbol}
    n_y::Int64
    n_x::Int64
    n_p::Int64
    n_ϵ::Int64
    n_z::Int64
    y::T1
    x::T2
    g_x::T3
    B::T5
    D::T6
    Q::T7  # can be nothing
    η::T8
    Γ::T9

    g_xx::T10
    g_σσ::T12
    A_0::T13
    A_1::T4
    A_2::T11

    C_1::T14
    C_0::T15
    C_2::T16
end

function SecondOrderPerturbationSolution(retcode, m::PerturbationModel, c::SolverCache)
    return SecondOrderPerturbationSolution(; retcode, m.mod.x_symbols, m.mod.y_symbols,
                                           m.mod.u_symbols, m.mod.p_symbols, c.p_d_symbols,
                                           m.n_x, m.n_y, m.n_p, m.n_ϵ, m.n_z, c.Q, c.η, c.y,
                                           c.x, c.B, D = maybe_diagonal(c.Ω), c.Γ, c.g_x,
                                           A_1 = c.h_x, c.g_xx, A_2 = 0.5 * c.h_xx, c.g_σσ,
                                           A_0 = 0.5 * c.h_σσ, c.C_1, c.C_0, c.C_2)
end


text:
default_lambda{T <: Float}(x::DenseMatrix{T}, y::DenseVector{T}) = sqrt( log(size(x,2)) / length(y)) :: T
default_lambda{T <: Float}(x::DenseMatrix{T}, y::DenseVector{T}, q::Int) = (ones(T, q) * sqrt( log(size(x,2)) / length(y))) :: Vector{T}
default_lambda(T::Type, n::Int, p::Int, q::Int) = (ones(T, q)) * sqrt( log(p / n)) :: Vector{T}
default_lambda(n::Int, p::Int, q::Int) = default_lambda(Float64, n, p, q)

"""
    fit_logistic(x, y, λ)

Refit a regularized logistic model using Newton's method.

Arguments:

- `x` is the `n` x `p` statistical model.
- `y` is the `n`-vector of binary responses.
- `λ` is the regularization parameter

Optional Arguments:

- `tol` is the convergence tolerance. Defaults to `1e-8`.
- `max_iter` sets the maximum number of backtracking steps. Defaults to `50`.
- `quiet` controls output. Defaults to `true` (no output).

Output:

- a vector of refit coefficients for `β`
- the number of backtracking steps taken.
"""
function fit_logistic!{T <: Float, V <: DenseVector}(
    v        :: IHTLogVariables{T, V},
    y        :: V,
    λ        :: T;
    tol      :: T    = convert(T, 1e-8),
    max_iter :: Int  = 50,
    quiet    :: Bool = true
)
    # number of cases?
    n = length(y)

    # if b is not warm-started, then ensure that it is not entirely zero
    if all(v.bk .== 0)
        v.bk[1] = logit(abs(mean(y)))
    end

    # initialize intermediate arrays for calculations
    #BLAS.gemv!('N', one(T), x, b, zero(T), xβ)
    A_mul_B!(v.xb, v.xk, v.bk)
    log2xb!(v.lxb, v.l2xb, v.xb)
    copy!(v.bk0, v.bk)
    fill!(v.db, zero(T))

    # track objective
    old_obj = convert(T, Inf)
    # new_obj = (sum(log(1 + exp(v.xb))) - dot(y, v.xb)) / n + λ*sumabs2(v.bk) / 2
    new_obj = logistic_loglik(v.xb, y, v.bk, λ)

    # output progress to console
    quiet || println("Iter\tHalves\tObjective")

    i = 0
    bktrk = 0
    # enter loop for Newton's method
    for i = 1:max_iter

        # db = (x'*(lxβ - y)) / n + λ*b
        BLAS.axpy!(-one(T), sdata(y), sdata(v.lxb))
        At_mul_B!(v.db, v.xk, v.lxb)
        scale!(v.db, 1/n)
        BLAS.axpy!(λ, v.bk, v.db)

        # d2b = (x'*diagm(l2xb)*x)/n + λ*I
        # note that log2xb!() already performs division by n on l2xb
        copy!(v.xk2, v.xk)
        scale!(v.l2xb, v.xk2)
        At_mul_B!(v.d2b, v.xk, v.xk2)
        v.d2b += λ*I

        # b = b0 - ntb = b0 - inv(d2b)*db
        #   = b0 - inv[ x' diagm(pi) diagm(1 - pi) x + λ*I] [x' (pi - y) + λ*b]
        try
            v.ntb = v.d2b \ v.db
        catch e
            warn("in fit_logistic, aborting after caught error: ", e)
            return i, div(bktrk, max_iter)
        end
        copy!(v.bk, v.bk0)
        BLAS.axpy!(-one(T), v.ntb, v.bk)

        # compute objective
        A_mul_B!(v.xb, v.xk, v.bk)
        # new_obj = (sum(log(1 + exp(v.xb))) - dot(y, v.xb)) / n + λ*sumabs2(v.bk) / 2
        new_obj = logistic_loglik(v.xb, y, v.bk, λ)

        # control against nonfinite objective
        isfinite(new_obj) || throw(error("in fit_logistic, objective is no longer finite"))

        # backtrack
        j = 0
        while (new_obj > old_obj + tol) && (j < 50)

            # increment iterator
            j += 1

            # b = b0 - 0.5*ntb
            copy!(v.bk, v.bk0)
            #BLAS.axpy!(p,-one(T) / (2^j),ntb,1,b,1)
            #BLAS.axpy!(-one(T) / (2^j),ntb,v.b)
            BLAS.axpy!(-1/(2^j), v.ntb, v.bk)

            # recalculate objective
            A_mul_B!(v.xb, v.xk, v.bk)
            new_obj = logistic_loglik(v.xb, y, v.bk, λ)

            # ensure a finite objective
            isfinite(new_obj) || throw(error("in fit_logistic, objective is no longer finite"))
        end

        # accumulate total backtracking steps
        bktrk += j

        # track distance between iterates
        dist = euclidean(v.bk, v.bk0) / (norm(v.bk0,2) + one(T))

        # track progress
        quiet || println(i, "\t", j, "\t", dist)

        # check for convergence
        # if converged, then return b
        dist < tol && return i, div(bktrk, i)

        # unconverged at this point, so update intermediate arrays
        #BLAS.gemv!('N', one(T), sdata(x), sdata(b), zero(T), sdata(xβ))
        A_mul_B!(v.xb, v.xk, v.bk)
        log2xb!(v.lxb, v.l2xb, v.xb)

        # save previous β
        copy!(v.bk0, v.bk)
        old_obj = new_obj
    end

    warn("fit_logistic failed to converge in $(max_iter) iterations, exiting...")
    return i, div(bktrk, max_iter)
end

"""
L0 PENALIZED LOGISTIC REGRESSION

    L0_log(x,y,k) -> IHTLogResults

This routine minimizes the loss function given by the negative logistic loglikelihood

    L(β) = sum(log(1 + exp(x*β))) - y'*x*β

subject to `β` lying in the set S_k = { x in R^p : || x ||_0 <= k }.
To ensure a stable model selection process, the optimization is performed over a Tikhonov-regularized copy of `L(β)`; the actual optimized objective is

    g(β) = L(β) + 0.5*lambda*sumabs2(β)

where `lambda` controls the strength of the L2 penalty.
This function extends the [MATLAB source code](http://users.ece.gatech.edu/sbahmani7/GraSP.html) for Sohail Bahmani's nonlinear hard thresholding pursuit framework [GraSP](http://jmlr.csail.mit.edu/papers/v14/bahmani13a.html).

Arguments:

- `x` is the `n` x `p` data matrix.
- `y` is the `n`-dimensional continuous response vector.
- `k` is the desired model size (support).

Optional Arguments:

- `lambda` is the strength of the regularization parameter. Defaults to `sqrt(log(p)/n)`.
- `mu` is the step size used in gradient descent. Defaults to `1.0`.
- `max_iter` is the maximum number of iterations for the algorithm. Defaults to `100`.
- `max_step` is the maximum number of backtracking steps for the step size calculation. Defaults to `100`.
- `tol` is the global tolerance. Defaults to `1e-6`.
- `tolG` is the tolerance for the gradient. Convergence occurs if the norm of the largest `3*k` elements of the gradient falls below `tolG`. Defaults to `1e-6`.
- `tolrefit` is the tolerance for Newton's method in the refitting routine. The refitting algorithm converges when the loss function calculated on the active set (that is, the refit coefficients) falls below `tolrefit`. Defaults to `1e-6`.
- `refit` is a `Bool` that controls refitting at every iterations. It is wise to refit the nonzero coefficients of `b` at every iteration as this may improve convergence behavior and estimation. Defaults to `true` (refit at every iteration).
- `quiet` is a `Bool` that controls algorithm output. Defaults to `true` (no output).

Outputs are wrapped into an `IHTLogResults` object with the following fields:

- 'time' is the compute time for the algorithm. Note that this does not account for time spent initializing optional argument defaults.
- 'iter' is the number of iterations that the algorithm took before converging.
- 'loss' is the optimal loss (half of residual sum of squares) at convergence.
- 'beta' is the final estimate of `b`.
- `active` is the final active set. The size of this set ranges from `k` to `2k`.
"""
function L0_log{T <: Float, V <: DenseVector}(
    x        :: DenseMatrix{T},
    y        :: V,
    k        :: Int;
    v        :: IHTLogVariables{T, V} = IHTLogVariables(sdata(x), sdata(y), k),
    lambda   :: T    = default_lambda(x, y),
    mu       :: T    = one(T),
    tol      :: T    = convert(T, 1e-6),
    tolG     :: T    = convert(T, 1e-3),
    tolrefit :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    max_step :: Int  = 100,
    refit    :: Bool = true,
    quiet    :: Bool = true
)

    # start timer
    tic()

    # problem dimensions?
    n,p = size(x)

    # check arguments
    length(y) == n      || throw(ArgumentError("Length n = $n of response vector y does not match number of rows = $(size(x,1)) in x"))
    k        >  p       && throw(ArgumentError("Value of argument k = $k exceeds number of predictors p = $p"))
    lambda   <  zero(T) && throw(ArgumentError("Value of argument lambda = $lambda must be nonnegative"))
    mu       >  zero(T) || throw(ArgumentError("Value of argument mu = $mu must be positive"))
    tol      >  eps(T)  || throw(ArgumentError("Value of argument tol must exceed machine precision"))
    tolG     >  eps(T)  || throw(ArgumentError("Value of argument tolG must exceed machine precision"))
    tolrefit >  eps(T)  || throw(ArgumentError("Value of argument tolrefit must exceed machine precision"))
    max_iter >= 0       || throw(ArgumentError("Value of max_iter must be nonnegative\n"))
    max_step >= 0       || throw(ArgumentError("Value of max_step must be nonnegative\n"))

    # initialize return values
    iter      = 0                 # number of iterations of L0_reg
    exec_time = zero(T)           # compute time *within* L0_reg
    loss      = convert(T, Inf)   # loss function value

    # initialize algorithm parameters
    num_df    = min(p, 3*k)       # largest 3*k active components of gradient
    short_df  = min(p, 2*k)       # largest 2*k active components of gradient
    converged = false             # is algorithm converged?
    stuck     = false             # is algorithm stuck in a cycle?
    normdf    = convert(T, Inf)   # norm of active portion of gradient
    loss0     = convert(T, Inf)   # penultimate loss function value
    loss00    = convert(T, Inf)   # antepenultimate loss function value
    bktrk     = 0                 # number of Newton backtracking steps
    nt_iter   = 0                 # number of Newton iterations
    lt        = length(v.active)  # size of current active set?

    # formatted output to monitor algorithm progress
    quiet || print_header_log()

    # main GraSP iterations
    for iter = 1:max_iter

        # notify and break if maximum iterations are reached
        # also break if algorithm is cycling
        if iter >= max_iter || stuck

            # warn about hitting maximum iterations
            iter >= max_iter && print_maxiter(max_iter, loss)

            # warn about cycling
            stuck && warn("L0_log appears to be cycling after $iter iterations, aborting early...\nCurrent loss: $loss\n")

            # send elements below tol to zero
            threshold!(v.b, tol)

            # if requested, apply final refit without regularization
            if refit
                copy!(v.idxs0, v.idxs)
                update_indices!(v.idxs, v.b)
                !isequal(v.idxs, v.idxs0) && copy!(v.xk, view(x, :, v.idxs))
                nt_iter, bktrk = fit_logistic!(v, y, zero(T), tol=tolrefit, max_iter=max_step, quiet=true)
                v.b[v.idxs] = v.bk
            end

            # stop timer
            exec_time = toq()

            # these are output variables for function
            # wrap them into a Dict and return
            return IHTLogResults(exec_time, iter, loss, copy(v.b), copy(v.active))
        end

        # save previous loss, iterate
        loss00 = loss0
        loss0  = loss
        copy!(v.b0, v.b)

        # size of current active set?
        if iter > 1
            v.active = union(v.dfidxs[1:short_df], v.bidxs[1:k])
        end
        lt = length(v.active)

        # update x*b
        # no need to compute anything if b = 0
        if all(v.b .== 0)
            fill!(v.xb, zero(T))
        else
            # A_mul_B!(v.xb, x, v.b)
            update_xb!(v.xb, x, v.b, v.active, lt)
        end

        # recompute active loss = (-dot(y,Xb) + sum(log(1.0 + exp(Xb)))) / n + 0.5*lambda*sumabs2(b[active])
        # special case: b = 0 --> Xb = 0 --> loss = n*log(1 + exp(0))/n + 0.5*lambda*norm(0)
        if iter < 2
            loss = 1 / eps(T)
        elseif all(v.xb .== 0)
            loss = log(2)
        else
            loss = logistic_loglik(v.xb, y, v.b, v.active, lambda, k)
        end

        # guard against numerical instabilities in loss function
        check_finiteness(loss)

        # recompute active gradient df[active] = (x[:,active]'*(logistic(Xb) - y)) / n + lambda*b[active]
        # arrange calculations differently if active set is entire support 1, 2, ..., p
        if lt == p
            logistic_grad!(v.df, v.lxb, x, y, v.b, v.xb, lambda)
        else
            logistic_grad!(v.df, v.lxb, x, y, v.b, v.xb, v.active, lt, lambda)
        end
#        logistic!(v.lxb, v.xb)
#        v.df = (x' * (v.lxb - y)) / n + lambda*v.b

        # need norm of top 3*k components of gradient
        normdf = df_norm(v.df, v.dfidxs, 1, num_df)

        # guard against numerical instabilities in gradient
        check_finiteness(normdf)

        # identify 2*k dominant directions in gradient
        selectperm!(v.dfidxs, v.df, 1:num_df, by=abs, rev=true, initialized=true)

        # clean b and fill, b[active] = b0[active] - μ*df[active]
        # note that sparsity level is size(active) which is one of [k, k+1, ..., 3*k]
        #update_x!(v.b, v.b0, v.df, v.active, mu)
        fill!(v.b, zero(T))
        v.b[v.active] = v.b0[v.active] .- mu.*v.df[v.active]

        # now apply hard threshold on model to original desired sparsity k
        project_k!(v.b, k)
        selectperm!(v.bidxs, v.b, 1:k, by=abs, rev=true, initialized=true)
        v.bk = v.b[view(v.bidxs, 1:k)]
        fill!(v.b, 0)
        v.b[view(v.bidxs, 1:k)] = v.bk
        update_xb!(v.xb, x, v.b, v.active, lt)
        loss = logistic_loglik(v.xb, y, v.b, v.active, lambda, k)

        # check for convergence
        converged_obj  = abs(loss - loss0) < tol
        converged_grad = normdf < tolG
        converged      = converged_obj || converged_grad
        stuck          = abs(loss - loss00) < tol
#        stuck          = abs(loss - loss00) < tol && converged_grad

        # output algorithm progress
        quiet || @printf("%d\t%d\t%d\t%3.7f\t%3.7f\n",iter,nt_iter,bktrk,loss,normdf)

        # check for convergence
        # if converged and in feasible set, then algorithm converged before maximum iteration
        # perform final computations and output return variables
        if converged

            # send elements below tol to zero
            threshold!(v.b, tol)

            # if requested, apply final refit without regularization
            # note that regularization is NOT set exactly to zero!
            # this aids numerical stability
            if refit
                copy!(v.idxs0, v.idxs)
                update_indices!(v.idxs, v.b)
#                !isequal(idxs,idxs0) && update_xk!(xk, x, idxs, k=k, n=n, p=p)
                copy!(v.xk, view(x, :, v.idxs))
                try
                    nt_iter, bktrk = fit_logistic!(v, y, 2*eps(T), tol=tolrefit, max_iter=max_step, quiet=true)
                    v.b[v.idxs] = v.bk
                catch e
#                    warn("in final refitting, caught error: ", e)
#                    warn("skipping final refit")
                    nt_iter = 0
                    bktrk   = 0
                end # end try-catch for refit
            end

            # stop time
            exec_time = toq()

            # announce convergence
            !quiet && print_log_convergence(iter, loss, exec_time, normdf)

            # return output
            return IHTLogResults(exec_time, iter, loss, copy(v.b), copy(v.active))
        end # end convergence check

        # unconverged at this point
        # refit nonzeroes in b?
        if refit

            # update boolean vector of nonzeroes
            copy!(v.idxs0, v.idxs)
            update_indices!(v.idxs, v.b)

            # update active set of x, if necessary
#            (iter == 1 || !isequal(idxs,idxs0)) && update_xk!(xk, x, idxs, k=k, n=n, p=p)
            copy!(v.xk, view(x, :, v.idxs))
            copy!(v.xk2, v.xk)

            # attempt refit but guard against possible instabilities or singularities
            # these are more frequent if lambda is small
            # if refitting destabilizes, then leave b alone
            try
                nt_iter, bktrk = fit_logistic!(v, y, lambda, tol=tolrefit, max_iter=max_step, quiet=true)
                v.b[v.idxs] = v.bk
            catch e
#                warn("in refitting, caught error: ", e)
#                warn("skipping refit")
                nt_iter = 0
                bktrk   = 0
            end # end try-catch for refit
        end # end refit

    end # end main GraSP iterations

    # return null result
    return IHTLogResults(zero(T), 0, zero(T), zeros(T, 1), zeros(Int, 1))
end # end L0_log

"""
COMPUTE AN IHT REGULARIZATION PATH FOR LOGISTIC REGRESSION

    iht_path_log(x,y,path) -> SparseCSCMatrix

This subroutine computes best logistic models for matrix `x` and response `y` by calling `L0_log` for each model over a regularization path denoted by `path`.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `path` is an `Int` vector that contains the model sizes to test.

Optional Arguments:

- `tol` is the global convergence tolerance for `L0_log`. Defaults to `1e-6`.
- `tolG` is the tolerance for the gradient. Convergence occurs if the norm of the largest `3*k` elements of the gradient falls below `tolG`. Defaults to `1e-6`.
- `tolrefit` is the tolerance for Newton's method in the refitting routine. The refitting algorithm converges when the loss function calculated on the active set (that is, the refit coefficients) falls below `tolrefit`. Defaults to `1e-6`.
- `max_iter` caps the number of iterations for the algorithm. Defaults to `100`.
- `max_step` caps the number of backtracking steps in the IHT kernel. Defaults to `100`.
- `refit` is a `Bool` that controls refitting at every iterations. It is wise to refit the nonzero coefficients of `b` at every iteration as this may improve convergence behavior and estimation. Defaults to `true` (refit at every iteration).
- `quiet` is a Boolean that controls the output. Defaults to `true` (no output).

Output:

- a sparse `p` x `length(path)` matrix where each column contains the computed model for each component of `path`.
"""
function iht_path_log{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    path     :: DenseVector{Int};
    lambdas  :: Vector{T} = default_lambda(x, y, length(path)),
    mu       :: T    = one(T),
    tol      :: T    = convert(T, 1e-6),
    tolG     :: T    = convert(T, 1e-3),
    tolrefit :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    max_step :: Int  = 50,
    refit    :: Bool = true,
    quiet    :: Bool = true,
)

    # dimensions of problems?
    n,p = size(x)

    # how many models will we compute?
    num_models = length(path)

    # preallocate temporary arrays
    v = IHTLogVariables(sdata(x), sdata(y), 1)

    # betas will be a sparse matrix with models
    betas = spzeros(T, p, num_models)

    # compute the path
    for i = 1:num_models

        # model size?
        q = path[i]
        quiet || println("Current model size: $q")

        # update temporary variables with new threshold size q
        update_variables!(v, x, q)

        # store projection of beta onto largest k nonzeroes in magnitude
        project_k!(v.b, q)

        # current regularization parameter?
        λ = lambdas[i]

        # now compute current model
        output = L0_log(x, y, q, v=v, lambda=λ, mu=mu, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, refit=refit, quiet=quiet)

        # extract and save model
        betas[:,i] = sparsevec(output.beta)

        # run garbage collector to clear temporary arrays for next step of path
        #gc()
    end

    # return a sparsified copy of the models
    return betas
end

"""
COMPUTE ONE FOLD IN A CROSSVALIDATION SCHEME FOR A PENALIZED LOGISTIC REGRESSION REGULARIZATION PATH

    one_fold_log(x,y,path,folds,fold) -> Vector{Float}

For a regularization path given by the `Int` vector `path`,
this function performs penalized logistic regression on `x` and `y` and computes an out-of-sample error based on the indices given in `folds`.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `path` is an `Int` vector that specifies which model sizes to include in the path, e.g. `path = collect(k0:increment:k_end)`.
- `folds` is an `Int` vector indicating which component of `y` goes to which fold, e.g. `folds = IHT.cv_get_folds(n,nfolds)`
- `fold` is the current fold to compute.

Optional Arguments:

- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-4`.
- `max_iter` caps the number of permissible iterations in the IHT algorithm. Defaults to `100`.
- `max_step` caps the number of permissible backtracking steps. Defaults to `100`.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).
- `deviance` is an `String` that specifies the crossvalidation criterion, either `"deviance"` for the logistic loglikelihood or `"class"` for misclassification error. Defaults to `"deviance"`.

Output:

- `errors` is a vector of out-of-sample errors (MSEs) for the current fold.
"""
function one_fold_log{T <: Float}(
    x         :: DenseMatrix{T},
    y         :: DenseVector{T},
    path      :: DenseVector{Int},
    folds     :: DenseVector{Int},
    fold      :: Int;
    lambdas   :: Vector{T} = default_lambda(x, y, length(path)),
    criterion :: String = "deviance",
    tol       :: T      = convert(T, 1e-6),
    tolG      :: T      = convert(T, 1e-3),
    tolrefit  :: T      = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    refit     :: Bool   = true,
    quiet     :: Bool   = true,
)
    # dimension of problem?
    n,p = size(x)
    nmodels = length(path)

    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))

    # make vector of indices for folds
    # train_idx is the vector that indexes the TRAINING set
    train_idx = folds .!= fold
    test_size = sum(!train_idx)

    # allocate the arrays for the training set
    x_train = x[train_idx,:]
    y_train = y[train_idx]

    # compute the regularization path on the training set
    betas = iht_path_log(x_train, y_train, path, tol=tol, tolG=tolG, tolrefit=tolrefit, max_iter=max_iter, max_step=max_step, refit=refit, quiet=quiet, lambdas=lambdas)

    # (masked) estimated responses
    # mask based on test_idx
    xβ  = x*betas
    fill!(view(xβ, train_idx, :), zero(T))
    lxb = zeros(T, size(xβ))
    logistic!(lxb, xβ)
    fill!(view(lxb, train_idx, :), zero(T))
    ytest  = copy(y)
    fill!(view(ytest, train_idx), zero(T))

    # compute the mean out-of-sample error for the TEST set
    errors = zeros(T, nmodels)
    if criterion == "class" # compute misclassification error
        errors .= vec(sumabs(lxb .- ytest, 1)) ./ sum(!train_idx)
    else # else -> criterion == "deviance" = 2*(neg logistic loss with no reg)
        #errors .= - one(T) ./ n .* ( vec( sum(xβ .* y, 1) .- sum(log(one(T) .+ exp(xβ)), 1) ) )
        errors .= - 2 ./ test_size .* ( vec( sum(xβ .* ytest, 1) .- sum(log(one(T) .+ exp(xβ)), 1) ) )
    end
    return errors :: Vector{T}
end

#"""
#    pfold_log(xfile, xtfile, x2file, yfile, meanfile, invstdfile, path, folds, numfolds[, pids=procs()])
#
#This function is the parallel execution kernel in `cv_log()`. It is not meant to be called outside of `cv_iht()`.
#It will distribute `numfolds` crossvalidation folds across the processes supplied by the optional argument `pids` and call `one_fold_log()` for each fold.
#Each fold will compute a regularization path given by `path`.
#`pfold()` collects the vectors of MSEs returned by calling `one_fold_log()` for each process, reduces them, and returns their average across all folds.
#"""
function pfold_log{T <: Float}(
    x         :: DenseMatrix{T},
    y         :: DenseVector{T},
    path      :: DenseVector{Int},
    folds     :: DenseVector{Int},
    q         :: Int;
    pids      :: Vector{Int} = procs(),
    lambdas   :: DenseVector{T} = default_lambda(x, y, length(path)),
    criterion :: String = "deviance",
    tol       :: T      = convert(T, 1e-6),
    tolG      :: T      = convert(T, 1e-3),
    tolrefit  :: T      = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    quiet     :: Bool   = true,
    refit     :: Bool   = true,
    header    :: Bool   = false
)

    # problem size?
    n,p = size(x)

    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))

    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = SharedArray(T, (length(path),q), pids=pids) :: SharedMatrix{T}

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > q && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        r = remotecall_fetch(worker) do
                                one_fold_log(x, y, path, folds, current_fold, max_iter=max_iter, max_step=max_step, quiet=quiet, tol=tol, tolG=tolG, tolrefit=tolrefit, refit=refit, criterion=criterion, lambdas=lambdas)
                        end # end remotecall_fetch()
                        setindex!(results, r, :, current_fold)
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (vec(sum(results, 2) ./ q)) :: Vector{T}
end

"""
PARALLEL CROSSVALIDATION ROUTINE FOR GRASP ALGORITHM FOR PENALIZED LOGISTIC REGRESSION

    cv_iht_log(x,y,path,q) -> Vector{Float}

This function will perform `q`-fold cross validation for the ideal model size in GraSP logistic regression.
It computes several paths as specified in the `path` argument using the design matrix `x` and the response vector `y`.
Each path is asynchronously spawned using any available processor.
For each path, one fold is held out of the analysis for testing, while the rest of the data are used for training.
The function to compute each path, `one_fold_log()`, will return a vector of out-of-sample errors (MCEs).

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `path` is an `Int` vector that specifies which model sizes to include in the path, e.g. `path = collect(k0:increment:k_end)`.
- `q` is the number of folds to compute.

Optional Arguments:

- `n` is the number of samples. Defaults to `length(y)`.
- `p` is the number of predictors. Defaults to `size(x,2)`.
- `folds` is the partition of the data. Defaults to `IHT.cv_get_folds(n,q)`.
- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-4`.
- `max_iter` caps the number of permissible iterations in the IHT algorithm. Defaults to `100`.
- `max_step` caps the number of permissible backtracking steps. Defaults to `100`.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).
   *NOTA BENE*: each processor outputs feed to the console without regard to the others,
   so setting `quiet = false` can yield very messy output!
- `refit` is a `Bool` to indicate whether or not to recompute the best model. Note that this argument affects both the continuous refitting in the GraSP algorithm in addition to a final refit after crossvalidation. Defaults to `true` (recompute).
- `deviance` is an `String` that specifies the crossvalidation criterion, either `"deviance"` for the logistic loglikelihood or `"class"` for misclassification error. Defaults to `"deviance"`.

Output:

- `errors` is the averaged MSE over all folds.

If called with `refit = true`, then the output also includes, for best model size `k_star`:

- `b`, a vector of `k_star` floats
- `bidx`, a vector of `k_star` indices indicating the support of the best model.
"""
function cv_log{T <: Float}(
    x         :: DenseMatrix{T},
    y         :: DenseVector{T};
    q         :: Int              = cv_get_num_folds(3, 5),
    path      :: DenseVector{Int} = collect(1:min(size(x,2),20)),
    folds     :: DenseVector{Int} = cv_get_folds(sdata(y),q),
    pids      :: Vector{Int}      = procs(),
    lambdas   :: Vector{T}        = default_lambda(x, y, length(path)),
    criterion :: String = "deviance",
    tol       :: T      = convert(T, 1e-6),
    tolG      :: T      = convert(T, 1e-3),
    tolrefit  :: T      = convert(T, 1e-6),
    max_iter  :: Int    = 100,
    max_step  :: Int    = 100,
    quiet     :: Bool   = true,
    refit     :: Bool   = true
)

    # problem size?
    n,p = size(x)
    @assert n == length(y)

    # ensure that crossvalidation criterion is valid
    criterion in ["deviance", "class"] || throw(ArgumentError("Argument criterion must be 'deviance' or 'class'"))

    # how many elements are in the path?
    num_models = length(path)

    # compute crossvalidation deviances
    errors = pfold_log(x, y, path, folds, q, max_iter=max_iter, max_step=max_step, quiet=quiet, pids=pids, tolG=tolG, tolrefit=tolrefit, refit=refit, criterion=criterion, lambdas=lambdas)

    # what is the best model size?
    k = path[indmin(errors)] :: Int

    # print results
    !quiet && print_cv_results(errors, path, k)

    # recompute ideal model
    # get λ value for best model
    λ = lambdas[path .== k][1]

    # use L0_log to extract model
    # with refit = true, L0_log will continuously refit predictors
    # no final refitting code necessary
    output = L0_log(x, sdata(y), k, max_iter=max_iter, max_step=max_step, quiet=quiet, tol=tol, tolG=tolG, tolrefit=tolrefit, refit=true, lambda=zero(T))#lambda=λ)

    # which components of beta are nonzero?
    bidx = find(output.beta)

    return IHTCrossvalidationResults(errors, sdata(path), output.beta[bidx], bidx, k)
end


text:
using Test
include("../src/classification.jl")
using Main.Classification
using Random

@testset "Gradient Descent" begin
    X = randn(10, 15)
    y = rand(0:1, 10)
    @testset "Context" begin
        gd = GDContext()
        gd = GDContext(learningrate = 0.1)
        @test true
    end
    @testset "Optimizer from context" begin
        gd = GDContext()
        o = GradientDescent(gd)
        o = GradientDescent()
        @test true
    end
    @testset "Optimizer methods" begin
        o = GradientDescent()
        addvalidationset!(o, X, y)
        @test_throws InexactError addvalidationset!(o, X, randn(10))
    end
    @testset "Fit" begin
        o = GradientDescent()
        β = zeros(15)
        fit!(o, β, X, y)
        @test length(o.loss) == o.iterations
        addvalidationset!(o, X, y)
        fit!(o, β, X, y)
        if o.converged
            @test length(o.validationloss) == o.iterations - 1
        else
            @test length(o.validationloss) == o.iterations
        end
    end
end

@testset "Stochastic Gradient Descent" begin
    X = randn(10, 15)
    y = rand(0:1, 10)
    @testset "Context" begin
        gd = SGDContext()
        gd = SGDContext(learningrate = 0.1)
        @test true
    end
    @testset "Optimizer from context" begin
        gd = SGDContext()
        o = StochasticGradientDescent(gd)
        o = StochasticGradientDescent()
        @test true
    end
    @testset "Optimizer methods" begin
        o = StochasticGradientDescent()
        addvalidationset!(o, X, y)
        @test_throws InexactError addvalidationset!(o, X, randn(10))
    end
    # Messed up types, causing this to fail
    # @testset "Fit" begin
        # o = StochasticGradientDescent()
        # β = zeros(15)
        # fit!(o, β, X, y)
        # @test length(o.loss) == o.iterations
        # addvalidationset!(o, X, y)
        # fit!(o, β, X, y)
        # @test true
    # end
end

@testset "Classifier" begin
    X = randn(10, 15)
    y = rand(0:1, 10)
    @testset "Gradient Descent Classifier" begin
        o = GradientDescent()
        clf = LogisticClassifier(o)
        add_designmatrix!(clf, X)
        fit!(clf, y)
        fit!(clf, X, y)
        @test true
    end
end


text:
"""
    simulate_Ising_parallel(
        algorithm ::Function,
        Ls        ::Array{Int, 1},
        βs        ::Union{Array{Float64}, AbstractRange{Float64}};
        n_sweeps  ::Int     =10^5,
        cutoff    ::Float64 =0.2,
        start_type::Symbol  =:cold,
        seed      ::Int     =8,
        )

Runs multiple Ising model simulations in parallel.

Arguments
* `algorithm`  : The algorithm with which to run the simulation
* `Ls`         : The side lengths for which to run the simulations with
* `βs`         : The inverse temperatures for which to run the simulations with
* `n_sweeps`   : The number of update sweeps for each simulation
* `cutoff`     : How much of the initial data to throw away (i.e. to ignore the observations
    when the system was thermalizing
* `start_type` : Whether or not to start the simulation with a cold or hot start, available
    options are `[:cold, :hot]`
* `seed`       : Random number generator seed value
Returns
* A dictionary with all of the simulations, where the keys are the specific `(L, β)` values
"""
function simulate_Ising_parallel(
    algorithm::Function,
    Ls::Array{Int,1},
    βs::Union{Array{Float64},AbstractRange{Float64}};
    n_sweeps::Int = 10^5,
    cutoff::Float64 = 0.2,
    start_type::Symbol = :cold,
    seed::Int = 8,
)

    out = Dict{Tuple,Ising}()
    Threads.@threads for β in βs
        Threads.@threads for L in Ls
            out[(L, β)] = algorithm(Ising(
                L,
                β,
                n_sweeps = n_sweeps,
                cutoff = cutoff,
                start_type = start_type,
                seed = seed,
            ))
        end
    end

    return out
end

"""
    simulate_Potts_parallel(
        algorithm ::Function,
        Ls        ::Array{Int, 1},
        βs        ::Union{Array{Float64}, AbstractRange{Float64}};
        n_sweeps  ::Int     =10^5,
        cutoff    ::Float64 =0.2,
        start_type::Symbol  =:cold,
        seed      ::Int     =8,
        )

Runs multiple Potts model simulations in parallel.

Arguments
* `q`          : The number of possible states
* `algorithm`  : The algorithm with which to run the simulation
* `Ls`         : The side lengths for which to run the simulations with
* `βs`         : The inverse temperatures for which to run the simulations with
* `n_sweeps`   : The number of update sweeps for each simulation
* `cutoff`     : How much of the initial data to throw away (i.e. to ignore the observations
    when the system was thermalizing
* `start_type` : Whether or not to start the simulation with a cold or hot start, available
    options are `[:cold, :hot]`
* `seed`       : Random number generator seed value
Returns
* A dictionary with all of the simulations, where the keys are the specific `(L, β)` values
"""
function simulate_Potts_parallel(
    q::Int,
    algorithm::Function,
    Ls::Array{Int,1},
    βs::Union{Array{Float64},AbstractRange{Float64}};
    n_sweeps::Int = 10^5,
    cutoff::Float64 = 0.2,
    start_type::Symbol = :cold,
    seed::Int = 8,
)

    out = Dict{Tuple,Potts}()
    Threads.@threads for β in βs
        Threads.@threads for L in Ls
            out[(L, β)] = algorithm(Potts(
                q,
                L,
                β,
                n_sweeps = n_sweeps,
                cutoff = cutoff,
                start_type = start_type,
                seed = seed,
            ))
        end
    end

    return out
end


text:
########## Collection of obseleted codes ##########

# function rolling_estim(func::Function, X::AbstractMatrix{T}, widx::AbstractArray{Int}, p::Int) where {T<:Real}
#     # wsize = w[end]-w[1]+1
#     # @assert wsize <= size(X,2)
#     return [func(X[:,n+widx]) for n=1:p:size(X,2) if n+widx[end]<size(X,2)]
#     # res = []
#     # for n=1:p:size(X,2)
#     #     idx = n+w
#     #     idx = idx[idx.<=size(X,2)]
#     #     push!(res, func(view(X, :,n+w)))
#     # end
# end

# function rolling_estim(fun::Function, X::AbstractVector{T}, wsize::Int, p::Int=1) where T
#     offset = wsize-1
#     res = [fun(view(X, (n*p):(n*p+wsize-1))) for n=1:p:N]
#     end

#     # y = fun(view(X, idx:idx+offset))  # test return type of fun
#     res = Vector{typeof(y)}(undef, div(length(data)-offset, p))
#     @inbounds for n=1:p:N
#         push!(res, fun(view(X, (idx*p):(idx*p+offset))))
#         end
#     @inbounds for n in eachindex(res)
#         res[n] = fun(hcat(X[idx*p:idx*p+offset]...))
#     end

#     return res
# end

# function scalogram_estim(Cxx::Matrix{Float64}, sclrng::AbstractArray{Int}, ρmax::Int=1)
#     nr, nc = size(Cxx)
#     @assert nr == nc == length(sclrng)
#     @assert ρmax >= 1
#     xvar = Float64[]
#     yvar = Float64[]
#     for ρ=1:ρmax
#         toto = [Cxx[j, ρ*j] for j in 1:nr if ρ*j<=nr]
#         xvar = vcat(xvar, sclrng[1:length(toto)])
#         yvar = vcat(yvar, abs.(toto))
#     end
#     df = DataFrames.DataFrame(
#         xvar=log2.(xvar),
#         yvar=log2.(yvar)
#     )
#     ols_hurst = GLM.lm(@GLM.formula(yvar~xvar), df)
#     hurst_estim = (GLM.coef(ols_hurst)[2]-1)/2
#     return hurst_estim, ols_hurst
# end

function powlaw_estim(X::Vector{Float64}, lags::AbstractArray{Int}, pows::AbstractArray{T}) where {T<:Real}
    @assert length(lags) > 1 && all(lags .> 1)

    # Define the function for computing the p-th moment of the increment
    moment_incr(X,d,p) = mean((abs.(X[d+1:end] - X[1:end-d])).^p)

    # Estimation of Hurst exponent and β
    H = zeros(Float64, length(pows))
    β = zeros(Float64, length(pows))
    C = zeros(Float64, length(pows))

    for (n,p) in enumerate(pows)
        powlaw_estim()
        C[n] = 2^(p/2) * gamma((p+1)/2)/sqrt(pi)

        yp = map(d -> log(moment_incr(X, d, p)), lags)
        xp = p * log.(lags)
        Ap = hcat(xp, ones(length(xp)))  # design matrix
        H[n], β[n] = Ap \ yp  # estimation of H and β

        # dg = DataFrames.DataFrame(xvar=Ap, yvar=yp)
        # ols = GLM.lm(@GLM.formula(yvar ~ xvar), dg)
        # β[n], H[n] = GLM.coef(ols)
    end

    Σ = exp.((β-log.(C))./pows)

    hurst = sum(H) / length(H)
    σ = sum(Σ) / length(Σ)

    return hurst, σ
end

function powlaw_estim_old(X::Vector{Float64}, lags::AbstractArray{Int}, pows::AbstractArray{T}) where {T<:Real}
    # Define the function for computing the p-th moment of the increment
    moment_incr(X,d,p) = mean((abs.(X[d+1:end] - X[1:end-d])).^p)

    # Estimation of Hurst exponent
    Y = zeros(Float64, (length(lags), length(pows)))
    for (n,p) in enumerate(pows)
        Y[:,n] = map(d -> log(moment_incr(X, d, p)), lags)
    end
    dY = diff(Y, 1)
    df = DataFrames.DataFrame(xvar=(diff(log.(lags)) * pows')[:],
                              yvar=dY[:])
    ols_hurst = GLM.lm(@GLM.formula(yvar ~ xvar), df)
    Hurst = GLM.coef(ols_hurst)[2]  # estimation of hurst exponent

    # Estimation of volatility
    # Constant in the p-th moment of normal distribution, see
    # https://en.wikipedia.org/wiki/Normal_distribution#Moments
    cps = [2^(p/2) * gamma((p+1)/2)/sqrt(pi) for p in pows]
    Z = Y -  Hurst * (log.(lags) * pows') - ones(lags) * log.(cps)'
    dg = DataFrames.DataFrame(xvar=(ones(lags) * pows')[:],
                              yvar=Z[:])
    ols_sigma = GLM.lm(@GLM.formula(yvar ~ xvar), dg)

    return Y, Dict('H'=>ols_hurst, 'σ'=>ols_sigma)
end

"""
Extract the estimates of Hurst and voaltility from the result of `powlaw_estim`.
"""
function powlaw_coeff(ols::Dict, h::Float64)
    H = GLM.coef(ols['H'])[2]
    σ = exp(GLM.coef(ols['σ'])[2] - H * log(h))

    return H, σ
end

##### Special functions #####

"""
Compute the continued fraction involved in the upper incomplete gamma function using the modified Lentz's method.
"""
function _uigamma_cf(s::Complex, z::Complex; N=100, epsilon=1e-20)
#     a::Complex = 0
#     b::Complex = 0
#     d::Complex = 0
    u::Complex = s
    v::Complex = 0
    p::Complex = 0

    for n=1:N
#         a, b = (n%2==1) ? ((-div(n-1,2)-s)*z, s+n) : (div(n,2)*z, s+n)
        a, b = (n%2==1) ? ((-div(n-1,2)-s), (s+n)/z) : (div(n,2), (s+n)/z)
        u = b + a / u
        v = 1/(b + a * v)
        d = log(u * v)
        (abs(d) < epsilon) ? break : (p += d)
#         println("$(a), $(b), $(u), $(v), $(d), $(p), $(exp(p))")
    end
    return s * exp(p)
end

doc"""
    uigamma0(z::Complex; N=100, epsilon=1e-20)

Upper incomplete gamma function with vanishing first argument:
$$ \Gamma(0,z) = \lim_{a\rightarrow 0} \Gamma(a,z) $$

Computed using the series expansion of the [exponential integral](https://en.wikipedia.org/wiki/Exponential_integral) $E_1(z)$.
"""
function uigamma0(z::Number; N=100, epsilon=1e-20)
    #     A::Vector{Complex} = [(-z)^k / k / exp(lgamma(k+1)) for k=1:N]
    #     s = sum(A[abs.(A)<epsilon])
    s::Complex = 0
    for k=1:N
        d = (-z)^k / k / exp(lgamma(k+1))
        (abs(d) < epsilon) ? break : (s += d)
    end
    r = -(eulergamma + log(z) + s)
    return (typeof(z) <: Real ? real(r) : r)
end

# """
# Upper incomplete gamma function.
# """
# function uigamma(a::Real, z::T; N=100, epsilon=1e-8) where {T<:Number}
#     z == 0 && return gamma(a)
#     u::T = z
#     v::T = 0
#     f::T = z
# #     f::Complex = log(z)
#     for n=1:N
#         an, bn = (n%2==1) ? (div(n+1,2)-a, z) : (div(n,2), 1)
#         u = bn + an / u
#         v = bn + an * v
#         f *= (u/v)
# #         f += (log(α) - log(β))
#         println("$(an), $(bn), $(u), $(v), $(f)")
#         if abs(u/v-1) < epsilon
#             break
#         end
#     end
#     return z^a * exp(-z) / f
# #     return z^a * exp(-z-f)
# end

doc"""
    uigamma(s::Complex, z::Complex; N=100, epsilon=1e-20)

Upper incomplete gamma function $\Gamma(s,z)$ with complex arguments.

Computed using the [continued fraction representation](http://functions.wolfram.com/06.06.10.0005.01).
The special case $\Gamma(0,z)$ is computed via the series expansion of the exponential integral $E_1(z)$.

# Reference
- [Upper incomplete gamma function](https://en.wikipedia.org/wiki/Incomplete_gamma_function)
- [Continued fraction representation](http://functions.wolfram.com/06.06.10.0005.01)
- [Exponential integral](https://en.wikipedia.org/wiki/Exponential_integral)
"""

function uigamma(s::Number, z::Number; N=100, epsilon=1e-20)
    if abs(s) == 0
        return uigamma0(z; N=N, epsilon=epsilon)
    end

    r = gamma(s) - z^s * exp(-z) / _uigamma_cf(Complex(s), Complex(z); N=N, epsilon=epsilon)
    return (typeof(s)<:Real && typeof(z)<:Real) ? real(r) : r
end

doc"""
    ligamma(s::Complex, z::Complex; N=100, epsilon=1e-20)

Lower incomplete gamma function $\gamma(s,z)$ with complex arguments.
"""
function ligamma(s::Number, z::Number; N=100, epsilon=1e-20)
    return gamma(s) - uigamma(s, z; N=N, epsilon=epsilon)
end

function wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractArray{Int}, v::Int, H::Real, σ::Real; mode::Symbol=:center)
    N, d = size(X)  # length and dim of X
    A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
    C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

    Σ = σ^2 * C1 .* A
    # Σ += Matrix(1.0I, size(Σ)) * max(1e-8, mean(abs.(Σ))*1e-5)
    # println("H=$H, σ=$σ, det(Σ)=$(det(Σ))")

    # method 1:
    # iX = Σ \ X'  # <- unstable!

    # method 2:
    # iΣ = pinv(Σ)  # regularization by pseudo-inverse
    # iX = iΣ * X'

    # method 3:
    iX = lsqr(Σ, X')

    return -1/2 * (tr(X*iX) + N*logdet(Σ) + N*d*log(2π))
end

function grad_wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractArray{Int}, v::Int, H::Real, σ::Real; mode::Symbol=:center)
    N, d = size(X)  # length and dim of X
    A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
    C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]
    dAda = [log(i*j) for i in sclrng, j in sclrng] .* A
    dC1da = [diff_C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

    Σ = σ^2 * C1 .* A
    # Σ += Matrix(1.0I, size(Σ)) * max(1e-8, mean(abs.(Σ))*1e-5)

    dΣda = σ^2 * (dC1da .* A + C1 .* dAda)
    dΣdb = 2σ * C1 .* A

    # method 1:
    # iX = Σ \ X'
    # da = N * tr(Σ \ dΣda) - tr(iX' * dΣda * iX)
    # db = N * tr(Σ \ dΣdb) - tr(iX' * dΣdb * iX)

    # method 2:
    # iΣ = pinv(Σ)  # regularization by pseudo-inverse
    # iX = iΣ * X'
    # da = N * tr(iΣ * dΣda) - tr(iX' * dΣda * iX)
    # db = N * tr(iΣ * dΣdb) - tr(iX' * dΣdb * iX)

    # method 3:
    iX = lsqr(Σ, X')
    da = N * tr(lsqr(Σ, dΣda)) - tr(iX' * dΣda * iX)
    db = N * tr(lsqr(Σ, dΣdb)) - tr(iX' * dΣdb * iX)

    return  -1/2 * [da, db]
end

function wavelet_MLE_estim(X::Matrix{Float64}, sclrng::AbstractArray{Int}, v::Int, d::Int=1, vars::Symbol=:all, init::Vector{Float64}=[0.5,1.]; cflag::Bool=false, dflag::Bool=false, mode::Symbol=:center)
    @assert size(X,2) == length(sclrng)
    @assert length(init) == 2
    @assert d >= 1

    func = x -> ()
    grad = x -> ()
    hurst_estim, σ_estim = init

    if vars == :all
        func = x -> -wavelet_MLE_obj(X[1:d:end, :], sclrng, v, x[1], x[2]; cflag=cflag, mode=mode)
        grad = x -> -grad_wavelet_MLE_obj(X[1:d:end,:], sclrng, v, x[1], x[2]; cflag=cflag, mode=mode)
    elseif vars == :hurst
        func = x -> -wavelet_MLE_obj(X[1:d:end, :], sclrng, v, x[1], σ_estim; cflag=cflag, mode=mode)
        grad = x -> -grad_wavelet_MLE_obj(X[1:d:end,:], sclrng, v, x[1], σ_estim; cflag=cflag, mode=mode)
    else
        func = x -> -wavelet_MLE_obj(X[1:d:end, :], sclrng, v, hurst_estim, x[2]; cflag=cflag, mode=mode)
        grad = x -> -grad_wavelet_MLE_obj(X[1:d:end,:], sclrng, v, hurst_estim, x[2]; cflag=cflag, mode=mode)
    end

    ε = 1e-10
    # opm = Optim.optimize(func, [ε, ε], [1-ε, 1/ε], init, Optim.Fminbox(Optim.GradientDescent()))
    opm = dflag ? Optim.optimize(func, grad, [ε, ε], [1-ε, 1/ε], init, Optim.Fminbox(Optim.GradientDescent()); inplace=false) : Optim.optimize(func, [ε, ε], [1-ε, 1/ε], init, Optim.Fminbox(Optim.GradientDescent()))
    # Optim.BFGS() or Optim.GradientDescent()
    res = Optim.minimizer(opm)

    if vars == :all
        hurst_estim = cflag ? sigmoid(res[1]) : res[1]
        σ_estim = cflag ? exp(res[2]) : res[2]
    elseif vars == :hurst
        hurst_estim = cflag ? sigmoid(res[1]) : res[1]
    else
        σ_estim = cflag ? exp(res[2]) : res[2]
    end

    return (hurst_estim, σ_estim), opm
end

# function partial_bspline_covmat(sclrng::AbstractVector{Int}, v::Int, H::Real, mode::Symbol)
#     return full_bspline_covmat(0, sclrng, v, H, mode)
# end

# function partial_bspline_log_likelihood_H(X::AbstractVecOrMat{T}, sclrng::AbstractVector{Int}, v::Int, H::Real; mode::Symbol=:center) where {T<:Real}
#     # @assert size(X,1) == length(sclrng)
#     Σ = partial_bspline_covmat(sclrng, v, H, mode)
#     # println(size(Σ))
#     # println(size(X))

#     return log_likelihood_H(Σ, X)
# end

# """
# B-Spline wavelet-MLE estimator with partial covariance matrix.
# """
# function partial_bspline_MLE_estim(X::AbstractVecOrMat{T}, sclrng::AbstractVector{Int}, v::Int; init::Real=0.5, mode::Symbol=:center) where {T<:Real}
#     @assert size(X,1) == length(sclrng)

#     func = h -> -partial_bspline_log_likelihood_H(X, sclrng, v, h; mode=mode)

#     ε = 1e-5
#     # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
#     # optimizer = Optim.BFGS()
#     # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
#     opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())

#     hurst = Optim.minimizer(opm)[1]

#     Σ = partial_bspline_covmat(sclrng, v, hurst, mode)
#     σ = sqrt(xiAx(Σ, X) / length(X))

#     return (hurst, σ), opm
# end

function partial_wavelet_log_likelihood_H(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int, H::Real; mode::Symbol=:center)
    N, J = size(X)  # length and dim of X

    A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
    Σ = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng] .* A

    iΣ = pinv(Σ)  # regularization by pseudo-inverse

    return -1/2 * (J*N*log(sum(X' .* (iΣ * X'))) + N*logdet(Σ))
end

function partial_wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int, H::Real, σ::Real; mode::Symbol=:center)
    N, d = size(X)  # length and dim of X

    A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
    C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

    Σ = σ^2 * C1 .* A
    # Σ += Matrix(1.0I, size(Σ)) * max(1e-10, mean(abs.(Σ))*1e-5)

    # println("H=$(H), σ=$(σ), mean(Σ)=$(mean(abs.(Σ)))")
    # println("logdet(Σ)=$(logdet(Σ))")

    # method 1:
    # iX = Σ \ X'

    # method 2:
    iΣ = pinv(Σ)  # regularization by pseudo-inverse
    iX = iΣ * X'  # regularization by pseudo-inverse

    # # method 3:
    # iX = lsqr(Σ, X')

    return -1/2 * (tr(X*iX) + N*logdet(Σ) + N*d*log(2π))
end

# function wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int, H::Real, σ::Real; mode::Symbol=:center)
#     N, d = size(X)  # length and dim of X

#     A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
#     C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

#     Σ = σ^2 * C1 .* A
#     # Σ += Matrix(1.0I, size(Σ)) * max(1e-10, mean(abs.(Σ))*1e-5)

#     # println("H=$(H), σ=$(σ), mean(Σ)=$(mean(abs.(Σ)))")
#     # println("logdet(Σ)=$(logdet(Σ))")

#     # method 1:
#     # iX = Σ \ X'

#     # method 2:
#     iΣ = pinv(Σ)  # regularization by pseudo-inverse
#     iX = iΣ * X'  # regularization by pseudo-inverse

#     # # method 3:
#     # iX = lsqr(Σ, X')

#     return -1/2 * (tr(X*iX) + N*logdet(Σ) + N*d*log(2π))
# end

# older version:
# function wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int, α::Real, β::Real; cflag::Bool=false, mode::Symbol=:center)
#     N, d = size(X)  # length and dim of X

#     H = cflag ? sigmoid(α) : α
#     σ = cflag ? exp(β) : β

#     A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
#     C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

#     Σ = σ^2 * C1 .* A
#     # Σ += Matrix(1.0I, size(Σ)) * max(1e-8, mean(abs.(Σ))*1e-5)

#     # println("H=$(H), σ=$(σ), α=$(α), β=$(β), mean(Σ)=$(mean(abs.(Σ)))")
#     # println("logdet(Σ)=$(logdet(Σ))")

#     # method 1:
#     # iX = Σ \ X'

#     # # method 2:
#     # iΣ = pinv(Σ)  # regularization by pseudo-inverse
#     # iX = iΣ * X'  # regularization by pseudo-inverse

#     # method 3:
#     iX = lsqr(Σ, X')

#     return -1/2 * (tr(X*iX) + N*log(abs(det(Σ))) + N*d*log(2π))
# end

# function grad_wavelet_MLE_obj(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int, α::Real, β::Real; cflag::Bool=false, mode::Symbol=:center)
#     N, d = size(X)  # length and dim of X

#     H = cflag ? sigmoid(α) : α
#     σ = cflag ? exp(β) : β

#     A = [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
#     C1 = [C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]
#     dAda = [log(i*j) for i in sclrng, j in sclrng] .* A
#     dC1da = [diff_C1rho(0, j/i, H, v, mode) for i in sclrng, j in sclrng]

#     if cflag
#         dAda *= diff_sigmoid(α)
#         dC1da *= diff_sigmoid(α)
#     end

#     Σ = σ^2 * C1 .* A
#     # Σ += Matrix(1.0I, size(Σ)) * max(1e-8, mean(abs.(Σ))*1e-5)
#     dΣda = σ^2 * (dC1da .* A + C1 .* dAda)
#     dΣdb = cflag ? 2*Σ : 2σ * C1 .* A

#     # method 1:
#     # iX = Σ \ X'
#     # da = N * tr(Σ \ dΣda) - tr(iX' * dΣda * iX)
#     # db = N * tr(Σ \ dΣdb) - tr(iX' * dΣdb * iX)

#     # method 2:
#     iΣ = pinv(Σ)  # regularization by pseudo-inverse
#     iX = iΣ * X'
#     da = N * tr(iΣ * dΣda) - tr(iX' * dΣda * iX)
#     db = N * tr(iΣ * dΣdb) - tr(iX' * dΣdb * iX)

#     # method 3:
#     # iX = lsqr(Σ, X')
#     # da = N * tr(lsqr(Σ, dΣda)) - tr(iX' * dΣda * iX)
#     # db = N * tr(lsqr(Σ, dΣdb)) - tr(iX' * dΣdb * iX)

#     return  -1/2 * [da, db]
# end

function wavelet_MLE_estim(X::Matrix{Float64}, sclrng::AbstractVector{Int}, v::Int; vars::Symbol=:all, init::Vector{Float64}=[0.5,1.], mode::Symbol=:center)
    @assert size(X,2) == length(sclrng)
    @assert length(init) == 2

    func = x -> ()
    hurst, σ = init
    # println(init)

    if vars == :all
        func = x -> -wavelet_MLE_obj(X, sclrng, v, x[1], x[2]; mode=mode)
    elseif vars == :hurst
        func = x -> -wavelet_MLE_obj(X, sclrng, v, x[1], σ; mode=mode)
    else
        func = x -> -wavelet_MLE_obj(X, sclrng, v, hurst, x[2]; mode=mode)
    end

    ε = 1e-8
    # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
    optimizer = Optim.BFGS()
    opm = Optim.optimize(func, [ε, ε], [1-ε, 1/ε], init, Optim.Fminbox(optimizer))
    # opm = Optim.optimize(func, [ε, ε], [1-ε, 1/ε], init, Optim.Fminbox(optimizer); autodiff=:forward)
    res = Optim.minimizer(opm)

    if vars == :all
        hurst, σ = res[1], res[2]
    elseif vars == :hurst
        hurst = res[1]
    else
        σ = res[2]
    end

    return (hurst, σ), opm
end

function split_by_day(data::TimeArray)
    res = []
    d0, d1 = Dates.Date(TimeSeries.timestamp(data[1])[1]), Dates.Date(TimeSeries.timestamp(data[end])[1])

    for d in d0:Dates.Day(1):d1
        m0 = Dates.DateTime(d)
        m1 = Dates.DateTime(d + Dates.Hour(23) + Dates.Minute(59))
        mtsp = m0:Dates.Minute(1):m1  # full timestamp
        stsp = intersect(mtsp, TimeSeries.timestamp(data))  # valid timestamp
        if length(stsp) > 0
            push!(res, data[mtsp])
        end
    end

    return res
end

function rolling_regress_predict(regressor::Function, predictor::Function, X0::AbstractVecOrMat{T}, Y0::AbstractVecOrMat{T}, λ::Int, (w,s,d)::Tuple{Int,Int,Int}, p::Int, trans::Function=(x->vec(x)); mode::Symbol=:causal) where {T<:Number}
    X = ndims(X0)>1 ? X0 : reshape(X0, 1, :)  # vec to matrix, create a reference not a copy
    Y = ndims(Y0)>1 ? Y0 : reshape(Y0, 1, :)
    @assert size(X,2) == size(Y,2)

    (any(isnan.(X)) || any(isnan.(Y))) && throw(ValueError("Inputs cannot contain nan values!"))
    L = size(X,2)  # total time
    @assert L >= w >= s

    res_time, res_reg, res_prd = [], [], []  # temporary lists

    if mode == :causal  # causal
        xf, yf = [], []  # future data (due to reversal of time) used for prediction

        for t = L:-p:w
            printfmtln("Processing time {}...\r", t)

            # xv = rolling_apply_hard(trans, view(X, :, max(1, t-w+1):t), s, d; mode=:causal)
            # yv0 = view(Y, :, max(1, t-w+1):t)[:, end:-d:1]  # time-reversed
            # yv = yv0[:,1:size(xv,2)][:,end:-1:1]  # reverse back
            xv = rolling_apply_hard(trans, view(X, :, t-w+1:t), s, d; mode=:causal)
            yv0 = view(Y, :, t-w+1:t)[:, end:-d:1]  # time-reversed
            yv = yv0[:,1:size(xv,2)][:,end:-1:1]  # reverse back

            (xs, ys) = if nan == :ignore
                # ignore columns containing nan values
                idx = findall(.!vec(any(isnan.(vcat(xv, yv)), dims=1)))  # without vec findall will return `CartesianIndex`
                if length(idx) > 0
                    (view(xv,:,idx), view(yv,:,idx))
                else
                    ([], [])
                end
            elseif nan == :zero
                # set all nan to zero
                xv[findall(isnan.(xv))] = 0
                yv[findall(isnan.(yv))] = 0
                (xv, yv)
            else
                (xv, yv)
            end

            if length(xs) > 0 && length(ys) > 0
                reg = regressor(ys, xs)  # regression
                prd = isempty(xf) ? [] : predictor(reg, xf, yf)  # prediction using future data
                pushfirst!(res_time, t)
                pushfirst!(res_reg, reg)
                pushfirst!(res_prd, prd)
            end
            (xf, yf) = if L-t < λ
                ([], [])
            else
                (xv[:,end], yv[:,end])  # yv[:,end]==Y[:,t] is the future truth
            end
        end
        # readjust the prediction
        n = sum([isempty(prd) for prd in res_prd])
        res_prd = circshift(res_prd, n)
    else  # anticausal: TODO
        for t = 1:p:L
            xv = rolling_apply_hard(trans, view(X, :, t:min(L, t+w-1)), s, d; mode=:anticausal)
            yv = rolling_apply_hard(trans, view(Y, :, t:min(L, t+w-1)), 1, d; mode=:anticausal)

            (xs, ys) = if nan == :ignore
                # ignore columns containing nan values
                idx = findall(.!vec(any(isnan.(xv), dims=1) + any(isnan.(yv), dims=1)))  # without vec findall return `CartesianIndex`

                if length(idx) > 0
                    (view(xv,:,idx), view(yv,:,idx))
                else
                    ([], [])
                end
            elseif nan == :zero
                # set all nan to zero
                xv[findall(isnan.(xv))] = 0
                yv[findall(isnan.(yv))] = 0
                (xv, yv)
            else
                (xv, yv)
            end

            if length(xs) > 0 && length(ys) > 0
                # push!(res, (t,regressor(ys, xs)))
                reg = regressor(ys', xs')
                prd = isempty(res_reg) ? [] : predictor(res_reg[end], xv, yv[:,end])
                pushfirst!(res_time, t)
                pushfirst!(res_reg, reg)
                pushfirst!(res_prd, prd)
            end
        end
    end
    res = [(time=t, regression=reg, prediction=prd) for (t, reg, prd) in zip(res_time, res_reg, res_prd)]  # final result is a named tuple
    return res
end

# function moment_incr(X::AbstractVector{T}, d::Integer, p::Real, w::StatsBase.AbstractWeights)
#     mean((abs.(X[d+1:end] - X[1:end-d])).^p, w)
# end

# moment_incr(X, d, p) = moment_incr(X, d, p, StatsBase.weights(ones(length(X)-d)))

"""
Compute the p-th moment of the increment of time-lag `d` of a 1d array.
"""
function moment_incr(X::AbstractVecOrMat{<:Real}, d::Integer, p::Real, k::Int=0)
    dX = X[d+1:end] - X[1:end-d]
    return if k==0  # q is the order of polynomial weight
        mean((abs.(dX)).^p)
    else
        w = StatsBase.weights(causal_weight(length(dX), k))
        mean((abs.(dX)).^p, w)
    end
    # median((abs.(dX)).^p)
end

function moment(X::AbstractVecOrMat{<:Real}, d::Integer, p::Real, k::Int=0)
    w = StatsBase.weights(causal_weight(length(dX), k))
    μX = mean(X, w, dims=1)
    return mean((abs.(X.-μX)).^p, w)
    # return median((abs.(X.-μX)).^p, w)
end

function powlaw_estim(X::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}, p::Real=2.; kt::Integer=1, ks::Integer=0, q::Real=2., ε::Real=1e-2)
    @assert length(lags) > 1 && all(lags .>= 1)
    @assert p > 0.

    cp = 2^(p/2) * gamma((p+1)/2)/sqrt(pi)  # constant depending on p

    # observation and explanatory vectors
    yp = map(d -> log(moment_incr(X, d, p, kt)), lags)
    xp = p * log.(lags)

    # weight for scales
    # wj = length(wj) < length(yp) ? StatsBase.weights(ones(length(yp))/length(yp)) : StatsBase.weights(wj/sum(wj))
    # # @assert all(wj .>= 0)
    w = StatsBase.weights(poly_weight(length(yp), ks))

    yc = yp .- mean(yp, w)
    xc = xp .- mean(xp, w)
    # func = h -> 1/2 * sum(w .* (yc - h*xc).^2)
    func = h -> 1/2 * mean((yc - h*xc).^2, w)
    # func = h -> 1/2 * mean(abs.(yc - h*xc).^q, W)

    # estimation of H and η
    # Gradient-free constrained optimization
    opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())
    # # Gradient-based optimization
    # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
    # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
    hurst = Optim.minimizer(opm)[1]
    η = mean(yp - hurst*xp, w)

    # # by manual inversion
    # Ap = hcat(xp, ones(length(xp))) # design matrix
    # hurst, β = Ap \ yp
    # # or by GLM
    # dg = DataFrames.DataFrame(xvar=xp, yvar=yp)
    # opm = GLM.lm(@GLM.formula(yvar~xvar), dg)
    # η, hurst = GLM.coef(opm)

    σ = exp((η-log(cp))/p)

    return (hurst, σ), opm
end
const fBm_powlaw_estim = powlaw_estim

"""
"""
function powlaw_estim_predict(X::AbstractVector{<:Real}, dlag::Integer, Δ::Integer, p::Real, plen::Integer, plag::Integer; kwargs...)
    (H, σ), opm = FracFin.powlaw_estim(X, dlag-Δ:dlag+Δ, p; kwargs...)
    dX0 = X[dlag+1:end] - X[1:end-dlag]
    # dX = view(dX0, length(dX0)-plen+1:length(dX0))
    pidx = reverse(length(dX0):-dlag:1)[end-plen+1:end]
    dX = dX0[pidx]
    # println(length(dX)+dlag)
    # println(length(X))
    μc, Σc = cond_mean_cov(FractionalGaussianNoise(H, dlag), length(dX)+1:length(dX)+plag, pidx, dX)
    return H, σ, μc, σ^2 * Σc
end

"""TODO
Multiscale fGn-MLE
"""
function ms_fGn_MLE_estim(X::AbstractVector{T}, lags::AbstractVector{Int}, w::Int) where {T<:Real}
    Hs = zeros(length(lags))
    Σs = zeros(length(lags))

    for (n,lag) in enumerate(lags)  # time lag for finite difference
        # vectorization with window size w
        dXo = rolling_vectorize(X[lag+1:end]-X[1:end-lag], w, 1, 1)
        # rolling mean with window size 2lag, then down-sample at step lag
        dX = rolling_mean(dXo, 2lag, lag; boundary=:hard)

        (hurst_estim, σ_estim), obj = fGn_MLE_estim(squeezedims(dX), lag)

        Hs[n] = hurst_estim
        Σs[n] = σ_estim
    end

    return Hs, Σs
end

##### B-Spline DCWT MLE (Not maintained) #####
# Implementation based on DCWT formulation, not working well in practice.

function fBm_bspline_covmat_lag(H::Real, v::Int, l::Int, sclrng::AbstractVector{Int}, mode::Symbol)
    return Amat_bspline(H, v, l, sclrng) .* [sqrt(i*j) for i in sclrng, j in sclrng].^(2H+1)
end

"""
Compute the covariance matrix of B-Spline DCWT coefficients of a pure fBm.

The full covariance matrix of `J`-scale transform and of time-lag `N` is a N*J-by-N*J symmetric matrix.

# Args
- l: maximum time-lag
- sclrng: scale range
- v: vanishing moments of B-Spline wavelet
- H: Hurst exponent
- mode: mode of convolution
"""
function fBm_bspline_covmat(l::Int, sclrng::AbstractVector{Int}, v::Int, H::Real, mode::Symbol)
    J = length(sclrng)
    Σ = zeros(((l+1)*J, (l+1)*J))
    Σs = [fBm_bspline_covmat_lag(H, v, d, sclrng, mode) for d = 0:l]

    for r = 0:l
        for c = 0:l
            Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (c>=r) ? Σs[c-r+1] : transpose(Σs[r-c+1])
        end
    end

    return Matrix(Symmetric(Σ))  #  forcing symmetry
    # return [(c>=r) ? Σs[c-r+1] : Σs[r-c+1]' for r=0:N-1, c=0:N-1]
end

"""
Evaluate the log-likelihood of B-Spline DCWT coefficients.
"""
function fBm_bspline_log_likelihood_H(X::AbstractVecOrMat{T}, sclrng::AbstractVector{Int}, v::Int, H::Real, mode::Symbol) where {T<:Real}
    @assert 0 < H < 1
    @assert size(X,1) % length(sclrng) == 0

    L = size(X,1) ÷ length(sclrng)  # integer division: \div
    # N = ndims(X)>1 ? size(X,2) : 1

    Σ = fBm_bspline_covmat(L-1, sclrng, v, H, mode)  # full covariance matrix

    # # strangely, the following does not work (logarithm of a negative value)
    # iΣ = pinv(Σ)  # regularization by pseudo-inverse
    # return -1/2 * (J*N*log(trace(X'*iΣ*X)) + logdet(Σ))

    return log_likelihood_H(Σ, X)
end

"""
B-Spline wavelet-MLE estimator.
"""
function fBm_bspline_DCWT_MLE_estim(X::AbstractVecOrMat{T}, sclrng::AbstractVector{Int}, v::Int, mode::Symbol; method::Symbol=:optim, ε::Real=1e-2) where {T<:Real}
    @assert size(X,1) % length(sclrng) == 0
    # number of wavelet coefficient vectors concatenated into one column of X
    L = size(X,1) ÷ length(sclrng)  # integer division: \div
    # N = ndims(X)>1 ? size(X,2) : 1

    func = x -> -fBm_bspline_log_likelihood_H(X, sclrng, v, x, mode)

    opm = nothing
    hurst = nothing

    if method == :optim
        # Gradient-free constrained optimization
        opm = Optim.optimize(func, ε, 1-ε, Optim.Brent())
        # # Gradient-based optimization
        # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
        # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
        hurst = Optim.minimizer(opm)[1]
    elseif method == :table
        Hs = collect(ε:ε:1-ε)
        hurst = Hs[argmin([func(h) for h in Hs])]
    else
        throw("Unknown method: ", method)
    end

    Σ = fBm_bspline_covmat(L-1, sclrng, v, hurst, mode)
    σ = sqrt(xiAx(Σ, X) / length(X))

    return (hurst, σ), opm
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminInviteRequestsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminInviteRequestsApprove(_api::AdminInviteRequestsApi, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.inviteRequests.approve", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "invite_request_id", invite_request_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Approve a workspace invite request.
Param: invite_request_id::String (required)
Param: team_id::String
Return: Dict{String, Any}
"""
function adminInviteRequestsApprove(_api::AdminInviteRequestsApi, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsApprove(_api, invite_request_id; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminInviteRequestsApprove(_api::AdminInviteRequestsApi, response_stream::Channel, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsApprove(_api, invite_request_id; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminInviteRequestsDeny(_api::AdminInviteRequestsApi, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.inviteRequests.deny", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "invite_request_id", invite_request_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Deny a workspace invite request.
Param: invite_request_id::String (required)
Param: team_id::String
Return: Dict{String, Any}
"""
function adminInviteRequestsDeny(_api::AdminInviteRequestsApi, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsDeny(_api, invite_request_id; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminInviteRequestsDeny(_api::AdminInviteRequestsApi, response_stream::Channel, invite_request_id::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsDeny(_api, invite_request_id; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminInviteRequestsList(_api::AdminInviteRequestsApi; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.inviteRequests.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

List all pending workspace invite requests.
Param: team_id::String
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function adminInviteRequestsList(_api::AdminInviteRequestsApi; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsList(_api; team_id=team_id, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminInviteRequestsList(_api::AdminInviteRequestsApi, response_stream::Channel; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsList(_api; team_id=team_id, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminInviteRequestsApprove, adminInviteRequestsDeny, adminInviteRequestsList