text:
function reduce_subtrees!(v,r,reduce)
    @assert length(r) >= length(v)
    n = length(v)
    if n == 0; return nothing; end

    @inbounds begin
        # Reduce no-child nodes
        for i = n:-1:heapparent(n)+1
            r[i] = reduce(v[i])
        end

        # Reduce one-child node, if present
        for i = heapparent(n+1) : heapparent(n)
            r[i] = reduce(v[i],r[n])
        end

        # Reduce two-children nodes
        for i = heapparent(n)-1:-1:1
            r[i] = reduce(v[i], r[heapleft(i)], r[heapright(i)])
        end
    end
    return nothing
end

function Reducer!(v,r,reduce)
    function ReducerF(i, _, nc)
        if nc == 0; r[i] = reduce(v[i]); end
        if nc == 1; r[i] = reduce(v[i], r[heapleft(i)]); end
        if nc == 2; r[i] = reduce(v[i], r[heapleft(i)], r[heapright(i)]); end
        return :continue
    end
    return ReducerF
end


text:
module HeatEqCalcOld

function laplacian4!(v, u, dx)
    v[begin] = (u[end] + u[begin+1] - 2u[begin])/dx^2
    @. @views v[begin+1:end-1] = (u[begin:end-2] + u[begin+2:end] - 2u[begin+1:end-1])/dx^2
    v[end] = (u[end-1] + u[begin] - 2u[1end])/dx^2
    return
end

function heateq4(u0, dx, tmax, N=200)
    t = 0:dx:tmax
    dt = step(t)
    u = similar(u0, length(u0), length(t)+1)
    u[:, 1] = u0
    v = similar(u0)
    N = 100
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian4!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

using StaticArrays

function calc_sol(
        n = 200,
        x = range(-π, π; length=n+1)[1:end-1],
        m = 200
    )
    dx = step(x)
    noise1, noise2 = 0.3randn(m), 0.3randn(m)
    f(x) = SVector{m}((1 .+ noise1)*sin(x) .- (1 .+ noise2)*cos(2x))
    u0 = f.(x)
    tmax = 1.0
    t, u = heateq4(u0, dx, tmax)
    (; t, x, u, m)
end

save_sol(sol) = write("sol.txt", repr(sol))

end # module

ext: jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
export solve_modelKTR
function solve_modelKTR(H, g, nlp_stop, X::PDataKTR, α::T) where {T}
    # target value should be close to satisfy α=||d||
    target = [(abs(α - X.norm_dirs[i])) for i = 1:X.nshifts]
    # pick the closest shift to the target within positive definite H+λI
    # before, check the shift = 0 for direction within the trust region

    positives = findall(X.positives)
    if (X.indmin == 0)  #  try Newton's direction
        X.indmin = 1
        if (positives[1] == 1) & (X.norm_dirs[1] <= α)
            X.d = X.xShift[1]
            X.λ = 0.0
            return X.d, X.λ
        end
    end
    start = X.indmin
    bidon, indmin = findmin(target[positives[start:end]])
    X.indmin = start + indmin - 1
    p_imin = positives[X.indmin]
    X.d = X.xShift[p_imin]
    # Refine this rough constrained direction
    X.λ = X.shifts[p_imin]

    return X.d, X.λ
end

# To replace.
#    1. Form the signed ( α - X.norm_dirs[i]) ) within positives of the form +++...----
#    2. If starts with a +, interpolate the interval where sign changes
#       else (hard case) obtain from Krylov shift solver negative curvature directions
#            and extrapolate to the boundary the short (longest) direction
#            in the descent negative curvature direction
#
# Keep on treating λ=0 special for Newton's direction (no extrapolation to boundary)


text:
using LinearAlgebra: checksquare

"""
    vtype(A) → T

Do some arithmetic to get a proper number type that a matrix A should operate
on.
"""
function vtype(A)
    T = eltype(A)
    typeof(zero(T)/sqrt(one(T)))
end

"""
```julia
partialschur(A; nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History
```

Find `nev` approximate eigenpairs of `A` with eigenvalues near a specified target.

The matrix `A` can be any linear map that implements `mul!(y, A, x)`, `eltype`
and `size`.

The method will run iteratively until the eigenpairs are approximated to
the prescribed tolerance or until `restarts` restarts have passed.

## Arguments

The most important keyword arguments:

| Keyword | Type | Default | Description |
|------:|:-----|:----|:------|
| `nev` | `Int` | `min(6, size(A, 1))` |Number of eigenvalues |
| `which` | `Target` | `LM()` | One of `LM()`, `LR()`, `SR()`, `LI()`, `SI()`, see below. |
| `tol` | `Real` | `√eps` | Tolerance for convergence: ‖Ax - xλ‖₂ < tol * ‖λ‖ |

The target `which` can be any of `subtypes(ArnoldiMethod.Target)`:

| Target | Description |
|------:|:-----|
| `LM()` | Largest magnitude: `abs(λ)` is largest |
| `LR()` | Largest real part: `real(λ)` is largest |
| `SR()` | Smallest real part: `real(λ)` is smallest |
| `LI()` | Largest imaginary part: `imag(λ)` is largest|
| `SI()` | Smallest imaginary part: `imag(λ)` is smallest|

!!! note

    The targets `LI()` and `SI()` only make sense in complex arithmetic. In real
    arithmetic `λ` is an eigenvalue iff `conj(λ)` is an eigenvalue and this
    conjugate pair converges simultaneously.

## Return values

The function returns a tuple

```julia
decomp, history = partialschur(A, ...)
```

where `decomp` is a [`PartialSchur`](@ref) struct which
forms a partial Schur decomposition of `A` to a prescribed tolerance:

```julia
> norm(A * decomp.Q - decomp.Q * decomp.R)
```

`history` is a [`History`](@ref) struct that holds some basic information about
convergence of the method:

```julia
> history.converged
true
> @show history
Converged after 359 matrix-vector products
```

## Advanced usage

Further there are advanced keyword arguments for tuning the algorithm:

| Keyword | Type | Default | Description |
|------:|:-----|:---|:------|
| `mindim` | `Int` | `min(max(10, nev), size(A,1))` | Minimum Krylov dimension (≥ nev) |
| `maxdim` | `Int` | `min(max(20, 2nev), size(A,1))` | Maximum Krylov dimension (≥ min) |
| `restarts` | `Int` | `200` | Maximum number of restarts |

When the algorithm does not converge, one can increase `restarts`. When the
algorithm converges too slowly, one can play with `mindim` and `maxdim`. It is
suggested to keep `mindim` equal to or slightly larger than `nev`, and `maxdim`
is usually about two times `mindim`.

"""
function partialschur(A;
                       nev::Int = min(6, size(A, 1)),
                       which::Target = LM(),
                       tol::Real = sqrt(eps(real(vtype(A)))),
                       mindim::Int = min(max(10, nev), size(A, 1)),
                       maxdim::Int = min(max(20, 2nev), size(A, 1)),
                       restarts::Int = 200)
    s = checksquare(A)
    if nev < 1
        throw(ArgumentError("nev cannot be less than 1"))
    end
    nev ≤ mindim ≤ maxdim ≤ s || throw(ArgumentError("nev ≤ mindim ≤ maxdim does not hold, got $nev ≤ $mindim ≤ $maxdim"))
    _partialschur(A, vtype(A), mindim, maxdim, nev, tol, restarts, which)
end

"""
    IsConverged(ritz, tol)

Functor to test whether Ritz values satisfy the convergence criterion. Current
convergence condition is ‖Ax - xλ‖₂ < max(ε‖H‖, tol * |λ|). This is supposed to
be scale invariant: the matrix `B = αA` for some constant `α` has the same
eigenvectors with eigenvalue λα, so this scaling with `α` cancels in the
inequality.
"""
struct IsConverged{RV<:RitzValues,T}
    ritz::RV
    tol::T
    H_frob_norm::RefValue{T}

    IsConverged(ritz::R, tol::T) where {R<:RitzValues,T} = new{R,T}(ritz, tol, RefValue(zero(T)))
end

(r::IsConverged{RV,T})(i::Integer) where {RV,T} =
    @inbounds return r.ritz.rs[i] < max(eps(T) * r.H_frob_norm[], r.tol * abs(r.ritz.λs[i]))

"""
    History(mvproducts, nconverged, converged, nev)

History shows whether the method has converged (when `nconverged` ≥ `nev`) and
how many matrix-vector products were necessary to do so.
"""
struct History
    mvproducts::Int
    nconverged::Int
    converged::Bool
    nev::Int
end

function _partialschur(A, ::Type{T}, mindim::Int, maxdim::Int, nev::Int, tol::Ttol, restarts::Int, which::Target) where {T,Ttol<:Real}
    n = size(A, 1)

    # Pre-allocated Arnoldi decomp
    arnoldi = Arnoldi{T}(n, maxdim)

    # Unpack for convenience
    H = arnoldi.H
    V = arnoldi.V

    # For a change of basis we have Vtmp as working space
    Vtmp = Matrix{T}(undef, n, maxdim)

    # Unitary matrix used for change of basis of V.
    Q = Matrix{T}(undef, maxdim, maxdim)

    # We only need to store one eignvector of the Hessenberg matrix
    x = zeros(complex(T), maxdim)

    # And we store the reflector to transform H back to Hessenberg separately
    G = Reflector{T}(maxdim)

    # Approximate residual norms for all Ritz values, and Ritz values
    ritz = RitzValues{T}(maxdim)
    isconverged = IsConverged(ritz, tol)
    ordering = get_order(which)
    groups = zeros(Int, maxdim)

    # V[:,1:active-1   ] ⸺ Locked vectors / approximately invariant subspace
    # V[:,active:maxdim] ⸺ Active part of decomposition.
    # V[:,active:k     ] ⸺ Smallest size of Arnoldi relation; k will increase over time.
    #                        We keep length(active:k) ≈ mindim, but we should also retain
    #                        space to add new directions to the Krylov subspace, so
    #                        length(k+1:maxdim) should not be too small either.
    active = 1
    k = mindim
    effective_nev = nev

    # Bookkeeping for number of mv-products
    prods = mindim

    # Initialize an Arnoldi relation of size `mindim`
    reinitialize!(arnoldi)
    iterate_arnoldi!(A, arnoldi, 1:mindim)

    for iter = 1 : restarts

        # Expand Krylov subspace dimension from `k` to `maxdim`.
        iterate_arnoldi!(A, arnoldi, k+1:maxdim)

        # Bookkeeping
        prods += length(k+1:maxdim)

        # Q accumulates the changes of basis via relfectors; initially it's just I.
        copyto!(Q, I)

        # Construct Schur decomposition of H[active:maxdim,active:maxdim] in-place
        local_schurfact!(view(H, OneTo(maxdim), :), active, maxdim, Q)

        # Update the Ritz values
        copyto!(ritz.ord, OneTo(maxdim))
        copy_eigenvalues!(ritz.λs, H)
        copy_residuals!(ritz.rs, H, Q, H[maxdim+1,maxdim], x, active:maxdim)

        # Create a permutation that sorts Ritz values from most wanted to least wanted
        sort!(ritz.ord, 1, maxdim, QuickSort, OrderPerm(ritz.λs, ordering))

        # Compute the Frobenius norm of H for the stopping criterion
        isconverged.H_frob_norm[] = norm(H)

        ### PARTITIONING OF SCHUR FORM IN [LOCKED | RETAINED | PURGED]

        # Plan de campagne: reorder the permutation `ritz.ord` such that
        # ritz.ord[1:nlock] gives the indices of the locked ritz values
        # ritz.ord[nlock+1:k] is the indices of the ritz values we wanna retain
        # ritz.ord[k+1:maxdim] is the indices of the ritz values we wanna truncate
        # Then we make ritz.groups[i] = {1,2,3} by iterating over them.

        # We keep at most `nev` or `nev+1` eigenvalues, depending on the split being
        # halfway a conjugate pair.
        effective_nev = include_conjugate_pair(T, ritz, nev)

        # Partition in converged & not converged.
        first_not_conv_idx = partition!(isconverged, ritz.ord, active:effective_nev)

        # Now ritz.ord[1:nlock] are converged eigenvalues that we want to lock, and
        # nlock ≤ effective_nev, so it's really just these that we are after!
        nlock = first_not_conv_idx === nothing ? effective_nev : first_not_conv_idx - 1

        # Next, purge the converged eigenvalues we do not want by moving them to the back.
        partition!(i -> !isconverged(i), ritz.ord, nlock+1:maxdim)

        # Determine the new length `k` of the truncated Krylov subspace:
        # 1. The dimension of the active part should be roughly `mindim`; so `k` will be
        #    larger than `mindim` when converged Ritz vectors have been locked.
        # 2. But `k` can't be so large that the expansion would barely give new information;
        #    hence we meet in the middle: `k` is at most halfway `mindim` and `maxdim`.
        # 3. If `k` ends up on the boundary of a conjugate pair, we increase `k` by 1.
        k = include_conjugate_pair(T, ritz, min(nlock + mindim, (mindim + maxdim) ÷ 2))

        # Locked ritz values
        @inbounds for i = 1:nlock
            groups[ritz.ord[i]] = 1
        end

        # Retained ritz values
        @inbounds for i = nlock+1:k
            groups[ritz.ord[i]] = 2
        end

        # Truncated ritz values
        @inbounds for i = k+1:maxdim
            groups[ritz.ord[i]] = 3
        end

        partition_schur_three_way!(H, Q, groups)

        # Restore the Hessenberg matrix via Householder reflections.
        # Note that we restore the new active part only -- Q[end, 1:nlock] is small enough
        # by convergence criterion.
        restore_arnoldi!(H, nlock + 1, k, Q, G)

        # Finally do the change of basis to get the length `k` Arnoldi relation.
        @views mul!(Vtmp[:,active:k], V[:,active:maxdim], Q[active:maxdim,active:k])
        @views copyto!(V[:,active:k], Vtmp[:,active:k])
        @views copyto!(V[:,k+1], V[:,maxdim+1])

        # The active part
        active = nlock + 1

        active > nev && break
    end

    nconverged = active - 1

    @views Vconverged = V[:,1:nconverged]
    @views Hconverged = H[1:nconverged,1:nconverged]

    # Sort the converged eigenvalues like the user wants them
    sortschur!(H, copyto!(Q, I), nconverged, ordering)

    # Change of basis
    @views mul!(Vtmp[:,1:nconverged], Vconverged, Q[1:nconverged,1:nconverged])
    @views copyto!(Vconverged, Vtmp[:,1:nconverged])

    # Copy the eigenvalues just one more time
    copy_eigenvalues!(ritz.λs, H, OneTo(nconverged))

    history = History(prods, nconverged, nconverged ≥ nev, nev)
    schur = PartialSchur(Vconverged, Hconverged, ritz.λs[1:nconverged])

    return schur, history
