text:
export build_many_body_op,
       build_sparse_many_body_op,
       build_spin1_many_body_op

"""
    build_many_body_op(L::Int, N::Int,
                       J::AbstractMatrix{T},
                       V::AbstractMatrix{T},
                       C::T=zero(T)) where T<:Number

Build many-body operator from the hopping matrix J with N particles.
"""
function build_many_body_op(L::Int, N::Int,
                            J::AbstractMatrix{T},
                            V::AbstractMatrix{T},
                            C::T=zero(T)) where T<:Number
    basis = make_LN_basis(L, N)
    return build_many_body_op(L, basis, J, V, C)
end

"""
    build_many_body_op(L::Int,
                       basis::Vector{Int},
                       J::AbstractMatrix{T},
                       V::AbstractMatrix{T},
                       C::T=zero(T)) where T<:Number

Build many-body operator from the hopping matrix J with a given basis of states.

# Arguments:
- L::Int: number of sites.
- basis::Vector{Int}: basis of Fock states for the mb operator.
- J::AbstractMatrix{T}: hopping matrix b^dagger_i*b_j.
- V::AbstractMatrix{T}: interaction matrix n_i*n_j.
- C::T=zero(T): constant term.
"""
function build_many_body_op(L::Int,
                            basis::Vector{Int},
                            J::AbstractMatrix{T},
                            V::AbstractMatrix{T},
                            C::T=zero(T)) where T<:Number
    # Dimension of the Hilbert space.
    dH = length(basis)

    # Use only upper triangular part of `V`.
    V = deepcopy(V)
    for i=1:L-1, j=i+1:L
        V[i, j] += V[j, i]
    end

    # Make the many-body operator of the same type as J.
    Op = zeros(T, (dH, dH))

    for s=1:dH
        state = basis[s]
        # Constant term.
        Op[s, s] += C

        # Diagonal terms.
        for i=0:L-1
            if (state>>i)&1 == 1
                Op[s, s] += J[i+1, i+1]
            end
            for j=i+1:L-1
                if !iszero(V[i+1, j+1]) && (state>>i)&1==1 && (state>>j)&1==1
                    Op[s, s] += V[i+1, j+1]
                end
            end
        end

        # Off-diagonal terms.
        for i=0:L-1
            for j=0:L-1
                if (!iszero(J[i+1, j+1])
                        && (state>>i)&1 == 0 && (state>>j)&1 == 1)
                    Jstate = state + 1<<i - 1<<j
                    js = searchsortedfirst(basis, Jstate)
                    Op[js, s] = J[i+1, j+1]
                end
            end
        end
    end
    return Op
end

"""
    build_sparse_many_body_op(L::Int,
                              basis::Vector{Int},
                              J::AbstractMatrix{T},
                              V::AbstractMatrix{T},
                              C::T=zero(T)) where {T<:Number}

Build sparse many-body operator from the matrix J with a N particles.

# Arguments:
- L::Int: number of sites.
- basis::Vector{Int}: basis of Fock states for the mb operator.
- J::AbstractMatrix{T}: hopping matrix b^dagger_i*b_j.
- V::AbstractMatrix{T}: interaction matrix n_i*n_j.
- C::T=zero(T): constant term.
"""
function build_sparse_many_body_op(L::Int,
                                   basis::Vector{Int},
                                   J::AbstractMatrix{T},
                                   V::AbstractMatrix{T},
                                   C::T=zero(T)) where {T<:Number}
    # Dimension of the Hilbert space.
    dH = length(basis)

    # Use only upper triangular part of `V`.
    V = deepcopy(V)
    for i=1:L-1, j=i+1:L
        V[i, j] += V[j, i]
    end

    # Number of non-zero off-diagonal elts in J.
    nnz_J = count(abs.(J) .> 1e-8) - count(abs.(diag(J)) .> 1e-8)
    # Number of nnz elts in Op due to off-diag terms and diag terms.
    nnz_in_Op = dH*(nnz_J + 1)
    rows = zeros(Int, nnz_in_Op)
    cols = zeros(Int, nnz_in_Op)
    vals = zeros(T, nnz_in_Op)

    cont = 1
    for s=1:dH
        state = basis[s]
        # Constant term.
        vals[cont] += C

        # Diagonal term.
        for i=0:L-1
            if (state>>i)&1 == 1
                vals[cont] += J[i+1, i+1]
            end
            for j=i+1:L-1
                if !iszero(V[i+1, j+1]) && (state>>i)&1==1 && (state>>j)&1==1
                    vals[cont] += V[i+1, j+1]
                end
            end
        end
        rows[cont] = s
        cols[cont] = s
        cont += 1

        # Off-diagonal terms.
        for i=0:L-1
            for j=0:L-1
                if (!iszero(J[i+1, j+1])
                        && (state>>i)&1 == 0 && (state>>j)&1 == 1)
                    Jstate = state + 1<<i - 1<<j
                    js = searchsortedfirst(basis, Jstate)
                    rows[cont] = js
                    cols[cont] = s
                    vals[cont] = J[i+1, j+1]
                    cont += 1
                end
            end
        end
    end
    # Remove zeros from rows, cols, vals.
    rows = rows[1:cont-1]
    cols = cols[1:cont-1]
    vals = vals[1:cont-1]
    Op = sparse(rows, cols, vals, dH, dH)
    return Op
end

"""
    build_spin1_many_body_op(L::Int,
                             Sz::Int,
                             J::AbstractMatrix{T},
                             W::AbstractMatrix{T},
                             Jz::Vector{T},
                             Jx::Vector{T},
                             C::T=zero(T)) where T<:Number

Build a spin 1 many-body operator from the hopping matrix J.
"""
function build_spin1_many_body_op(L::Int,
                                  Sz::Int,
                                  J::AbstractMatrix{T},
                                  W::AbstractMatrix{T},
                                  Jz::Vector{T},
                                  Jx::Vector{T},
                                  C::T=zero(T)) where T<:Number
    basis = make_spin1_LSz_basis(L, Sz)
    return build_spin1_many_body_op(L, basis, J, W, Jz, Jx, C)
end

"""
    build_spin1_many_body_op(L::Int,
                             basis::Vector{Int},
                             J::AbstractMatrix{T},
                             W::AbstractMatrix{T},
                             Jz::Vector{T},
                             Jx::Vector{T},
                             C::T=zero(T)) where T<:Number

Build a spin 1 many-body operator from the hopping matrix J.

# Arguments:
- L::Int: number of sites.
- basis::Vector{Int}: basis of Fock states for the mb operator.
- J::AbstractMatrix{T}: hopping matrix, ignore diagonal values S^+_i*S^-_j.
- W::AbstractMatrix{T}: interaction matrix S^z_i*S^z_j.
- Jz::Vector{T}: local Sz terms S^z_i.
- Jx::Vector{T}: local Sx terms S^x_i.
- C::T=zero(T): constant term.
"""
function build_spin1_many_body_op(L::Int,
                                  basis::Vector{Int},
                                  J::AbstractMatrix{T},
                                  W::AbstractMatrix{T},
                                  Jz::Vector{T},
                                  Jx::Vector{T},
                                  C::T=zero(T)) where T<:Number
    # Dimension of the Hilbert space.
    dH = length(basis)

    # Make the many-body operator of the same type as J.
    Op = zeros(T, (dH, dH))

    for s=1:dH
        state = basis[s]
        # Constant term.
        Op[s, s] += C

        # Sz and Sz^2 terms.
        for i=0:L-1
            Op[s, s] += Jz[i+1]*((state>>2i)&1 - (state>>(2i+1))&1)
            Op[s, s] += W[i+1, i+1]*((state>>2i)&1 + (state>>(2i+1))&1)
        end

        # S^z_i S^z_j terms, i < j.
        for j=1:L-1
            for i=0:j-1
                iszero(W[i+1, j+1]) && continue
                # Four possibilities to find the state. In positions i and j,
                # respectively: (+, +), (-, -), (-, +), (+, -).
                if (((state>>2i)&3 == 1 && (state>>2j)&3 == 1) ||
                    ((state>>2i)&3 == 2 && (state>>2j)&3 == 2))
                    Op[s, s] += W[i+1, j+1]
                elseif (((state>>2i)&3 == 2 && (state>>2j)&3 == 1) ||
                        ((state>>2i)&3 == 1 && (state>>2j)&3 == 2))
                    Op[s, s] -= W[i+1, j+1]
                end
            end
        end

        # S^+_i S^-_j, with i != j, terms.
        for i=0:L-1
            for j=0:L-1
                i == j && continue
                iszero(J[i+1, j+1]) && continue
                # Four possibilities to find the state. In positions i and j,
                # respectively: (0, +), (0, 0), (-, +), (-, 0).
                if (state>>2i)&3 == 0 && (state>>2j)&3 == 1
                    Jstate = state + 1<<2i - 1<<2j
                    js = searchsortedfirst(basis, Jstate)
                    Op[js, s] += 2*J[i+1, j+1]
                elseif (state>>2i)&3 == 0 && (state>>2j)&3 == 0
                    Jstate = state + 1<<2i + 2<<2j
                    js = searchsortedfirst(basis, Jstate)
                    Op[js, s] += 2*J[i+1, j+1]
                elseif (state>>2i)&3 == 2 && (state>>2j)&3 == 1
                    Jstate = state - 2<<2i - 1<<2j
                    js = searchsortedfirst(basis, Jstate)
                    Op[js, s] += 2*J[i+1, j+1]
                elseif (state>>2i)&3 == 2 && (state>>2j)&3 == 0
                    Jstate = state - 2<<2i + 2<<2j
                    js = searchsortedfirst(basis, Jstate)
                    Op[js, s] += 2*J[i+1, j+1]
                end
            end
        end
    end

    # Sx terms.
    if norm(Jx) > 1e-8
        for s=1:dH
            state = basis[s]

            # S^x_i terms.
            for i=0:L-1
                iszero(Jx[i+1]) && continue
                if (state>>2i)&3 == 0
                    # state_i = |1, 0>.
                    Jstate1 = state + 1<<2i # New state: |1, 1>.
                    js1 = searchsortedfirst(basis, Jstate1)
                    # If new state is in base.
                    if js1 <= dH && Jstate1 == basis[js1]
                        Op[js1, s] += Jx[i+1]/sqrt(2)
                    end

                    Jstate2 = state + 2<<2i # New state: |1, -1>.
                    js2 = searchsortedfirst(basis, Jstate2)
                    if js2 <= dH && Jstate2 == basis[js2]
                        Op[js2, s] += Jx[i+1]/sqrt(2)
                    end
                elseif (state>>2i)&3 == 1
                    # state_i = |1, 1>.
                    Jstate = state - 1<<2i # New state: |1, 0>.
                    js = searchsortedfirst(basis, Jstate)
                    if js <= dH && Jstate == basis[js]
                        Op[js, s] += Jx[i+1]/sqrt(2)
                    end
                elseif (state>>2i)&3 == 2
                    # state_i = |1, -1>.
                    Jstate = state - 2<<2i # New state: |1, 0>.
                    js = searchsortedfirst(basis, Jstate)
                    if js <= dH && Jstate == basis[js]
                        Op[js, s] += Jx[i+1]/sqrt(2)
                    end
                end
            end

        end
    end
    return Op
end


text:
using PrettyPlots, Plots, DataFrames
#gr();
#pgfplots();
pyplot();

# TODO consider moving this to PrettyPlots
default(guidefont = font(17), tickfont = font(15), legendfont = font(12), titlefont = font(20))
c=defineCase(;(:mode=>:caseStudyPath));

plotSettings(;(:simulate=>true),(:plant=>true),(:plantOnly=>false));