end

function partition_schur_three_way!(R, Q, groups::AbstractVector{Int})
    # Partitioning goes like this:
    # |2|33|1|2|33|2|1|
    #  h
    #  m
    #  l

    # |2|33|1|2|33|2|1|  -> |2|33|1|2|33|2|1|
    #    h               ->       h
    #  m                 ->    m
    #  l                 ->  l

    # |2|33|1|2|33|2|1|  -> |1|2|33|2|33|2|1|
    #       h            ->         h
    #    m               ->      m
    #  l                 ->    l

    # |1|2|33|2|33|2|1|  -> |1|2|2|33|33|2|1|
    #         h          ->           h
    #      m             ->        m
    #    l               ->    l

    # |1|2|2|33|33|2|1|  -> |1|2|2|33|33|2|1|
    #           h        ->              h
    #        m           ->        m
    #    l               ->    l

    # |1|2|2|33|33|2|1|  -> |1|2|2|2|33|33|1|
    #              h     ->                h
    #        m           ->          m
    #    l               ->    l

    # |1|2|2|2|33|33|1|  -> |1|1|2|2|2|33|33|
    #                h   ->                  h
    #          m         ->            m
    #    l               ->      l

    hi = 1
    mi = 1
    lo = 1

    # Be very scared whenever eigenvalues are too close!
    while hi ≤ length(groups)
        group = groups[hi]
        blocksize = is_start_of_11_block(R, hi) ? 1 : 2

        if group == 3
            hi += blocksize
        elseif group == 2
            rotate_right!(R, mi, hi, Q)
            hi += blocksize
            mi += blocksize
        else # group == 1
            rotate_right!(R, lo, hi, Q)
            hi += blocksize
            mi += blocksize
            lo += blocksize
        end
    end

    nothing
end

"""
    sortschur!(R, Q, to::Int, ordering::Ordering)

Reorder the elements on the diagonal of R as indicates by `ordering`.
Eigenvalues are computed from R.
"""
function sortschur!(R, Q, to::Int, ordering::Ordering)
    # Some convoluted insertion sort algorithm
    # `next` is the index of the eigenvalue to be put in position
    # `curr` is the eigenvalue as it is being moved to the front of R
    # `prev` is the index of the block prior to `curr`.

    to ≤ 1 && return

    next_idx = 1

    while next_idx ≤ to
        # Points to the start of the block we're swapping
        curr_idx = next_idx

        # Get the current eigenvalue and remember whether we're on a block
        # so that we skip either 1 or 2 steps.
        curr_size = is_start_of_11_block(R, curr_idx) ? 1 : 2
        curr_λ = eigenvalue(R, curr_idx)

        # Insertion part of the sort
        while curr_idx > 1

            # Identify previous block size (1x1 or 2x2)
            prev_size = is_end_of_11_block(R, curr_idx - 1) ? 1 : 2
            prev_idx = curr_idx - prev_size
            prev_λ = eigenvalue(R, prev_idx)

            # Maybe we're done.
            lt(ordering, curr_λ, prev_λ) === false && break

            # Maybe swap.
            swap!(R, prev_idx, prev_size == 1, curr_size == 1, Q)
            curr_idx -= prev_size
        end

        next_idx += curr_size
    end
end

"""
    include_conjugate_pair(T, ritz, i) → {i, i + 1}

Returns i or i + 1 depending on whether Ritz value i and i + 1 form a conjugate pair
together
"""
@inline function include_conjugate_pair(::Type{<:Real}, ritz::RitzValues, i)
    i >= length(ritz.ord) && return i
    @inbounds λ1 = ritz.λs[ritz.ord[i+0]]
    @inbounds λ2 = ritz.λs[ritz.ord[i+1]]
    return imag(λ1) != 0 && λ1' == λ2 ? i + 1 : i
end

@inline include_conjugate_pair(::Type{T}, ritz::RitzValues, i) where {T} = i

"""
    update_residual_norms!(rs, H, Q, hₖ₊₁ₖ) → rs

Computes the Ritz residuals ‖Ax - λx‖₂ = |yₖ| * |hₖ₊₁ₖ| for each eigenvalue
"""
function copy_residuals!(rs::AbstractVector{T}, H, Q, hₖ₊₁ₖ, x::AbstractVector, range::AbstractRange) where {T<:Real}
    fill!(rs, 0)
    m = size(H, 2)
    @inbounds for i = range
        fill!(x, zero(T))
        len = collect_eigen!(x, H, i)
        tmp = zero(complex(T))
        for j = 1 : len
            tmp += Q[m, j] * x[j]
        end
        rs[i] = abs(tmp * hₖ₊₁ₖ)
    end

    rs
end


text:
@testset "test_synthesis" begin
P = tf(1.,[1.,1])
C = tf([1.,1],[1.,0])
L = P*C
Lsys = ss(L)

B = [1]
A = [1,1]
R = [1,1]
S = [1]
T = [1]
@testset "minreal + feedback" begin
@test isapprox(minreal(feedback(P,C),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5)
@test isapprox(numpoly(minreal(feedback(L),1e-5))[1].coeffs, numpoly(tf(1,[1,1]))[1].coeffs)# This test is ugly, but numerical stability is poor for minreal
@test feedback2dof(B,A,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
@test feedback2dof(P,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)])
@test isapprox(poles(minreal(tf(feedback(Lsys)),1e-5)) , poles(minreal(feedback(L),1e-5)), atol=1e-5)

Pint = tf(1,[1,1])
Cint = tf([1,1],[1,0])
Lint = P*C

@test isapprox(minreal(feedback(Pint,Cint),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5) # TODO consider keeping minreal of Int system Int
@test isapprox(numpoly(minreal(feedback(Lint),1e-5))[1].coeffs, numpoly(tf(1,[1,1]))[1].coeffs)# This test is ugly, but numerical stability is poor for minreal
@test isapprox(poles(minreal(tf(feedback(Lsys)),1e-5)) , poles(minreal(feedback(L),1e-5)), atol=1e-5)

@test feedback(ss(1),ss(1)) == ss(0.5)
@test_throws ErrorException feedback(ss([1 0; 0 1], ones(2,2), ones(1,2),0))

# Test Feedback Issue: 163
g = tf([1],[1,1])
gfb = feedback(g)
gfb2 = tf(feedback(ss(g)))
@test hinfnorm(gfb - gfb2)[1] <= 1e-14

# Test more feedback
s = tf("s")
ftf = 1.0*(2s+3)/((5s+7)*(11s+13))
fzpk = zpk(ftf)

ffb = feedback(fzpk)            # Zpk feedback
ffb2 = zpk(feedback(ss(fzpk)))  # ss feedback + ss conversion
ffb3 = feedback(ftf)            # tf feedback
ffb4 = feedback(ss(ftf))        # ss feedback
ffb5 = minreal(fzpk/(1+fzpk))   # Zpk feedback manual

z1,p1,k1 = zpkdata(ffb)
z2,p2,k2 = zpkdata(ffb2)
z3,p3,k3 = zpkdata(ffb3)
z4,p4,k4 = zpkdata(ffb4)
z5,p5,k5 = zpkdata(ffb5)
@test sort(real.(z1[1])) ≈ sort(real.(z2[1])) ≈ sort(real.(z3[1])) ≈ sort(real.(z4[1])) ≈ sort(real.(z5[1]))
@test sort(real.(p1[1])) ≈ sort(real.(p2[1])) ≈ sort(real.(p3[1])) ≈ sort(real.(p4[1])) ≈ sort(real.(p5[1]))
@test k1 ≈ k2 ≈ k3 ≈ k4 ≈ k5
end

@testset "place" begin
    sys = ss(-4, 2, 3, 0)
    A, B, C, _ = ssdata(sys)

    @test place(A, B, [-10]) == [3][:,:]
    @test place(A, B, [-10], :c) == [3][:,:]
    @test place(A, C, [-10], :o) == [2][:,:]

    A = [0 1; 0 0]
    B = [0; 1]
    C = [1 0]
    sys = ss(A, B, C, 0)

    @test place(A, B, [-1.0, -1]) ≈ [1 2]
    @test place(sys, [-1.0, -1]) ≈ [1 2]
    @test place(A, B, [-1.0, -1], :c) ≈ [1 2]
    @test place(sys, [-1.0, -1], :c) ≈ [1 2]
    @test place(A, C, [-2.0, -2], :o) ≈ [4; 4]
    @test place(sys, [-2.0, -2], :o) ≈ [4; 4]

    @test place(A, B, [-2 + im, -2 - im]) ≈ [5 4]
    @test place(A, C, [-4 + 2im, -4 - 2im], :o) ≈ [8; 20]

    A = ones(3,3) - diagm([3, 4, 5])
    B = [1; 0; 2]
    C = [1 1 0]
    @test place(A, B, [-2 + 2im, -2 - 2im, -4]) ≈ [-2.6 5.2 0.8]
    @test place(A, C, [-2 + 3im, -2 - 3im, -4], :o) ≈ [11; -12; 1]
end

@testset "acker" begin
Random.seed!(0)
A = randn(3,3)
B = randn(3,1)
p = [3.0,2,1]
K = ControlSystems.acker(A,B,p)
@test ControlSystems.eigvalsnosort(A-B*K) ≈ p

p = [-1+im, -1-im, -1]
K = ControlSystems.acker(A,B,p)
@test ControlSystems.eigvalsnosort(A-B*K) ≈ p
end

@testset "LQR" begin
    Ts = 0.1
    A = [1 Ts; 0 1]
    B = [0;1] # Note B is vector, B'B is scalar, but compatible with I
    C = [1 0]
    Q = I
    R = I
    L = dlqr(A,B,Q,R)
    @test L ≈ [0.5890881713787511 0.7118839434795103]
    sys = ss(A,B,C,0,Ts)
    L = lqr(sys, Q, R)
    @test L ≈ [0.5890881713787511 0.7118839434795103]

    L = dlqr(sys, Q, R)
    @test L ≈ [0.5890881713787511 0.7118839434795103]

    B = reshape(B,2,1)  # Note B is matrix, B'B is compatible with I
    L = dlqr(A,B,Q,R)
    @test L ≈ [0.5890881713787511 0.7118839434795103]

    Q = eye_(2)
    R = eye_(1)
    L = dlqr(A,B,Q,R)
    @test L ≈ [0.5890881713787511 0.7118839434795103]

    B = [0;1]   # Note B is vector, B'B is scalar
    Q = eye_(2)
    R = eye_(1)
    L ≈ dlqr(A,B,Q,R)
    #L ≈ [0.5890881713787511 0.7118839434795103]
end

end


text:
@testset "util" begin
    import DiffLinearAlgebra: importable
    @test Expr(:import, importable(:Foo)...) == :(import Foo)
    @test Expr(:import, importable(:(Package.Foo))...) == :(import Package.Foo)
    @test Expr(:import, importable(:(Package.Subpackage.Foo))...) ==
        :(import Package.Subpackage.Foo)

    @test import_expr(DLA.DiffOp(:Foo, :(Tuple{Foo}), [true])) ==
        :(import Foo)
    @test import_expr(DLA.DiffOp(:(Package.Foo), :(Tuple{Foo}), [true])) ==
        :(import Package.Foo)
    @test import_expr(DLA.DiffOp(:(Package.Subpackage.Foo), :(Tuple{Foo}), [true])) ==
        :(import Package.Subpackage.Foo)
end


text:
load("core.jl")

srand(1)
means = [0.1, 0.1, 0.1, 0.1, 0.9]
n_arms = length(means)
means = means[randperm(n_arms)]
arms = map(mu -> BernoulliArm(mu), means)
println("Best arm is $(ind_max(means))")

total_results = zeros(0, 5)

algo = AnnealingEpsilonGreedy(zeros(Int64, n_arms), zeros(n_arms))
results = test_algorithm(algo, arms, 5000, 250)
total_results = vcat(total_results, results)

csvwrite("algorithms/epsilon_greedy/annealing_results.csv", total_results)

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
# Julia wrappers around the NumPy API, as part of the PyCall package

#########################################################################
# Initialization (UGLY)

# Linking NumPy's C API from Julia requires some serious hackery,
# because NumPy does not export its symbols in the usual way for
# shared libraries.  Instead, it provides a Python variable
# numpy.core.multiarray._ARRAY_API that points to a lookup table of
# pointers to the API functions and global variables.  Moreover, the
# meaning of this table, along with a lot of other important
# constants, is defined in a C header file which changes between NumPy
# versions, so we need to do some regex parsing of this header file in
# order to extract the necessary information to call NumPy.  Ugly and
# mildly insane, but I don't see much alternative (at least to do this
# purely in Julia).
#
# The result of npy_api_initialize, below, is to produce the following
# tables of API pointers:
const npy_api = Dict{Symbol, Ptr{Cvoid}}()

# need a global to cache pyimport("numpy.core.multiarray"), in order
# to ensure the module is not garbage-collected as long as we are using it
# for the npy_api pointers.
const npy_multiarray = PyNULL()

npy_initialized = false # global to prevent multiple initializations

# Macro version of npyinitialize() to inline npy_initialized? check
macro npyinitialize()
    :(npy_initialized::Bool ? nothing : npyinitialize())
end

# numpy.number types, used to detect scalars for conversion routines
const npy_number = PyNULL()
const npy_integer = PyNULL()
const npy_floating = PyNULL()
const npy_complexfloating = PyNULL()
const npy_bool = PyNULL()

import LinearAlgebra # for MKL check

function npyinitialize()
    global npy_initialized

    if npy_initialized::Bool
        return
    end
    try
        copy!(npy_multiarray, pyimport("numpy.core.multiarray"))
    catch e
        error("numpy.core.multiarray required for multidimensional Array conversions - ", e)
    end
    if pyversion.major < 3
        PyArray_API = @pycheck ccall((@pysym :PyCObject_AsVoidPtr),
                                     Ptr{Ptr{Cvoid}}, (PyPtr,),
                                     npy_multiarray."_ARRAY_API")
    else
        PyArray_API = @pycheck ccall((@pysym :PyCapsule_GetPointer),
                                     Ptr{Ptr{Cvoid}}, (PyPtr,Ptr{Cvoid}),
                                     npy_multiarray."_ARRAY_API", C_NULL)
    end

    numpy = pyimport("numpy")

    # emit a warning if both Julia and NumPy are linked with MKL (#433)
    if LinearAlgebra.BLAS.vendor() === :mkl &&
       LinearAlgebra.BLAS.BlasInt === Int64 && hasproperty(numpy, "__config__")
        config = numpy."__config__"
        if hasproperty(config, "blas_opt_info")
            blaslibs = get(config."blas_opt_info", Vector{String}, "libraries", String[])
            if any(s -> occursin("mkl", lowercase(s)), blaslibs)
                @warn "both Julia and NumPy are linked with MKL, which may cause conflicts and crashes (#433)."
            end
        end
    end

    # directory for numpy include files to parse
    inc = pycall(numpy."get_include", AbstractString)

    # numpy.number types
    copy!(npy_number, numpy."number")
    copy!(npy_integer, numpy."integer")
    copy!(npy_floating, numpy."floating")
    copy!(npy_complexfloating, numpy."complexfloating")
    copy!(npy_bool, numpy."bool_")

    # Parse __multiarray_api.h to obtain length and meaning of PyArray_API
    try
        hdrfile = open(joinpath(inc, "numpy", "__multiarray_api.h"))
        hdr = read(hdrfile, String)
        close(hdrfile)
    catch e
        error("could not read __multiarray_api.h to parse PyArray_API ", e)
    end
    hdr = replace(hdr, r"\\\s*\n"=>" "); # rm backslashed newlines
    r = r"^#define\s+([A-Za-z]\w*)\s+\(.*\bPyArray_API\s*\[\s*([0-9]+)\s*\]\s*\)\s*$"m # regex to match #define PyFoo (... PyArray_API[nnn])
    PyArray_API_length = 0
    for m in eachmatch(r, hdr) # search for max index into PyArray_API
        PyArray_API_length = max(PyArray_API_length, parse(Int, m.captures[2])+1)
    end
    API = unsafe_wrap(Array, PyArray_API, (PyArray_API_length,))
    for m in eachmatch(r, hdr) # build npy_api table
        npy_api[Symbol(m.captures[1])] = API[parse(Int, m.captures[2])+1]
    end
    if !haskey(npy_api, :PyArray_New)
        error("failure parsing NumPy PyArray_API symbol table")
    end

    npy_initialized = true
    return
end

#########################################################################
# Hard-coded constant values, copied from numpy/ndarraytypes.h ...
# the values of these seem to have been stable for some time, and
# the NumPy developers seem to have some awareness of binary compatibility

# NumPy Types:

const NPY_BOOL = Int32(0)
const NPY_BYTE = Int32(1)
const NPY_UBYTE = Int32(2)
const NPY_SHORT = Int32(3)
const NPY_USHORT = Int32(4)
const NPY_INT = Int32(5)
const NPY_UINT = Int32(6)
const NPY_LONG = Int32(7)
const NPY_ULONG = Int32(8)
const NPY_LONGLONG = Int32(9)
const NPY_ULONGLONG = Int32(10)
const NPY_FLOAT = Int32(11)
const NPY_DOUBLE = Int32(12)
const NPY_LONGDOUBLE = Int32(13)
const NPY_CFLOAT = Int32(14)
const NPY_CDOUBLE = Int32(15)
const NPY_CLONGDOUBLE = Int32(16)
const NPY_OBJECT = Int32(17)
const NPY_STRING = Int32(18)
const NPY_UNICODE = Int32(19)
const NPY_VOID = Int32(20)
const NPY_HALF = Int32(23)

#########################################################################
# conversion from Julia types to NumPy types

npy_type(::Type{Bool}) = NPY_BOOL
npy_type(::Type{Int8}) = NPY_BYTE
npy_type(::Type{UInt8}) = NPY_UBYTE
npy_type(::Type{Int16}) = NPY_SHORT
npy_type(::Type{UInt16}) = NPY_USHORT
npy_type(::Type{Int32}) = NPY_INT
npy_type(::Type{UInt32}) = NPY_UINT
npy_type(::Type{Int64}) = NPY_LONGLONG
npy_type(::Type{UInt64}) = NPY_ULONGLONG
npy_type(::Type{Float16}) = NPY_HALF
npy_type(::Type{Float32}) = NPY_FLOAT
npy_type(::Type{Float64}) = NPY_DOUBLE
npy_type(::Type{ComplexF32}) = NPY_CFLOAT
npy_type(::Type{ComplexF64}) = NPY_CDOUBLE
npy_type(::Type{PyPtr}) = NPY_OBJECT

# flags:
const NPY_ARRAY_ALIGNED = Int32(0x0100)
const NPY_ARRAY_WRITEABLE = Int32(0x0400)

#########################################################################
# no-copy conversion of Julia arrays to NumPy arrays.

# Julia arrays are in column-major order, but in some cases it is useful
# to pass them to Python as row-major arrays simply by reversing the
# dimensions. For example, although NumPy works with both row-major and
# column-major data, some Python libraries like OpenCV seem to require
# row-major data (the default in NumPy). In such cases, use PyReverseDims(array)
function NpyArray(a::StridedArray{T}, revdims::Bool) where T<:PYARR_TYPES
    @npyinitialize
    size_a = revdims ? reverse(size(a)) : size(a)
    strides_a = revdims ? reverse(strides(a)) : strides(a)
    p = @pycheck ccall(npy_api[:PyArray_New], PyPtr,
          (PyPtr,Cint,Ptr{Int},Cint, Ptr{Int},Ptr{T}, Cint,Cint,PyPtr),
          npy_api[:PyArray_Type],
          ndims(a), Int[size_a...], npy_type(T),
          Int[strides_a...] * sizeof(eltype(a)), a, sizeof(eltype(a)),
          NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE,
          C_NULL)
    return PyObject(p, a)
end

# function PyObject(a::StridedArray{T}) where T<:PYARR_TYPES
#     try
#         return NpyArray(a, false)
#     catch
#         return array2py(a) # fallback to non-NumPy version
#     end
# end

function PyReverseDims(a::StridedArray{T,N}) where {T<:PYARR_TYPES,N}
    try
        return NpyArray(a, true)
    catch
        return array2py(permutedims(a, N:-1:1)) # fallback to non-NumPy version
    end
end
PyReverseDims(a::BitArray) = PyReverseDims(Array(a))

# fallback to physically transposing the array
PyReverseDims(a::AbstractArray{<:Any,N}) where {N} = PyObject(permutedims(a, N:-1:1))
PyReverseDims(a::AbstractMatrix) = PyObject(permutedims(a))

"""
    PyReverseDims(array)

Passes a Julia `array` to Python as a NumPy row-major array
(rather than Julia's native column-major order) with the
dimensions reversed (e.g. a 2×3×4 Julia array is passed as
a 4×3×2 NumPy row-major array).  This is useful for Python
libraries that expect row-major data.
"""
PyReverseDims(a::AbstractArray)

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

# # transposed arrays can be passed to NumPy without copying
# PyObject(a::Union{LinearAlgebra.Adjoint{<:Real},LinearAlgebra.Transpose}) =
#     PyReverseDims(a.parent)

# PyObject(a::LinearAlgebra.Adjoint) = PyObject(Matrix(a)) # non-real arrays require a copy


text:
module TestAdd

using CxxCall
dir = mktempdir()
lib = joinpath(dir, "libHelloWorld")
filepath = joinpath(dir, "HelloWorld.cxx")
eval(cxxsetup())
eval(cxxnewfile(filepath,
"""
#include <iostream>
"""))

@cxx lib function hello(name::Cstring)::Cvoid
    """
    std::cout << "Hello " << name << "!" << std::endl;
    """
end

using Test
using Libdl

@testset "hello_world" begin
    @test !ispath(filepath)
    cxx_write_code!()
    @test isfile(filepath)
    libpath = "$(lib).$dlext"
    run(`g++ -shared -fPIC $(filepath) -o $libpath`)
    @test isfile(libpath)
    dlopen(libpath)
    hello("world")
end

end#module


text:
## cholesky.jl : Gaussian random field generator based on the Cholesky decomposition of the covariance matrix

"""
    Cholesky <: GaussianRandomFieldGenerator

A [`GaussiandRandomFieldGenerator`](@ref) based on a Cholesky factorization of the covariance matrix.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts1 = 0:0.02:1; pts2 = 0:0.02:1
0.0:0.02:1.0

julia> grf = GaussianRandomField(c,Cholesky(),pts1,pts2)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51x51 structured grid, using a Cholesky decomposition

julia> plot(grf)
[...]

```
Note that the Cholesky factorization requires the covariance matrix to be `Symmetric` and positive definite. If the covariance matrix is not `Symmetric`, a warning will be thrown but the method will try to continue with an approximate symmetric matrix. If the covariance matrix is not positive definite, an error will be thrown. Try using the `Spectral` method in that case.

See also: [`Spectral`](@ref), [`KarhunenLoeve`](@ref), [`CirculantEmbedding`](@ref)
"""
struct Cholesky <: GaussianRandomFieldGenerator end

const CholeskyGRF = GaussianRandomField{C,Cholesky} where {C}

# compose a GaussianRandomField using Cholesky factorization
function _GaussianRandomField(mean,cov,method::Cholesky,pts...)
    C = apply(cov,pts,pts)

    # error checking
    issymmetric(C) || warn("to use a Cholesky factorization, the covariance matrix must be symmetric/hermitian")
    isposdef(C) || throw(ArgumentError("to use a Cholesky factorization, the covariance matrix must be positive definite"))

    # compute Cholesky factorization
    L = (cholesky(Symmetric(C))).U'

    GaussianRandomField{typeof(cov),Cholesky,typeof(pts)}(mean,cov,pts,L)
end

# returns the required dimension of the random points
randdim(grf::CholeskyGRF) = size(grf.data,1)

# sample from the GaussianRandomField using Cholesky factorization
function _sample(grf::CholeskyGRF, xi)
	grf.mean + std(grf.cov)*reshape(grf.data*xi,size(grf.mean))
end

show(io::IO,::Cholesky) = print(io,"Cholesky decomposition")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule cilantro_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("cilantro")
JLLWrappers.@generate_main_file("cilantro", UUID("a3b35723-f448-5ea8-82d7-f10131ef8160"))
end  # module cilantro_jll


text:
#using TensorToolbox, Test, LinearAlgebra

println("\n\n****Testing ttensor.jl")

X=rand(20,10,50,5)

println("\n...Test tensor X of size: ", size(X))
println("\n...Testing hosvd.")
println("\nCreating exact decomposition with rank = size(X):")
T=hosvd(X)
println("Results:")
println("Type of output T: ", typeof(T))
@test isa(T,ttensor)
println("Core tensor size: ", coresize(T))
@test coresize(T) == size(X)
print("Factor matrices sizes: ")
for n=1:ndims(T)
	print(size(T.fmat[n])," ")
  @test size(T.fmat[n]) == (size(X,n),size(X,n))
end

err=norm(full(T) - X)
println("\n\n...Testing function full, i.e. n-mode multiplication (ttm): norm(full(T)-X) = ", err)
@test err ≈ 0 atol=1e-10

A=MatrixCell(undef,ndims(T))
for n=1:ndims(T)
  A[n]=rand(rand(1:10),size(T,n))
end
err=norm(full(ttm(T,A)) - ttm(full(T),A))
println("\n...Testing ttm for ttensor T and array of matrices A : norm(full(ttm(T,A))-ttm(full(T),A)) = ", err)
@test err ≈ 0 atol=1e-10

R=[5,5,5,5]
println("\n...Testing hosvd with smaller multilinear rank: ", R)
T=hosvd(X,reqrank=R)
println("Results:")
println("Type of output T: ", typeof(T))
@test isa(T,ttensor)
println("Core tensor size: ", coresize(T))
@test coresize(T) == tuple(R...)
print("Factor matrices sizes: ")
for n=1:ndims(T)
  print(size(T.fmat[n])," ")
  @test size(T.fmat[n]) == (size(X,n),R[n])
end
println("\n\n...Testing size of ttensor T : ", size(T))
@test size(T) == size(X)
println("\n...Testing ndims of ttensor T : ", ndims(T))
@test ndims(T) == ndims(X)
println("\n...Testing nrank of ttensor T for mode 1: ", nrank(T,1))
@test nrank(T,1) == R[1]
println("\n...Testing mrank of ttensor T: ", mrank(T))
@test mrank(T) == tuple(R...)

println("\n...Testing functions matten and tenmat (by mode).")
for n=1:ndims(T)
  Tn=tenmat(T,n)
  @test norm(Tn-tenmat(full(T),n)) ≈ 0 atol=1e-10
end

R=[5,5,5];
T=hosvd(rand(60,50,40),reqrank=R);
sz = size(T);
X=full(T);
println("\n\n...Testing hosvd for tensor with noise.")
println("For ttensor T, X=full(T) tensor of size ",size(X)," and rank ",R," , N noise tensor and S=hosvd(X+N,",R,").")
rsz = tuple([[sz...]' 1 1]...)
Rdm = reshape(randn(rsz), sz);
N=1e-3 * norm(X) * Rdm / norm(Rdm);
Y = X + N;
S=hosvd(Y,reqrank=R)
err=norm(T-S)
noise=norm(N)
println("Error( norm(T-S) ): ",err,". Noise norm: ",noise,".")

f(x,y,z)=1/(x+y+z)
dx=dy=dz=[n*0.1 for n=1:20]
X=[ f(x,y,z) for x=dx, y=dy, z=dz ]

println("\n\n...Testing hosvd with eps_abs=1e-5 on function defined tensor X of size ", size(X), " and multlinear rank", mrank(X))
T=hosvd(X,eps_abs=1e-5)
println("Results:")
println("Type of output T: ", typeof(T))
@test isa(T,ttensor)
println("Core tensor size: ",  coresize(T))
print("Factor matrices sizes: ")
for fmat in T.fmat
	print(size(fmat)," ")
end
@test norm(full(T)-X) ≈ 0 atol=1e-5

println("\n\n...Testing hosvd with eps_rel=1e-5 on function defined tensor X of size ", size(X))
T=hosvd(float(X),eps_rel=1e-5)
println("Results:")
println("Type of output T: ", typeof(T))
@test isa(T,ttensor)
println("Core tensor size: ", coresize(T))
print("Factor matrices sizes: ")
for fmat in T.fmat
	print(size(fmat)," ")
end
println("\n\n...Testing if factor matrices of ttensor T are orthogonal: ", T.isorth)
@test T.isorth

R=[3,3,3]
println("\n\n...Recompress Tucker tensor T to smaller rank: ", R)
S=hosvd(T,reqrank=R)
println("Results:")
println("Type of output T: ", typeof(S))
@test isa(S,ttensor)
println("Core tensor size: ", coresize(S))
@test coresize(S) == tuple(R...)
print("Factor matrices sizes: ")
for n=1:ndims(S)
  print(size(S.fmat[n])," ")
  @test size(S.fmat[n]) == (size(T,n),R[n])
end

S=randttensor([4,5,2],[3,3,3])
println("\n\n...Testing orthogonal flag and reorthogonalization for random ttensor S of size ", size(S),".")
println("Has orthogonal factor matrices: ", S.isorth)
println("After reorthogonalization: ", reorth!(S).isorth)
@test S.isorth

println("\n\n...Testing norm of ttensor T.")
err=abs(norm(T) - norm(full(T)))
println("|norm(T) - norm(full(T))| = ", err)
@test err ≈ 0 atol=1e-10

println("\n\n...Testing scalar multiplication 3*T.")
err=norm(full(3*T) - 3*full(T))
println("norm(full(3*T) - 3*full(T)) = ", err )
@test err ≈ 0 atol=1e-10

X=randttensor([6,8,2,5,4],[4,3,2,2,3]);
Y=randttensor([6,8,2,5,4],[3,6,3,4,3]);
println("\n\n...Creating two random ttensors X and Y of size ", size(X),".")
println("\n...Testing addition.")
Z=X+Y;
F=full(X)+full(Y);
err=norm(full(Z) - F)
println("norm(full(X+Y) - (full(X)+full(Y))) = ",err )
@test err ≈ 0 atol=1e-10

println("\n\n...Testing inner product.")
Z=innerprod(X,Y)
err=abs(Z - innerprod(full(X),full(Y)))
println("|innerprod(X,Y) - innerprod(full(X),full(Y))| = ",err )
@test err ≈ 0 atol=1e-10

println("\n\n...Testing Hadamard product.")
err=norm(full(ewprod(X,Y)) - full(X).*full(Y))
println("norm(full(ewprod(X,Y)) - full(X).*full(Y)) = ", err)
@test  err ≈ 0 atol=1e-10

println("\n\n...Testing singular values of matricizations of Tucker Tensor.")
R=[3,3,3]
T=randttensor([10,9,8],R)
println("\nSingular values of matricizations of random Tucker tensor of size ", size(T), ", rank ",R," and norm ",norm(T),".")
for n=1:ndims(T)
  sv = msvdvals(T,n)
  Tn=tenmat(T,n);
  println("Mode-$n singular values error: ",norm(sv-svd!(tenmat(full(T),n)).S[1:length(sv)]))
  @test norm(sv-svd!(tenmat(full(T),n)).S[1:length(sv)]) ≈ 0 atol=1e-10
end

R=[5,5,5,5]
T=randttensor([20,20,20,20],R);
println("\nSingular values of matricizations of random Tucker tensor of size ", size(T), ", rank ",R," and norm ",norm(T),".")
for n=1:ndims(T)
  sv = msvdvals(T,n)
  Tn=tenmat(T,n)
  println("Mode-$n singular values error: ",norm(sv-svd!(tenmat(full(T),n)).S[1:length(sv)]))
  @test norm(sv-svd!(tenmat(full(T),n)).S[1:length(sv)]) ≈ 0 atol=1e-10
end

println("\n...Testing function mhadtv.")
X=randttensor([5,4,3],[3,3,3])
Y=randttensor([5,4,3],[2,2,2])
println("\n\n...Creating two random ttensors X and Y of size ", size(X),".")
v=rand(12)
mode=1
println("Multiplying mode-$mode matricized ewprod(X,Y) by a random vector.")
Z=mhadtv(X,Y,v,mode,'n')
Hn=tenmat(ewprod(X,Y),mode)
err = norm(Z-Hn*v)
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10
v=rand(5)
Z=mhadtv(X,Y,v,mode,'t')
err = norm(Z-Hn'*v)
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10
v=rand(5)
Z=mhadtv(X,Y,v,mode,'b')
err = norm(Z-Hn*Hn'*v)
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing function mttkrp.")
X=randttensor([5,4,3],[3,3,3])
mode=1
M1=rand(5,5);
M2=rand(4,5);
M3=rand(3,5);
M=[M1,M2,M3]
println("Multiplying mode-$mode matricized tensor X by Khatri-Rao product of matrices.")
Z=mttkrp(X,M,mode)
err = norm(Z-tenmat(X,mode)*khatrirao(M3,M2))
println("Multiplication error: ",err)
@test err ≈ 0 atol=1e-10

println("\n...Testing cp_als on ttensor.")
I=5;J=4;K=3;
a=rand(I);
b=rand(J);
c=rand(K);
X=zeros(I,J,K);
for i=1:I,j=1:J,k=1:K
    X[i,j,k]=a[i]*b[j]*c[k];
end
T=hosvd(X)
println("Test rank-1 ttensor T of size: ", size(T))
Z=cp_als(T,1)
err=norm(full(Z)-full(T))
println("norm(full(Tcp)-T) = ",err)
@test err ≈ 0 atol=1e-10


text:
"""
DioOption represents a set of options that is used in DIO, a fundamental
frequency analysis.

**Fields**

$(FIELDS)
"""
mutable struct DioOption
    f0floor::Cdouble
    f0ceil::Cdouble
    channels_in_octave::Cdouble
    "frame period in ms"
    period::Cdouble
    speed::Cint
    "added in v0.2.1-2 (WORLD 0.2.0_2)"
    allowed_range::Cdouble

    function DioOption(f0floor, f0ceil, channels_in_octave, period, speed,
                       allowed_range = 0.1)
        f0floor > f0ceil && throw(ArgumentError("F0floor must be larger than F0ceil"))
        f0floor < 0 && throw(ArgumentError("f0floor must be positive"))
        if channels_in_octave < 0
            throw(ArgumentError("channels_in_octave must be positive"))
        end
        period <= 0 && throw(ArgumentError("period must be positive"))
        (1 <= speed <= 12) || throw(ArgumentError("1 ≤ speed ≤ 12 is supprted"))
        allowed_range >= 0 || throw(ArgumentError("allowed_range >= 0 is supported"))
        new(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
    end
end

# Note that the default options assume that the sampling frequency of a input
# speech signal is 44.1 kHz.
function DioOption(;
                   f0floor::Float64=71.0,
                   f0ceil::Float64=800.0,
                   channels_in_octave::Float64=2.0,
                   period::Float64=5.0,
                   speed::Integer=1,
                   allowed_range::Float64=0.1
                   )
    DioOption(f0floor, f0ceil, channels_in_octave, period, speed, allowed_range)
end

"""
HarvestOption represents a set of options that is used in Harvest, a fundamental
frequency analysis.

**Fields**

$(FIELDS)
"""
mutable struct HarvestOption
    f0floor::Cdouble
    f0ceil::Cdouble
    "frame period in ms"
    period::Cdouble
    function HarvestOption(f0floor, f0ceil, period)
        f0floor > f0ceil && throw(ArgumentError("F0floor must be larger than F0ceil"))
        f0floor < 0 && throw(ArgumentError("f0floor must be positive"))
        period <= 0 && throw(ArgumentError("period must be positive"))
        new(f0floor, f0ceil, period)
    end
end

function HarvestOption(;
                   f0floor::Float64=71.0,
                   f0ceil::Float64=800.0,
                   period::Float64=5.0)
   HarvestOption(f0floor, f0ceil, period)
end

"""
CheapTrick options

**Fields**

$(FIELDS)
"""
mutable struct CheapTrickOption
    q1::Cdouble
    f0floor::Cdouble
    fftsize::Cint
    function CheapTrickOption(fs; q1=-0.15, f0floor=71.0)
        p = new()
        ccall((:InitializeCheapTrickOption, libworld), Cvoid,
              (Cint, Ptr{CheapTrickOption},), fs, Ref(p))
        p.q1 = q1
        p.f0floor = f0floor
        p.fftsize = get_fftsize_for_cheaptrick(Int(fs), p)
        p
    end
end

"""
D4C options (nothing for now, but for future changes)

**Fields**

$(FIELDS)
"""
mutable struct D4COption
    threshold::Cdouble
    D4COption(threshold=0.85) = new(threshold)
end

function get_samples_for_dio(fs::Real, len::Integer, period::Real)
    ccall((:GetSamplesForDIO, libworld), Cint,
          (Cint, Cint, Cdouble), fs, len, period)
end

function get_samples_for_harvest(fs::Real, len::Integer, period::Real)
    ccall((:GetSamplesForHarvest, libworld), Cint,
          (Cint, Cint, Cdouble), fs, len, period)
end

"""
$(SIGNATURES)

Dio estimates F0 trajectory given a monoral input signal.

**Paremters**

- `x`  : Input signal
- `fs` : Sampling frequency
- `opt` : DioOption

**Returns**

- `time_axis`  : Temporal positions.
- `f0`         : F0 contour.

"""
function dio(x::StridedVector{Cdouble}, fs::Real, opt::DioOption=DioOption())
    expectedlen = get_samples_for_dio(fs, length(x), opt.period)
    f0 = Array{Cdouble}(undef, expectedlen)
    timeaxis = Array{Cdouble}(undef, expectedlen)
    ccall((:Dio, libworld),  Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ref{DioOption}, Ptr{Cdouble}, Ptr{Cdouble}),
          x, length(x), fs, opt, timeaxis, f0)
    f0, timeaxis
end

"""
$(SIGNATURES)

Harvest estimates F0 trajectory given a monoral input signal.

**Paremters**

- `x`  : Input signal
- `fs` : Sampling frequency
- `opt` : HarvestOption

**Returns**

- `time_axis`  : Temporal positions.
- `f0`         : F0 contour.
"""
function harvest(x::StridedVector{Cdouble}, fs::Real, opt::HarvestOption=HarvestOption())
    expectedlen = get_samples_for_harvest(fs, length(x), opt.period)
    f0 = Array{Cdouble}(undef, expectedlen)
    timeaxis = Array{Cdouble}(undef, expectedlen)
    ccall((:Harvest, libworld),  Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ref{HarvestOption}, Ptr{Cdouble}, Ptr{Cdouble}),
          x, length(x), fs, opt, timeaxis, f0)
    f0, timeaxis
end

"""
$(SIGNATURES)

StoneMask refines the estimated F0 by Dio,

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Temporal information
- `f0` : f0 contour

**Returns**

- `refined_f0` : Refined F0
"""
function stonemask(x::StridedVector{Cdouble}, fs::Integer,
                   timeaxis::StridedVector{Cdouble},
                   f0::StridedVector{Cdouble})
    refinedF0 = Array{Cdouble}(undef, length(f0))
    ccall((:StoneMask, libworld),  Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{Cdouble}),
          x, length(x), fs, timeaxis, f0, length(f0), refinedF0)
    refinedF0
end

# Array{T,2} -> Array{Ptr{T}}
function ptrarray2d!(dst::Array{Ptr{T},1}, src::Array{T,2}) where {T<:Real}
    for i=1:size(src, 2)
        @inbounds dst[i] = pointer(view(src, 1:size(src, 1), i), 1)
    end
end

"""
$(SIGNATURES)

GetFFTSizeForCheapTrick calculates the FFT size based on the sampling
frequency and the lower limit of f0 (It is defined in world.h).

**Parameters**

- `fs`: Sampling frequency
- `opt`: CheapTrickOption

**Returns**

- `fftsize` : FFT size
"""
function get_fftsize_for_cheaptrick(fs::Integer,
                                    opt::CheapTrickOption=CheapTrickOption(fs))
    fftsize = ccall((:GetFFTSizeForCheapTrick, libworld), Cint,
                    (Cint,Ptr{CheapTrickOption}), fs, Ref(opt))
    convert(Int, fftsize)
end

"""
$(SIGNATURES)

CheapTrick calculates the spectrogram that consists of spectral envelopes
estimated by CheapTrick.

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Time axis
- `f0` : F0 contour
- `opt` : CheapTrick option

**Returns**

- `spectrogram`  : Spectrogram estimated by CheapTrick.
"""
function cheaptrick(x::StridedVector{Cdouble}, fs::Integer,
                    timeaxis::StridedVector{Cdouble},
                    f0::StridedVector{Cdouble};
                    opt::CheapTrickOption=CheapTrickOption(fs)
    )
    freqbins = get_fftsize_for_cheaptrick(fs, opt)>>1 + 1
    spectrogram = Array{Cdouble}(undef, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array{Ptr{Cdouble}}(undef, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    ccall((:CheapTrick, libworld), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Ptr{CheapTrickOption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), Ref(opt), cspectrogram)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        @inbounds spectrogram[j,i] = unsafe_load(cspectrogram[i], j)
    end

    spectrogram
end

"""
$(SIGNATURES)

D4C calculates the aperiodicity estimated by D4C.

**Parameters**

- `x` : Input signal
- `fs` : Sampling frequency
- `time_axis` : Time axis
- `f0` : F0 contour

**Returns**

- `aperiodicity` : Aperiodicity estimated by D4C.
"""
function d4c(x::StridedVector{Cdouble}, fs::Integer,
             timeaxis::StridedVector{Cdouble},
             f0::StridedVector{Cdouble};
             opt::D4COption=D4COption())
    fftsize = get_fftsize_for_cheaptrick(fs, CheapTrickOption(fs))
    freqbins = fftsize>>1 + 1
    aperiodicity = zeros(Cdouble, freqbins, length(f0))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    caperiodicity = Array{Ptr{Cdouble}}(undef, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    ccall((:D4C, libworld), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint,
           Ptr{D4COption}, Ptr{Ptr{Cdouble}}),
          x, length(x), fs, timeaxis, f0, length(f0), fftsize, Ref(opt),
          caperiodicity)

    # Array{Cdouble,2} <- Array{Ptr{Cdouble}}
    for i=1:length(f0), j=1:freqbins
        aperiodicity[j,i] = unsafe_load(caperiodicity[i], j)
    end

    aperiodicity
end

"""
$(SIGNATURES)

Synthesis synthesize the voice based on f0, spectrogram and
aperiodicity (not excitation signal.

**Parameters**

- `f0` : f0 contour
- `spectrogram` : Spectrogram estimated by CheapTrick
- `aperiodicity` : Aperiodicity spectrogram based on D4C
- `period` :  Temporal period used for the analysis
- `fs` : Sampling frequency
- `len` : Length of the output signal

**Returns**

- `y` : Calculated speech
"""
function synthesis(f0::StridedVector{Cdouble},
                   spectrogram::StridedMatrix{Cdouble},
                   aperiodicity::StridedMatrix{Cdouble},
                   period::Real, fs::Integer, len::Integer)
    fftsize = get_fftsize_for_cheaptrick(fs, CheapTrickOption(fs))

    # Array{Cdouble,2} -> Array{Ptr{Cdouble}}
    cspectrogram = Array{Ptr{Cdouble}}(undef, size(spectrogram, 2))
    ptrarray2d!(cspectrogram, spectrogram)

    caperiodicity = Array{Ptr{Cdouble}}(undef, size(aperiodicity, 2))
    ptrarray2d!(caperiodicity, aperiodicity)

    synthesized = Array{Cdouble}(undef, len)
    ccall((:Synthesis, libworld), Cvoid,
          (Ptr{Cdouble}, Cint, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}},
           Cint, Cdouble, Cint, Cint, Ptr{Cdouble}),
          f0, length(f0), cspectrogram, caperiodicity, fftsize, period, fs, len,
          synthesized)

    synthesized
end

# matlabfunctions

"""
$(SIGNATURES)

inplace version of interp1

**Parameters**

- `x` : Input vector (Time axis)
- `y` : Values at x[n]
- `xi`: Required vector
- `yi` : Interpolated vector
"""
function interp1!(x::StridedVector{Cdouble},
                  y::StridedVector{Cdouble},
                  xi::StridedVector{Cdouble},
                  yi::StridedVector{Cdouble})
    @assert length(x) == length(y)
    @assert length(xi) == length(yi)
    ccall((:interp1, libworld), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, y, length(x), xi, length(xi), yi)
    yi
end

"""
$(SIGNATURES)

interp1 interpolates to find yi, the values of the underlying function Y
at the points in the vector or array xi. x must be a vector.
http://www.mathworks.co.jp/help/techdoc/ref/interp1.html

**Parameters**

- `x` : Input vector (Time axis)
- `y` : Values at x[n]
- `xi`: Required vector

**Returns**

- `yi` : Interpolated vector
"""
function interp1(x::StridedVector{Cdouble},
                 y::StridedVector{Cdouble},
                 xi::StridedVector{Cdouble})
    yi = similar(xi)
    interp1!(x, y, xi, yi)
end


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

"""
    find_marginalmodes(marg::MarginalDist)

*BAT-internal, not part of stable public API.*

Find the modes of a MarginalDist.
Returns a vector of the bin-centers of the bin(s) with the heighest weight.
"""
function find_marginalmodes(marg::MarginalDist)
    hist = convert(Histogram, marg.dist.dist)
    dims = ndims(hist.weights)

    max = maximum(hist.weights)
    maxima_idx = findall(x->x==max, hist.weights)

    bin_centers = get_bin_centers(marg)

    return [[bin_centers[d][maxima_idx[i][d]] for d in 1:dims] for i in 1:length(maxima_idx) ]
end

"""
    get_bin_centers(marg::MarginalDist)

*BAT-internal, not part of stable public API.*

Returns a vector of the bin-centers.
"""
function get_bin_centers(marg::MarginalDist)
    hist = convert(Histogram, marg.dist.dist)
    edges = hist.edges
    dims = ndims(hist.weights)

    centers = [[edges[d][i]+0.5*(edges[d][i+1]-edges[d][i]) for i in 1:length(edges[d])-1] for d in 1:dims]

    return centers
end

function islower(weights, idx)
    if idx==1 && weights[idx]>0
        return true
    elseif weights[idx]>0 && weights[idx-1]==0 && idx < length(weights)
        return true
    else
        return false
    end
end

function isupper(weights, idx)
    if idx==length(weights) && weights[idx-1]>0
        return true
    elseif weights[idx]==0 && weights[idx-1]>0
        return true
    else
        return false
    end
end

# return the lower and upper edges for clusters in which the bincontent is non-zero for all dimensions of a StatsBase.Histogram
# clusters that are seperated <= atol are combined
function get_interval_edges(h::StatsBase.Histogram; atol::Real = 0)
    weights = h.weights
    len = length(weights)

    lower = [h.edges[1][i] for i in 1:len if islower(weights, i)]
    upper = [h.edges[1][i] for i in 2:len if isupper(weights, i)]

    if atol != 0
        idxs = [i for i in 1:length(upper)-1 if lower[i+1]-upper[i] <= atol]
        deleteat!(upper, idxs)
        deleteat!(lower, idxs.+1)
    end

    return lower, upper
end

# for 1d and 2d histogramsm
function get_smallest_intervals(
    histogram::StatsBase.Histogram,
    intervals::Array{Float64, 1}
)
    intervals = sort(intervals)

    hist = deepcopy(histogram)
    dims = size(hist.weights)
    weights = vec(hist.weights)
    totalweight = sum(weights)
    rel_weights = weights/totalweight

    hists_weights = [zeros(length(weights)) for i in 1:length(intervals)]
    weight_ids = sortperm(weights, rev=true) # starting with highest weight

    for (i, intv) in enumerate(intervals)
        sum_weights = 0.
        for w in weight_ids
            if sum_weights < intv
                hists_weights[i][w] = weights[w]
                sum_weights += rel_weights[w]
            else
                break
            end
        end
    end

    hists = Array{StatsBase.Histogram}(undef, length(intervals))

    for i in 1:length(intervals)
        hists[i] = deepcopy(hist)
        hists[i].weights = reshape(hists_weights[i], dims)
    end

    realintervals = get_probability_content(hist, hists)

    return reverse(hists), reverse(realintervals)
end

# for 1d histograms
function split_central(
    histogram::StatsBase.Histogram,
    intervals::Array{Float64, 1}
)
    intervals = sort(intervals)
    intervals = (1 .-intervals)/2

    hist = deepcopy(histogram)
    hists = Array{StatsBase.Histogram}(undef, length(intervals))

    for i in 1:length(intervals)
        hists[i] = deepcopy(hist)
    end

    weights = vec(hist.weights)
    totalweight = sum(weights)
    rel_weights = weights/totalweight

    for (i, intv) in enumerate(intervals)
        sum_left = 0.
        sum_right = 0.

        for l in 1:length(weights)
           if sum_left + rel_weights[l] < intv
                sum_left = sum_left + rel_weights[l]
                hists[i].weights[l] = 0
           else
                break
           end
        end

        for r in length(weights):-1:1
           if sum_right + rel_weights[r] < intv
                sum_right = sum_right + rel_weights[r]
                hists[i].weights[r] = 0
           else
                break
           end
        end
    end

    realintervals = get_probability_content(hist, hists)

    return reverse(hists), reverse(realintervals)
end

# calculate probability percentage enclosed inside the intervals of hists
function get_probability_content(
    hist::StatsBase.Histogram,
    hists::Array{StatsBase.Histogram, 1}
)
    totalweight = sum(hist.weights)
    return [sum(hists[i].weights)/totalweight for i in 1:length(hists)]
end

function calculate_levels(
    hist::StatsBase.Histogram,
    intervals::Array{<:Real, 1}
)
    intervals = sort(intervals)
    levels = Vector{Real}(undef, length(intervals)+1)

    weights = sort(vec(hist.weights), rev=true)

    weight_ids = sortperm(weights, rev=true);
    sum_of_weights = sum(weights)

    sum_w = 0.0

    for w_id in weight_ids
        if(sum_w <= intervals[end])
            i = findfirst(x -> sum_w <= x, intervals)
            levels[i] = weights[w_id]
            sum_w += weights[w_id]/sum_of_weights
        end
    end

    levels[end] = 1.1*sum_of_weights
    return sort(levels)
end


text:
using RandomMatrices
using Base.Test

srand(1)
include("FormalPowerSeries.jl")
include("Haar.jl")


text:
module InterfaceTests

using Test
using Gridap.Arrays

a = rand(20,12)

test_array(a,a)
test_array(a,a,≈)

# @test array_caches() == ()
# @test getitems!((),(),1) == ()

end # module


text:
struct DynamicBudgetedBipartiteMatchingSolution{T, U} <: BudgetedBipartiteMatchingSolution{T, U}
  instance::BudgetedBipartiteMatchingInstance{T, U}
  solution::Vector{Edge{T}}
  value::Float64
  state::Dict{Tuple{Int, Int, Int, Int}, Float64} # Data structure built by the dynamic-programming recursion.
  solutions::Dict{Tuple{Int, Int, Int, Int}, Vector{Edge{T}}} # From the indices of state to the corresponding solution.
end

function matching_dp_budgeted(instance::BudgetedBipartiteMatchingInstance{T, U}) where {T, U}
  # First, determine the size of a perfect matching.
  nr = Dict(k => 1.0 for (k, v) in instance.matching.reward)
  ni = BipartiteMatchingInstance(instance.matching, nr)
  ns = matching_dp_imperfect(ni)
  m = length(ns.solution)

  # Solve this budgeted matching to ensure the constraints are met.

  # Four indices:
  # - two vertices, like matching_dp_imperfect;
  # - the number of edges still to add to have a perfect matching, like matching_dp;
  # - the minimum budget to respect.
  state = Dict{Tuple{Int, Int, Int, Int}, Float64}((i, j, k, l) => -Inf for i in instance.matching.vertex_left, j in instance.matching.vertex_right, k in 0:m, l in 0:instance.budget)
  solutions = Dict{Tuple{Int, Int, Int, Int}, Vector{Edge{T}}}((i, j, m, zero(U)) => Edge{T}[] for i in instance.matching.vertex_left, j in instance.matching.vertex_right)

  # Initialise: k=m-1, l=0, i.e. take at most one edge, don't consider the budget yet.
  v1 = first(instance.matching.vertex_right)
  for i in instance.matching.vertex_left
    if has_edge(instance.matching.graph, i, v1)
      state[i, v1, m - 1, 0] = instance.matching.reward[Edge(i, v1)]
      solutions[i, v1, m - 1, 0] = [Edge(i, v1)]
    end
  end

  v1 = first(instance.matching.vertex_left)
  for j in instance.matching.vertex_right
    if has_edge(instance.matching.graph, v1, j)
      state[v1, j, m - 1, 0] = instance.matching.reward[Edge(v1, j)]
      solutions[v1, j, m - 1, 0] = [Edge(v1, j)]
    end
  end

  for (i, prev_i) in zip(instance.matching.vertex_left[2:end], instance.matching.vertex_left[1:end-1])
    for (j, prev_j) in zip(instance.matching.vertex_right[2:end], instance.matching.vertex_right[1:end-1])
      # Take (i, j).
      if has_edge(instance.matching.graph, i, j)
        state[i, j, m - 1, 0] = instance.matching.reward[Edge(i, j)]
        solutions[i, j, m - 1, 0] = [Edge(i, j)]
      end

      # Don't take (i, j): copy a previous solution. No mixing possible, as only one edge is considered.
      if state[prev_i, j, m - 1, 0] > state[i, j, m - 1, 0]
        state[i, j, m - 1, 0] = state[prev_i, j, m - 1, 0]
        solutions[i, j, m - 1, 0] = solutions[prev_i, j, m - 1, 0]
      end
      if state[i, prev_j, m - 1, 0] > state[i, j, m - 1, 0]
        state[i, j, m - 1, 0] = state[i, prev_j, m - 1, 0]
        solutions[i, j, m - 1, 0] = solutions[i, prev_j, m - 1, 0]
      end
    end
  end

  # Initialisation for l=0, k < m-1: create perfect matchings, regardless of budget.
  for k in (m - 2):-1:0
    # Don't handle the cases with just one vertex in one part of the graph:
    # this loop requires at least two edges to be taken, but the considered
    # subset of edges can lead to at most only one being taken.

    # Consider the other vertices.
    for (i, prev_i) in zip(instance.matching.vertex_left[2:end], instance.matching.vertex_left[1:end-1])
      for (j, prev_j) in zip(instance.matching.vertex_right[2:end], instance.matching.vertex_right[1:end-1])
        # Take (i, j).
        if has_edge(instance.matching.graph, i, j) && state[prev_i, prev_j, k + 1, 0] > -Inf
          state[i, j, k, 0] = instance.matching.reward[Edge(i, j)] + state[prev_i, prev_j, k + 1, 0]
          solutions[i, j, k, 0] = solutions[prev_i, prev_j, k + 1, 0]
          push!(solutions[i, j, k, 0], Edge(i, j))
        end

        # Don't take (i, j): mix previous solutions (as a special case, maybe take a previous solution completely).
        for k2 in 0:m
          for k3 in 0:m
            if state[prev_i, j, k2, 0] + state[i, prev_j, k3, 0] > state[i, j, k, 0]
              potential_sol = union(Set{Edge{T}}(solutions[prev_i, j, k2, 0]), Set{Edge{T}}(solutions[i, prev_j, k3, 0]))
              val = sum(instance.matching.reward[Edge(src(e), dst(e))] for e in potential_sol)

              if length(potential_sol) == m - k && val > state[i, j, k, 0]
                state[i, j, k, 0] = val
                solutions[i, j, k, 0] = collect(potential_sol)
              end
            end
          end
        end
      end
    end
  end

  # Main iteration, considering the first vertex of each side is already done.
  for l in 1:instance.budget
    # Repeat the loop as often as necessary to handle the corner case with zero-weight edges.
    while true
      changes = false

      # Deal with the first two vertices (one on each side). All constraints are imposed.
      # Just one edge to be taken, as one vertex is always fixed on one side:
      # no feasible solution for k < m-1.
      vl = first(instance.matching.vertex_left)
      vr = first(instance.matching.vertex_right)
      if has_edge(instance.matching.graph, vl, vr) && instance.weight[Edge(vl, vr)] >= l
        state[vl, vr, m - 1, l] = instance.matching.reward[Edge(vl, vr)]
        solutions[vl, vr, m - 1, l] = [Edge(vl, vr)]
      end

      v1 = first(instance.matching.vertex_right)
      for i in instance.matching.vertex_left
        if has_edge(instance.matching.graph, i, v1) && instance.weight[Edge(i, v1)] >= l
          state[i, v1, m - 1, l] = instance.matching.reward[Edge(i, v1)]
          solutions[i, v1, m - 1, l] = [Edge(i, v1)]
        end
      end

      v1 = first(instance.matching.vertex_left)
      for j in instance.matching.vertex_right
        if has_edge(instance.matching.graph, v1, j) && instance.weight[Edge(v1, j)] >= l
          state[v1, j, m - 1, l] = instance.matching.reward[Edge(v1, j)]
          solutions[v1, j, m - 1, l] = [Edge(v1, j)]
        end
      end

      # Do the rest.
      for k in (m - 1):-1:0
        for (i, prev_i) in zip(instance.matching.vertex_left[2:end], instance.matching.vertex_left[1:end-1])
          for (j, prev_j) in zip(instance.matching.vertex_right[2:end], instance.matching.vertex_right[1:end-1])

            # Take (i, j) if possible.
            if has_edge(instance.matching.graph, i, j)
              remaining_budget = max(0, l - instance.weight[Edge(i, j)])

              if state[prev_i, prev_j, k + 1, remaining_budget] > -Inf
                state[i, j, k, l] = instance.matching.reward[Edge(i, j)] + state[prev_i, prev_j, k + 1, remaining_budget]
                solutions[i, j, k, l] = solutions[prev_i, prev_j, k + 1, remaining_budget]
                push!(solutions[i, j, k, l], Edge(i, j))

                changes = true
              end
            end

            # Don't take (i, j): mix previous solutions (as a special case, maybe take a previous solution completely).
            for k2 in 0:m
              for k3 in 0:m
                if state[prev_i, j, k2, l] + state[i, prev_j, k3, l] > state[i, j, k, l]
                  potential_sol = union(Set{Edge{T}}(solutions[prev_i, j, k2, l]), Set{Edge{T}}(solutions[i, prev_j, k3, l]))
                  val = sum(instance.matching.reward[Edge(src(e), dst(e))] for e in potential_sol)
                  wei = sum(instance.weight[Edge(src(e), dst(e))] for e in potential_sol)

                  if length(potential_sol) == m - k && val > state[i, j, k, l] && wei >= l
                    state[i, j, k, l] = val
                    solutions[i, j, k, l] = collect(potential_sol)

                    changes = true
                  end
                end
              end
            end
          end
        end
      end

      # If this iteration made no change, stop and go on to the next budget.
      if ! changes
        break
      end
    end
  end

  val = state[last(instance.matching.vertex_left), last(instance.matching.vertex_right), 0, instance.budget]
  sol = solutions[last(instance.matching.vertex_left), last(instance.matching.vertex_right), 0, instance.budget]

  return DynamicBudgetedBipartiteMatchingSolution(instance, sol, val, state, solutions)
  # return SimpleBudgetedBipartiteMatchingSolution(instance, sol)
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Simulation
s = get_simulation("quadruped", "flat_2D_lc", "flat");
model = s.model
env = s.env

# ## Reference Trajectory
ref_traj = deepcopy(get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
    load_type = :split_traj_alt))
update_friction_coefficient!(ref_traj, model, env);
H = ref_traj.H
h = ref_traj.h

# ## MPC setup
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 1000
κ_mpc = 2.0e-4

obj = TrackingObjective(model, env, H_mpc,
    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.nq-3)]) for t = 1:H_mpc],
    u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