description = string(
"In this test: \n",c.m.name,"\n
* m.Ni=",c.m.Ni," \n
* m.Nck=",c.m.Nck,"\n
* m.tp=",c.m.tp," \n
* m.tex=",c.m.tex,"\n
* m.max_cpu_time=",c.m.max_cpu_time," \n
* m.Nck=",c.m.Nck,"\n
* m.PredictX0=",c.m.PredictX0," \n
* m.FixedTp=",c.m.FixedTp,"\n
")
results_dir=string("test1_10_obstacles")
resultsDir!(r,results_dir;description=description);
savePlantData(n,r)

if _pretty_defaults[:simulate];
  println("Plotting the Final Results!")
  mainSim(n,r,c,pa;(:mode=>:path3))
end

# TODO make this an example
if n.r.status==:Infeasible
  #Settings(;evalConstraints=true,save=true,MPC=false);
  n.s.evalConstraints=true;n.s.save=true;n.s.MPC=false;
  postProcess!(n)
  # trouble getting a feasible solution? -> look at the constraints
  print(n.r.constraint.value)
end # then consider relaxing tolerances etc. to make it work
# then looking at the output
# These are the dual infeasibilities
#=
│ Row │ step │ x0_con   │
├─────┼──────┼──────────┤
│ 1   │ 1    │ -267.777 │
│ 2   │ 2    │ 64.6356  │
│ 3   │ 3    │ 271.738  │
│ 4   │ 4    │ 2.44213  │
│ 5   │ 5    │ 3496.5   │
=#
# These are the dual infeasibilities

# we cam try to relax the constraint on the

#=
# define tolerances
XF_tol=[NaN,NaN,NaN,NaN,NaN];
X0_tol=[0.05,0.05,0.05,0.05,0.01];
defineTolerances(n;X0_tol=X0_tol,XF_tol=XF_tol);
=#

#using PrettyPlots
#s=Settings(;save=true,MPC=false,simulate=false,format=:png);
#pSimPath(n,r,s,c,r.eval_num)

ext: jl


text:
# Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
# Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
# pp. 829–858 (2017).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    DavisYinIteration(; <keyword-arguments>)

Instantiate the Davis-Yin splitting algorithm (see [1]) for solving
convex optimization problems of the form

    minimize f(x) + g(x) + h(A x),

where `h` is smooth and `A` is a linear mapping (for example, a matrix).

# Arguments
- `x0`: initial point.
- `f=Zero()`: proximable objective term.
- `g=Zero()`: proximable objective term.
- `h=Zero()`: smooth objective term.
- `A=I`: linear operator (e.g. a matrix).
- `Lh=nothing`: Lipschitz constant of the gradient of x ↦ h(Ax).
- `gamma=nothing`: stepsize to use, defaults to `1/Lh` if not set (but `Lh` is).

# References
- [1] Davis, Yin. "A Three-Operator Splitting Scheme and its Optimization
Applications", Set-Valued and Variational Analysis, vol. 25, no. 4,
pp. 829–858 (2017).
"""

@Base.kwdef struct DavisYinIteration{R,C<:Union{R,Complex{R}},T<:AbstractArray{C},Tf,Tg,Th,TA}
    f::Tf = Zero()
    g::Tg = Zero()
    h::Th = Zero()
    A::TA = I
    x0::T
    lambda::R = real(eltype(x0))(1)
    Lh::Maybe{R} = nothing
    gamma::Maybe{R} = Lh !== nothing ? (1 / Lh) : error("You must specify either Lh or gamma")
end

Base.IteratorSize(::Type{<:DavisYinIteration}) = Base.IsInfinite()

struct DavisYinState{T,S}
    z::T
    xg::T
    y::S
    grad_h_y::S
    z_half::T
    xf::T
    res::T
end

function Base.iterate(iter::DavisYinIteration)
    z = copy(iter.x0)
    xg, = prox(iter.g, z, iter.gamma)
    y = iter.A * xg
    grad_h_y, = gradient(iter.h, y)
    z_half = 2 .* xg .- z .- iter.gamma .* (iter.A' * grad_h_y)

    # TODO: check the following
    xf = similar(z_half)
    prox!(xf, iter.f, z_half, iter.gamma)
    # xf, = prox(iter.f, z_half, iter.gamma)  # this fails the tests
    ############################################################################

    res = xf - xg
    z .+= iter.lambda .* res
    state = DavisYinState(z, xg, y, grad_h_y, z_half, xf, res)
    return state, state
end

function Base.iterate(iter::DavisYinIteration, state::DavisYinState)
    prox!(state.xg, iter.g, state.z, iter.gamma)
    mul!(state.y, iter.A, state.xg)
    gradient!(state.grad_h_y, iter.h, state.y)
    mul!(state.z_half, iter.A', state.grad_h_y)
    state.z_half .*= -iter.gamma
    state.z_half .+= 2 .* state.xg .- state.z
    prox!(state.xf, iter.f, state.z_half, iter.gamma)
    state.res .= state.xf .- state.xg
    state.z .+= iter.lambda .* state.res
    return state, state
end

# Solver

struct DavisYin{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::DavisYin)(x0; kwargs...)
    stop(state::DavisYinState) = norm(state.res, Inf) <= solver.tol
    disp((it, state)) = @printf("%5d | %.3e\n", it, norm(state.res, Inf))
    iter = DavisYinIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.xf, state_final.xg, num_iters
end

DavisYin(; maxit=10_000, tol=1e-8, verbose=false, freq=100, kwargs...) =
    DavisYin(maxit, tol, verbose, freq, kwargs)


text:
using Random, LinearAlgebra

# Standardised tests for the eigen decomposition of a square kronecker product
function eigen_tests(rng, C::AbstractKroneckerProduct)

    # Check approximate correctness of decomposition
    λ, Γ = eigen(C)
    @test Γ * Diagonal(λ) * Γ' ≈ C

    # Check approximate correctness of shifted decomposition
    σ² = abs2(randn(rng))
    λ′, Γ′ = (eigen(C) + σ² * I)
    @test Γ′ * Diagonal(λ′) * Γ′' ≈ Matrix(C + σ² * I)
    @test σ² * I + C == C + σ² * I

    # Test with decomposition
    v = randn(rng, length(λ))
    @test eigen(C) \ v ≈ Float64.(Matrix(C)) \ v

    # Test various linear algebra operations with decomposition
    C_dense = Float64.(Matrix(C))
    @test det(eigen(C)) ≈ det(C_dense)
    @test logdet(eigen(C)) ≈ logdet(C_dense)
    @test Matrix(inv(eigen(C))) ≈ inv(C_dense)
end

@testset "eigen" begin
    rng, P, Q, R = MersenneTwister(123456), 3, 5, 7

    # Generate some positive definite matrices so that logdet can be tested.
    A_, B_, C_ = randn(rng, P, P), randn(rng, Q, Q), randn(rng, R, R)
    A, B, C = Symmetric(A_ * A_' + I), Symmetric(B_ * B_' + I), Symmetric(C_ * C_' + I)

    D = kronecker(A, B)
    eigen_tests(rng, D)
    eigen_tests(rng, kronecker(D, C))
    eigen_tests(rng, kronecker(C, D))
    eigen_tests(rng, kronecker(kronecker(A, B), kronecker(B, A)))
    eigen_tests(rng, kronecker(A, 4))

    function test_non_square_extension()
        local n, m, A, B, K, M
        n, m = 3, 5
        # 1. K is square, while A, B aren't
        A = randn(n, m)
        B = randn(m, n)
        K = A ⊗ B
        M = Matrix(K)
        E = eigen(K)
        F = eigen(M)
        @test E.values ≈ F.values
        @test det(E) ≈ det(F)
        @test Matrix(E) ≈ K

        # 2. K is not square
        A = randn(n, m)
        B = randn(m, m)
        K = A ⊗ B
        @test_throws DimensionMismatch eigen(K)
    end
    test_non_square_extension()
end


text:
struct RandomForest
  forests::Vector
end

function train(d::Data{Normalized}, ::Val{:randomforest})
  M = length(d.y_T)

  forests = map(1:M) do m
    X_m = d.X_T[m]
    y_m = d.y_T[m]
    n_m = length(y_m)

    rf = ensemble.RandomForestRegressor(n_estimators=100)
    ScikitLearn.fit!(rf, X_m, y_m)

    rf
  end

  RandomForest(forests)
end

function getoutcomes(new_data::Data{Normalized}, rf::RandomForest)
  M = length(rf.forests)
  n = length(new_data.y)
  outcomes = Matrix{Float64}(undef, n, M)

  for m = 1:M
    outcomes[:, m] = ScikitLearn.predict(rf.forests[m], new_data.X)
  end

  outcomes
end


text:
#= test/pqr.jl
=#

  @testset "pqr" begin
    m = 128
    n =  64
    M = matrixlib(:fourier, rand(m), rand(n))
    opts = LRAOptions(maxdet_tol=0., sketch_randn_niter=1)

    @time for (t, s) in ((:none,                 :none ),
                         (:RandomGaussian,       :randn),
                         (:RandomSubset,         :sub  ),
                         (:SRFT,                 :srft ),
                         (:SparseRandomGaussian, :sprn ))
      opts.sketch = s
      for T in (Float32, Float64, ComplexF32, ComplexF64)
        let A
            println("  $t/$T")

            rtol = 5*eps(real(T))
            approx_rtol = 100*rtol
            opts.rtol = rtol

            A = convert(Array{T}, T <: Real ? real(M) : M)
            nrm = norm(A)

            F  = pqrfact(:n, A, opts)
            Fc = pqrfact(:c, A, opts)
            @test norm(A  - Matrix(F )) < approx_rtol*nrm
            @test norm(A' - Matrix(Fc)) < approx_rtol*nrm

            xm = rand(T, m)
            xn = rand(T, n)
            y = A  *xn; @test norm(y - F  *xn) < approx_rtol*norm(y)
            y = A' *xm; @test norm(y - F' *xm) < approx_rtol*norm(y)
            y = transpose(A)*xm; @test norm(y - transpose(F)*xm) < approx_rtol*norm(y)

            Bm = rand(T, m, m)
            Bn = rand(T, n, n)
            C = A   *Bn  ; @test norm(C - F   *Bn  ) < approx_rtol*norm(C)
            C = A   *Bn' ; @test norm(C - F   *Bn' ) < approx_rtol*norm(C)
            C = A * transpose(Bn); @test norm(C - F   *transpose(Bn)) < approx_rtol*norm(C)
            C = A'  *Bm  ; @test norm(C - F'  *Bm  ) < approx_rtol*norm(C)
            C = A'  *Bm' ; @test norm(C - F'  *Bm' ) < approx_rtol*norm(C)
            C = transpose(A) *Bm  ; @test norm(C - transpose(F) *Bm  ) < approx_rtol*norm(C)
            C = transpose(A) *transpose(Bm); @test norm(C - transpose(F) *transpose(Bm)) < approx_rtol*norm(C)
            C = Bm  *A   ; @test norm(C - Bm  *F   ) < approx_rtol*norm(C)
            C = Bn  *A'  ; @test norm(C - Bn  *F'  ) < approx_rtol*norm(C)
            C = Bn  *transpose(A) ; @test norm(C - Bn  *transpose(F) ) < approx_rtol*norm(C)
            C = Bm' *A   ; @test norm(C - Bm' *F   ) < approx_rtol*norm(C)
            C = Bn' *A'  ; @test norm(C - Bn' *F'  ) < approx_rtol*norm(C)
            C = transpose(Bm)*A   ; @test norm(C - transpose(Bm)*F   ) < approx_rtol*norm(C)
            C = transpose(Bn)*transpose(A) ; @test norm(C - transpose(Bn)*transpose(F) ) < approx_rtol*norm(C)

            x = A*rand(T, n)
            y = F\x; @test norm(x - A*y) < approx_rtol*norm(x)

            X = A*rand(T, n, n)
            C = F\X; @test norm(X - A*C) < approx_rtol*norm(X)
        end
      end
    end
  end


text:
import IterativeSolvers

struct GMRESSolver{L,T} <: LinearMap{T}
    linear_operator::L
    maxiter::Int
    restart::Int
    tol::T
    verbose::Bool
end

Base.axes(A::GMRESSolver) = reverse(axes(A.linear_operator))

function GMRESSolver(op; maxiter=0, restart=0, tol=sqrt(eps(real(eltype(op)))), verbose=true)

    m, n = size(op)
    @assert m == n

    maxiter == 0 && (maxiter = div(n, 5))
    restart == 0 && (restart = n)

    GMRESSolver(op, maxiter, restart, tol, verbose)
end

operator(solver::GMRESSolver) = solver.linear_operator

function solve(solver::GMRESSolver, b)

    T = promote_type(eltype(solver), eltype(b))
    # x = PseudoBlockVector{T}(undef, BlockArrays.blocksizes(solver)[1])
    x = similar(Array{T}, axes(solver)[2])
    fill!(x,0)
    x, ch = solve!(x, solver, b)

    # op = operator(solver)
    # x, ch = IterativeSolvers.gmres(op, b, log=true,  maxiter=solver.maxiter,
    #     restart=solver.restart, reltol=solver.tol, verbose=solver.verbose)
    # return x, ch
end

function solve!(x, solver::GMRESSolver, b)
    op = operator(solver)
    x, ch = IterativeSolvers.gmres!(x, op, b, log=true,  maxiter=solver.maxiter,
        restart=solver.restart, reltol=solver.tol, verbose=solver.verbose)
    return x, ch
end

function Base.:*(A::GMRESSolver, b::AbstractVector)

    T = promote_type(eltype(A), eltype(b))
    y = PseudoBlockVector{T}(undef, BlockArrays.blocksizes(A)[1])

    mul!(y, A, b)

    # x, ch = solve(solver, b)
    # println("Number of iterations: ", ch.iters)
    # ch.isconverged || error("Iterative solver did not converge.")
    # return x
end

Base.size(solver::GMRESSolver) = reverse(size(solver.linear_operator))

function LinearAlgebra.mul!(y::AbstractVecOrMat, solver::GMRESSolver, x::AbstractVector)
    fill!(y,0)
    y, ch = solve!(y, solver, x)
    solver.verbose && println("Number of iterations: ", ch.iters)
    ch.isconverged || error("Iterative solver did not converge.")
    return y
end

LinearAlgebra.adjoint(A::GMRESSolver) = GMRESSolver(adjoint(A.linear_operator), A.maxiter, A.restart, A.tol, A.verbose)
LinearAlgebra.transpose(A::GMRESSolver) = GMRESSolver(transpose(A.linear_operator), A.maxiter, A.restart, A.tol, A.verbose)

# function LinearAlgebra.mul!(
#     y::AbstractVecOrMat,
#     transA::LinearMaps.AdjointMap{<:Any,<:GMRESSolver},
#     x::AbstractVector)

#     LinearMaps.check_dim_mul(y, transA, x)
#     A = transA.lmap
#     B = GMRESSolver(adjoint(A.linear_operator), A.maxiter, A.restart, A.tol, A.verbose)
#     return mul!(y, B, x)
# end

# function LinearAlgebra.mul!(
#     y::AbstractVecOrMat,
#     transA::LinearMaps.TransposeMap{<:Any,<:GMRESSolver},
#     x::AbstractVector)

#     LinearMaps.check_dim_mul(y, transA, x)
#     A = transA.lmap
#     B = GMRESSolver(transpose(A.linear_operator), A.maxiter, A.restart, A.tol, A.verbose)
#     return mul!(y, B, x)
# end

function gmres(eq::DiscreteEquation; maxiter=0, restart=0, tol=0)

    test_space_dict  = eq.test_space_dict
    trial_space_dict = eq.trial_space_dict

    lhs = eq.equation.lhs
    rhs = eq.equation.rhs

    b = assemble(rhs, test_space_dict)
    Z = assemble(lhs, test_space_dict, trial_space_dict)

    block_sizes = zeros(Int, length(trial_space_dict))
    for (p,x) in eq.trial_space_dict
        block_sizes[p] = numfunctions(x)
    end

    T = promote_type(eltype(Z), eltype(b))
    x = PseudoBlockVector{T}(undef, block_sizes)
    fill!(x, 0)

    if tol == 0
        invZ = GMRESSolver(Z, maxiter=maxiter, restart=restart)
    else
        invZ = GMRESSolver(Z, maxiter=maxiter, restart=restart, tol=tol)
    end
    # x = invZ * b
    mul!(x, invZ, b)
end


text:
import Base: put!, take!, pathof, readdir, log

abstract type AbstractModelConfig end

"""
    struct ModelConfig <: AbstractModelConfig

```
model :: Union{Function,String,Pkg.Types.PackageSpec} = "anonymous"
configuration :: Union{Function,String} = "anonymous"
options :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
channel :: Channel{Any} = Channel{Any}(10)
folder :: String = tempdir()
ID :: UUID = UUIDs.uuid4()
```
"""
Base.@kwdef struct ModelConfig <: AbstractModelConfig
    model :: Union{Function,String,Pkg.Types.PackageSpec} = "anonymous"
    configuration :: Union{Function,String} = "anonymous"
    options :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    inputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    outputs :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    status :: OrderedDict{Any,Any} = OrderedDict{Any,Any}()
    channel :: Channel{Any} = Channel{Any}(10)
    folder :: String = tempdir()
    ID :: UUID = UUIDs.uuid4()
end

"""
    pathof(x::AbstractModelConfig)

Returns the run directory path for x ; i.e. joinpath(x.folder,string(x.ID))
"""
pathof(x::AbstractModelConfig) = joinpath(x.folder,string(x.ID))

"""
    readdir(x::AbstractModelConfig)

Same as readdir(pathof(x)).
"""
readdir(x::AbstractModelConfig) = readdir(pathof(x))

"""
    setup(x)

Defaults to `default_ClimateModelSetup(x)`. Can be expected to be
specialized for most concrete types of `AbstractModelConfig`

```
f=ClimateModels.RandomWalker
tmp=ModelConfig(model=f)
setup(tmp)
```
"""
setup(x :: AbstractModelConfig) = default_ClimateModelSetup(x)

function default_ClimateModelSetup(x::AbstractModelConfig)
    !isdir(joinpath(x.folder)) ? mkdir(joinpath(x.folder)) : nothing
    pth=pathof(x); !isdir(pth) ? mkdir(pth) : nothing

    if isa(x.model,Pkg.Types.PackageSpec)
        hasfield(Pkg.Types.PackageSpec,:url) ? url=x.model.url : url=x.model.repo.source
        @suppress Pkg.develop(url=url)
        if x.configuration=="anonymous"
            put!(x.channel,run_the_tests)
        else
            put!(x.channel,x.configuration)
        end
    elseif isa(x.model,Function)
        put!(x.channel,x.model)
    elseif isa(x.configuration,Function)
        put!(x.channel,x.configuration)
    else
        nothing
    end

    !isdir(joinpath(pth,"log")) ? git_log_init(x) : nothing
    git_log_prm(x)

    fil_in=ClimateModels.Pkg.project().path
    fil_out=joinpath(pth,"log","Project.toml")
    cp(fil_in,fil_out)
    git_log_fil(x,fil_out,"add Project.toml to log")
    fil_in=joinpath(dirname(fil_in),"Manifest.toml")
    fil_out=joinpath(pth,"log","Manifest.toml")
    cp(fil_in,fil_out)
    git_log_fil(x,fil_out,"add Manifest.toml to log")

    return x
end

"""
    git_log_init(x :: AbstractModelConfig)

Create `log` subfolder, initialize git, and commit initial README.md
"""
function git_log_init(x :: AbstractModelConfig)
    !isdir(joinpath(x.folder)) ? mkdir(joinpath(x.folder)) : nothing
    p=pathof(x)
    !isdir(p) ? mkdir(p) : nothing
    p=joinpath(pathof(x),"log")
    !isdir(p) ? mkdir(p) : nothing

    f=joinpath(p,"README.md")
    q=pwd()
    cd(p)

    msg=("## Initial Setup\n\n",
    "ID            = ```"*string(x.ID)*"```\n\n",
    "model         = ```"*string(x.model)*"```\n\n",
    "configuration = ```"*string(x.configuration)*"```\n\n")
    open(f, "w") do io
        write(io, msg...)
    end

    try
        @suppress run(`$(git()) init -b main`)
    catch e
        @suppress run(`$(git()) init`)
    end
    run(`$(git()) add README.md`)
    try
        @suppress run(`$(git()) commit README.md -m "initial setup"`)
    catch e
        run(`$(git()) config user.email "you@example.com"`)
        run(`$(git()) config user.name "Your Name"`)
        @suppress run(`$(git()) commit README.md -m "initial setup"`)
    end

    cd(q)
end

"""
    run_the_tests(x)

Default for launching model when it is a cloned julia package

```jldoctest; output = false
using ClimateModels, Pkg
tmp0=PackageSpec(url="https://github.com/JuliaOcean/AirSeaFluxes.jl")
tmp1=ModelConfig(model=tmp0)
setup(tmp1)
build(tmp1)
launch(tmp1)

clean(tmp1)=="no task left in pipeline"

# output

true
```
"""
function run_the_tests(x)
    hasfield(Pkg.Types.PackageSpec,:url) ? url=x.model.url : url=x.model.repo.source
    try
        @suppress Pkg.test(split(url,"/")[end][1:end-3])
    catch e
        txt=split(url,"/")[end][1:end-3]
        println("could not run Pkg.test($txt)")
    end
end

"""
    build(x)

Defaults to `default_ClimateModelBuild(x)`. Can be expected to be
specialized for most concrete types of `AbstractModelConfig`

```jldoctest; output = false
using ClimateModels
tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
build(tmp)

isa(tmp,AbstractModelConfig)

# output

true
```
"""
build(x :: AbstractModelConfig) = default_ClimateModelBuild(x)

"""
    compile(x)

Defaults to `default_ClimateModelBuild(x)`. Can be expected to be
specialized for most concrete types of `AbstractModelConfig`

```jldoctest; output = false
using ClimateModels, Pkg
tmp0=PackageSpec(url="https://github.com/JuliaOcean/AirSeaFluxes.jl")
tmp=ModelConfig(model=tmp0)
setup(tmp)
compile(tmp)

isa(tmp,AbstractModelConfig)

# output

true
```
"""
compile(x :: AbstractModelConfig) = default_ClimateModelBuild(x)

function default_ClimateModelBuild(x::AbstractModelConfig)
    isa(x.model,Pkg.Types.PackageSpec) ? build_the_pkg(x) : nothing
end

"""
   build_the_pkg(x)

Default for building/compiling model when it is a cloned julia package
"""
function build_the_pkg(x)
    @suppress Pkg.build()
end

"""
    launch(x)

Defaults to `default_ClimateModelLaunch(x)` which consists in `take!(x)`
for `AbstractModelConfig`. Can be expected to be specialized for most
concrete types of `AbstractModelConfig`

```
f=ClimateModels.RandomWalker
tmp=ModelConfig(model=f)
setup(tmp)
build(tmp)
launch(tmp)
```
"""
launch(x :: AbstractModelConfig) = default_ClimateModelLaunch(x)

function default_ClimateModelLaunch(x::AbstractModelConfig)
    !isempty(x.channel) ? take!(x) : "no task left in pipeline"
end

"""
    clean(x :: AbstractModelConfig)

Cancel any remaining task (x.channel) and rm the run directory (pathof(x))

```
tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
clean(tmp)
```
"""
function clean(x :: AbstractModelConfig)
    #cancel any remaining task
    while !isempty(x.channel)
        take!(x.channel)
    end
    #clean up run directory
    isdir(pathof(x)) ? rm(pathof(x),recursive=true) : nothing
    #
    return "no task left in pipeline"
end

"""
    monitor(x)

Show `x.status[end]` by default.

```
tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
monitor(tmp)
```
"""
function monitor(x :: AbstractModelConfig)
    k=keys(x.status)
    n=length(k)
    if n>0
        for i in k
            j=x.status[i]
            println(i*" = $j")
        end
    else
        println("no status information")
    end
end

help(x :: AbstractModelConfig) = println("Please consider using relevant github issue trackers for questions")

"""
    show(io::IO, z::AbstractModelConfig)

```
tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
show(tmp)
```
"""
function Base.show(io::IO, z::AbstractModelConfig)
    printstyled(io, "  ID            = ",color=:normal)
    printstyled(io, "$(z.ID)\n",color=:slateblue1)
    printstyled(io, "  model         = ",color=:normal)
    if isa(z.model,Pkg.Types.PackageSpec)
        hasfield(Pkg.Types.PackageSpec,:url) ? url=z.model.url : url=z.model.repo.source
        printstyled(io, "$(url)\n",color=:slateblue1)
    else
        printstyled(io, "$(z.model)\n",color=:slateblue1)
    end
    printstyled(io, "  configuration = ",color=:normal)
    printstyled(io, "$(z.configuration)\n",color=:slateblue1)
#    printstyled(io, "  status        = ",color=:normal)
#    printstyled(io, "$(z.status)\n",color=:slateblue1)
    printstyled(io, "  run folder    = ",color=:normal)
    rundir=pathof(z)
    printstyled(io, "$(rundir)\n",color=:slateblue1)
    logdir=joinpath(pathof(z),"log")
    printstyled(io, "  log subfolder = ",color=:normal)
    printstyled(io, "$(logdir)\n",color=:slateblue1)
    for i in z.channel.data
        printstyled(io, "  task(s)       = ",color=:normal)
        printstyled(io, "$(i)\n",color=:slateblue1)
    end
end

"""
    put!(x :: AbstractModelConfig,v)

Adds `v` to x.channel (i.e. `put!(x.channel,v)`)

```jldoctest; output = false
using ClimateModels, Suppressor
tmp=ModelConfig()
setup(tmp)
put!(tmp,ClimateModels.RandomWalker)
ClimateModels.pause(tmp)
@suppress ClimateModels.monitor(tmp)
@suppress ClimateModels.help(tmp)
launch(tmp)

isa(tmp,AbstractModelConfig)

# output

true
```
"""
put!(x :: AbstractModelConfig,v) = put!(x.channel,v)
pause(x :: AbstractModelConfig) = put!(x.channel,"pausing now")

"""
    take!(x :: AbstractModelConfig)

Takes command `v` from x.channel (i.e. `take!(x.channel)`) and execute `v(x)`
(if a Function) or return `v` (if not a Function, e.g. a String).

```
tmp=ModelConfig()
put!(tmp,ClimateModels.RandomWalker)
take!(tmp)
```
"""
function take!(x :: AbstractModelConfig)
    tmp=take!(x.channel)
    taskID=string(UUIDs.uuid4())

    msg=("## Task Started \n\n name = ```"*string(tmp)*
                    "```  \n\n ID = "*taskID*" \n\n")
    git_log_msg(x,msg,"task started ["*taskID*"]")

    if isa(tmp,Function)
        tmp(x)
    else
        tmp
    end

    #rewrite current parameters to git log (if changed)
    git_log_prm(x)

    msg=("## Task Ended   \n\n name = ```"*string(tmp)*
                    "```  \n\n ID = "*taskID*" \n\n")
    git_log_msg(x,msg,"task ended   ["*taskID*"]")
end

"""
    git_log_msg(x :: AbstractModelConfig,msg,commit_msg)

Add message `msg` to the `log/README.md` file and git commit.
"""
function git_log_msg(x :: AbstractModelConfig,msg,commit_msg)
    p=joinpath(pathof(x),"log")
    f=joinpath(p,"README.md")
    if isfile(f)
        q=pwd()
        cd(p)
        open(f, "a") do io
            write(io, msg...)
        end
        @suppress run(`$(git()) commit README.md -m "$commit_msg"`)
        cd(q)
    end
end

"""
    git_log_fil(x :: AbstractModelConfig,fil,commit_msg)

Commit changes to file `log/fil` with message `commit_msg`. If `log/fil` is
unknown to git (i.e. commit errors out) then try adding `log/fil` first.
"""
function git_log_fil(x :: AbstractModelConfig,fil,commit_msg)
    p=joinpath(pathof(x),"log")
    f=joinpath(p,fil)
    if isfile(f)
        q=pwd()
        cd(p)
        try
            @suppress run(`$(git()) commit $f -m "$commit_msg"`)
        catch
            try
                run(`$(git()) add $f`)
                @suppress run(`$(git()) commit $f -m "$commit_msg"`)
            catch
                @suppress println("no change to file -> skipping git commit")
            end
        end
        cd(q)
    end
end

"""
    git_log_prm(x :: AbstractModelConfig)

Add files found in `tracked_parameters/` (if any) to git log.
"""
function git_log_prm(x :: AbstractModelConfig)
    p=joinpath(pathof(x),"log")

    if !isempty(x.inputs)
        fil=joinpath(p,"tracked_parameters.toml")
        isfile(fil) ? txt="modify" : txt="initial"
        open(fil, "w") do io
            TOML.print(io, x.inputs)
        end
        q=pwd()
        cd(p)
        try
            run(`$(git()) add tracked_parameters.toml`)
            @suppress run(`$(git()) commit tracked_parameters.toml -m "$(txt) tracked_parameters.toml"`)
        catch e
            #should be skipped when no modification
        end
        cd(q)
    end

    if isdir(joinpath(p,"tracked_parameters"))
        q=pwd()
        cd(p)
        try
            commit_msg="add files in `tracked_parameters/` to git"
            tmp1=readdir("tracked_parameters")
            [run(`$(git()) add tracked_parameters/$i`) for i in tmp1]
            @suppress run(`$(git()) commit -m "$commit_msg"`)
        catch e
            #not sure why this would fail
        end
    cd(q)
    end
end

"""
    git_log_show(x :: AbstractModelConfig)

Show the record of git commits that have taken place in the `log` folder.
"""
function git_log_show(x :: AbstractModelConfig)
    p=joinpath(pathof(x),"log")
    q=pwd()
    cd(p)
    stdout=joinpath(pathof(x),"tmp.txt")
    @suppress run(pipeline(`$(git()) log --decorate --oneline --reverse`,stdout))
    cd(q)
    return readlines(stdout)
end

"""
    log(x :: AbstractModelConfig)

Show the record of git commits that have taken place in the `log` folder.
"""
log(x :: AbstractModelConfig) = git_log_show(x)

"""
    log(x :: AbstractModelConfig, commit_msg :: String;
                 fil="", msg="", init=false, prm=false)

- init=true : create `log` subfolder, initialize git, and commit initial README.md
- prm=true  : add files found in `input` or `tracked_parameters/` (if any) to git log
- !isempty(fil) : commit changes to file `log/fil` with message `commit_msg`.
  If `log/fil` is unknown to git (i.e. commit errors out) then try adding `log/fil` first.

```jldoctest; output = false
using ClimateModels

f=ClimateModels.RandomWalker
i=ClimateModels.OrderedDict(); i["NS"]=100

tmp=ModelConfig(model=f,inputs=i)
setup(tmp)
build(tmp)
launch(tmp)

log(tmp,
    "update tracked_parameters.toml (or skip)",
    fil="tracked_parameters.toml")
ClimateModels.@suppress log(tmp)
isa(tmp,AbstractModelConfig)

# output

true
```
"""
function log(x :: AbstractModelConfig, commit_msg :: String; fil="", msg="", init=false, prm=false)
    init ? git_log_init(x) : nothing
    prm ? git_log_prm(x) : nothing
    !isempty(fil) ? git_log_fil(x :: AbstractModelConfig,fil,commit_msg) : nothing
    !isempty(msg) ? git_log_msg(x :: AbstractModelConfig,msg,commit_msg) : nothing
end

#train(x :: AbstractModelConfig,y) = missing
#compare(x :: AbstractModelConfig,y) = missing
#analyze(x :: AbstractModelConfig,y) = missing


text:
#!/usr/bin/env julia --project=.

InFile1="test-dataset/scenario_1/case.con" # Contingency information
InFile2="test-dataset/scenario_1/case.inl" # Unit Inertia and Governor Response Data File
InFile3="test-dataset/scenario_1/case.raw" # Power Flow Raw Data File
InFile4="test-dataset/scenario_1/case.rop" # Generator Cost Data File
TimeLimitInSeconds=600 # 10 min.
ScoringMethod=2
NetworkModel="IEEE 14"

println("  $(InFile1)") # Print the contents of specified file.
println("  $(InFile2)")
println("  $(InFile3)")
println("  $(InFile4)")
println("  $(TimeLimitInSeconds)")
println("  $(ScoringMethod)")
println("  $(NetworkModel)")

include("MyJulia2.jl") # Insert the contents of MyJulia2.jl here.

compute_solution2(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds, ScoringMethod, NetworkModel, output_dir="test-dataset/scenario_1")
 # Calls the function named "compute_solution2"


text:
### Param Server State ###

const DISPLAY_FILTERS=false

if DISPLAY_FILTERS
	using Images, ImageView
end

@with_kw type ParamServerState
	params
	master_mailbox
	shared_pserver_mailbox
	pserver_mailbox
	accumulated_gradients=SimpleGradient(Any[PABASTO.dummy_weights1,PABASTO.dummy_biases1])
	n_accumulated_gradients=0
	index::Int
	time_var=now()
	param_request_pending::Bool=false
	tau=0.5
	exit::Bool=false
end

function paramserver(master_mailbox, shared_pserver_mailbox,pserver_mailbox, index,starting_params)
	if DISPLAY_FILTERS
		global c = canvasgrid(4,5)
	end

	state = ParamServerState(
	params=starting_params,
	master_mailbox=master_mailbox,
	shared_pserver_mailbox=shared_pserver_mailbox,
	pserver_mailbox=pserver_mailbox,
	index=index)
	println("[PARAM SERVER] initialized")

	while !state.exit
		msg=take!(state.pserver_mailbox)
		if msg==nothing
			msg=take!(state.shared_pserver_mailbox,state.index)
		end
		#handle(state,take!(state.pserver_mailbox))
		handle(state,msg)
		if state.index>=2
			low=max(state.index/4,1)
			high=min(state.index/2,state.index-1)
			@assert high >= low
			time_elapsed = Int(now() - state.time_var)
			if time_elapsed >= state.tau * 1000 && !state.param_request_pending
				println("[PARAM SERVER] Requesting parameter value updates")
				println("[PARAM SERVER] Time last parameter update request: $(time_elapsed)ms");
				put!(state.shared_pserver_mailbox, ParameterUpdateRequestMessage(state.pserver_mailbox),low,high);
				state.param_request_pending = true;
			end
		end
	end

	put!(state.master_mailbox, CeasedOperationMessage(myid()))
	println("[PARAM SERVER] Shutting down")
end

function handle(state:: ParamServerState, message::Void)
	#println("[PARAM SERVER] Spinning")
	#sleep(0.01)
end

function handle(state::ParamServerState, message::ParameterUpdateMessage)
	state.params.data = message.parameters.data
	state.time_var = now()
	state.param_request_pending=false
	println("[PARAM SERVER] Processed parameter value update")
end

function handle(state::ParamServerState, message::ParameterUpdateRequestMessage)
	println("[PARAM SERVER] Reading params")
	put!(message.worker_mailbox,ParameterUpdateMessage(state.params))
end

function updateview(state)
	for i in 1:10
		view(c[i],f(state.params.data[1][i,:]),interactive=false)
	end
end
function f(x)
	a=x-minimum(x)
	a/=maximum(a)
	grayim(transpose(reshape(a,(28,28))))
end
function handle(state::ParamServerState,message::GradientUpdateMessage)
	if state.index==1
		println("[PARAM SERVER] Committing gradients")
		update(state.params, fetch(message.gradient))
	else
		println("[PARAM SERVER] Accumulating gradients")
		state.accumulated_gradients.data+=fetch(message.gradient).data
		state.n_accumulated_gradients+=1
		if state.n_accumulated_gradients > 2
			low=max(state.index/4,1)
			high=min(state.index/2,state.index-1)
			@assert high >= low
			println("[PARAM SERVER] Pushing gradients")
			put!(state.shared_pserver_mailbox,GradientUpdateMessage(state.accumulated_gradients),low,high)
			state.accumulated_gradients.data*=0
			state.n_accumulated_gradients=0
		end
	end
	if DISPLAY_FILTERS
		updateview(state)
	end
end

function handle(state::ParamServerState,message::CeaseOperationMessage)
	state.exit = true
end


text:
__precompile__()

module OSQP

export OSQPMathProgBaseInterface

# Compatibility stuff
using Compat
using Compat.SparseArrays
using Compat.Iterators

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include("../deps/deps.jl")
else
    error("OSQP not properly installed. Please run Pkg.build(\"OSQP\")")
end

function __init__()
    # Get version
    ver_array = split(version(), ".")
    ver_string = string(ver_array[1], ".", ver_array[2], ".", ver_array[3])  # Get string without dev vers
    vnum = VersionNumber(ver_string)

    depsdir = realpath(joinpath(dirname(@__FILE__), "..", "deps"))
    if (vnum.major != 0 && vnum.minor != 2)
        error("Current OSQP version installed is $(osqp_version()), but we require version 0.2.*. Delete the contents of the `$depsdir` directory except for the files `build.jl` and `.gitignore`, then rerun Pkg.build(\"OSQP\").")
    end
end

include("constants.jl")
include("types.jl")
include("interface.jl")
include("mpbinterface.jl")
include("MathOptInterfaceOSQP.jl")

end # module


text:
"""
  lagrangesolve(graph)
  solve graph using lagrange decomposition
"""
function lagrangesolve(graph;
  max_iterations=10,
  update_method=:subgradient, #probingsubgradient
  ϵ=0.001, # ϵ-convergence tolerance
  timelimit=3600,
  α=2, # default subgradient step
  lagrangeheuristic=fixbinaries, # function to calculate the upper bound
  initialmultipliers=:zero, # :relaxation for LP relaxation
  δ = 0.5, # Factor to shrink step when subgradient stuck
  maxnoimprove = 3,
  cpbound=1e6) # Amount of iterations that no improvement is allowed before shrinking step

  ### INITIALIZATION ###
  lgprepare(graph,δ,maxnoimprove,cpbound)
  n = getattribute(graph,:normalized)

  if initialmultipliers == :relaxation
    initialrelaxation(graph)
  end

  starttime = time()
  s = Solution(method=:dual_decomposition)
  λ = getattribute(graph,:λ)[end]
  x = getattribute(graph,:x)[end]
  res = getattribute(graph,:res)[end]
  nmult = getattribute(graph,:numlinks)
  nodes = [node for node in values(getnodes(graph))]
  setattribute(graph,:α, [1.0α])
  iterval = 0

  ### ITERATIONS ###
  for iter in 1:max_iterations
    variant = iter == 1 ? :default : update_method # Use default version in the first iteration

    iterstart = time()
    # Solve subproblems
    Zk = 0
    for node in nodes
       (x,Zkn) = solvenode(node,λ,x,variant)
       Zk += Zkn
    end
    setattribute(graph, :steptaken, true)

    # If no improvement, increase counter
    if Zk < getattribute(graph, :Zk)[end]
      setattribute(graph, :noimprove, getattribute(graph, :noimprove) + 1)
    end
    # If too many iterations without improvement, decrease :α
    if getattribute(graph, :noimprove) >= getattribute(graph, :maxnoimprove)
      setattribute(graph, :noimprove, 0)
      getattribute(graph, :α)[end] *= getattribute(graph, :δ)
    end
    # Save info
    push!(getattribute(graph, :Zk),Zk)
    push!(getattribute(graph, :x),x)
    α = getattribute(graph, :α)[end]
    push!(getattribute(graph, :α), α)

    # Update residuals
    res = x[:,1] - x[:,2]
    push!(getattribute(graph, :res),res)

    # Save iteration data
    itertime = time() - iterstart
    tstamp = time() - starttime
    saveiteration(s,tstamp,[n*iterval,n*Zk,itertime,tstamp],n)
    printiterationsummary(s,singleline=false)

    # Check convergence
    if norm(res) < ϵ
      s.termination = "Optimal"
      return s
    end

    # Check time limit
    if tstamp > timelimit
      s.termination = "Time Limit"
      return s
    end

    # Update multipliers
    println("α = $α")
    (λ, iterval) = updatemultipliers(graph,λ,res,update_method,lagrangeheuristic)
    push!(getattribute(graph, :λ), λ)
    # Update iteration time
    s.itertime[end] = time() - iterstart
  end
  s.termination = "Max Iterations"
  return s
end

# Preprocess function
"""
  lgprepare(graph::PlasmoGraph)

  Prepares the graph to apply lagrange decomposition algorithm
"""
function lgprepare(graph::ModelGraph, δ=0.5, maxnoimprove=3,cpbound=nothing)
  if hasattribute(graph,:preprocessed)
    return true
  end
  n = normalizegraph(graph)
  links = getlinkconstraints(graph)
  nmult = length(links) # Number of multipliers
  setattribute(graph, :numlinks, nmult)
  setattribute(graph, :λ, [zeros(nmult)]) # Array{Float64}(nmult)
  setattribute(graph, :x, [zeros(nmult,2)]) # Linking variables values
  setattribute(graph, :res, [zeros(nmult)]) # Residuals
  setattribute(graph, :Zk, [0.0]) # Bounds
  setattribute(graph, :cuts, [])
  setattribute(graph, :δ, δ)
  setattribute(graph, :noimprove, 0)
  setattribute(graph, :maxnoimprove, maxnoimprove)
  setattribute(graph, :explore, [])
  setattribute(graph, :steptaken, false)

  # Create Lagrange Master
  ms = Model(solver=getsolver(graph))
  @variable(ms, η, upperbound=cpbound)
  @variable(ms, λ[1:nmult])
  @objective(ms, Max, η)

  setattribute(graph, :lgmaster, ms)

  # Each node save its initial objective and set a solver if they don't have one
  for n in values(getnodes(graph))
    mn = getmodel(n)
    if mn.solver == JuMP.UnsetSolver()
      JuMP.setsolver(mn, getsolver(graph))
    end
    mn.ext[:preobj] = mn.obj
    mn.ext[:multmap] = Dict()
    mn.ext[:varmap] = Dict()
  end

  # Maps
  # Multiplier map to know which component of λ to take
  # Varmap knows what values to post where
  for (i,lc) in enumerate(links)
    for j in 1:length(lc.terms.vars)
      var = lc.terms.vars[j]
      var.m.ext[:multmap][i] = (lc.terms.coeffs[j],lc.terms.vars[j])
      var.m.ext[:varmap][var] = (i,j)
    end
  end

  setattribute(graph, :preprocessed, true)
end

# Solve a single subproblem
function solvenode(node,λ,x,variant=:default)
  m = getmodel(node)
  # Restore objective function
  m.obj = m.ext[:preobj]
  m.ext[:lgobj] = m.ext[:preobj]
  # Add dualized part to objective function
  for k in keys(m.ext[:multmap])
    coef = m.ext[:multmap][k][1]
    var = m.ext[:multmap][k][2]
    m.ext[:lgobj] += λ[k]*coef*var
    m.obj += λ[k]*coef*var
    if variant == :ADMM
      j = 3 - m.ext[:varmap][var][2]
      m.obj += 1/2*(coef*var - coef*x[k,j])^2
    end
  end
  # Solve
  solve(m)
  # Pass output
  for v in keys(m.ext[:varmap])
    val = JuMP.getvalue(v)
    x[m.ext[:varmap][v]...] = val
  end
  objval = JuMP.getvalue(m.ext[:lgobj])
  setattribute(node, :objective, objval)
  setattribute(node, :solvetime, getsolvetime(m))

  return x, objval
end

# Multiplier Initialization
function initialrelaxation(graph)
  if !hasattribute(graph,:mflat)
    setattribute(graph, :mflat, create_jump_graph_model(graph))
    getattribute(graph, :mflat).solver = getsolver(graph)
  end
  n = getattribute(graph , :normalized)
  nmult = getattribute(graph , :numlinks)
  mf = getattribute(graph , :mflat)
  solve(mf,relaxation=true)
  getattribute(graph , :λ)[end] = n*mf.linconstrDuals[end-nmult+1:end]
  return getobjectivevalue(mf)
end

function updatemultipliers(graph,λ,res,method,lagrangeheuristic=nothing)
  if method == :subgradient
    subgradient(graph,λ,res,lagrangeheuristic)
  elseif method == :probingsubgradient
    probingsubgradient(graph,λ,res,lagrangeheuristic)
  elseif method == :cuttingplanes
    cuttingplanes(graph,λ,res)
  elseif  method == :interactive
    interactive(graph,λ,res,lagrangeheuristic)
  end
end

# Update functions
function subgradient(graph,λ,res,lagrangeheuristic)
  α = getattribute(graph , :α)[end]
  n = getattribute(graph , :normalized)
  bound = n*lagrangeheuristic(graph)
  Zk = getattribute(graph , :Zk)[end]
  step = α*abs(Zk-bound)/(norm(res)^2)
  λ += step*res
  return λ,bound
end

function αeval(αv,graph,bound)
  xv = deepcopy(getattribute(graph , :x)[end])
  res = getattribute(graph , :res)[end]
  Zk = getattribute(graph , :Zk)[end]
  λ = getattribute(graph , :λ)[end]
  nodes = [node for node in values(getnodes(graph))]
  step = abs(Zk-bound)/(norm(res)^2)
  zk = 0
  for node in nodes
     (xv,Zkn) = solvenode(node,λ+αv*step*res,xv,:default)
     zk += Zkn
  end
  return zk
end

function αexplore(graph,bound)
  df = getattribute(graph , :explore)
  n = getattribute(graph , :normalized)
  z = Float64[]
  for α in 0:0.1:2
    push!(z,n*αeval(α,graph,bound))
  end
  push!(df,z)
end

function probingsubgradient(graph,λ,res,lagrangeheuristic,α=getattribute(graph , :α)[end],Δ=0.01;exhaustive=false)
  res = getattribute(graph , :res)[end]
  Zk = getattribute(graph , :Zk)[end]
  n = getattribute(graph , :normalized)
  bound = n*lagrangeheuristic(graph)
  step = abs(Zk-bound)/(norm(res)^2)
  # First point
  α1 = 0
  z1 = Zk

  # Second point
  α2 = α
  z2 = αeval(α2,graph,bound)

  if (z2 - z1)/abs(z1) >= Δ
    # If second point gives an increase of more than Δ%, take it.
    return (λ += α2*step*res), bound
  elseif abs((z2 - z1)/z1) < Δ
    # If second point is similar to the first one, take the midpoint
    return (λ += α2/2*step*res), bound
  else
    # If second point is larger than first, take quarter-step
    if exhaustive
      return probingsubgradient(graph,λ,res,lagrangeheuristic,α/4,Δ;exhaustive=true)
    else
      return (λ += α2/4*step*res), bound
    end
  end
end

@require Plots begin
function interactive(graph,λ,res,lagrangeheuristic)
  α = getattribute(graph , :α)[end]
  n = getattribute(graph , :normalized)
  bound = n*lagrangeheuristic(graph)
  Zk = getattribute(graph , :Zk)[end]
  αexplore(graph,bound)
  plot(0:0.1:2,getattribute(graph , :explore)[end])
  print("α = ")
  α = parse(Float64,readline(STDIN))
  step = α*abs(Zk-bound)/(norm(res)^2)
  λ += step*res
  return λ,bound
end
end

function cuttingplanes(graph,λ,res)
  ms = getattribute(graph , :lgmaster)
  Zk = getattribute(graph , :Zk)[end]
  nmult = getattribute(graph , :numlinks)

  λvar = getindex(ms, :λ)
  η = getindex(ms,:η)

  cut = @constraint(ms, η <= Zk + sum(λvar[j]*res[j] for j in 1:nmult))
  push!(getattribute(graph , :cuts), cut)

  solve(ms)
  return getvalue(λvar), getobjectivevalue(ms)
end

# Standard Lagrangean Heuristics
function fixbinaries(graph::ModelGraph,cat=[:Bin])
  if !haskey(graph.attributes,:mflat)
    setattribute(graph , :mflat, create_jump_graph_model(graph))
    getattribute(graph, :mflat).solver = getsolver(graph)
  end
  n = getattribute(graph , :normalized)
  mflat = getattribute(graph , :mflat)
  mflat.colVal = vcat([getmodel(n).colVal for n in values(getnodes(graph))]...)
  for j in 1:mflat.numCols
    if mflat.colCat[j] in cat
      mflat.colUpper[j] = mflat.colVal[j]
      mflat.colLower[j] = mflat.colVal[j]
    end
  end
  status = solve(mflat)
  if status == :Optimal
    return n*getobjectivevalue(mflat)
  else
    error("Heuristic model not infeasible or unbounded")
  end
end

function fixintegers(graph::ModelGraph)
  fixbinaries(graph,[:Bin,:Int])
end


text: -"foo"

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
################################################################################
#
# (q/p)adic integers
#
# complete enough to support hnf
################################################################################
# CHECK precision!!!

mutable struct QadicRing{S, T} <: Generic.Ring
  Q::S #The corresponding local field
  basis::Vector{T} #The OK-basis of the ring, where OK is
                   #the maximal order of the base field of Q
  function QadicRing{S, T}(x::S) where {S <: Union{LocalField, FlintQadicField, FlintPadicField}, T}
    z = new{S, T}()
    z.Q = x
    return z
  end

end

function Base.show(io::IO, Q::QadicRing)
  println("Integers of ", Q.Q)
end

function MaximalOrder(Q::FlintQadicField)
  return QadicRing{FlintQadicField, qadic}(Q)
end

function MaximalOrder(Q::FlintPadicField)
  return QadicRing{FlintPadicField, padic}(Q)
end
#integers(Q::FlintQadicField) = ring_of_integers(Q)
function MaximalOrder(Q::LocalField{S, T}) where {S, T <: Union{EisensteinLocalField, UnramifiedLocalField}}
  return QadicRing{LocalField{S, T}, LocalFieldElem{S, T}}(Q)
end
#integers(Q::FlintPadicField) = ring_of_integers(Q)

mutable struct QadicRingElem{S, T} <: RingElem
  P::QadicRing{S, T}
  x::T
  function QadicRingElem(P::QadicRing{S, T}, a::T) where {S, T}
    r = new{S, T}(P, a)
  end
end

function Base.show(io::IO, a::QadicRingElem)
  print(io, a.x)
end

*(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.P, a.x*b.x)
+(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.P, a.x+b.x)
-(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.P, a.x-b.x)
-(a::QadicRingElem) = QadicRingElem(a.P, -a.x)
^(a::QadicRingElem, b::QadicRingElem) = QadicRingElem(a.P, a.x^b.x)
^(a::T, b::QadicRingElem{S, T}) where {S, T} = a^b.x

function inv(a::QadicRingElem)
  valuation(a.x) == 0 || error("The element is not invertible!")
  return QadicRingElem(a.P, inv(a.x))
end

==(a::QadicRingElem, b::QadicRingElem) = a.x == b.x

function divexact(a::QadicRingElem, b::QadicRingElem)
  @assert !iszero(b.x)
  iszero(a) && return a
  valuation(a.x) >= valuation(b.x) || error("division not exact")
  return QadicRingElem(a.P, a.x//b.x)
end

function Base.divrem(a::QadicRingElem, b::QadicRingElem)
  if valuation(a.x) < valuation(b.x)
    return setprecision(a.P(0), precision(a)), a
  end
  q = divexact(a, b)
  return q, a-q*b
end

function Base.div(a::QadicRingElem, b::QadicRingElem)
  if valuation(a.x) < valuation(b.x)
    return setprecision(a.P(0), precision(a))
  end
  q = divexact(a, b)
  return q
end

parent(a::QadicRingElem) = a.P
elem_type(::Type{QadicRing{S, T}}) where {S, T} = QadicRingElem{S, T}
parent_type(::Type{QadicRingElem{S, T}}) where {S, T} = QadicRing{S, T}

zero(Q::QadicRing) = QadicRingElem(Q, Q.Q(0))
one(Q::QadicRing) = QadicRingElem(Q, Q.Q(1))

(Q::QadicRing{S, T})(a::T) where {S, T} = QadicRingElem(Q, a)
(Q::QadicRing)(a::QadicRingElem) = QadicRingElem(a.P, a.x)
(Q::QadicRing)(a::Integer) = QadicRingElem(Q, Q.Q(a))
(Q::QadicRing)(a::fmpz) = QadicRingElem(Q, Q.Q(a))

function (Q::QadicRing)(a::fmpq)
  p = prime(Q.Q)
  if iszero(mod(denominator(a), p))
    error("The element is not in the ring!")
  end
  return QadicRingElem(Q, Q.Q(a))
end

(Q::QadicRing)() = QadicRingElem(Q, Q.Q())

(Q::LocalField{S, T})(a::QadicRingElem{S, T}) where {S, T} = a.x
(Q::FlintPadicField)(a::QadicRingElem{FlintPadicField, padic}) = a.x
(Q::FlintQadicField)(a::QadicRingElem{FlintQadicField, qadic}) = a.x

valuation(a::QadicRingElem) = valuation(a.x)
isunit(a::QadicRingElem) = !iszero(a) && valuation(a) == 0
(Q::FlintQadicField)(a::padic) =  _map(Q, a) #TODO: do properly

function _map(Q::FlintQadicField, a::padic)
  K = parent(a)
  v = valuation(a)
  if v >= 0
    q = Q(lift(a))
    return setprecision(q, a.N)
  else
    d = uniformizer(K)^-v
    n = a*d
    n1 = divexact(Q(lift(n)), Q(lift(d)))
    return n1
  end
end

function Base.deepcopy_internal(a::QadicRingElem, dict::IdDict)
  return QadicRingElem(a.P, a.x)
end

function canonical_unit(a::QadicRingElem)
  iszero(a.x) && return setprecision(a.P(1), precision(a))
  v = valuation(a.x)
  return QadicRingElem(a.P, inv(a.x//prime(a.P.Q)^v))
end

function gcdx(a::QadicRingElem, b::QadicRingElem)
  if iszero(a)
    c = canonical_unit(b)
    return b*c, a, c
  end
  if iszero(b)
    c = canonical_unit(a)
    return a*c, c, b
  end
  if valuation(a.x) < valuation(b.x)
    c = canonical_unit(a)
    return a*c, c, setprecision(a.P(0), precision(a))
  else
    c = canonical_unit(b)
    return b*c, setprecision(b.P(0), precision(b)), c
  end
end

function mul_red!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem, f::Bool = false)
  return b*c
end

function mul!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem)
  return b*c
end

function add!(a::QadicRingElem, b::QadicRingElem, c::QadicRingElem)
  return b+c
end

function addeq!(a::QadicRingElem, b::QadicRingElem)
  return a+b
end

Base.iszero(a::QadicRingElem) = iszero(a.x)
Base.isone(a::QadicRingElem) = isone(a.x)

Base.precision(Q::QadicRing) = precision(Q.Q)
Base.precision(a::QadicRingElem) = precision(a.x)

function setprecision!(Q::QadicRing, n::Int)
  setprecision!(Q.Q, n)
end

function Base.setprecision(a::QadicRingElem, n::Int)
  return a.P(setprecision(a.x, n))
end

function setprecision!(a::QadicRingElem, n::Int)
  a.x = setprecision!(a.x, n)
end

function Base.setprecision(a::Generic.MatSpaceElem{QadicRingElem{qadic}}, n::Int)
  return map_entries(x -> setprecision(x, n), a)
end

coefficient_ring(Q::QadicRing) = integers(coefficient_ring(Q.Q))


text:
"""
    documentation(::Type{QE}, searchstring::AbstractString)

Searches through the description of all the known flags that can be set in QuantumEspresso,
returns the flags where the description contains the `searchstring`.

    documentation(::Type{Elk}, flag::Symbol)

Returns the documentation for a given flag.
"""
function documentation(::Type{QE}, searchstring::AbstractString)
    found = Pair{String, Vector{QEFlagInfo}}[]
    for inputinfo in QEInputInfos
        foundflags = QEFlagInfo[]
        for fi in allflags(inputinfo)
            if occursin(searchstring, fi.description)
                push!(foundflags, fi)
            end
        end
        if !isempty(foundflags)
            push!(found, inputinfo.exec => foundflags)
        end
    end
    return found
end

function documentation(::Type{QE}, flagsymbol::Symbol)
    flag = qe_flaginfo(flagsymbol)
    if flag == nothing
        error("No documentation found for flag $flagsymbol, are you sure it is a valid flag for any QE executables?")
    end
    return flag
end

"""
    documentation(::Type{Elk}, searchstring::AbstractString)

Searches through the description of all the known flags that can be set in Elk,
returns the flags where the description contains the `searchstring`.

    documentation(::Type{Elk}, flag::Symbol)

Returns the documentation for a given flag.
"""
function documentation(::Type{Elk}, searchstring::AbstractString)
    found = Pair{ElkControlBlockInfo, Vector{ElkFlagInfo}}[]
    for b in ELK_CONTROLBLOCKS
        found_flags = ElkFlagInfo[]
        for f in b.flags
            if occursin(searchstring, f.description)
                push!(found_flags, f)
            end
        end
        if !isempty(found_flags) || occursin(searchstring, b.description)
            push!(found, b => found_flags)
        end
    end
    return found
end

function documentation(::Type{Elk}, flagsymbol::Symbol)
    flag = elk_flaginfo(flagsymbol)
    if flag == nothing
        block = getfirst(x->x.name == flagsymbol, ELK_CONTROLBLOCKS)
        if block == nothing
            error("No documentation found for flag $flagsymbol, are you sure it is a valid flag for Elk?")
        else
            return block
        end
    else
        return flag
    end
end


text:
using NURBS

include("bezier.jl")
include("knot.jl")
include("basis.jl")
include("bspline.jl")
include("rbspline.jl")
include("beziersurfaces.jl")
include("bsurfaces.jl")
include("rbsurfaces.jl")
include("raise.jl")

include("bezsurftest.jl")


text:
#This has very limited methods created to use the PDFIO functionality for
#testing and quick prototyping. While you can use them to get some of the
#regular testing these are not currently part of the core library.

using PDFIO
using PDFIO.Common
using PDFIO.Cos
using PDFIO.PD

function pdfhlp_extract_doc_content_to_dir(filename,dir=tempdir())
  file=rsplit(filename, '/',limit=2)
  filenm=file[end]
  dirpath=joinpath(dir,filenm)
  if isdir(dirpath)
    rm(dirpath; force=true, recursive=true)
  end
  mkdir(dirpath)
  doc=pdDocOpen(filename)
  try
    npage= pdDocGetPageCount(doc)
    for i=1:npage
      page = pdDocGetPage(doc, i)
      if pdPageIsEmpty(page)==false
        contents=pdPageGetContents(page)
        bufstm = get(contents)
        buf = read(bufstm)
        close(bufstm)
        path=joinpath(dirpath,string(i)*".txt")
        write(path, buf)
      end
    end
  finally
    pdDocClose(doc)
  end
end

function pdfhlp_extract_doc_embedded_files(filename,dir=tempdir())
  file=rsplit(filename, '/',limit=2)
  filenm=file[end]
  dirpath=joinpath(dir,filenm)
  if isdir(dirpath)
    rm(dirpath; force=true, recursive=true)
  end
  mkdir(dirpath)
  doc=pdDocOpen(filename)
  try
    catalog = pdDocGetCatalog(doc)
    names = get(catalog, CosName("Names"))
    cosDoc = pdDocGetCosDoc(doc)
    nmdict = cosDocGetObject(cosDoc, names)
    println(nmdict)
    if nmdict !== CosNull
      efref  = get(nmdict, CosName("EmbeddedFiles"))
      efroot = cosDocGetObject(cosDoc, efref)
      #simple case no tree just a few files attached in the root node.
      #A proper implementation needs full names tree traversal.
      efarr = get(efroot, CosName("Names"))
      data = get(efarr)
      len=length(data)
      println(len)
      for i=1:len:2
        key=data[i]
        println(key)
        val=data[i+1]
        filespec=cosDocGetObject(cosDoc, val)
        ef=get(filespec, CosName("EF"))
        filename=get(filespec,CosName("F")) #UF could be there as well.
        stmref=get(ef, CosName("F"))
        stm=cosDocGetObject(cosDoc,stmref)
        bufstm=decode(stm)
        buf=read(bufstm)
        close(bufstm)
        path=joinpath(dirpath,get(filename))
        write(path,buf)
      end
    end
  finally
    pdDocClose(doc)
  end
end

function pdfhlp_extract_doc_attachment_files(filename,dir=tempdir())
  file=rsplit(filename, '/',limit=2)
  filenm=file[end]
  dirpath=joinpath(dir,filenm)
  if isdir(dirpath)
    rm(dirpath; force=true, recursive=true)
  end
  mkdir(dirpath)
  doc=pdDocOpen(filename)
  cosDoc=pdDocGetCosDoc(doc)
  try
    npage= pdDocGetPageCount(doc)
    for i=1:npage
      page = pdDocGetPage(doc, i)
      cospage = pdPageGetCosObject(page)
      annots=cosDocGetObject(cosDoc, get(cospage, CosName("Annots")))
      if (annots === CosNull)
        continue
      end
      annotsarr=get(annots)
      for annot in annotsarr
        annotdict = cosDocGetObject(cosDoc, annot)
        subtype = get(annotdict,CosName("Subtype"))
        if (subtype == CosName("FileAttachment"))
          filespec=cosDocGetObject(cosDoc, get(annotdict,CosName("FS")))
          ef=get(filespec, CosName("EF"))
          filename=get(filespec,CosName("F")) #UF could be there as well.
          stmref=get(ef, CosName("F"))
          stm=cosDocGetObject(cosDoc,stmref)
          bufstm=decode(stm)
          buf=read(bufstm)
          close(bufstm)
          path=joinpath(dirpath,get(filename))
          write(path,buf)
        end
      end
    end
  finally
    pdDocClose(doc)
  end
end


text:
# --- calculates moment's tensor
"""

  momel(X::Matrix{Float}, ind::Tuple)

Returns Float, an element of moment's tensor at ind multiindex

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> momel(M, (1,1,1,1))
0.4318298020613279

```
"""
@inline momel(X::Matrix{T}, multind::Tuple) where T<: AbstractFloat = blockel(X, multind, multind, 0)

"""
  naivemoment(data::Matrix{Float}, m::Int)

Returns Array{Float, m} the m'th moment tensor

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> naivemoment(M, 3)
2×2×2 Array{Float64,3}:
[:, :, 1] =
  -0.523092   0.142552
  0.142552  -0.0407653

[:, :, 2] =
  0.142552   -0.0407653
  -0.0407653   0.0120729

```
"""
function naivemoment(X::Matrix{T}, m::Int = 4) where T<: AbstractFloat
  n = size(X, 2)
  moment = zeros(T, fill(n, m)...,)
  for i = 1:(n^m)
    ind = Tuple(CartesianIndices((fill(n, m)...,))[i])
    @inbounds moment[ind...] = momel(X, ind)
  end
  moment
end

# --- uses the naive method to calculate cumulants 2 - 6
"""

  mixel(X::Matrix{T}, ind::Tuple)

Returns Float, mixed element for cumulants 4-6 at ind multi-index

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> mixel(M, (1,1,1,1))
-1.232956812564408

mixel(M, (1,1,1,1,1,1))
1.015431116914347
```
"""
function mixel(X::Matrix{T}, i::Tuple) where T<: AbstractFloat
  a = zero(T)
  if length(i) == 4
    a -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[4]))
    a -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[4])) + momel(X, (i[1],i[4]))*momel(X, (i[2],i[3]))
  elseif length(i) == 5
    a -= momel(X, (i[1],i[2],i[3]))*momel(X, (i[4],i[5])) + momel(X, (i[1],i[2],i[4]))*momel(X, (i[3],i[5]))
    a -= momel(X, (i[1],i[2],i[5]))*momel(X, (i[3],i[4])) + momel(X, (i[2],i[3],i[4]))*momel(X, (i[1],i[5]))
    a -= momel(X, (i[2],i[3],i[5]))*momel(X, (i[1],i[4])) + momel(X, (i[1],i[3],i[4]))*momel(X, (i[2],i[5]))
    a -= momel(X, (i[1],i[3],i[5]))*momel(X, (i[2],i[4])) + momel(X, (i[3],i[4],i[5]))*momel(X, (i[1],i[2]))
    a -= momel(X, (i[2],i[4],i[5]))*momel(X, (i[1],i[3])) + momel(X, (i[1],i[4],i[5]))*momel(X, (i[2],i[3]))
  elseif length(i) == 6
    a1 = -momel(X, (i[1],i[2],i[3]))*momel(X, (i[4],i[5], i[6])) - momel(X, (i[1],i[2],i[4]))*momel(X, (i[3],i[5], i[6]))
    a1 -= momel(X, (i[1],i[2],i[5]))*momel(X, (i[3],i[4], i[6])) + momel(X, (i[1],i[2],i[6]))*momel(X, (i[3],i[4], i[5]))
    a1 -= momel(X, (i[1],i[3],i[4]))*momel(X, (i[2],i[5], i[6])) + momel(X, (i[1],i[3],i[5]))*momel(X, (i[2],i[4], i[6]))
    a1 -= momel(X, (i[1],i[3],i[6]))*momel(X, (i[2],i[4], i[5])) + momel(X, (i[1],i[4],i[5]))*momel(X, (i[2],i[3], i[6]))
    a1 -= momel(X, (i[1],i[4],i[6]))*momel(X, (i[2],i[3], i[5])) + momel(X, (i[1],i[5],i[6]))*momel(X, (i[2],i[3], i[4]))
    a2 = -momel(X, (i[1],i[2],i[3],i[4]))*momel(X, (i[5], i[6])) - momel(X, (i[1],i[2],i[3],i[5]))*momel(X, (i[4],i[6]))
    a2 -= momel(X, (i[1],i[2],i[3],i[6]))*momel(X, (i[4], i[5])) + momel(X, (i[1],i[2],i[4],i[5]))*momel(X, (i[3], i[6]))
    a2 -= momel(X, (i[1],i[2],i[4],i[6]))*momel(X, (i[3], i[5])) + momel(X, (i[1],i[2],i[5],i[6]))*momel(X, (i[3], i[4]))
    a2 -= momel(X, (i[3],i[4],i[5],i[6]))*momel(X, (i[1], i[2])) + momel(X, (i[1],i[3],i[4],i[5]))*momel(X, (i[2], i[6]))
    a2 -= momel(X, (i[1],i[3],i[4],i[6]))*momel(X, (i[2], i[5])) + momel(X, (i[1],i[3],i[5],i[6]))*momel(X, (i[2], i[4]))
    a2 -= momel(X, (i[2],i[4],i[5],i[6]))*momel(X, (i[1], i[3])) + momel(X, (i[1],i[4],i[5],i[6]))*momel(X, (i[2], i[3]))
    a2 -= momel(X, (i[2],i[3],i[5],i[6]))*momel(X, (i[1], i[4])) + momel(X, (i[2],i[3],i[4],i[6]))*momel(X, (i[1], i[5]))
    a2 -= momel(X, (i[2],i[3],i[4],i[5]))*momel(X, (i[1], i[6]))
    a3 = -momel(X, (i[1],i[2]))*momel(X, (i[3],i[4]))*momel(X, (i[5], i[6]))
    a3 -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[5]))*momel(X, (i[4], i[6]))
    a3 -= momel(X, (i[1],i[2]))*momel(X, (i[3],i[6]))*momel(X, (i[4], i[5]))
    a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[4]))*momel(X, (i[5], i[6]))
    a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[5]))*momel(X, (i[4], i[6]))
    a3 -= momel(X, (i[1],i[3]))*momel(X, (i[2],i[6]))*momel(X, (i[4], i[5]))
    a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[3]))*momel(X, (i[5], i[6]))
    a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[3]))*momel(X, (i[4], i[6]))
    a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[3]))*momel(X, (i[4], i[5]))
    a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[5]))*momel(X, (i[3], i[6]))
    a3 -= momel(X, (i[1],i[4]))*momel(X, (i[2],i[6]))*momel(X, (i[3], i[5]))
    a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[4]))*momel(X, (i[3], i[6]))
    a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[4]))*momel(X, (i[3], i[5]))
    a3 -= momel(X, (i[1],i[5]))*momel(X, (i[2],i[6]))*momel(X, (i[3], i[4]))
    a3 -= momel(X, (i[1],i[6]))*momel(X, (i[2],i[5]))*momel(X, (i[3], i[4]))
    a += a1+a2-2*a3
  end
  a
end

"""
  naivecumulant(data::Matrix, m::Int)

Returns Array{Float, m} the m'th cumulant tensor

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> naivecumulant(M, 3)
2×2×2 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

```
"""
function naivecumulant(X::Matrix{T}, m::Int = 4) where T<: AbstractFloat
  m < 7 || throw(AssertionError("naive implementation of $m cumulant not supported"))
  if m == 1
    return naivemoment(X,m)
  end
  X = X .- mean(X, dims=1)
  ret = naivemoment(X,m)
  if m in [4,5,6]
    n = size(X, 2)
    for i = 1:(n^m)
      ind = Tuple(CartesianIndices((fill(n, m)...,))[i])
      @inbounds ret[ind...] += mixel(X, ind)
    end
  end
    return ret
end


text:
module JuliaLSP

using LanguageServer, LanguageServer.SymbolServer

function julia_main()::Cint
    server = LanguageServer.LanguageServerInstance(stdin, stdout)
    server.runlinter = true
    run(server)
    return 0 # if things finished successfully
end

end # module


text:
using RDatasets, KernelDensity
mlmf = dataset("mlmRev", "Gcsemv");

df = mlmf[complete_cases(mlmf[[:Written, :Course]]), :];

dc = array(df[:Course]);
dcc = zeros(0);
for i in 1:length(dc)
  if !isnan(dc[i])
    push!(dcc,dc[i])
  end
end
kdc = kde(dcc);

dw = array(df[:Written]);
dww = zeros(0);
for i in 1:length(dw)
  if !isnan(dw[i])
    push!(dww,dw[i])
  end
end
kdw = kde(dww);

using Winston
kdc = kde(dcc); kdw = kde(dww);
plot(kdc.x, kdc.density, "r--", kdw.x, kdw.density, "b;")

max_stars_repo_path: Module 3/Chapter05/kdensity.jl

max_issues_repo_path: Module 3/Chapter05/kdensity.jl

max_forks_repo_path: Module 3/Chapter05/kdensity.jl


text:
"""
Reference:
1. Structure Discovery in Nonparametric Regression through Compositional Kernel Search
		David Duvenaud, James Robert Lloyd, Roger Grosse, Joshua B. Tenenbaum, Zoubin Ghahramani (2013)

2. https://www.cs.toronto.edu/~duvenaud/cookbook/

3. GaussianProcesses.jl: A Nonparametric Bayes package for the Julia Language
"""

"""
RBF Kernel
"""
struct ArdGaussKernel{T, VT<:AbstractVector{T}}
	l::VT
	σ::VT
end
@functor ArdGaussKernel
function (GK::ArdGaussKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	length(GK.l) == size(x, 1) || throw(DimensionMismatch("size of length scale parameter should be the same as number of features"))
	n = size(x, 2)
	scaled_x = x ./ GK.l
	d = square_metric(scaled_x)
	σ_square = GK.σ[1]*GK.σ[1]
	(@. σ_square*exp(-0.5*d)) + Diagonal(λ*ones(n))
end
function (GK::ArdGaussKernel{T, VT})(x, xo) where {T, VT}
	σ_square = GK.σ[1]*GK.σ[1]
	scaled_x = x ./ GK.l
	scaled_xo = xo ./ GK.l
	d = square_metric(scaled_x, scaled_xo)
	@. σ_square*exp(-0.5*d)
end

function ArdGaussKernel(n_features::Int)
	l = rand(n_features)
	σ = rand(1)
	GaussKernel(l, σ)
end

reset(GK::ArdGaussKernel) = ArdGaussKernel(length(GK.l))

## derived kernel
struct IsoGaussKernel{T, VT<:AbstractVector{T}}
	l::VT
	σ::VT
end
@functor IsoGaussKernel
function (IsoGK::IsoGaussKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x ./ IsoGK.l[1]
	d = square_metric(scaled_x)
	σ_square = IsoGK.σ[1]*IsoGK.σ[1]
	(@. σ_square*exp(-0.5*d)) + Diagonal(λ*ones(n))
end
function (IsoGK::IsoGaussKernel{T, VT})(x, xo) where {T, VT}
	σ_square = IsoGK.σ[1]*IsoGK.σ[1]
	scaled_x = x ./ IsoGK.l[1]
	scaled_xo = xo ./ IsoGK.l[1]
	d = square_metric(scaled_x, scaled_xo)
	@. σ_square*exp(-0.5*d)
end

function IsoGaussKernel()
	l = rand(1)
	σ = rand(1)
	IsoGaussKernel(l, σ)
end

reset(IsoGK::IsoGaussKernel) = IsoGaussKernel()

"""
Rational quardratic kernel
"""
struct ArdRQKernel{T, VT<:AbstractVector{T}}
	l::VT
	α::VT
	σ::VT
end
@functor ArdRQKernel
function (RQK::ArdRQKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	length(RQK.l) == size(x, 1) || throw(DimensionMismatch("size of length scale parameter should be the same as number of features"))
	n = size(x, 2)
	scaled_x = x ./ RQK.l
	d = square_metric(scaled_x)
	σ_square = RQK.σ[1]*RQK.σ[1]
	(@. σ_square*exp((-RQK.α[1])*log(1.0 + 0.5*d/RQK.α[1]))) + Diagonal(λ*ones(n))
end
function (RQK::ArdRQKernel{T, VT})(x, xo) where {T, VT}
	σ_square = RQK.σ[1]*RQK.σ[1]
	scaled_x = x ./ RQK.l
	scaled_xo = xo ./ RQK.l
	d = square_metric(scaled_x, scaled_xo)
	@. σ_square*exp((-RQK.α[1])*log(1.0 + 0.5*d/RQK.α[1]))
end

function ArdRQKernel(n_features::Int)
	l = rand(n_features)
	α = rand(1)
	σ = rand(1)
	RQKernel(l, α, σ)
end

reset(RQK::ArdRQKernel) = ArdRQKernel(length(RQK.l))

## derived kernel
struct IsoRQKernel{T, VT<:AbstractVector{T}}
	l::VT
	α::VT
	σ::VT
end
@functor IsoRQKernel
function (IsoRQK::IsoRQKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x ./ IsoRQK.l[1]
	d = square_metric(scaled_x)
	σ_square = IsoRQK.σ[1]*IsoRQK.σ[1]
	(@. σ_square*exp((-IsoRQK.α[1])*log(1.0 + 0.5*d/IsoRQK.α[1]))) + Diagonal(λ*ones(n))
end
function (IsoRQK::IsoRQKernel{T, VT})(x, xo) where {T, VT}
	σ_square = IsoRQK.σ[1]*IsoRQK.σ[1]
	scaled_x = x ./ IsoRQK.l[1]
	scaled_xo = xo ./ IsoRQK.l[1]
	d = square_metric(scaled_x, scaled_xo)
	@. σ_square*exp((-IsoRQK.α[1])*log(1.0 + 0.5*d/IsoRQK.α[1]))
end

function IsoRQKernel()
	l = rand(1)
	α = rand(1)
	σ = rand(1)
	IsoRQKernel(l, α, σ)
end

reset(IsoRQK::IsoRQKernel) = IsoRQKernel()

"""
Iso Periodic kernel
"""
struct IsoPeriodKernel{T, VT<:AbstractVector{T}}
	p::VT
	l::VT
	σ::VT
end
@functor IsoPeriodKernel
function (IsoPK::IsoPeriodKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	n = size(x, 2)
	scaled_x = x ./ IsoPK.p[1]
	d = norm2_metric(scaled_x)
	σ_square = IsoPK.σ[1]*IsoPK.σ[1]
	l_square = IsoPK.l[1]*IsoPK.l[1]
	(@. σ_square*exp(-2.0*sinpi(d)*sinpi(d)/l_square)) + Diagonal(λ*ones(n))
end
function (IsoPK::IsoPeriodKernel{T, VT})(x, xo) where {T, VT}
	scaled_x = x ./ IsoPK.p[1]
	scaled_xo = xo ./ IsoPK.p[1]
	d = norm2_metric(scaled_x, scaled_xo)
	σ_square = IsoPK.σ[1]*IsoPK.σ[1]
	l_square = IsoPK.l[1]*IsoPK.l[1]
	@. σ_square*exp(-2.0*sinpi(d)*sinpi(d)/l_square)
end

function IsoPeriodKernel()
	p = 1.0 ./ rand(1)
	l = rand(1)
	σ = rand(1)
	IsoPeriodKernel(p, l, σ)
end

reset(IsoPK::IsoPeriodKernel) = IsoPeriodKernel()

"""
Linear Kernel
"""
struct ArdLinearKernel{T, VT<:AbstractVector{T}}
	σb::VT
	σv::VT
	c::VT
end
@functor ArdLinearKernel
function (LK::ArdLinearKernel{T, VT})(x; λ=T(1e-6)) where {T, VT}
	length(LK.c) == size(x, 1) && length(Lk.σv) == size(x, 1) || throw(DimensionMismatch("size of offset parameter should be the same as number of features"))
	n = size(x, 2)
	shifted_x = x .- LK.c
	scaled_shifted_x = shifted_x ./ LK.σv
	d = inner_prod_metric(scaled_shifted_x)
	σb_square = LK.σb[1]*LK.σb[1]
	(@. d+σb_square) + Diagonal(λ*ones(n))
end
function (LK::ArdLinearKernel{T, VT})(x, xo) where {T, VT}
	shifted_x = x .- LK.c
	scaled_shifted_x = shifted_x ./ LK.σv
	shifted_xo = xo .- LK.c
	scaled_shifted_xo = shifted_xo ./ LK.σv
	d = inner_prod_metric(scaled_shifted_x, scaled_shifted_xo)
	σb_square = LK.σb[1]*LK.σb[1]
	d+σb_square
end

function ArdLinearKernel(n_features::Int)
	σb = rand(1)
	σv = rand(n_features)
	c  = rand(n_features)
	LinearKernel(σb, σv, c)
end

reset(LK::ArdLinearKernel) = ArdLinearKernel(length(LK.c))

"""
White noise kernel
"""
struct WhiteNoiseKernel{T, VT<:AbstractVector{T}}
	σ::VT
end
@functor WhiteNoiseKernel
function (WNK::WhiteNoiseKernel{T, VT})(x) where {T, VT}
	n = size(x, 2)
	σ_square = WNK.σ[1]*WNK.σ[1]
	Diagonal(σ_square*ones(n))
end

function WhiteNoiseKernel()
	σ = rand(1)
	WhiteNoiseKernel(σ)
end

reset(WNK::WhiteNoiseKernel) = WhiteNoiseKernel()

"""
Constant Kernel (suffers PSD problem !!!)
"""
struct ConstantKernel{T, VT<:AbstractVector{T}}
	σ::VT
end
@functor ConstantKernel
function (CK::ConstantKernel{T, VT})(x) where {T, VT}
	n = size(x, 2)
	CK.σ[1]*ones(n, n)
end

function ConstantKernel()
	σ = rand(1)
	ConstantKernel(σ)
end

reset(CK::ConstantKernel) = ConstantKernel()

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
"""
Typed wrapper for `Vector{String}`.
"""
struct Path{T}
    steps::Vector{String}

    # Hooray for inner constructors!
    Path{T}() where {T} = new{T}(String[])
    function Path{T}(strings::Vector{<:AbstractString}) where {T}
        return new{T}(String.(strings))
    end
end

function Path{T}(string::AbstractString) where {T}
    return isempty(string) ? Path{T}() : Path{T}(split(string, "."))
end
Path{T}(strings::AbstractString...) where {T} = Path{T}(collect(strings))

# only equal if types T are equal and steps are equal
Base.:(==)(::Path, ::Path) = false
Base.:(==)(a::Path{T}, b::Path{T}) where {T} = a.steps == b.steps

Base.hash(c::Path{T}, u::UInt64) where {T} = hash(c.steps, hash(T, u))

Base.length(p::Path) = length(p.steps)
Base.first(p::Path) = first(p.steps)
Base.last(p::Path) = last(p.steps)

# Methods for growing a path.
path_promote(::Type{T}, ::Type{T}) where {T} = T
path_demote(::Type{T}) where {T} = T

function catpath(a::Path{T}, b::Path{U}) where {T,U}
    return Path{path_promote(T, U)}(vcat(a.steps, b.steps))
end

catpath(a::String, b::Path{T}) where {T} = Path{T}(vcat(a, b.steps))

function striplast(a::Path{T}) where {T}
    length(a) == 0 && throw(BoundsError(a, 0))
    return Path{path_demote(T)}(a.steps[1:(end - 1)])
end

function stripfirst(a::Path{T}) where {T}
    length(a) == 0 && throw(BoundsError(a, 0))
    return Path{T}(a.steps[2:end])
end

function splitpath(a::Path{T}, i::Integer) where {T}
    0 <= i <= length(a) || throw(BoundsError())

    pre = Path{path_demote(T)}(a.steps[1:(end - i)])
    post = Path{T}(a.steps[(end + 1 - i):end])

    return pre, post
end

Base.show(io::IO, p::Path) = print(io, string(p))
Base.string(c::Path{T}) where {T} = "Path{$(string(T))} $(join(c.steps, "."))"


text:
using InteractionNetworkModels, StructuredDistances, BenchmarkTools
using Plots
mode = [[1,2,1],[1,1,1]]
mode = [[1,2,1],[1]]
γ = 2.9
d_S = EditDistance(LCS())
d_E = MatchingDistance(LCS())
d_E = AfpMatchingDistance(FastLCS(100), 3.0)
V = 1:10

model = SIM(mode, γ, d_E, V, 10, 20)
β = 0.7
mcmc_move = InvMcmcMixtureMove(
    (
        EditAllocationMove(ν=1),
        InsertDeleteMove(ν=1, len_dist=TrGeometric(0.8, 1, model.K_inner.u))
    ),
    (β, 1-β)
)

mcmc = InvMcmcSampler(
    mcmc_move,
    burn_in=2000, lag=75
    )

x = mcmc(model, desired_samples=10000, lag=1, burn_in=0)
plot(x)
acceptance_prob(mcmc)
summaryplot(x)

x.sample[findall(length.(x.sample) .> 5)]

@btime draw_sample(mcmc, model, desired_samples=1000, burn_in=0, lag=1)

acceptance_prob(mcmc)

mcmc_old = SimMcmcInsertDelete(
    ν_ed=1, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.8, 1, model.K_inner.u),
    burn_in=2000, lag=75
)

@btime draw_sample(mcmc_old, model, desired_samples=1000, burn_in=0, lag=1)

x_old = mcmc_old(model, desired_samples=4000)
plot!(x_old)
summaryplot!(x_old)

x.sample
x_old.sample

x_old.sample

x.sample


text:
using DungBase
using Test


text:
n = 3
m = 2
p = 5

@testset "Testing $(MLK.dotvectors!)" begin
    for T in FloatingPointTypes
        Set_X = [rand(T, p) for i = 1:n]
        Set_Y = [rand(T,p)  for i = 1:m]

        X = permutedims(hcat(Set_X...))
        Y = permutedims(hcat(Set_Y...))

        @test isapprox(MLK.dotvectors(Val(:row),    X), vec(sum((X.*X), dims = 2)))
        @test isapprox(MLK.dotvectors(Val(:col), X), vec(sum((X.*X), dims = 1)))

        @test isapprox(MLK.dotvectors!(Val(:row),    Vector{T}(undef, n), X), vec(sum((X.*X), dims = 2)))
        @test isapprox(MLK.dotvectors!(Val(:col), Vector{T}(undef, p), X), vec(sum((X.*X), dims = 1)))

        @test_throws DimensionMismatch MLK.dotvectors!(Val(:row),    Vector{T}(undef, 2), Array{T}(undef, 3,2))
        @test_throws DimensionMismatch MLK.dotvectors!(Val(:row),    Vector{T}(undef, 4), Array{T}(undef, 3,4))
        @test_throws DimensionMismatch MLK.dotvectors!(Val(:col), Vector{T}(undef, 2), Array{T}(undef, 2,3))
    end
end

@testset "Testing $(MLK.gramian!)" begin
    for T in FloatingPointTypes
        Set_X = [rand(T, p) for i = 1:n]
        Set_Y = [rand(T,p)  for i = 1:m]

        X = permutedims(hcat(Set_X...))
        Y = permutedims(hcat(Set_Y...))

        P = [LinearAlgebra.dot(x,y) for x in Set_X, y in Set_X]

        @test isapprox(MLK.gramian!(Val(:row),    Array{T}(undef, n,n), X,  true), P)
        @test isapprox(MLK.gramian!(Val(:col), Array{T}(undef, n,n), permutedims(X), true), P)

        P = [LinearAlgebra.dot(x,y) for x in Set_X, y in Set_Y]
        @test isapprox(MLK.gramian!(Val(:row),    Array{T}(undef, n,m), X,  Y),  P)
        @test isapprox(MLK.gramian!(Val(:col), Array{T}(undef, n,m), permutedims(X), permutedims(Y)), P)

        @test_throws DimensionMismatch MLK.gramian!(Val(:row), Array{T}(undef, n+1,n),   X, true)
        @test_throws DimensionMismatch MLK.gramian!(Val(:row), Array{T}(undef, n,n+1),   X, true)
        @test_throws DimensionMismatch MLK.gramian!(Val(:row), Array{T}(undef, n+1,n+1), X, true)

        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, n+1,n),   permutedims(X), true)
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, n,n+1),   permutedims(X), true)
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, n+1,n+1), permutedims(X), true)

        @test_throws DimensionMismatch MLK.gramian!(Val(:row),    Array{T}(undef, n+1,m), X,  Y)
        @test_throws DimensionMismatch MLK.gramian!(Val(:row),    Array{T}(undef, n,m+1), X,  Y)
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, n+1,m), permutedims(X), permutedims(Y))
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, n,m+1), permutedims(X), permutedims(Y))

        @test_throws DimensionMismatch MLK.gramian!(Val(:row),    Array{T}(undef, m+1,n), Y,  X)
        @test_throws DimensionMismatch MLK.gramian!(Val(:row),    Array{T}(undef, m,n+1), Y,  X)
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, m+1,n), permutedims(Y), permutedims(X))
        @test_throws DimensionMismatch MLK.gramian!(Val(:col), Array{T}(undef, m,n+1), permutedims(Y), permutedims(X))
    end
end


text:
function estimate_return(::UnweightedIS, π, τ::BanditExperience)
    return exp(logpdf(π, τ.action) - τ.blogp) * τ.reward
end

function estimate_return(::WeightedIS, π, τ::BanditExperience)
    ρ = exp(logpdf(π, τ.action) - τ.blogp)
    return ρ * τ.reward, ρ
end

function estimate_return_withentropy(::UnweightedIS, π, τ::BanditExperience, α)
    logp = logpdf(π, τ.action)
    return exp(logp - τ.blogp) * (τ.reward - α * logp)
end

function estimate_return_withentropy(::WeightedIS, π, τ::BanditExperience, α)
    logp = logpdf(π, τ.action)
    ρ = exp(logp - τ.blogp)
    return ρ * (τ.reward - α * logp), ρ
end