p = ci_mpc_policy(ref_traj, s, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
	mode = :configuration,
    n_opts = NewtonOptions(
		solver = :lu_solver,
		r_tol = 3e-4,
		max_iter = 5,
		),
    mpc_opts = CIMPCOptions(),
	ip_opts = InteriorPointOptions(
					undercut = 5.0,
					γ_reg = 0.1,
					κ_tol = κ_mpc,
					r_tol = 1.0e-8,
					diff_sol = true,
					solver = :empty_solver,
					# max_time = 1000.0,
					),
    )

# ## Initial conditions
q1_sim, v1_sim = initial_conditions(ref_traj);

# ## Simulator
sim = simulator(s, H_sim, h=h_sim, policy=p);

# ## Simulate
simulate!(sim, q1_sim, v1_sim);

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
anim = visualize_meshrobot!(vis, model, sim.traj, h=h_sim * 5, sample=5);

# ## Timing result
# Julia is [JIT-ed](https://en.wikipedia.org/wiki/Just-in-time_compilation) so re-run the MPC setup through Simulate for correct timing results.
process!(sim.stats, N_sample) # Time budget
H_sim * h_sim / sum(sim.stats.policy_time) # Speed ratio


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("datasets", "iris"), x=:SepalLength, y=:SepalWidth,
     color=:Species, Geom.point)


text:
@testset "Fitting" begin
  @testset "Block" begin
    gev_seed = StableRNG(12345)
    dist1    = GeneralizedExtremeValue(0.0, 1.0, 0.0)
    dist2    = GeneralizedExtremeValue(0.0, 1.0, 0.5)
    dist3    = GeneralizedExtremeValue(1.0, 1.0, 0.5) # The problem distribution...
    r_gev_1  = readdlm(joinpath(datadir,"R_gev_dist1_fits.csv"), ',')
    r_gev_2  = readdlm(joinpath(datadir,"R_gev_dist2_fits.csv"), ',')
    r_gev_3  = readdlm(joinpath(datadir,"R_gev_dist3_fits.csv"), ',')
    r_result = zip(eachrow(r_gev_1), eachrow(r_gev_2), eachrow(r_gev_3))
    for (test_ix, r_gev_ix) in enumerate(r_result)
      for (r_mle, dist) in zip(r_gev_ix, (dist1, dist2, dist3))
        sample = rand(gev_seed, dist, 1000)
        try
          julia_fit = fit(GeneralizedExtremeValue, BlockMaxima(sample, 1))
          j_mle = [julia_fit.μ, julia_fit.σ, julia_fit.ξ]
          @test isapprox(j_mle, r_mle, atol=1e-3)
        catch er
          @info "Optimization failed for $dist on trial $test_ix with error $er."
        end
      end
    end
  end

  @testset "PoT" begin
    gpd_seed = StableRNG(12345)
    dist1    = GeneralizedPareto(0.0, 1.0, 0.0)
    dist2    = GeneralizedPareto(0.0, 1.0, 0.5)
    dist3    = GeneralizedPareto(1.0, 1.0, 0.5)
    r_gpd_1  = readdlm(joinpath(datadir,"R_gpd_dist1_fits.csv"), ',')
    r_gpd_2  = readdlm(joinpath(datadir,"R_gpd_dist2_fits.csv"), ',')
    r_gpd_3  = readdlm(joinpath(datadir,"R_gpd_dist3_fits.csv"), ',')
    r_result = zip(eachrow(r_gpd_1), eachrow(r_gpd_2), eachrow(r_gpd_3))
    for (test_ix, r_gpd_ix) in enumerate(r_result)
      for (r_mle, dist) in zip(r_gpd_ix, (dist1, dist2, dist3))
        sample = rand(gpd_seed, dist, 1000)
        try
          julia_fit = fit(GeneralizedPareto, PeakOverThreshold(sample, dist.μ))
          j_mle     = [julia_fit.σ, julia_fit.ξ]
          @test isapprox(j_mle, r_mle, atol=1e-3)
        catch er
          @info "Optimization failed for $dist on trial $test_ix with error $er."
        end
      end
    end
  end
end

#= Code to generate the R test files:
using RCall, StableRNGs, Distributions
R"library(ismev)"
gev_R_results = NTuple{3, Vector{Float64}}[]
gev_seed = StableRNG(12345)
dist1    = GeneralizedExtremeValue(0.0, 1.0, 0.0)
dist2    = GeneralizedExtremeValue(0.0, 1.0, 0.5)
dist3    = GeneralizedExtremeValue(1.0, 1.0, 0.5)
for test_ix in 1:10
  mles_trialj = map((dist1, dist2, dist3)) do dist
    sample = rand(gev_seed, dist, 1000)
    R"""
      r_fit = gev.fit($sample, show=FALSE)
      r_mle = r_fit$mle
    """
    @rget r_mle
  end
  push!(gev_R_results, mles_trialj)
end

gpd_R_results = NTuple{3, Vector{Float64}}[]
gpd_seed = StableRNG(12345)
dist1    = GeneralizedPareto(0.0, 1.0, 0.0)
dist2    = GeneralizedPareto(0.0, 1.0, 0.5)
dist3    = GeneralizedPareto(1.0, 1.0, 0.5)
for test_ix in 1:10
  mles_trialj = map((dist1, dist2, dist3)) do dist
    sample = rand(gpd_seed, dist, 1000)
    R"""
      r_fit = gpd.fit($sample, threshold=$(dist.μ), show=FALSE)
      r_mle = r_fit$mle
    """
    @rget r_mle
  end
  push!(gpd_R_results, mles_trialj)
end

using DelimitedFiles
writedlm("data/R_gev_dist1_fits.csv", getindex.(gev_R_results, 1), ',')
writedlm("data/R_gev_dist2_fits.csv", getindex.(gev_R_results, 2), ',')
writedlm("data/R_gev_dist3_fits.csv", getindex.(gev_R_results, 3), ',')
writedlm("data/R_gpd_dist1_fits.csv", getindex.(gpd_R_results, 1), ',')
writedlm("data/R_gpd_dist2_fits.csv", getindex.(gpd_R_results, 2), ',')
writedlm("data/R_gpd_dist3_fits.csv", getindex.(gpd_R_results, 3), ',')
=#


text:
module AD4SM

include("adiff.jl")
include("materials.jl")
include("elements.jl")

# using .adiff, .Materials, .Elements
export adiff, Materials, Elements

end # module


text:
module MyPlots

import Plots as plt
plt.pyplot()

function f(x, y)
    return (1.5-x+x*y)^2+(2.25-x+x*y^2)^2+(2.625-x+x*y^3)^2
end

const plot = function ()
    x = -5.0:0.1:5.0
    y = -5.0:0.1:5.0
    # plt.heatmap(x, y, f)
    # plt.png("MyPlots/heat.png")
    plt.plot(x, y, f, st=:wireframe)
    plt.png("MyPlots/wireframe.png")
end

end # module

max_stars_repo_name: applied-math-coding/article-snippets

max_issues_repo_name: applied-math-coding/article-snippets

max_forks_repo_name: applied-math-coding/article-snippets


text:
const orb_sampling_pattern = [[8, -3, 9, 5], [4, 2, 7, -12], [-11, 9, -8, 2], [7, -12, 12, -13], [2, -13, 2, 12], [1, -7, 1, 6],
                        [-2, -10, -2, -4], [-13, -13, -11, -8], [-13, -3, -12, -9], [10, 4, 11, 9], [-13, -8, -8, -9], [-11, 7, -9, 12],
                        [7, 7, 12, 6], [-4, -5, -3, 0], [-13, 2, -12, -3], [-9, 0, -7, 5], [12, -6, 12, -1], [-3, 6, -2, 12], [-6, -13, -4, -8],
                        [11, -13, 12, -8], [4, 7, 5, 1], [5, -3, 10, -3], [3, -7, 6, 12], [-8, -7, -6, -2], [-2, 11, -1, -10], [-13, 12, -8, 10],
                        [-7, 3, -5, -3], [-4, 2, -3, 7], [-10, -12, -6, 11], [5, -12, 6, -7], [5, -6, 7, -1], [1, 0, 4, -5], [9, 11, 11, -13],
                        [4, 7, 4, 12], [2, -1, 4, 4], [-4, -12, -2, 7], [-8, -5, -7, -10], [4, 11, 9, 12], [0, -8, 1, -13], [-13, -2, -8, 2],
                        [-3, -2, -2, 3], [-6, 9, -4, -9], [8, 12, 10, 7], [0, 9, 1, 3], [7, -5, 11, -10], [-13, -6, -11, 0], [10, 7, 12, 1],
                        [-6, -3, -6, 12], [10, -9, 12, -4], [-13, 8, -8, -12], [-13, 0, -8, -4], [3, 3, 7, 8], [5, 7, 10, -7], [-1, 7, 1, -12],
                        [3, -10, 5, 6], [2, -4, 3, -10], [-13, 0, -13, 5], [-13, -7, -12, 12], [-13, 3, -11, 8], [-7, 12, -4, 7], [6, -10, 12, 8],
                        [-9, -1, -7, -6], [-2, -5, 0, 12], [-12, 5, -7, 5], [3, -10, 8, -13], [-7, -7, -4, 5], [-3, -2, -1, -7], [2, 9, 5, -11],
                        [-11, -13, -5, -13], [-1, 6, 0, -1], [5, -3, 5, 2], [-4, -13, -4, 12], [-9, -6, -9, 6], [-12, -10, -8, -4], [10, 2, 12, -3],
                        [7, 12, 12, 12], [-7, -13, -6, 5], [-4, 9, -3, 4], [7, -1, 12, 2], [-7, 6, -5, 1], [-13, 11, -12, 5], [-3, 7, -2, -6],
                        [7, -8, 12, -7], [-13, -7, -11, -12], [1, -3, 12, 12], [2, -6, 3, 0], [-4, 3, -2, -13], [-1, -13, 1, 9], [7, 1, 8, -6],
                        [1, -1, 3, 12], [9, 1, 12, 6], [-1, -9, -1, 3], [-13, -13, -10, 5], [7, 7, 10, 12], [12, -5, 12, 9], [6, 3, 7, 11],
                        [5, -13, 6, 10], [2, -12, 2, 3], [3, 8, 4, -6], [2, 6, 12, -13], [9, -12, 10, 3], [-8, 4, -7, 9], [-11, 12, -4, -6],
                        [1, 12, 2, -8], [6, -9, 7, -4], [2, 3, 3, -2], [6, 3, 11, 0], [3, -3, 8, -8], [7, 8, 9, 3], [-11, -5, -6, -4],
                        [-10, 11, -5, 10], [-5, -8, -3, 12], [-10, 5, -9, 0], [8, -1, 12, -6], [4, -6, 6, -11], [-10, 12, -8, 7], [4, -2, 6, 7],
                        [-2, 0, -2, 12], [-5, -8, -5, 2], [7, -6, 10, 12], [-9, -13, -8, -8], [-5, -13, -5, -2], [8, -8, 9, -13], [-9, -11, -9, 0],
                        [1, -8, 1, -2], [7, -4, 9, 1], [-2, 1, -1, -4], [11, -6, 12, -11], [-12, -9, -6, 4], [3, 7, 7, 12], [5, 5, 10, 8],
                        [0, -4, 2, 8], [-9, 12, -5, -13], [0, 7, 2, 12], [-1, 2, 1, 7], [5, 11, 7, -9], [3, 5, 6, -8], [-13, -4, -8, 9],
                        [-5, 9, -3, -3], [-4, -7, -3, -12], [6, 5, 8, 0], [-7, 6, -6, 12], [-13, 6, -5, -2], [1, -10, 3, 10], [4, 1, 8, -4],
                        [-2, -2, 2, -13], [2, -12, 12, 12], [-2, -13, 0, -6], [4, 1, 9, 3], [-6, -10, -3, -5], [-3, -13, -1, 1], [7, 5, 12, -11],
                        [4, -2, 5, -7], [-13, 9, -9, -5], [7, 1, 8, 6], [7, -8, 7, 6], [-7, -4, -7, 1], [-8, 11, -7, -8], [-13, 6, -12, -8],
                        [2, 4, 3, 9], [10, -5, 12, 3], [-6, -5, -6, 7], [8, -3, 9, -8], [2, -12, 2, 8], [-11, -2, -10, 3], [-12, -13, -7, -9],
                        [-11, 0, -10, -5], [5, -3, 11, 8], [-2, -13, -1, 12], [-1, -8, 0, 9], [-13, -11, -12, -5], [-10, -2, -10, 11], [-3, 9, -2, -13],
                        [2, -3, 3, 2], [-9, -13, -4, 0], [-4, 6, -3, -10], [-4, 12, -2, -7], [-6, -11, -4, 9], [6, -3, 6, 11], [-13, 11, -5, 5],
                        [11, 11, 12, 6], [7, -5, 12, -2], [-1, 12, 0, 7], [-4, -8, -3, -2], [-7, 1, -6, 7], [-13, -12, -8, -13], [-7, -2, -6, -8],
                        [-8, 5, -6, -9], [-5, -1, -4, 5], [-13, 7, -8, 10], [1, 5, 5, -13], [1, 0, 10, -13], [9, 12, 10, -1], [5, -8, 10, -9],
                        [-1, 11, 1, -13], [-9, -3, -6, 2], [-1, -10, 1, 12], [-13, 1, -8, -10], [8, -11, 10, -6], [2, -13, 3, -6], [7, -13, 12, -9],
                        [-10, -10, -5, -7], [-10, -8, -8, -13], [4, -6, 8, 5], [3, 12, 8, -13], [-4, 2, -3, -3], [5, -13, 10, -12], [4, -13, 5, -1],
                        [-9, 9, -4, 3], [0, 3, 3, -9], [-12, 1, -6, 1], [3, 2, 4, -8], [-10, -10, -10, 9], [8, -13, 12, 12], [-8, -12, -6, -5],
                        [2, 2, 3, 7], [10, 6, 11, -8], [6, 8, 8, -12], [-7, 10, -6, 5], [-3, -9, -3, 9], [-1, -13, -1, 5], [-3, -7, -3, 4],
                        [-8, -2, -8, 3], [4, 2, 12, 12], [2, -5, 3, 11], [6, -9, 11, -13], [3, -1, 7, 12], [11, -1, 12, 4], [-3, 0, -3, 6],
                        [4, -11, 4, 12], [2, -4, 2, 1], [-10, -6, -8, 1], [-13, 7, -11, 1], [-13, 12, -11, -13], [6, 0, 11, -13], [0, -1, 1, 4],
                        [-13, 3, -9, -2], [-9, 8, -6, -3], [-13, -6, -8, -2], [5, -9, 8, 10], [2, 7, 3, -9], [-1, -6, -1, -1], [9, 5, 11, -2],
                        [11, -3, 12, -8], [3, 0, 3, 5], [-1, 4, 0, 10], [3, -6, 4, 5], [-13, 0, -10, 5], [5, 8, 12, 11], [8, 9, 9, -6], [7, -4, 8, -12],
                        [-10, 4, -10, 9], [7, 3, 12, 4], [9, -7, 10, -2], [7, 0, 12, -2], [-1, -6, 0, -11]]

const freak_sampling_pattern = [[29, 27], [30, 26], [41, 39], [33, 16], [20, 11], [11, 8], [26, 12], [43, 14], [40, 34], [34, 16], [39, 17], [22, 21],
                                [30, 12], [21, 16], [6, 2], [34, 33], [18, 14], [24, 13], [10, 4], [26, 7], [19, 13], [42, 38], [23, 20], [5, 3], [12, 7],
                                [7, 1], [12, 1], [39, 15], [10, 9], [30, 7], [35, 32], [17, 15], [40, 21], [38, 30], [37, 31], [41, 3], [36, 31], [32, 23],
                                [34, 22], [33, 21], [24, 19], [37, 36], [40, 22], [29, 15], [20, 8], [41, 17], [40, 33], [9, 4], [38, 7], [42, 26], [34, 21],
                                [41, 15], [38, 12], [37, 25], [23, 11], [42, 12], [22, 16], [23, 8], [42, 30], [14, 6], [38, 26], [32, 20], [42, 7], [39, 5],
                                [17, 3], [33, 22], [35, 23], [39, 29], [27, 4], [27, 15], [35, 11], [32, 11], [29, 17], [32, 8], [35, 20], [29, 4], [18, 2],
                                [15, 3], [29, 3], [27, 3], [36, 25], [39, 27], [41, 29], [22, 9], [22, 4], [31, 25], [18, 6], [41, 27], [27, 17], [15, 5],
                                [29, 5], [35, 8], [17, 5], [40, 10], [21, 4], [40, 9], [39, 4], [42, 24], [14, 2], [21, 10], [26, 18], [27, 5], [40, 28],
                                [21, 9], [15, 10], [32, 2], [36, 19], [14, 1], [35, 2], [41, 4], [38, 24], [18, 1], [42, 19], [15, 4], [18, 7], [19, 12],
                                [38, 14], [22, 10], [30, 14], [28, 16], [30, 13], [30, 18], [13, 12], [14, 12], [33, 28], [24, 7], [38, 18], [26, 13], [34, 28],
                                [7, 6], [32, 6], [38, 19], [17, 4], [31, 24], [31, 19], [42, 14], [13, 7], [16, 9], [24, 1], [16, 4], [26, 14], [8, 5], [19, 1],
                                [37, 19], [30, 6], [36, 13], [36, 24], [17, 9], [26, 6], [23, 2], [26, 2], [14, 7], [12, 2], [24, 12], [28, 21], [11, 2], [30, 2],
                                [28, 10], [10, 3], [15, 9], [20, 6], [39, 8], [16, 10], [35, 18], [13, 1], [28, 22], [37, 24], [17, 10], [28, 9], [30, 19],
                                [19, 7], [42, 2], [6, 1], [26, 24], [29, 16], [4, 3], [27, 16], [43, 31], [5, 4], [41, 11], [11, 3], [2, 1], [20, 5], [27, 22],
                                [23, 6], [8, 6], [20, 2], [18, 12], [35, 14], [35, 1], [9, 5], [25, 19], [30, 24], [39, 22], [40, 15], [29, 21], [41, 8],
                                [25, 13], [25, 24], [40, 29], [35, 5], [40, 27], [39, 20], [40, 17], [22, 15], [32, 1], [39, 21], [33, 3], [32, 3], [23, 5],
                                [17, 8], [41, 23], [31, 12], [23, 3], [34, 17], [32, 18], [41, 21], [26, 19], [7, 2], [34, 5], [41, 22], [31, 7], [34, 27],
                                [20, 3], [12, 6], [34, 3], [39, 11], [32, 5], [9, 3], [22, 3], [37, 1], [43, 39], [19, 18], [24, 14], [41, 20], [14, 11], [11, 5],
                                [33, 17], [29, 22], [8, 3], [33, 5], [16, 3], [33, 29], [14, 13], [21, 5], [18, 8], [17, 16], [21, 3], [21, 17], [27, 21],
                                [18, 13], [13, 6], [16, 15], [15, 11], [35, 3], [43, 10], [19, 6], [24, 2], [22, 5], [25, 1], [31, 30], [27, 8], [39, 23],
                                [20, 18], [29, 8], [8, 4], [32, 12], [10, 5], [23, 1], [36, 26], [20, 1], [24, 6], [13, 2], [16, 5], [42, 25], [23, 14], [20, 17],
                                [29, 11], [38, 36], [38, 25], [11, 6], [5, 2], [42, 31], [32, 15], [11, 4], [33, 27], [8, 2], [19, 2], [6, 3], [37, 30], [39, 34],
                                [37, 26], [35, 17], [39, 35], [42, 11], [41, 32], [20, 4], [8, 1], [23, 7], [27, 11], [28, 3], [23, 4], [11, 9], [24, 18],
                                [12, 8], [20, 12], [11, 1], [20, 7], [23, 15], [41, 33], [38, 8], [28, 5], [34, 29], [36, 18], [19, 14], [21, 15], [11, 10],
                                [23, 18], [28, 15], [22, 17], [23, 12], [31, 18], [38, 35], [38, 20], [35, 13], [6, 5], [42, 23], [32, 17], [32, 4], [41, 2],
                                [35, 4], [30, 11], [20, 9], [29, 28], [9, 8], [41, 40], [36, 14], [36, 30], [28, 17], [15, 8], [37, 14], [31, 14], [3, 2],
                                [31, 26], [38, 11], [31, 6], [23, 10], [17, 11], [18, 11], [20, 14], [20, 15], [30, 25], [12, 11], [41, 28], [28, 27], [14, 5],
                                [17, 2], [39, 6], [23, 9], [26, 11], [18, 3], [36, 2], [39, 28], [35, 15], [32, 13], [42, 32], [41, 18], [23, 17], [8, 7],
                                [14, 8], [35, 26], [26, 8], [40, 39], [27, 2], [32, 30], [42, 37], [26, 25], [42, 8], [15, 6], [40, 8], [6, 4], [29, 20], [17, 6],
                                [38, 23], [38, 37], [15, 2], [25, 18], [39, 14], [25, 2], [18, 5], [16, 8], [4, 2], [14, 3], [42, 36], [42, 20], [3, 1], [30, 8],
                                [27, 23], [35, 24], [32, 29], [25, 14], [38, 31], [7, 5], [40, 11], [29, 2], [13, 11], [13, 8], [16, 11], [5, 1], [35, 30],
                                [20, 13], [39, 32], [21, 8], [32, 19], [22, 8], [27, 20], [27, 6], [29, 23], [20, 16], [23, 13], [35, 27], [22, 11], [41, 1],
                                [21, 11], [41, 14], [23, 16], [42, 5], [10, 6], [39, 33], [35, 19], [14, 4], [39, 1], [26, 23], [40, 20], [39, 18], [18, 4],
                                [32, 16], [32, 26], [15, 1], [41, 35], [41, 34], [19, 11], [18, 17], [32, 24], [40, 23], [24, 11], [30, 20], [30, 5], [35, 16],
                                [17, 1], [26, 20], [38, 3], [18, 15], [19, 8], [26, 3], [32, 27], [33, 20], [24, 8], [18, 10], [30, 23], [15, 7], [23, 21],
                                [29, 18], [24, 20], [23, 22], [7, 4], [16, 2], [28, 8], [42, 35], [32, 21], [18, 9], [34, 23], [17, 12], [17, 14], [29, 1],
                                [13, 3], [15, 12], [15, 14], [12, 4], [28, 11], [24, 23], [21, 20], [38, 32], [35, 22], [27, 14], [40, 2], [35, 29], [27, 1],
                                [30, 3], [22, 20], [9, 1], [26, 5], [21, 2], [20, 19], [16, 6], [17, 7], [32, 22], [22, 2], [34, 20], [34, 6], [13, 5], [22, 6],
                                [38, 5], [23, 19], [10, 1], [35, 21], [29, 12], [36, 11], [27, 7], [33, 23], [31, 8], [42, 3]]

const freak_orientation_sampling_pattern = [[1, 4], [2, 5], [3, 6], [1, 3], [2, 4], [3, 5], [4, 6], [5, 1], [6, 2], [7, 10], [8, 11], [9, 12], [7, 9], [8, 10],
                                            [9, 11], [10, 12], [11, 7], [12, 8], [13, 16], [14, 17], [15, 18], [13, 15], [14, 16], [15, 17], [16, 18], [17, 13],
                                            [18, 14], [19, 22], [20, 23], [21, 24], [19, 21], [20, 22], [21, 23], [22, 24], [23, 19], [24, 20], [25, 28], [26, 29],
                                            [27, 30], [31, 34], [32, 35], [33, 36], [37, 40], [38, 41], [39, 42]]

const freak_num_circular_pattern = [6, 6, 6, 6, 6, 6, 6, 1]

const freak_largest_radius = 2 / 3

const freak_smallest_radius = 2 / 24

const freak_circular_spacing = (freak_largest_radius - freak_smallest_radius) / 21

const freak_radii = [freak_largest_radius, freak_largest_radius - 6 * freak_circular_spacing, freak_largest_radius - 11 * freak_circular_spacing,
                     freak_largest_radius - 15 * freak_circular_spacing, freak_largest_radius - 18 * freak_circular_spacing,
                     freak_largest_radius - 20 * freak_circular_spacing, freak_smallest_radius, 0]

const freak_sigma = vcat(freak_radii[1:end - 1] / 2, freak_radii[end - 1] / 2)

const freak_orientation_steps = 256;

const brisk_radii = [0, 2.9, 4.9, 7.4, 10.8]

const brisk_sigma = [0.65, 1.16, 1.41, 2.00, 2.19]

const brisk_num_circular_pattern = [1, 10, 14, 15, 20]

const brisk_orientation_steps = 1024

const brisk_dmin = 8.2

const brisk_dmax = 5.85

const brisk_points = 60


text:
const dec_algs = Dict(
    :cholesky => choleskydec,
    :EG => EG,
    :MD => MD,
    :LEX_M => LEX_M,
    :MCS_M => MCS_M,
    :minLC => minLC,
    :maxLC => maxLC,
    :minTrig => minTrig,
    :maxTrig => maxTrig
)

function decompose(g::AbstractGraph, dec_type=:cholesky; kwargs...)
    cg, data = dec_algs[dec_type](g; kwargs...)
    cliques, cliquetree, nb_lc = get_decomposition_data(cg)
    return cliques, cliquetree, cg, nb_lc, data
end

function get_decomposition_data(g::AbstractGraph)
    cliques = maximal_cliques(g)
    cliquetree, nb_lc = get_cliquetree(cliques)
    return cliques, cliquetree, nb_lc
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using NeuronBuilder
using Test

@testset "NeuronBuilder.jl" begin
    # Write your tests here.
end

max_stars_repo_name: sg-s/NeuronBuilder

max_issues_repo_name: sg-s/NeuronBuilder


text:
import Compat: is_unix, is_windows, is_apple

depsfile = joinpath(dirname(@__FILE__),"deps.jl")
if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    f = open(depsfile,"w")
    println(f,"const libgurobi = \"$path\"")
    close(f)
end

aliases = ["gurobi65","gurobi60","gurobi56","gurobi55"]

paths_to_try = copy(aliases)

for a in aliases
    if haskey(ENV, "GUROBI_HOME")
        if is_unix()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"],"lib",string("lib",a,".so")))
        end
        if is_windows()
            push!(paths_to_try, joinpath(ENV["GUROBI_HOME"],"bin",string(a,".",Libdl.dlext)))
        end
    end
    # gurobi uses .so on OS X for some reason
    if is_apple()
        push!(paths_to_try, string("lib$a.so"))
    end
end

found = false
for l in paths_to_try
    d = Libdl.dlopen_e(l)
    if d != C_NULL
        found = true
        write_depsfile(l)
        break
    end
end

if !found
    error("Unable to locate Gurobi installation. Note that this must be downloaded separately from gurobi.com")
end


text:
using RoME
using Test
using DelimitedFiles
using Statistics

# TODO this initial test is likely obsolete
@testset "Test syntax for null hypothesis matrix substitution" begin

function ff(X::Array{Float64,2})
  return X.+1.0
end

global c = Categorical([0.5;0.5])
global A = zeros(3,5)

global skipdos = rand(c, 5)
global dos = skipdos.==2

global B = zeros(3,5)
B[:,dos] = ff(A[:,dos])

@test sum(sum(B[:,dos],dims=2) .== sum(dos),dims=1)[1,1] == 3

end

@testset "Test if null hypothesis occurs as expected..." begin

global N = 100
global fg = initfg()

global initCov = diagm( [1;1;1; 0.01;0.01;0.01].^2 )
global odoCov = deepcopy(initCov)

end

@testset "Adding PriorPose3 to graph..." begin

global v1 = addVariable!(fg, :x1, Pose3, N=N)
global initPosePrior = PriorPose3( MvNormal( zeros(6), initCov) )
global f1  = addFactor!(fg,[v1], initPosePrior, graphinit=true)

# ls(fg, :x1)

ensureAllInitialized!(fg)

end

@testset "Ensure variables initialized properly" begin

# start with to tight an initialization
global muX1 = Statistics.mean(getVal(fg,:x1),dims=2)
global stdX1 = Statistics.std(getVal(fg,:x1),dims=2)
@test sum(map(Int,abs.(muX1[1:3]) .< 0.4)) == 3
@test sum(map(Int,abs.(muX1[4:6]) .< 0.04)) == 3
@test sum(map(Int,abs.(1.0 .- stdX1[1:3]) .< 0.3)) == 3
@test sum(map(Int,abs.(0.01 .- stdX1[4:6]) .< 0.1)) == 3

end

@testset "Testing PriorPose3 evaluation..." begin

global priorpts = approxConv(fg, :x1f1, :x1)
# priorpts = evalFactor(fg, fg.g.vertices[2], 1)
global means = Statistics.mean(priorpts,dims=2)
@test sum(map(Int,abs.(means[1:3]) .> 0.5)) == 0
@test sum(map(Int,abs.(means[4:6]) .> 0.05)) == 0

# global v2, f2 = addOdoFG!(fg, Pose3Pose3( MvNormal([25;0;0;0;0;0.0], odoCov)) )
# global v3, f3 = addOdoFG!(fg, Pose3Pose3( MvNormal([25;0;0;0;0;0.0], odoCov)) )

addVariable!(fg, :x2, Pose3)
addFactor!(fg, [:x1;:x2], Pose3Pose3( MvNormal([25;0;0;0;0;0.0], odoCov)) )

addVariable!(fg, :x3, Pose3)
addFactor!(fg, [:x2;:x3], Pose3Pose3( MvNormal([25;0;0;0;0;0.0], odoCov)) )

end

@testset "Testing Pose3Pose3 evaluation..." begin

global X1pts = getVal(fg, :x1)
global X2pts = approxConv(fg, :x2x3f1, :x2, N=N)
# X2pts = evalFactor(fg, fg.g.vertices[6], 3, N=N)
global X3pts = approxConv(fg, :x2x3f1, :x3)
global X2ptsMean = Statistics.mean(X2pts,dims=2)
global X3ptsMean = Statistics.mean(X3pts,dims=2)

@test  sum(map(Int, abs.(X2ptsMean) - [25.0;0;0;0;0;0] .< 5.0 ))  == 6
@test  sum(map(Int, abs.(X3ptsMean -  [50.0;0;0;0;0;0]) .< 5.0 )) == 6

tree,smt,hist = solveTree!(fg)

end

@testset "Adding Pose3Pose3 w nullhypo to graph..." begin

global odo3 = SE3([-35;0;0], Quaternion(0))
global odoc3 = Pose3Pose3( MvNormal(veeEuler(odo3), odoCov) )

@test length(intersect(DFG.listVariables(fg), [:x1;:x2;:x3])) == 3
# define 50/50% hypothesis
addFactor!(fg,[:x3;:x1],odoc3, nullhypo=0.5)

global X1pts = approxConv(fg, :x3x1f1, :x1, N=N)
global X2pts = approxConv(fg, :x3x1f1, :x3, N=N)

global p1 = manikde!(X1pts, Pose3)
global p2 = manikde!(X2pts, Pose3)

end

@testset "loading validation data for testing." begin

global tstdtdir = dirname(@__FILE__)
global X1ptst = readdlm(joinpath(tstdtdir, "X1ptst.csv"),',')
global X2ptst = readdlm(joinpath(tstdtdir, "X2ptst.csv"),',')

global p1t = manikde!(X1ptst, Pose3)
global p2t = manikde!(X2ptst, Pose3)

# plotKDE([p2t;p2],c=["red";"blue"],dims=[1;2],levels=3)
# kld(marginal(p1,[2]), marginal(p1t,[2]), method=:unscented)

global t1 = mmd(X1ptst, X1pts, AMP.SE3_Manifold, size(X1ptst,2), size(X1pts,2), bw=[0.001;])
global t2 = mmd(X2ptst, X2pts, AMP.SE3_Manifold, size(X1ptst,2), size(X1pts,2), bw=[0.001;])
# TODO Change to mmd
# global t1 = minimum([abs(kld(p1, p1t)[1]) ; abs(kld(p1t, p1)[1])])
# global t2 = minimum([abs(kld(p2, p2t)[1]) ; abs(kld(p2t, p2)[1])])

@test t1 < 0.5
@test t2 < 0.5

end

# using KernelDensityEstimatePlotting
# plotKDE([p1; p1t], dims=[1], c=["black";"red"])