text:
"""
CPLEX basis, proposed by Bixby in DOI: 10.1287/ijoc.4.3.267.

Original form:

minimize c x
subject to
-Inf <= A1 x <= b1
  b2 <= A2 x <= Inf
  bl <= A3 x <= bu
  b4 <= A4 x <= b4
  xl <= x <= xu

Canonical form:

minimize c x + 0 s1 + 0 s2 + 0 s3
subjec to
  A1 x + s1 == b1
  A2 x - s2 == b2
  A3 x - s3 == bl
  Ab x      == b4
  xl <= x <= xu
  0 <= s1 <= Inf
  0 <= s2 <= Inf
  0 <= s3 <= bu - bl,

where we treat the constraint bl <= A3 x <= bu as the form A x >= b.
"""

module Cplex

using MatrixOptInterface
using SparseArrays
import ..PhaseOne
import ..Artificial
import Simplex
import Simplex: nrows, ncols

const MatOI = MatrixOptInterface
const MOI = MatOI.MOI

function run(lp::MatOI.LPSolverForm, Tv::Type; kwargs...)::Simplex.SpxData

    if !haskey(kwargs, :original_lp)
        @error "Argument :original_lp is not provided."
    end
    if !haskey(kwargs, :pivot_rule)
        @error "Argument :pivot_rule is not provided."
    end
    original_lp = kwargs[:original_lp]
    pivot_rule = kwargs[:pivot_rule]

    # get cplex basis
    B = cplex_basis(original_lp)

    # add artificial variables that are identified in basis
    cpxlp, cpxartif = Artificial.reformulate(lp, basis = B)

    # compute solution x (maybe infeasible)
    x, is_feasible = compute_x(cpxlp, B)

    # number of canonical variables and artificial variables
    num_orgvars = ncols(lp)
    num_artif = ncols(cpxlp) - num_orgvars
    # @show num_orgvars, num_artif

    # If x is infeasible, we need to solve a piecewise LP.
    p1artif = Int[]
    if !is_feasible
        # @warn("CPLEX initial basis solution is infeasible.")
        # convert to phase-one form
        cpxlp, B, p1artif = reformulate(cpxlp, x, num_orgvars, num_artif, B)
        @assert size(cpxlp.A,1) == length(B)
    end

    # load the problem
    spx = Simplex.SpxData(cpxlp, Tv)
    spx.pivot_rule = deepcopy(pivot_rule)

    # set basis
    Simplex.set_basis(spx, B)

    # Run simplex method
    Simplex.run_core(spx)

    # restore the original column bounds
    # for j in 1:ncols(cpxlp)
    #     spx.lpdata.v_lb[j] = cpxlp.v_lb[j]
    #     spx.lpdata.v_ub[j] = cpxlp.v_ub[j]
    # end

    for j in [cpxartif; p1artif]
        spx.lpdata.v_lb[j] = 0.0
        spx.lpdata.v_ub[j] = 0.0
    end

    return spx
end

function cplex_basis(lp::MatOI.LPForm)::Array{Int64}

    m1_idx = Int[]
    m2_idx = Int[]
    m3_idx = Int[]

    for i = 1:nrows(lp)
        if lp.c_lb[i] == lp.c_ub[i]
            push!(m3_idx, i)
        elseif lp.c_lb[i] > -Inf
            push!(m2_idx, i)
        elseif lp.c_ub[i] < Inf
            push!(m1_idx, i)
        else
            @error "Constraint type is not supported"
        end
    end
    m1 = length(m1_idx) # number of <= constraints
    m2 = length(m2_idx) # number of >= constraints
    m3 = length(m3_idx) # number of == constraints
    # It is important to have the constraints in the following order: <=, >=, and ==.
    # So we use a permuatation vector.
    perm_A = [m1_idx; m2_idx; m3_idx]
    # @show m1, m2, m3

    C2 = Int[] # free
    C3 = Int[] # single bounded
    C4 = Int[] # bounded

    # penalty
    q = Array{Float64}(undef, ncols(lp))

    for j = 1:ncols(lp)
        if lp.v_lb[j] > -Inf
            if lp.v_ub[j] < Inf
                push!(C4, j)
                q[j] = lp.v_lb[j] - lp.v_ub[j]
            else
                push!(C3, j)
                q[j] = lp.v_lb[j]
            end
        elseif lp.v_ub[j] < Inf
            push!(C3, j)
            q[j] = -lp.v_ub[j]
        else
            push!(C2, j)
            q[j] = 0.0
        end
    end

    # @show lp.c
    gamma = maximum(abs.(lp.c))
    cmax = gamma > 0 ? 1000 * gamma : 1.0
    q .+= lp.c / cmax

    # sort C2, C3, C4 by q
    sort!(C2, by = x -> q[x])
    sort!(C3, by = x -> q[x])
    sort!(C4, by = x -> q[x])

    # ordered column indices
    C = [C2; C3; C4]

    # Step 1
    I = [ones(Int64, m1 + m2); zeros(Int64, nrows(lp) - m1 - m2)]
    r = [ones(Int64, m1 + m2); zeros(Int64, nrows(lp) - m1 - m2)]
    v = fill(Inf, nrows(lp))
    B = collect(1:(m1+m2)) .+ ncols(lp) # {n+1, ..., n+m1+m2}

    # Step 2
    for k = 1:ncols(lp)
        # Step 2a
        alpha = 0.
        max_l = -1
        for i = 1:nrows(lp)
            if r[i] == 0
                absA = abs(lp.A[perm_A[i],C[k]])
                if alpha < absA
                    alpha = absA
                    max_l = i
                end
            end
        end
        if alpha >= 0.99
            push!(B, C[k])
            # println("Step 2a: Adding $(C[k]) to basis with alpha ($alpha) >= 0.99")
            I[max_l] = 1
            v[max_l] = alpha
            for i = 1:nrows(lp)
                if abs(lp.A[perm_A[i],C[k]]) > 0
                    r[i] += 1
                end
            end
            continue
        end
        # Step 2b
        alpha = 0.
        max_l = -1
        next_k = false
        for i = 1:nrows(lp)
            if abs(lp.A[perm_A[i],C[k]]) > 0.01 * v[i]
                next_k = true
                break
            end
            if I[i] == 0
                absA = abs(lp.A[perm_A[i],C[k]])
                if alpha < absA
                    alpha = absA
                    max_l = i
                end
            end
        end
        if next_k
            continue
        end
        if alpha > 0.0
            push!(B, C[k])
            # println("Step 2b: Adding $(C[k]) to basis with alpha ($alpha) > 0.0")
            I[max_l] = 1
            v[max_l] = alpha
            for i = 1:nrows(lp)
                if abs(lp.A[perm_A[i],C[k]]) > 0
                    r[i] += 1
                end
            end
        end
    end

    # Step 3
    for i = (m1+m2+1):nrows(lp)
        if I[i] == 0
            # Add an artificial variable to cover each remaining uncovered row
            push!(B, ncols(lp) + m1 + m2 + perm_A[i]) # B \cup {n+m1+m2+i}
            # println("Step 3: Adding $(lp.ncols + m1 + m2 + perm_A[i]) to basis")
        end
    end
    # @show length(B), nrows(lp), B
    @assert length(B) == nrows(lp)

    return sort(B)
end

"""
    Input:
    minimize a
    subject to
        A x + d a == b
        xl <= x <= xu,
    where
        d_i = +1 if b_i - A_i xN > 0,
        d_i = -1 otherwise.

    Output:
    minimize a + sl + su
    subject to
        A x + d a == b
        x + sl - sl2 == xl (nxl)
        x - su + su2 == xu (nxu)
        a, sl, sl2, su, su2 >= 0,

...
# Arguments
- `lp::MatOI.LPSolverForm{T,AT,VT}`: LP problem object in the canonical form
- `x::VT`: the solution vector in the dimension of the number of columns in lp
- `num_orgvars::Int`: the number of variables in the original problem (not in lp)
- `num_artif::Int`: the number of artificial variables introduced to lp
- `basis::Vector{Int}`: basis vector identified for lp
...
"""
function reformulate(lp::MatOI.LPSolverForm{T,AT,VT}, x::VT, num_orgvars::Int, num_artif::Int, basis::Vector{Int}) where {T, AT, VT}
    @assert length(basis) == nrows(lp)
    @assert length(x) == ncols(lp)

    # compute on cpu
    c = Array(lp.c); xl = Array(lp.v_lb); xu = Array(lp.v_ub); b = Array(lp.b);
    # x = Array(lp.x)

    # count slacks for column bounds
    xl_idx = Int[]
    xu_idx = Int[]
    sizehint!(xl_idx, num_orgvars)
    sizehint!(xu_idx, num_orgvars)
    for j in 1:num_orgvars
        if x[j] < xl[j]
            push!(xl_idx, j)
        end
        if x[j] > xu[j]
            push!(xu_idx, j)
        end
    end
    nxl = length(xl_idx)
    nxu = length(xu_idx)

    # Piecewise function is replaced by a set of linear inequalities.
    num_auxvars = nxl + nxu # number of variables to replace the piecewise functions
    num_slacks = nxl + nxu # slacks for the linear inequalities
    lp_nrows = nrows(lp) + num_auxvars # number of rows
    lp_ncols = ncols(lp) + num_auxvars + num_slacks # number of columns

    # row bounds
    b = VT([b; xl[xl_idx]; xu[xu_idx]])

    # TODO: Do we need to compute x?

    # solution x
    x = Array([x; zeros(num_auxvars + num_slacks)])

    # basis construction
    for (i,j) in enumerate(xl_idx)
        if x[j] < xl[j]
            push!(basis, ncols(lp) + i)
            x[ncols(lp) + i] = xl[j] - x[j]
        else
            push!(basis, ncols(lp) + num_auxvars + i)
            x[ncols(lp) + num_auxvars + i] = x[j] - xl[j]
        end
    end
    for (i,j) in enumerate(xu_idx)
        if x[j] > xu[j]
            push!(basis, ncols(lp) + nxl + i)
            x[ncols(lp) + nxl + i] = x[j] - xu[j]
        else
            push!(basis, ncols(lp) + num_auxvars + nxl + i)
            x[ncols(lp) + num_auxvars + nxl + i] = xu[j] - x[j]
        end
    end

    # objective coefficient
    c = VT([c; ones(num_auxvars); zeros(num_slacks)])

    # column bounds
    xl = VT([xl; zeros(num_auxvars + num_slacks)])
    xu = VT([xu; fill(Inf, num_auxvars + num_slacks)])
    xl[xl_idx] = fill(-Inf,nxl)
    xu[xu_idx] = fill(+Inf,nxu)

    artif_col_idx = Int[]

    # create the submatrix for auxiliary and slack variables
    I = Array{Int64}(undef, num_auxvars + num_slacks)
    J = collect(1:(num_auxvars + num_slacks))
    V = Array{Float64}(undef, num_auxvars + num_slacks)
    pos = 1
    # sl
    for i = 1:nxl
        I[pos] = nrows(lp) + i
        V[pos] = 1
        push!(artif_col_idx, ncols(lp) + J[pos])
        pos += 1
    end
    # su
    for i = 1:nxu
        I[pos] = nrows(lp) + nxl + i
        V[pos] = -1
        push!(artif_col_idx, ncols(lp) + J[pos])
        pos += 1
    end
    # sl2
    for i = 1:nxl
        I[pos] = nrows(lp) + i
        V[pos] = -1
        pos += 1
    end
    # su2
    for i = 1:nxu
        I[pos] = nrows(lp) + nxl + i
        V[pos] = 1
        pos += 1
    end
    S = sparse(I, J, V, lp_nrows, num_auxvars + num_slacks)

    # create the matrix in canonical form
    A = AT([[sparse(Matrix(lp.A));
        sparse(collect(1:nxl), xl_idx, ones(nxl), nxl, ncols(lp));
        sparse(collect(1:nxu), xu_idx, ones(nxu), nxu, ncols(lp))] S])

    senses = fill(MatOI.EQUAL_TO, lp_nrows)
    p1lp = MatOI.LPSolverForm{T,typeof(A),typeof(c)}(lp.direction, c, A, b, senses, xl, xu)

    return p1lp, sort(basis), artif_col_idx
end

"""
    compute solution x based on basis information
    and return whether x is feasible or not
"""
function compute_x(lp::MatOI.LPSolverForm{T,AT,VT}, basis::Vector{Int}) where {T, AT, VT}
    @assert nrows(lp) == length(basis)

    nonbasis = Int[]
    sizehint!(nonbasis, ncols(lp) - length(basis))
    x = VT(undef, ncols(lp))

    # construct nonbasic indices
    for j in 1:ncols(lp)
        if !in(j, basis)
            push!(nonbasis, j)

            # set x values
            axl = abs(lp.v_lb[j])
            axu = abs(lp.v_ub[j])
            if lp.v_lb[j] > -Inf && axl <= axu
                x[j] = lp.v_lb[j]
            elseif lp.v_ub[j] < Inf && axl > axu
                x[j] = lp.v_ub[j]
            else
                x[j] = 0.0
            end
        end
    end

    x[basis] .= lp.A[:,basis] \ (lp.b .- lp.A[:,nonbasis] * x[nonbasis])

    is_feasible = true
    for j = 1:ncols(lp)
        if x[j] < lp.v_lb[j] || x[j] > lp.v_ub[j]
            is_feasible = false
            # @show j, x[j], lp.v_lb[j], lp.v_ub[j]
            break
        end
    end

    return x, is_feasible
end

end


text:
using Metropolis
using Documenter

makedocs(;
    modules=[Metropolis],
    authors="Edwin Bedolla",
    repo="https://github.com/edwinb-ai/Metropolis.jl/blob/{commit}{path}#L{line}",
    sitename="Metropolis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://edwinb-ai.github.io/Metropolis.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/edwinb-ai/Metropolis.jl")

ext: jl


text:
export iauRxp
"""
Multiply a p-vector by an r-matrix.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
    r        double[3][3]    r-matrix
    p        double[3]       p-vector

Returned:
    rp       double[3]       r * p

Note:
    It is permissible for p and rp to be the same array.

Called:
    iauCp        copy p-vector

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRxp(double r[3][3], double p[3], double rp[3])
function iauRxp(r::AbstractMatrix{<:Real}, p::AbstractVector{<:Real})
    rp = zeros(Float64, 3)
    ccall((:iauRxp, libsofa_c), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Array{Float64, 2}, r'),
          convert(Array{Float64, 1}, p),
          rp)

    return SVector{3}(rp)
 end


text:
using Oceananigans.Operators

function solve_for_pressure!(pressure, solver, arch, grid, U, G, Δt)
    if solver.type isa Channel && arch isa GPU
        ϕ = RHS = solver.storage.storage1
    else
        ϕ = RHS = solver.storage
    end

    @launch(device(arch), config=launch_config(grid, :xyz),
            calculate_pressure_right_hand_side!(RHS, solver.type, arch, grid, U, G, Δt))

    solve_poisson_equation!(solver, grid)

    @launch(device(arch), config=launch_config(grid, :xyz),
            copy_pressure!(pressure, ϕ, solver.type, arch, grid))

    return nothing
end

"""
Calculate the right-hand-side of the Poisson equation for the non-hydrostatic
pressure and in the process apply the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along any direction we need to perform a GPU fast cosine transform algorithm.
"""
function calculate_pressure_right_hand_side!(RHS, solver_type, arch, grid, U, G, Δt)
    @loop_xyz i j k grid begin
        i′, j′, k′ = permute_index(solver_type, arch, i, j, k, grid.Nx, grid.Ny, grid.Nz)

        @inbounds RHS[i′, j′, k′] = divᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w) / Δt +
                                    divᶜᶜᶜ(i, j, k, grid, G.u, G.v, G.w)
    end
    return nothing
end

"""
Copy the non-hydrostatic pressure into `p` from `solver.storage` and
undo the permutation

    [a, b, c, d, e, f, g, h] -> [a, c, e, g, h, f, d, b]

along any dimensions in which a GPU fast cosine transform was performed.
"""
function copy_pressure!(p, ϕ, solver_type, arch, grid)
    @loop_xyz i j k grid begin
        i′, j′, k′ = unpermute_index(solver_type, arch, i, j, k, grid.Nx, grid.Ny, grid.Nz)
        @inbounds p[i′, j′, k′] = real(ϕ[i, j, k])
    end
    return nothing
end


text:
using AlgebraicDynamics.UWDDynam
using Catlab.WiringDiagrams
using Test

const UWD = UndirectedWiringDiagram

@testset "UWDDynam" begin
  dx(x, p, t) = [x[1]^2, 2*x[1]-x[2]]
  dy(y, p, t) = [1 - y[1]^2]

  r = ContinuousResourceSharer{Float64}(2, dx)
  @test eltype(r) == Float64
  #identity
  d = UWD(2)
  add_box!(d, 2)
  add_junctions!(d, 2)
  set_junction!(d, [1,2])
  set_junction!(d, [1,2], outer=true)

  r2 = oapply(d, [r])
  @test nstates(r) == nstates(r2)
  @test nports(r) == nports(r2)
  @test portmap(r) == portmap(r2)
  x0 = [10.0, 7.5]
  @test eval_dynamics(r, x0) == eval_dynamics(r2, x0)
  @test exposed_states(r, x0) == exposed_states(r2, x0)

  h = 0.1
  drs = oapply(d, [euler_approx(r, h)])
  drs2 = euler_approx(r2, h)
  @test eval_dynamics(drs, x0) == eval_dynamics(drs2, x0)
  drs3 = euler_approx(r2)
  drs4 = oapply(d, [euler_approx(r)])
  @test eval_dynamics(drs, x0) == eval_dynamics(drs3, x0, [h])
  @test eval_dynamics(drs, x0) == eval_dynamics(drs4, x0, [h])

  # merge
  d = UWD(1)
  add_box!(d, 2)
  add_junctions!(d, 1)
  set_junction!(d, [1,1])
  set_junction!(d, [1], outer=true)

  r2 = oapply(d, [r])
  @test nstates(r2) == 1
  @test nports(r2) == 1
  @test portmap(r2) == [1]
  @test eval_dynamics(r2, [5.0]) == [30.0]
  @test exposed_states(r2, [5.0]) == [5.0]

  r = ContinuousResourceSharer{Float64}(2, (u,p,t) -> [u[1]*p[1], u[2]*p[2] + t])
  r2 = oapply(d, [r])
  @test eval_dynamics(r2, [5.0], [1.0, 2.0], 10.0) == [25.0]

  # copy
  r = ContinuousResourceSharer{Float64}(1, 2, dx, [2])
  d = UWD(2)
  add_box!(d, 1)
  add_junctions!(d, 1)
  set_junction!(d, [1])
  set_junction!(d, [1,1], outer = true)
  r2 = oapply(d, r)
  @test nstates(r2) == 2
  @test nports(r2) == 2
  @test portmap(r2) == [2,2]
  @test eval_dynamics(r2, x0) == eval_dynamics(r, x0)
  @test exposed_states(r2, x0) == [x0[2], x0[2]]

  # copy states and merge back together
  r = ContinuousResourceSharer{Float64}(2, 2, dx, [1,1])
  d = UWD(1)
  add_box!(d, 2)
  add_junctions!(d, 1)
  set_junction!(d, [1,1])
  set_junction!(d, [1], outer = true)
  r2 = oapply(d, r)
  @test nstates(r2) == 2
  @test nports(r2) == 1
  @test portmap(r2) == [1]
  @test eval_dynamics(r2, x0) == eval_dynamics(r, x0)
  @test exposed_states(r2, x0) == [x0[1]]

  # copy states and merge with otherwise
  r = ContinuousResourceSharer{Float64}(1, dy)
  rcopy = ContinuousResourceSharer{Float64}(2, 1, dy, [1,1])
  d = HypergraphDiagram{Nothing, Symbol}(2)
  add_box!(d, 1, name = :r); add_box!(d, 2, name = :copy); add_box!(d, 1, name = :r)
  add_junctions!(d, 2)
  set_junction!(d, [1,1,2,2])
  set_junction!(d, [1,2], outer = true)
  xs = Dict(:r => r, :copy => rcopy)
  r2 = oapply(d, xs)
  @test nstates(r2) == 1
  @test nports(r2) == 2
  @test portmap(r2) == [1,1]
  @test eval_dynamics(r2, [7.0]) == [3 *(1 - 7.0^2)]
  @test exposed_states(r2, [7.0]) == [7.0, 7.0]

  # add a state
  d = UWD(2)
  add_box!(d, 1)
  add_junctions!(d, 2)
  set_junction!(d, [1])
  set_junction!(d, [1,2], outer = true)
  r2 = oapply(d, [r])
  @test nstates(r2) == 2
  @test nports(r2) == 2
  @test portmap(r2) == [1,2]
  @test eval_dynamics(r2, [7.0, 11.0]) == [-48.0, 0.0]
  @test exposed_states(r2, [7.0, 11.0]) == [7.0, 11.0]

  # lots of boxes
  r = ContinuousResourceSharer{Float64}(2, dx)
  s = ContinuousResourceSharer{Float64}(1, dy)
  xs = Dict(:r => r, :s => s)
  d = HypergraphDiagram{Nothing, Symbol}(5)
  add_box!(d, 2, name = :r); add_box!(d, 1, name = :s); add_box!(d, 2, name = :r)
  add_junctions!(d, 4)
  set_junction!(d, [1,1,1,4,2])
  set_junction!(d, [1,1,2,3,3], outer = true)
  r2 = oapply(d, xs)
  @test nstates(r2) == 4
  @test nports(r2) == 5
  @test portmap(r2) == [1,1, 3, 4,4]
  x0 = [2.0, 7.0, 3.0, 5.0]
  @test eval_dynamics(r2, x0) == [3.0, 49.0, 11.0, 0.0]
  @test exposed_states(r2, x0) == [2.0, 2.0, 3.0, 5.0, 5.0]

  h = 0.1
  dr = oapply(d, euler_approx([r,s,r], h))
  dr2 = euler_approx(r2, h)
  dr3 = oapply(d, euler_approx([r,s,r]))
  dr4 = euler_approx(r2)
  @test eval_dynamics(dr, x0) == [2.3, 11.9, 4.1, 5.0]
  @test eval_dynamics(dr, x0) == eval_dynamics(dr2, x0)
  @test eval_dynamics(dr, x0) == eval_dynamics(dr3, x0, [h], 0)
  @test eval_dynamics(dr, x0) == eval_dynamics(dr4, x0, [h], 0)

  h = 0.25
  dr = oapply(d, euler_approx(xs, h))
  dr2 = euler_approx(r2, h)
  dr3 = oapply(d, euler_approx(xs))
  dr4 = euler_approx(r2)
  @test eval_dynamics(dr, x0) == eval_dynamics(dr2, x0)
  @test eval_dynamics(dr, x0) == eval_dynamics(dr3, x0, [h], 0)
  @test eval_dynamics(dr, x0) == eval_dynamics(dr4, x0, [h], 0)

  # substitute and oapply commute
  d = UWD(4)
  add_box!(d, 1); add_box!(d, 2)
  add_junctions!(d, 3)
  set_junction!(d, [1,2,1])
  set_junction!(d, [1,2,3,3], outer = true)

  din = UWD(2)
  add_box!(din, 1); add_box!(din, 2)
  add_junctions!(din, 2)
  set_junction!(din, [1,1,2])
  set_junction!(din, [1,2], outer = true)

  dtot = ocompose(d, 2, din)

  s = ContinuousResourceSharer{Float64}(1, 2, dx, [2])

  r1 = oapply(dtot, [s,s,r])
  r2 = oapply(d, [s, oapply(din, [s,r])])

  @test nstates(r1) == nstates(r2)
  @test nports(r1) == nports(r2)
  @test portmap(r1) == portmap(r2)
  x0 = [2.0, 3.0, 5.0, 7.0, 11.0]
  @test eval_dynamics(r1, x0) == eval_dynamics(r2, x0)

end # test set


text:
import Arpack: eigs
using LinearAlgebra

function pcca(T::Matrix, n::Integer; optimize = false, solver=BaseSolver())
    israte = isratematrix(T)
    pi     = stationarydensity(T, israte)
    X, λ   = schurvectors(T, pi, n, israte, solver)
	assertstructure(X, pi)
    chi    = makeprobabilistic(X, optimize)
end

function isratematrix(T::Matrix)
    s = sum(T[1,:])
    isapprox(s, 0) && return true
    isapprox(s, 1) && return false
    error("given matrix is neither a rate nor a probability matrix")
end

function stationarydensity(T, israte=isratematrix(T))
    which = israte ? :SM : :LM
    pi = eigs(T', nev=1, which=which)[2] |> vec
	#pi = eigvecs(T, sortby=real)[:,end]
    @assert isreal(pi)
    pi = abs.(pi)
    pi = pi / sum(pi)
end

function makeprobabilistic(X::Matrix, optimize::Bool)
    n = size(X, 2)
    A = innersimplexalgorithm(X)
    if n > 2 && optimize
        A = opt(A, X)
    end
    chi = X*A
end

function crispassignments(chi)
     assignments = mapslices(argmax, chi, dims=2) |> vec
end

function assertstructure(X, pi)
	@assert X' * Diagonal(pi) * X ≈ I
	@assert X[:,1] ≈ ones(length(pi))
end

# compute initial guess based on indexmap
innersimplexalgorithm(X) = feasiblize!(inv(X[indexmap(X), :]), X)

function indexmap(X)
    # get indices of rows of X to span the largest simplex
    rnorm(x) = sqrt.(sum(abs2.(x), dims=2)) |> vec
    ind=zeros(Int, size(X,2))
    for j in 1:length(ind)
        rownorm=rnorm(X)
        # store largest row index
        ind[j]=argmax(rownorm)
        if j == 1
            # translate to origin
            X = X .- X[ind[1],:]'
        else
            # remove subspace
            X=X/rownorm[ind[j]]
            vt=X[ind[j],:]'
            X=X-X*vt'*vt
        end
    end
    return ind
end

function feasiblize!(A,X)
    A[:,1] = -sum(A[:,2:end], dims=2)
    A[1,:] = -minimum(X[:,2:end] * A[2:end,:], dims=1)
    A / sum(A[1,:])
end

# crispness criterion, cf. Roeblitz (2013)
# only applies if X is normalized (eq. 8)
# TODO: check does this apply?
function roeblitzcrit(A)
    n = size(A,1)
    trace = 0
    for i = 1:n, j=1:n
        trace += A[i,j]^2 / A[1,j]
    end
    return trace
end

using Optim
function opt(A0, X)
    A = copy(A0)
    Av = view(A, 2:size(A,1), 2:size(A,2)) # view on the variable part

    function obj(a)
        Av[:] = a
        -roeblitzcrit(feasiblize!(A, X))
    end

    result = optimize(obj, Av[:], NelderMead())
    Av[:] = result.minimizer
    return feasiblize!(A, X)
end

function isreversible(P, pi=stationarydensity(P))
    db = [pi[i] * P[i,j] - pi[j] * P[j,i] for i =1:5, j=1:5]
    all(isapprox.(db, 0, atol=1e-8))
end


text:
@testset "1002.find-common-characters.jl" begin
    @test common_chars(["bella","label","roller"]) == ['l', 'l', 'e']
    @test common_chars(["cool","lock","cook"]) == ['c', 'o']
end


text:
function install_projs_args()

    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "--deep", "-d"
            help = "define how many extra version of each pkg to install, " *
                    "from newer to older."
            arg_type = Int
            default = 0
        "--root", "-r"
            help = "The root folder for searching projects."
            arg_type = String
            default = "."
        "--dry-run", "-y"
            help = "Run without consequences."
            action = :store_true
    end

    parsed_args = ArgParse.parse_args(ARGS, argset)
    deep = parsed_args["deep"]
    dry_run = parsed_args["dry-run"]
    root = abspath(parsed_args["root"])

    install_projs(root, deep, dry_run)
end


text:
function river_length(River)
    x, z, dx = River.x, River.z, River.Δx^2
    L = zero(x)
    @tturbo for i in 2:length(x)
        L[i] = √(dx  + (z[i]-z[i-1])^2)
    end
    cumsum(L)
end

function lazycumsum(A::Vector{T}, idx) where T
    s = zero(T)
    @turbo for i in 1:idx
        s += A[i]
    end
    s
end


text:
module ImageUtils

using ImageMagick
using Reexport
using ColorVectorSpace
using AxisArrays
const axes = Base.axes
using Random
using Statistics
using LinearAlgebra
using Interpolations
using Unitful
using FFTW

@reexport using Colors
@reexport using FixedPointNumbers
@reexport using Images

include("Arrays.jl")
include("Coloring.jl")
include("Slicing.jl")
include("Transformation.jl")
include("PerfusionMaps.jl")
include("FileExport.jl")
include("Permutation.jl")
include("AnalyzeNifti.jl")

end # module


text:
module Switching

using MAT
using DataFrames
using Dates
using OrderedCollections
using ShiftedArrays
#using StaticArrays

include("raw_data_paths.jl")
include("utilities.jl")
include("process_pokes.jl")
include("process_streaks.jl")
include("process_rawdata.jl")
include("drug_calendar.jl")

export matread, DataFrame, by
export get_rawdata, get_rawfile, read_rawfile
export process_pokes, process_session, process_all
export get_alldata, read_alldata_format, process_alldata

end # module


text:
using DataStructures

"""
	type SimpleNetwork <: Network
		time::Int64
		layers::DataStructures.OrderedDict{Symbol, Layer}
		plasticconnections::Array{Connection, 1}
		layeriterator::Base.ValueIterator{DataStructures.OrderedDict{Symbol,Simulator.Layer}}
"""

type SimpleNetwork <: Network
	time::Int64
	layers::DataStructures.OrderedDict{Symbol, Layer}
	plasticconnections::Array{Connection, 1}
	# currently iterating a OrderedDict causes memory allocation, see
	# https://github.com/JuliaLang/DataStructures.jl/issues/93
	# shouldn't be a big deal for this application
	layeriterator::Base.ValueIterator{DataStructures.OrderedDict{Symbol,Layer}}
end

"""
	SimpleNetwork()

Create an empty network
"""
function SimpleNetwork()
	layers = DataStructures.OrderedDict{Symbol, Layer}()
	it = values(layers)
	SimpleNetwork(0, layers, Connection[], it)
end

export SimpleNetwork

import Base.reverse!

function reverse!(d::DataStructures.OrderedDict)
    reverse!(d.keys)
    reverse!(d.vals)
    DataStructures.rehash!(d)
end

"""
	reverselayerorder!(net::SimpleNetwork)

"""
function reverselayerorder!(net::SimpleNetwork)
	reverse!(net.layers)
	net.layeriterator = values(net.layers)
end
export reverselayerorder!

"""
    connect!(net::SimpleNetwork, from::Symbol, to::Symbol, conorkind; label = :default)
"""
function connect!(net::SimpleNetwork, from::Symbol, to::Symbol, conorkind;
				  label = :default)
	contype = connect!(net.layers[from], net.layers[to], conorkind, label = label)
	if ismatch(r"Plastic", contype)
		push!(net.plasticconnections,
		      net.layers[to].inputconnections[label][from])
	end
end

"""
	addlayer!(net::SimpleNetwork, layer::Layer)
"""
function addlayer!(net::SimpleNetwork, layer::Layer)
	net.layers[layer.name] = layer
	net.layeriterator = values(net.layers)
end
export addlayer!

"""
	updatenet!(net::SimpleNetwork)

Call updateneuron! for all layers in `net.layers`
and updateconnection! for all connections in `net.plasticconnections`
and increment `net.time`.
"""
function updatenet!(net::SimpleNetwork)
	net.time += 1
	for l in net.layeriterator
		collectmessages!(l.neurons, l.inputconnections)
	end
	for l in net.layeriterator
		updateneuron!(l.neurons)
	end
	for c in net.plasticconnections
		updateconnection!(c)
	end
end
export updatenet!

function updatenetasync!(net::SimpleNetwork)
	net.time += 1
	for l in net.layeriterator
		collectmessages!(l.neurons, l.inputconnections)
		updateneuron!(l.neurons)
	end
	for c in net.plasticconnections
		updateconnection!(c)
	end
end
export updatenetasync!

max_stars_repo_name: jbrea/NetSim

max_issues_repo_name: jbrea/NetSim

max_forks_repo_name: jbrea/NetSim


text:
abstract type AbstractOptions end
mutable struct Options <: AbstractOptions
    n::Int64
    v::Bool

    Options(n::Int64 = 0, v::Bool = false) = new(n, v)
end

function __repr__(self::AbstractOptions)::String
    return AbstractOptions(self.n, self.v)
end

function __eq__(self::AbstractOptions, other::AbstractOptions)::Bool
    return __key(self) == __key(other)
end

function __key(self::AbstractOptions)
    (self.n, self.v)
end

function fib(i::Int64)::Int64
    if i == 0 || i == 1
        return 1
    end
    return fib(i - 1) + fib(i - 2)
end

if abspath(PROGRAM_FILE) == @__FILE__
    args = parse_args(Options)
    if args.n == 0
        args.n = 5
    end
    println(fib(args.n))
end


text:
abstract type OrdinaryDiffEqAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
abstract type OrdinaryDiffEqAdaptiveAlgorithm <: OrdinaryDiffEqAlgorithm end
abstract type OrdinaryDiffEqCompositeAlgorithm <: OrdinaryDiffEqAlgorithm end

abstract type OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} <: OrdinaryDiffEqAdaptiveAlgorithm end
abstract type OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} <: OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} end
abstract type OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} <: OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} end

abstract type OrdinaryDiffEqImplicitAlgorithm{CS,AD} <: OrdinaryDiffEqAlgorithm end
abstract type OrdinaryDiffEqNewtonAlgorithm{CS,AD} <:  OrdinaryDiffEqImplicitAlgorithm{CS,AD} end
abstract type OrdinaryDiffEqRosenbrockAlgorithm{CS,AD} <:  OrdinaryDiffEqImplicitAlgorithm{CS,AD} end

abstract type OrdinaryDiffEqExponentialAlgorithm <: OrdinaryDiffEqAlgorithm end
abstract type OrdinaryDiffEqAdaptiveExponentialAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm end
const ExponentialAlgorithm = Union{OrdinaryDiffEqExponentialAlgorithm,OrdinaryDiffEqAdaptiveExponentialAlgorithm}

abstract type OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm <: OrdinaryDiffEqAdaptiveAlgorithm end

struct FunctionMap{scale_by_time} <: OrdinaryDiffEqAlgorithm end
FunctionMap(;scale_by_time=false) = FunctionMap{scale_by_time}()

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

# RK methods

struct ExplicitRK{TabType} <: OrdinaryDiffEqAdaptiveAlgorithm
  tableau::TabType
end
ExplicitRK(;tableau=ODE_DEFAULT_TABLEAU) = ExplicitRK(tableau)

@inline trivial_limiter!(u, f, t) = nothing

struct Euler <: OrdinaryDiffEqAlgorithm end
struct Heun <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Ralston <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Midpoint <: OrdinaryDiffEqAdaptiveAlgorithm end
struct RK4 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Anas5{T} <: OrdinaryDiffEqAlgorithm
  w::T
end
Anas5(; w=1) = Anas5(w)

struct OwrenZen3 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct OwrenZen4 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct OwrenZen5 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct CarpenterKennedy2N54 <: OrdinaryDiffEqAlgorithm end
struct SSPRK22{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK22(stage_limiter! = trivial_limiter!) = SSPRK22(stage_limiter!, trivial_limiter!)
struct SSPRK33{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK33(stage_limiter! = trivial_limiter!) = SSPRK33(stage_limiter!, trivial_limiter!)
struct SSPRK53{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK53(stage_limiter! = trivial_limiter!) = SSPRK53(stage_limiter!, trivial_limiter!)
struct SSPRK63{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK63(stage_limiter! = trivial_limiter!) = SSPRK63(stage_limiter!, trivial_limiter!)
struct SSPRK73{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK73(stage_limiter! = trivial_limiter!) = SSPRK73(stage_limiter!, trivial_limiter!)
struct SSPRK83{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK83(stage_limiter! = trivial_limiter!) = SSPRK83(stage_limiter!, trivial_limiter!)
struct SSPRK432{StageLimiter,StepLimiter} <: OrdinaryDiffEqAdaptiveAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK432(stage_limiter! = trivial_limiter!) = SSPRK432(stage_limiter!, trivial_limiter!)
struct SSPRK932{StageLimiter,StepLimiter} <: OrdinaryDiffEqAdaptiveAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK932(stage_limiter! = trivial_limiter!) = SSPRK932(stage_limiter!, trivial_limiter!)
struct SSPRK54{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK54(stage_limiter! = trivial_limiter!) = SSPRK54(stage_limiter!, trivial_limiter!)
struct SSPRK104{StageLimiter,StepLimiter} <: OrdinaryDiffEqAlgorithm
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end
SSPRK104(stage_limiter! = trivial_limiter!) = SSPRK104(stage_limiter!, trivial_limiter!)
struct BS3 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct DP5 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct DP5Threaded <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Tsit5 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct DP8 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct TanYam7 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct TsitPap8 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Feagin10 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Feagin12 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct Feagin14 <: OrdinaryDiffEqAdaptiveAlgorithm end

struct BS5 <: OrdinaryDiffEqAdaptiveAlgorithm
  lazy::Bool
  BS5(;lazy=true) = new(lazy)
end
struct Vern6 <: OrdinaryDiffEqAdaptiveAlgorithm
  lazy::Bool
  Vern6(;lazy=true) = new(lazy)
end
struct Vern7 <: OrdinaryDiffEqAdaptiveAlgorithm
  lazy::Bool
  Vern7(;lazy=true) = new(lazy)
end
struct Vern8 <: OrdinaryDiffEqAdaptiveAlgorithm
  lazy::Bool
  Vern8(;lazy=true) = new(lazy)
end
struct Vern9 <: OrdinaryDiffEqAdaptiveAlgorithm
  lazy::Bool
  Vern9(;lazy=true) = new(lazy)
end

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

# Symplectic methods

struct SymplecticEuler <: OrdinaryDiffEqAlgorithm end
struct VelocityVerlet <: OrdinaryDiffEqAlgorithm end
struct VerletLeapfrog <: OrdinaryDiffEqAlgorithm end
struct PseudoVerletLeapfrog <: OrdinaryDiffEqAlgorithm end
struct McAte2 <: OrdinaryDiffEqAlgorithm end
struct Ruth3 <: OrdinaryDiffEqAlgorithm end
struct McAte3 <: OrdinaryDiffEqAlgorithm end
struct CandyRoz4 <: OrdinaryDiffEqAlgorithm end
struct McAte4 <: OrdinaryDiffEqAlgorithm end
struct CalvoSanz4 <: OrdinaryDiffEqAlgorithm end
struct McAte42 <: OrdinaryDiffEqAlgorithm end
struct McAte5 <: OrdinaryDiffEqAlgorithm end
struct Yoshida6 <: OrdinaryDiffEqAlgorithm end
struct KahanLi6 <: OrdinaryDiffEqAlgorithm end
struct McAte8 <: OrdinaryDiffEqAlgorithm end
struct KahanLi8 <: OrdinaryDiffEqAlgorithm end
struct SofSpa10 <: OrdinaryDiffEqAlgorithm end

# Nyström methods

struct IRKN3 <: OrdinaryDiffEqAlgorithm end
struct Nystrom4 <: OrdinaryDiffEqAlgorithm end
struct Nystrom4VelocityIndependent <: OrdinaryDiffEqAlgorithm end
struct IRKN4 <: OrdinaryDiffEqAlgorithm end
struct Nystrom5VelocityIndependent <: OrdinaryDiffEqAlgorithm end
struct DPRKN6 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct DPRKN8 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct DPRKN12 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct ERKN4 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct ERKN5 <: OrdinaryDiffEqAdaptiveAlgorithm end

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

# Adams Bashforth and Adams moulton methods

struct AB3 <: OrdinaryDiffEqAlgorithm end
struct AB4 <: OrdinaryDiffEqAlgorithm end
struct AB5 <: OrdinaryDiffEqAlgorithm end
struct ABM32 <: OrdinaryDiffEqAlgorithm end
struct ABM43 <: OrdinaryDiffEqAlgorithm end
struct ABM54 <: OrdinaryDiffEqAlgorithm end

# Variable Step Size Adams methods

struct VCAB3 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct VCAB4 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct VCAB5 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct VCABM3 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct VCABM4 <: OrdinaryDiffEqAdaptiveAlgorithm end
struct VCABM5 <: OrdinaryDiffEqAdaptiveAlgorithm end

# Variable Order and Variable Step Size Adams methods

struct VCABM <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm end

# IMEX Multistep methods

struct CNAB2{CS,AD,F,F2,FDT,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
CNAB2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                      extrapolant=:linear,new_jac_conv_bound = 1e-3) =
                      CNAB2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                      typeof(new_jac_conv_bound)}(
                      linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound)

struct CNLF2{CS,AD,F,F2,FDT,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
CNLF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                      extrapolant=:linear,new_jac_conv_bound = 1e-3) =
                      CNLF2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                      typeof(new_jac_conv_bound)}(
                      linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound)

struct QNDF1{CS,AD,F,F2,FDT,T2,κType,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
  kappa::κType
end
QNDF1(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                 linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                                  extrapolant=:linear,new_jac_conv_bound = 1e-3,kappa = -0.1850,
                 controller = :Predictive) =
                 QNDF1{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                 typeof(new_jac_conv_bound),typeof(kappa),controller}(
                 linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound,kappa)

QBDF1(;kwargs...) = QNDF1(;kappa=0,kwargs...)

struct QNDF2{CS,AD,F,F2,FDT,T2,κType,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
  kappa::κType
end
QNDF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                 linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                                  extrapolant=:linear,new_jac_conv_bound = 1e-3,kappa = -1//9,
                 controller = :Predictive) =
                 QNDF2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                 typeof(new_jac_conv_bound),typeof(kappa),controller}(
                 linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound,kappa)

QBDF2(;kwargs...) = QNDF2(;kappa=0,kwargs...)

struct QNDF{CS,AD,F,F2,FDT,K,T,T2,κType,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  κ::K
  tol::T
  extrapolant::Symbol
  new_jac_conv_bound::T2
  kappa::κType
end
Base.@pure QNDF(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),κ=nothing,tol=nothing,
                extrapolant=:linear,new_jac_conv_bound = 1e-3,kappa=promote(-0.1850,-1//9,-0.0823,-0.0415,0),
                controller = :Predictive) =
                QNDF{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                typeof(κ),typeof(tol),typeof(new_jac_conv_bound),typeof(kappa),controller}(
                linsolve,nlsolve,diff_type,κ,tol,extrapolant,new_jac_conv_bound,kappa)

Base.@pure QBDF(;kwargs...) = QNDF(;kappa=tuple(0//1,0//1,0//1,0//1,0//1),kwargs...)

struct SBDF{CS,AD,F,F2,FDT,K,T,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  κ::K
  tol::T
  extrapolant::Symbol
  new_jac_conv_bound::T2
  order::Int
end
SBDF(order;chunk_size=0,autodiff=true,diff_type=Val{:central},
     linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),κ=nothing,tol=nothing,
     extrapolant=:linear,new_jac_conv_bound = 1e-3) =
     SBDF{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
     typeof(κ),typeof(tol),typeof(new_jac_conv_bound)}(
     linsolve,nlsolve,diff_type,κ,tol,extrapolant,new_jac_conv_bound,order)
IMEXEuler(;kwargs...) = SBDF(1;kwargs...)
SBDF2(;kwargs...) = SBDF(2;kwargs...)
SBDF3(;kwargs...) = SBDF(3;kwargs...)
SBDF4(;kwargs...) = SBDF(4;kwargs...)

# Adams/BDF methods in Nordsieck forms
struct AN5   <: OrdinaryDiffEqAdaptiveAlgorithm end
struct JVODE{bType,aType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm
  algorithm::Symbol
  bias1::bType
  bias2::bType
  bias3::bType
  addon::aType
end

JVODE(algorithm=:Adams;bias1=6, bias2=6,bias3=10,
                 addon=1//10^6) = JVODE(algorithm,bias1,bias2,bias3,addon)
JVODE_Adams(;kwargs...) = JVODE(:Adams;kwargs...)
JVODE_BDF(;kwargs...) = JVODE(:BDF;kwargs...)

# ROCK methods
struct ROCK2 <: OrdinaryDiffEqAdaptiveAlgorithm end

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

# Generic implicit methods

struct GenericImplicitEuler{F} <: OrdinaryDiffEqAdaptiveAlgorithm
  nlsolve::F
  extrapolant::Symbol
end
GenericImplicitEuler(;
            nlsolve=NLSOLVEJL_SETUP(),extrapolant=:linear) =
            GenericImplicitEuler{typeof(nlsolve)}(nlsolve,extrapolant)

struct GenericTrapezoid{F} <: OrdinaryDiffEqAdaptiveAlgorithm
  nlsolve::F
  extrapolant::Symbol
end
GenericTrapezoid(;
            nlsolve=NLSOLVEJL_SETUP(),extrapolant=:linear) =
            GenericTrapezoid{typeof(nlsolve)}(nlsolve,extrapolant)

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

# Linear Methods

struct MidpointSplitting <: OrdinaryDiffEqAlgorithm end

struct LinearExponential <: OrdinaryDiffEqExponentialAlgorithm
  krylov::Symbol
  m::Int
  iop::Int
end
LinearExponential(;krylov=:off, m=10, iop=0) = LinearExponential(krylov, m, iop)

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

# SDIRK Methods

struct ImplicitEuler{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
ImplicitEuler(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                          extrapolant=:constant,new_jac_conv_bound=1e-3,
                          controller=:Predictive) =
                          ImplicitEuler{chunk_size,autodiff,typeof(linsolve),
                          typeof(nlsolve),typeof(diff_type),
                          typeof(new_jac_conv_bound),controller}(linsolve,
                          nlsolve,diff_type,extrapolant,new_jac_conv_bound)

struct ImplicitMidpoint{CS,AD,F,F2,FDT,T2} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
ImplicitMidpoint(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                      extrapolant=:linear,new_jac_conv_bound=1e-3) =
                      ImplicitMidpoint{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                      typeof(new_jac_conv_bound)}(
                      linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound)

struct Trapezoid{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Trapezoid(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                                            extrapolant=:linear,new_jac_conv_bound = 1e-3,
                      controller = :PI) =
                      Trapezoid{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                      typeof(new_jac_conv_bound),controller}(
                      linsolve,nlsolve,diff_type,extrapolant,new_jac_conv_bound)

struct TRBDF2{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
TRBDF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                 linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                 smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                 controller = :Predictive) =
TRBDF2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
      typeof(new_jac_conv_bound),controller}(
      linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct SDIRK2{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
SDIRK2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 SDIRK2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct SSPSDIRK2{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD} # Not adaptive
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
SSPSDIRK2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:constant,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 SSPSDIRK2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct Kvaerno3{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Kvaerno3(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 Kvaerno3{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct KenCarp3{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
KenCarp3(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 KenCarp3{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct Cash4{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
  embedding::Int
end
Cash4(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive,embedding=3) =
 Cash4{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound,embedding)

struct Hairer4{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Hairer4(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 Hairer4{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct Hairer42{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Hairer42(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 Hairer42{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct Kvaerno4{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Kvaerno4(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 Kvaerno4{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct Kvaerno5{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
Kvaerno5(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 Kvaerno5{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct KenCarp4{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
KenCarp4(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 KenCarp4{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

struct KenCarp5{CS,AD,F,F2,FDT,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
KenCarp5(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:linear,new_jac_conv_bound = 1e-3,
                   controller = :Predictive) =
 KenCarp5{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
        typeof(new_jac_conv_bound),controller}(
        linsolve,nlsolve,diff_type,smooth_est,extrapolant,new_jac_conv_bound)

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

# Rosenbrock Methods

struct Rosenbrock23{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rosenbrock23(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rosenbrock23{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rosenbrock32{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rosenbrock32(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rosenbrock32{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct ROS3P{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
ROS3P(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = ROS3P{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rodas3{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rodas3(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rodas3{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct RosShamp4{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
RosShamp4(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = RosShamp4{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Veldd4{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Veldd4(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Veldd4{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Velds4{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Velds4(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Velds4{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct GRK4T{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
GRK4T(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = GRK4T{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct GRK4A{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
GRK4A(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = GRK4A{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Ros4LStab{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Ros4LStab(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Ros4LStab{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rodas4{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rodas4(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rodas4{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rodas42{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rodas42(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rodas42{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rodas4P{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rodas4P(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rodas4P{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct Rodas5{CS,AD,F,FDT} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  linsolve::F
  diff_type::FDT
end
Rodas5(;chunk_size=0,autodiff=true,diff_type=Val{:central},linsolve=DEFAULT_LINSOLVE) = Rodas5{chunk_size,autodiff,typeof(linsolve),typeof(diff_type)}(linsolve,diff_type)

struct GeneralRosenbrock{CS,AD,F,TabType} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD}
  tableau::TabType
  factorization::F
end

GeneralRosenbrock(;chunk_size=0,autodiff=true,
                    factorization=lu!,tableau=ROSENBROCK_DEFAULT_TABLEAU) =
                    GeneralRosenbrock{chunk_size,autodiff,typeof(factorization),typeof(tableau)}(tableau,factorization)

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

struct GenericIIF1{F} <: OrdinaryDiffEqExponentialAlgorithm
  nlsolve::F
end
GenericIIF1(;nlsolve=NLSOLVEJL_SETUP()) = GenericIIF1{typeof(nlsolve)}(nlsolve)

struct GenericIIF2{F} <: OrdinaryDiffEqExponentialAlgorithm
  nlsolve::F
end
GenericIIF2(;nlsolve=NLSOLVEJL_SETUP()) = GenericIIF2{typeof(nlsolve)}(nlsolve)

for Alg in [:LawsonEuler, :NorsettEuler, :ETDRK2, :ETDRK3, :ETDRK4, :HochOst4]
  @eval struct $Alg <: OrdinaryDiffEqExponentialAlgorithm
    krylov::Bool
    m::Int
    iop::Int
    autodiff::Bool
    chunksize::Int
  end
  @eval $Alg(;krylov=false, m=30, iop=0, autodiff=true, chunksize=0) = $Alg(krylov, m, iop, autodiff, chunksize)
end
ETD1 = NorsettEuler # alias
for Alg in [:Exprb32, :Exprb43]
  @eval struct $Alg <: OrdinaryDiffEqAdaptiveExponentialAlgorithm
    m::Int
    iop::Int
    autodiff::Bool
    chunksize::Int
  end
  @eval $Alg(;m=30, iop=0, autodiff=true, chunksize=0) = $Alg(m, iop, autodiff, chunksize)
end
for Alg in [:Exp4, :EPIRK4s3A, :EPIRK4s3B, :EPIRK5s3, :EXPRB53s3, :EPIRK5P1, :EPIRK5P2]
  @eval struct $Alg <: OrdinaryDiffEqExponentialAlgorithm
    adaptive_krylov::Bool
    m::Int
    iop::Int
    autodiff::Bool
    chunksize::Int
  end
  @eval $Alg(;adaptive_krylov=true, m=30, iop=0, autodiff=true, chunksize=0) = $Alg(adaptive_krylov, m, iop, autodiff, chunksize)
end
struct SplitEuler <: OrdinaryDiffEqExponentialAlgorithm end
struct ETD2 <: OrdinaryDiffEqExponentialAlgorithm end

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

struct ABDF2{CS,AD,F,F2,FDT,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::FDT
  κ::K
  tol::T
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
end
ABDF2(;chunk_size=0,autodiff=true,diff_type=Val{:central},
      κ=nothing,tol=nothing,linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
      smooth_est=true,extrapolant=:linear,new_jac_conv_bound=1e-3,
      controller=:Predictive) =
ABDF2{chunk_size,autodiff,typeof(linsolve),typeof(nlsolve),typeof(diff_type),
      typeof(κ),typeof(tol),typeof(new_jac_conv_bound),controller}(
      linsolve,nlsolve,diff_type,κ,tol,smooth_est,extrapolant,new_jac_conv_bound)

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

struct CompositeAlgorithm{T,F} <: OrdinaryDiffEqCompositeAlgorithm
  algs::T
  choice_function::F
end

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

### Algorithm Groups

const MassMatrixAlgorithms = Union{OrdinaryDiffEqRosenbrockAlgorithm,
                                   OrdinaryDiffEqRosenbrockAdaptiveAlgorithm,
                                   ImplicitEuler,ImplicitMidpoint}

const MultistepAlgorithms = Union{IRKN3,IRKN4,
                                  ABDF2,
                                  AB3,AB4,AB5,ABM32,ABM43,ABM54}

const SplitAlgorithms = Union{CNAB2,CNLF2,SBDF,
                              GenericIIF1,GenericIIF2,
                              KenCarp3,KenCarp4,KenCarp5}


text:
#
# Test inference setup/interface
#

let
    bn = rand_discrete_bn()

    @test_throws ArgumentError InferenceState(bn, [:waldo, :N3])
    @test_throws ArgumentError InferenceState(bn, [:N2, :N3], Assignment(:N1 => 1, :N3 => 3, :N7 => 2016))
end

let
    bn = rand_discrete_bn()
    inf = InferenceState(bn, [:N3, :N5, :N2, :N3])

    @test inf.query == [:N3, :N5, :N2]
end

let
    bn = rand_discrete_bn()
    inf = InferenceState(bn, :N6)

    @test inf.query == [:N6]
end

let
    bn = rand_discrete_bn()

    @test_throws ArgumentError InferenceState(bn, [:N3, :waldo, :N5])
    @test_throws ArgumentError InferenceState(bn, :N1, Assignment(:N1=>2))
end

let
    # A → C ← B
    bn = DiscreteBayesNet()
    push!(bn, DiscreteCPD(:a, [1.0,0.0]))
    push!(bn, DiscreteCPD(:b, [0.0,1.0]))
    push!(bn, DiscreteCPD(:c, [:a, :b], [2,2], [Categorical([0.1,0.9]),
                                                Categorical([0.2,0.8]),
                                                Categorical([1.0,0.0]),
                                                Categorical([0.4,0.6]),
                                                ]))

    for im in [ExactInference(), LikelihoodWeightingInference(), LoopyBelief(), GibbsSamplingNodewise(), GibbsSamplingFull()]
        srand(0)
        ϕ = infer(im, bn, :a)::Factor
        @test length(ϕ) == 2
        f = ϕ[:a=>1]::Factor
        @test length(f) == 1
        @test isapprox(f.potential[1], 1.0, atol=0.02)

        f = ϕ[:a=>2]::Factor
        @test length(f) == 1
        @test isapprox(f.potential[1], 0.0, atol=0.02)
    end

    for im in [ExactInference(), LikelihoodWeightingInference(), LoopyBelief(), GibbsSamplingNodewise(), GibbsSamplingFull()]
        srand(0)
        ϕ = infer(im, bn, :c)::Factor
        @test length(ϕ) == 2
        @test isapprox(ϕ[:c=>1].potential[1], 1.0, atol=0.02)
        @test isapprox(ϕ[:c=>2].potential[1], 0.0, atol=0.02)
    end

    for im in [ExactInference(), LikelihoodWeightingInference(), GibbsSamplingNodewise(), GibbsSamplingFull()] # LoopyBelief(),
        srand(0)
        ϕ = infer(im, bn, [:b, :c])
        @test size(ϕ) == (2,2)
        @test isapprox(ϕ[:b=>1, :c=>1].potential[1], 0.0, atol=0.02)
        @test isapprox(ϕ[:b=>2, :c=>1].potential[1], 1.0, atol=0.02)
        @test isapprox(ϕ[:b=>1, :c=>2].potential[1], 0.0, atol=0.02)
        @test isapprox(ϕ[:b=>2, :c=>2].potential[1], 0.0, atol=0.02)
    end

    # Student Example (PGM page 53)
    bn = DiscreteBayesNet()
    push!(bn, DiscreteCPD(:D, [0.6,0.4])) # difficulty
    push!(bn, DiscreteCPD(:I, [0.7,0.3])) # intelligence
    push!(bn, DiscreteCPD(:G, [:D, :I], [2,2], [Categorical([0.3,0.4,0.3]), # grade
                                                Categorical([0.9,0.08,0.02]),
                                                Categorical([0.05,0.25,0.7]),
                                                Categorical([0.5,0.3,0.2]),
                                                ]))
    push!(bn, DiscreteCPD(:L, [:G], [3], [Categorical([0.1,0.9]), # letter
                                          Categorical([0.4,0.6]),
                                          Categorical([0.99,0.01]),
                                         ]))
    push!(bn, DiscreteCPD(:S, [:I], [2], [Categorical([0.95,0.05]), # SAT
                                          Categorical([0.2,0.8]),
                                         ]))

    # P(D) = [0.6, 0.4]
    for im in [ExactInference(), LikelihoodWeightingInference(), LoopyBelief(), GibbsSamplingNodewise(), GibbsSamplingFull()]
        srand(0)
        ϕ = infer(im, bn, :D)
        @test size(ϕ) == (2,)
        @test isapprox(ϕ[:D=>1].potential[1], 0.6, atol=0.05)
        @test isapprox(ϕ[:D=>2].potential[1], 0.4, atol=0.05)
    end

    # P(G|d₁, i₁) = [0.3, 0.4, 0.3]
    for im in [ExactInference(), LikelihoodWeightingInference(), LoopyBelief(), GibbsSamplingNodewise(), GibbsSamplingFull()]
        srand(0)
        ϕ = infer(im, bn, :G, evidence=Assignment(:D=>1, :I=>1))
        @test size(ϕ) == (3,)
        @test isapprox(ϕ.potential[1], 0.3, atol=0.05)
        @test isapprox(ϕ.potential[2], 0.4, atol=0.05)
        @test isapprox(ϕ.potential[1], 0.3, atol=0.05)
    end
end


text:
# Note: This code assumes the presence of a global variable `project_root_dir`.
import BSON

function clear_cache!()
    empty!(results_cache)
end

function unload_cache!()
    global results_cache = nothing
end

function save_cache!(result_group)
    save_path = joinpath(project_root_dir, "data/$result_group/cache.bson")
    @assert !isfile(save_path)
    BSON.bson(save_path, results_cache)
end

function load_cache(result_group)
    result_cache_file = joinpath(project_root_dir, "data/$result_group/cache.bson")
    if !isfile(result_cache_file)
        nothing
    else
        BSON.load(result_cache_file)
    end
end

macro run_cached(assigned_computation_expr)
    @assert assigned_computation_expr.head == :(=)
    var, fun = assigned_computation_expr.args

    quote
        $(esc(var)) = run_cached!(result_group, $(Meta.quot(var))) do
            $fun
        end
    end
end

function run_cached!(f, result_group, key; force_run = false)
    @assert !isempty(result_group)
    prefixed_key = Symbol(result_group, :., key)
    result = force_run ? f() : get(f, results_cache, prefixed_key)
    results_cache[prefixed_key] = result
    result
end

function load_cache_if_not_defined!(the_result_group; filly_emtpy = true)
    global results_cache = if !isdefined(Main, :results_cache) || isnothing(results_cache)
        loaded_cache = load_cache(the_result_group)
        if isnothing(loaded_cache)
            @info "No persisted results cache file found. Resuming with an empty cache."
            Dict()
        else
            @info "Loaded cached results for group \"$the_result_group\" from file!"
            loaded_cache
        end
    elseif isdefined(Main, :results_cache) &&
           any(!startswith(string(k), "$the_result_group.") for k in keys(results_cache))
        error("Skipping. Cache contains results for another group. In order to proceed, call
              `unload_cache()`. If you want to keep the cached results, make sure to call
              `save_cache!()` first.")
    elseif results_cache isa Dict
        @info "Cache for group \"$the_result_group\" present. No additional cache loaded."
        results_cache
    else
        error("Unknown cache type.")
    end

    global result_group = the_result_group
end


text:
using Test
using DataFrames, TableOperations, Tables, Random
using Recommenders:
    Movielens100k,
    load_dataset,
    ratio_split,
    ImplicitMF,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG,
    EvaluateValidData

ml100k = Movielens100k()
download(ml100k)
rating, _, _ = load_dataset(ml100k)

Random.seed!(1234);
train_table, test_table = ratio_split(rating, 0.8)

prec10 = MeanPrecision(10)
recall10 = MeanRecall(10)
ndcg10 = MeanNDCG(10)
metrics = [prec10, recall10, ndcg10]

cb = EvaluateValidData(ndcg10, test_table, 1)

model = ImplicitMF(16, true, 0.01)
result = evaluate_u2i(
    model,
    train_table,
    test_table,
    metrics,
    10,
    callbacks = [cb],
    col_item = :movieid,
    n_epochs = 20,
    n_negatives = 1,
    learning_rate = 0.01,
    drop_history = true,
    verbose = 1,
)

@test !(model.user_embedding === nothing)
@test !(model.item_embedding === nothing)
@test !(model.user_bias === nothing)
@test !(model.item_bias === nothing)
@test !(model.μ === nothing)

@test result[:ndcg10] >= 0
@test result[:precision10] >= 0
@test result[:recall10] >= 0

tempfile = tempname()
save_model(model, tempfile)
@test isfile(tempfile)
model = load_model(tempfile)
@test model.user_embedding !== nothing
@test model.item_embedding !== nothing


text:
function test_io(AT)
    @testset "input/output" begin
        @testset "showing" begin
          io = IOBuffer()
          A = AT(Int64[1])
          B = AT(Int64[1 2;3 4]) # vectors and non-vector arrays showing
                                 # are handled differently in base/arrayshow.jl

          show(io, MIME("text/plain"), A)
          seekstart(io)
          @test occursin(Regex("^1-element $AT{Int64,1.*}:\n 1\$"), String(take!(io)))

          show(io, A)
          seekstart(io)
          msg = String(take!(io)) # result of e.g. `print` differs on 32bit and 64bit machines
          # due to different definition of `Int` type
          # print([1]) shows as [1] on 64bit but Int64[1] on 32bit
          @test msg == "[1]" || msg == "Int64[1]"

          show(io, MIME("text/plain"), B)
          seekstart(io)
          @test occursin(Regex("^2×2 $AT{Int64,2.*}:\n 1  2\n 3  4\$"), String(take!(io)))

          show(io, B)
          seekstart(io)
          msg = String(take!(io))
          @test msg == "[1 2; 3 4]" || msg == "Int64[1 2; 3 4]"

          show(io, MIME("text/plain"), A')
          seekstart(io)
          msg = String(take!(io)) # the printing of Adjoint depends on global state
          @test occursin(Regex("^1×1 Adjoint{Int64,$AT{Int64,1.*}}:\n 1\$"), msg) ||
                occursin(Regex("^1×1 LinearAlgebra.Adjoint{Int64,$AT{Int64,1.*}}:\n 1\$"), msg)
        end
    end
end


text:
using DelimitedFiles
using OffsetArrays
using PyPlot
using ImageFiltering

export readImage, readSimImage, writeImage, saveImage

function readImage(path::String)
        image = readdlm(path, Float64)
        z = reshape(image[:, 3], (301, 301))
        z = centered(z) # creates an OffsetArray where the x and y indices are from -150 to 150 instead of 1 to 301
        return z
end

function readSimImage(path::String, column_nr::Int64)
        image = readdlm(path, Float64)
        z = reshape(image[:, column_nr], (301, 301))
        z = copy(transpose(centered(z)))
        return z
end

function writeImage(path::String, image::OffsetArray{Float64,2,Array{Float64,2}})
        x = reshape(repeat(transpose(collect(axes(image)[1])), 301, 1), (90601, 1))
        y = reshape(repeat(collect(reverse(axes(image)[2])), 1, 301), (90601, 1))
        z = reshape(image, (90601, 1))
        open(path, "w") do io
           writedlm(io, [x y z])
        end
end

function saveImage(path::String, image::Array{Float64,2})
# save current figure as pdf (path = name.pdf) or as jpg (path = name.jpg )
        savefig(path, bbox_inches="tight")
end

function saveImage(path::String, image::OffsetArray{Float64,2,Array{Float64,2}})
# save current figure as pdf (path = name.pdf) or as jpg (path = name.jpg )
        savefig(path, bbox_inches="tight")
end


text:
__precompile__()

module SmartParser

using DataStructures

export TPattern, SIMILARITY_LEVEL
export RTYPE, __DEFAULT__R__
include("types_and_settings.jl")

include("__REFDIC__.jl")

export MASK_RULES, MASK_RULES_DIC_INV
include("rules.jl")

export mask, tokenize, tokenize0, load_file, encode_line, revert
export preprocess_raw_input
include("tokenize.jl")

export Block, Singleline, MultiDict, tree_print
export khash, copy, isequal
export DFS
export children, label, num_nodes
export max_depth, min_depth
export collect_action_dfs, collect_action_bfs
export is_multi, is_single
export Block
include("Tree.jl")

export increaseindex!
include("special_dict_op.jl")

export similarity
export MostFreqSimilarSubsq
include("block_similarity.jl")

export find_block_MostFreqSimilarSubsq
export loop_until_stable
export verify_block, is_valid_block, is_valid_x, is_valid_C
export build_block_init_by_linebreak
export build_block_init, typical_blocks
export merge_children
include("find_block.jl")

export block_print, treep
include("block_print.jl")

export parse_file!
export extract_DATA
include("parse_file.jl")

export StructuredOutputFile
export structurize_file
include("structurize_file.jl")

export lookup_codes,            lookup_code
export get_DATA_by_codes,       get_DATA_by_code
export get_n_blocks_by_codes,   get_n_blocks_by_code
export next_block_by_codes,     next_block_by_code
export get_blocks_max_by_codes, get_blocks_max_by_code
export select_block_by_patt
export get_DATA_by_patt
include("search.jl")

end # module

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
@testset "MTL" begin

@testset "devices" begin

devs = devices()
@test length(devs) > 0

dev = first(devs)
@test dev == devs[1]

if length(devs) > 1
    @test dev != devs[2]
end

compact_str = sprint(io->show(io, dev))
full_str = sprint(io->show(io, MIME"text/plain"(), dev))

@test dev.name isa String
@test dev.isLowPower isa Bool
@test dev.isRemovable isa Bool
@test dev.hasUnifiedMemory isa Bool
@test dev.registryID isa Integer
@test dev.maxTransferRate isa Integer

@test dev.recommendedMaxWorkingSetSize isa Integer
@test dev.maxThreadgroupMemoryLength isa Integer
@test dev.maxThreadsPerThreadgroup isa MTL.MtSize
@test dev.maxBufferLength isa Integer

@test dev.currentAllocatedSize isa Integer

end

@testset "compile options" begin

opts = MtlCompileOptions()

compact_str = sprint(io->show(io, opts))
full_str = sprint(io->show(io, MIME"text/plain"(), opts))

@test opts.fastMathEnabled isa Bool
val = !opts.fastMathEnabled
opts.fastMathEnabled = val
@test opts.fastMathEnabled == val

@test opts.languageVersion isa VersionNumber
opts.languageVersion = v"1.0"
@test opts.languageVersion == v"1.0"

end

@testset "libraries" begin

dev = first(devices())
opts = MtlCompileOptions()

let lib = MtlLibrary(dev, "", opts)
    @test lib.device == dev
    @test lib.label == nothing
    @test isempty(lib.functionNames)
end

metal_code = read(joinpath(@__DIR__, "dummy.metal"), String)
let lib = MtlLibrary(dev, metal_code, opts)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns
end

binary_path = joinpath(@__DIR__, "dummy.metallib")
let lib = MtlLibraryFromFile(dev, binary_path)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns
end

binary_code = read(binary_path)
let lib = MtlLibraryFromData(dev, binary_code)
    @test lib.device == dev
    @test lib.label == nothing
    fns = lib.functionNames
    @test length(fns) == 2
    @test "kernel_1" in fns
    @test "kernel_2" in fns

    compact_str = sprint(io->show(io, lib))
    full_str = sprint(io->show(io, MIME"text/plain"(), lib))
end

end

@testset "functions" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "dummy.metallib"))
fun = MtlFunction(lib, "kernel_1")

compact_str = sprint(io->show(io, fun))
full_str = sprint(io->show(io, MIME"text/plain"(), fun))

@test fun.device == dev
@test fun.label == nothing
@test fun.name == "kernel_1"
@test fun.functionType == MTL.MtFunctionTypeKernel

end

@testset "events" begin

dev = first(devices())

let ev = MtlEvent(dev)
    @test ev.device == dev
    @test ev.label == nothing
end

let ev = MtlSharedEvent(dev)
    @test ev.device == dev
    @test ev.label == nothing
    @test ev.signaledValue == 0
end

end

@testset "fences" begin

dev = first(devices())

let fen = MtlFence(dev)
    @test fen.device == dev
end

end

@testset "heap" begin

dev = first(devices())

let desc = MtlHeapDescriptor()
    @test desc.type == MTL.MtHeapTypeAutomatic
    desc.type = MTL.MtHeapTypePlacement
    @test desc.type == MTL.MtHeapTypePlacement

    @test desc.size == 0
    desc.size = 1024
    @test desc.size == 1024

    @test desc.storageMode == MTL.MtStorageModePrivate
    desc.storageMode = MTL.MtStorageModeShared
    @test desc.storageMode == MTL.MtStorageModeShared

    @test desc.cpuCacheMode == MTL.MtCPUCacheModeDefaultCache
    desc.cpuCacheMode = MTL.MtCPUCacheModeWriteCombined
    @test desc.cpuCacheMode == MTL.MtCPUCacheModeWriteCombined

    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeDefault
    desc.hazardTrackingMode = MTL.MtHazardTrackingModeUntracked
    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeUntracked

    @test desc.resourceOptions == MTL.MtResourceStorageModeShared |
                                  MTL.MtResourceCPUCacheModeWriteCombined |
                                  MTL.MtResourceHazardTrackingModeUntracked
    desc.resourceOptions = MTL.MtResourceStorageModePrivate |
                           MTL.MtResourceCPUCacheModeDefaultCache |
                           MTL.MtResourceHazardTrackingModeDefault
    @test desc.resourceOptions == MTL.MtResourceStorageModePrivate |
                                  MTL.MtResourceCPUCacheModeDefaultCache |
                                  MTL.MtResourceHazardTrackingModeDefault

    # setting resource options should be reflected in individual fields
    @test desc.storageMode == MTL.MtStorageModePrivate
    @test desc.cpuCacheMode == MTL.MtCPUCacheModeDefaultCache
    @test desc.hazardTrackingMode == MTL.MtHazardTrackingModeDefault
end

desc = MtlHeapDescriptor()
desc.size = 0x4000 # TODO: use heapBufferSizeAndAlign
let heap = MtlHeap(dev, desc)
    @test heap.label == nothing

    @test heap.type == desc.type

    @test heap.size == desc.size

    # NOTE: these checks are fragile, as the heap options seems to depend on the requested size

    #=@test=# heap.storageMode == desc.storageMode

    #=@test=# heap.cpuCacheMode == desc.cpuCacheMode

    #=@test=# heap.hazardTrackingMode == desc.hazardTrackingMode

    #=@test=# heap.resourceOptions == desc.resourceOptions
end

end

@testset "buffers" begin

dev = first(devices())

buf = MtlBuffer{Int}(dev, 1)

@test buf.length == 8

# MtlResource properties
@test buf.device == dev
@test buf.label == nothing

@test sizeof(buf) == 8

free(buf)

end

@testset "command queue" begin

dev = first(devices())

cmdq = MtlCommandQueue(dev)

@test cmdq.device == dev
@test cmdq.label == nothing

end

@testset "command buffer" begin

dev = first(devices())
cmdq = MtlCommandQueue(dev)

cmdbuf = MtlCommandBuffer(cmdq)

@test cmdbuf.device == dev
@test cmdbuf.commandQueue == cmdq
@test cmdbuf.label == nothing
@test cmdbuf.error === nothing
@test cmdbuf.status == MTL.MtCommandBufferStatusNotEnqueued
@test cmdbuf.kernelStartTime == 0
@test cmdbuf.kernelEndTime == 0
@test cmdbuf.gpuStartTime == 0
@test cmdbuf.gpuEndTime == 0

let ev = MtlSharedEvent(dev)
    @test ev.signaledValue == 0
    encode_signal!(cmdbuf, ev, 42)
    encode_wait!(cmdbuf, ev, 21)
    commit!(cmdbuf)
    wait_completed(cmdbuf)
    @test ev.signaledValue == 42
end

cmdbuf = MtlCommandBuffer(cmdq)
@test cmdbuf.status == MTL.MtCommandBufferStatusNotEnqueued
enqueue!(cmdbuf)
@test cmdbuf.status == MTL.MtCommandBufferStatusEnqueued
commit!(cmdbuf)
@test cmdbuf.status == MTL.MtCommandBufferStatusCommitted
# Completion happens too quickly to test for committed status to be checked
wait_completed(cmdbuf) == MTL.MtCommandBufferStatusCompleted
@test cmdbuf.status == MTL.MtCommandBufferStatusCompleted

# CommandBufferDescriptor tests
desc = MTL.mtNewCommandBufferDescriptor()
@test MTL.mtCommandBufferDescriptorRetainedReferences(desc) == true
MTL.mtCommandBufferDescriptorRetainedReferencesSet(desc,false)
@test MTL.mtCommandBufferDescriptorRetainedReferences(desc) == false

@test MTL.mtCommandBufferDescriptorErrorOptions(desc) == MTL.MtCommandBufferErrorOptionNone
MTL.mtCommandBufferDescriptorErrorOptionsSet(desc,MTL.MtCommandBufferErrorOptionEncoderExecutionStatus)
@test MTL.mtCommandBufferDescriptorErrorOptions(desc) == MTL.MtCommandBufferErrorOptionEncoderExecutionStatus

cmq = MtlCommandQueue(current_device())
cmdbuf = MtlCommandBuffer(cmq; retainReferences=false, errorOption=MTL.MtCommandBufferErrorOptionEncoderExecutionStatus)
@test cmdbuf.retainedReferences == false
@test cmdbuf.errorOptions == MTL.MtCommandBufferErrorOptionEncoderExecutionStatus

end

@testset "compute pipeline" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "dummy.metallib"))
fun = MtlFunction(lib, "kernel_1")

pipeline = MtlComputePipelineState(dev, fun)

@test pipeline.device == dev
@test pipeline.label === nothing

@test pipeline.maxTotalThreadsPerThreadgroup isa Integer
@test pipeline.threadExecutionWidth isa Integer
@test pipeline.staticThreadgroupMemoryLength == 0

end

@testset "argument encoder" begin

dev = first(devices())
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "vadd.metallib"))
fun = MtlFunction(lib, "vadd")

encoder = MtlArgumentEncoder(fun, 1)

@test encoder.encodedLength == 0
@test encoder.alignment == 1

# TODO: actually encode arguments

end

# TODO: continue adding tests

end


text:
@doc raw"""
    Moments <: AbstractImageBinarizationAlgorithm
    Moments()

    binarize([T,] img, f::Moments)
    binarize!([out,] img, f::Moments)

The following rule determines the binarization threshold:  if one assigns all
observations below the threshold to a value z₀ and all observations above the
threshold to a value z₁, then the first three moments of the original histogram
must match the moments of this specially constructed bilevel histogram.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Details

Let ``f_i`` ``(i=1 \ldots I)`` denote the number of observations in the
``i``th bin of the histogram and ``z_i`` ``(i=1 \ldots I)`` the observed value
associated with the ``i``th bin.  Then the probability that an observation ``z_i``
belongs to the ``i``th bin is given by  ``p_i = \frac{f_i}{N}`` (``i = 1,
\ldots, I``), where ``N = \sum_{i=1}^{I}f_i``.

Moments can be computed from the histogram ``f`` in the following way:

```math
m_k = \frac{1}{N} \sum_i p_i (z_i)^k \quad k = 0,1,2,3, \ldots.
```
The principle of moment-preserving thresholding is to select a threshold value,
as well as two representative values ``z_0`` and ``z_1`` (``z_0 < z_1``),
such that if all below-threshold values in ``f`` are replaced by ``z_0`` and
all above-threshold values are replaced by ``z_1``, then this specially constructed
bilevel histogram ``g`` will have the same first three moments as ``f``.

Concretely, let ``q_0`` and ``q_1`` denote the fractions of observations below
and above the threshold in ``f``, respectively. The constraint that the first
three moments in ``g`` must equal the first three moments in ``f`` can be
expressed by the following system of four equations

```math
\begin{aligned}
   q_0 (z_0)^0 + q_1 (z_1)^0   & = m_0 \\
   q_0 (z_0)^1 + q_1 (z_1)^1   & = m_1 \\
   q_0 (z_0)^2 + q_1 (z_1)^2   & = m_2 \\
   q_0 (z_0)^3 + q_1 (z_1)^3   & = m_3 \\
\end{aligned}
```
where the left-hand side represents the moments of ``g`` and the right-hand side
represents the moments of ``f``. To find the desired treshold value, one first solves
the four equations to obtain ``q_0`` and ``q_1``, and then chooses the threshold
``t`` such that ``q_0 = \sum_{z_i \le t} p_i``.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, Moments())
```

# Reference

[1] W.-H. Tsai, “Moment-preserving thresolding: A new approach,” Computer Vision, Graphics, and Image Processing, vol. 29, no. 3, pp. 377–393, Mar. 1985. [doi:10.1016/0734-189x(85)90133-1](https://doi.org/10.1016/0734-189x%2885%2990133-1)
"""
struct Moments <: AbstractImageBinarizationAlgorithm end

function (f::Moments)(out::GenericGrayImage, img::GenericGrayImage)
    edges, counts = build_histogram(img,  256)
    t = find_threshold(HistogramThresholding.Moments(), counts[1:end], edges)
    for i in CartesianIndices(img)
        out[i] = img[i] < t ? 0 : 1
    end
    out
end

(f::Moments)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
using DelayEmbeddings
using Test, Random, DelimitedFiles

println("\nTesting pecuzal_method.jl...")
@testset "PECUZAL" begin

## Test Lorenz example
# For comparison reasons using Travis CI we carry out the integration on a UNIX
# OS and save the resulting time series
# See https://github.com/JuliaDynamics/JuliaDynamics for the storage of
# the time series used for testing
#
# u0 = [0, 10.0, 0.0]
# lo = Systems.lorenz(u0; σ=10, ρ=28, β=8/3)
# tr = trajectory(lo, 100; Δt = 0.01, Ttr = 100)
tr = readdlm(joinpath(tsfolder, "test_time_series_lorenz_standard_N_10000_multivariate.csv"))
tr = Dataset(tr)

@testset "Univariate example" begin

    s = vec(tr[:, 1]) # input timeseries = x component of lorenz
    w = estimate_delay(s, "mi_min")
    Tmax = 100

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding(s[1:5000];
                                        τs = 0:Tmax , w = w, L_threshold = 0.05, econ=true)

    @test -1.08 < Ls[1] < -1.04
    @test -0.51 < Ls[2] < -0.47
    @test -0.18 < Ls[3] < -0.14

    @test τ_vals[2] == 19
    @test τ_vals[3] == 96
    @test τ_vals[4] == 80
    @test length(ts_vals) == 4

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding(s;
                                        τs = 0:Tmax , w = w, L_threshold = 0.05)
    @test -0.95 < Ls[1] < -0.91
    @test -0.48 < Ls[2] < -0.44

    @test τ_vals[2] == 18
    @test τ_vals[3] == 9
    @test length(ts_vals) == 3

    YY1 = DelayEmbeddings.hcat_lagged_values(Dataset(s), s, 21)
    YY2 = DelayEmbeddings.hcat_lagged_values(YY1, s, 13)
    @test length(YY1) == length(YY2)
    @test YY1 == YY2[:,1:2]

end

@testset "Multivariate example" begin
    ## Test of the proposed PECUZAL-embedding-method (multivariate case)

    w1 = estimate_delay(tr[:,1], "mi_min")
    w2 = estimate_delay(tr[:,2], "mi_min")
    w3 = estimate_delay(tr[:,3], "mi_min")
    w = w1
    Tmax = 100

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding(tr[1:5000,:];
                                        τs = 0:Tmax , w = w, econ = true)

    @test length(ts_vals) == 4
    @test ts_vals[2] == ts_vals[3] == ts_vals[4] == 1
    @test ts_vals[1] == 3
    @test τ_vals[1] == 0
    @test τ_vals[2] == 9
    @test τ_vals[3] == 64
    @test τ_vals[4] == 53
    @test -1.40 < Ls[1] < -1.36
    @test -0.76 < Ls[2] < -0.72
    @test -0.1 < Ls[3] < -0.06

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding(tr[1:5000,:];
                                        τs = 0:Tmax , w = w, econ = true, L_threshold = 0.2)
    @test -1.40 < Ls[1] < -1.36
    @test -0.76 < Ls[2] < -0.72
    @test size(Y,2) == 3

end

@testset "Dummy example" begin
    # Dummy input
    Random.seed!(1234)
    d1 = randn(1000)
    d2 = rand(1000)
    Tmax = 50
    dummy_set = Dataset(hcat(d1,d2))

    w1 = estimate_delay(d1, "mi_min")
    w2 = estimate_delay(d2, "mi_min")
    w = w1

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding(dummy_set;
                                    τs = 0:Tmax , w = w, econ = true)

    @test size(Y,2) == 1
end

end


text:
using ArgParse

s=ArgParseSettings()
@add_arg_table s begin
    "--part2"
        help="Do part 2"
        action=:store_true
    "--test"
        help="Run test data"
        action=:store_true
end
parsed_args = parse_args(ARGS, s)

function parseLine(line)
    parse(Int32, line)
end

function growcounts!(counts)
    n = counts[1]
    for i in 1:8
        counts[i] = counts[i+1]
    end
    counts[7] += n
    counts[9] = n
end

test_suffix = parsed_args["test"] ? "_test" : ""
lines = readlines("input/day_6$test_suffix.txt")
fish = [parse(Int8, f) for f in split(lines[1], ",")]
census = zeros(UInt64, 9)
function countfish(age)
    ageindex = age + 1
    census[ageindex] += 1
end
foreach(countfish, fish)
for _ in 1:(parsed_args["part2"] ? 256 : 80)
    growcounts!(census)
end
answer = sum(census)
println("The answer is: $answer")


text:
"""
AffineDiffEqOperator{T} <: AbstractDiffEqOperator{T}

`Ex: (A₁(t) + ... + Aₙ(t))*u + B₁(t) + ... + Bₙ(t)`

AffineDiffEqOperator{T}(As,Bs,du_cache=nothing)

Takes in two tuples for split Affine DiffEqs

1. update_coefficients! works by updating the coefficients of the component
   operators.
2. Function calls L(u, p, t) and L(du, u, p, t) are fallbacks interpretted in this form.
   This will allow them to work directly in the nonlinear ODE solvers without
   modification.
3. f(du, u, p, t) is only allowed if a du_cache is given
4. B(t) can be Union{Number,AbstractArray}, in which case they are constants.
   Otherwise they are interpreted they are functions v=B(t) and B(v,t)

Solvers will see this operator from integrator.f and can interpret it by
checking the internals of As and Bs. For example, it can check isconstant(As[1])
etc.
"""
struct AffineDiffEqOperator{T,T1,T2,U} <: AbstractDiffEqOperator{T}
    As::T1
    Bs::T2
    du_cache::U
    function AffineDiffEqOperator{T}(As,Bs,du_cache=nothing) where T
        all([size(a) == size(As[1])
             for a in As]) || error("Operator sizes do not agree")
        new{T,typeof(As),typeof(Bs),typeof(du_cache)}(As,Bs,du_cache)
    end
end

Base.size(L::AffineDiffEqOperator) = size(L.As[1])

function (L::AffineDiffEqOperator)(u,p,t::Number)
    update_coefficients!(L,u,p,t)
    du = sum(A*u for A in L.As)
    for B in L.Bs
        if typeof(B) <: Union{Number,AbstractArray}
            du .+= B
        else
            du .+= B(t)
        end
    end
    du
end

function (L::AffineDiffEqOperator)(du,u,p,t::Number)
    update_coefficients!(L,u,p,t)
    L.du_cache === nothing && error("Can only use inplace AffineDiffEqOperator if du_cache is given.")
    du_cache = L.du_cache
    fill!(du,zero(first(du)))
    # TODO: Make type-stable via recursion
    for A in L.As
        mul!(du_cache,A,u)
        du .+= du_cache
    end
    for B in L.Bs
        if typeof(B) <: Union{Number,AbstractArray}
            du .+= B
        else
            B(du_cache,t)
            du .+= du_cache
        end
    end
end

function update_coefficients!(L::AffineDiffEqOperator,u,p,t)
    # TODO: Make type-stable via recursion
    for A in L.As; update_coefficients!(A,u,p,t); end
    for B in L.Bs; update_coefficients!(B,u,p,t); end
end

@deprecate is_constant(L::AbstractDiffEqOperator) isconstant(L)


text:
using SEDF
using Test


text:
@doc raw"""
    gradient_descent(M, F, gradF, x)

perform a gradient descent

```math
x_{k+1} = \operatorname{retr}_{x_k}\bigl( s_k\operatorname{grad}f(x_k) \bigr)
```

with different choices of ``s_k`` available (see `stepsize` option below).

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F: \mathcal M→ℝ`` to minimize
* `gradF` – the gradient ``\operatorname{grad}F: \mathcal M → T\mathcal M`` of F
* `x` – an initial value ``x ∈ \mathcal M``

# Optional
* `evaluation` – ([`AllocatingEvaluation`](@ref)) specify whether the gradient works by allocation (default) form `gradF(M, x)`
  or [`MutatingEvaluation`](@ref) in place, i.e. is of the form `gradF!(M, X, x)`.
* `retraction_method` – (`default_retraction_method(M)`) a `retraction(M,x,ξ)` to use.
* `return_options` – (`false`) – if activated, the extended result, i.e. the
    complete [`Options`](@ref) are returned. This can be used to access recorded values.
    If set to false (default) just the optimal value `x_opt` if returned
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`) specify a [`Stepsize`](@ref)
  functor.
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(200), `[`StopWhenGradientNormLess`](@ref)`(10.0^-8))`)
  a functor inheriting from [`StoppingCriterion`](@ref) indicating when to stop.
...
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `x_opt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function gradient_descent(
    M::AbstractManifold, F::TF, gradF::TDF, x; kwargs...
) where {TF,TDF}
    x_res = allocate(x)
    copyto!(M, x_res, x)
    return gradient_descent!(M, F, gradF, x_res; kwargs...)
end
@doc raw"""
    gradient_descent!(M, F, gradF, x)

perform a gradient_descent ``x_{k+1} = \mathrm{retr}_{x_k} s_k\operatorname{grad}f(x_k)`` in place of `x`
with different choices of ``s_k`` available.

# Input
* `M` – a manifold ``\mathcal M``
* `F` – a cost function ``F:\mathcal M→ℝ`` to minimize
* `gradF` – the gradient ``\operatorname{grad}F:\mathcal M→ T\mathcal M`` of F
* `x` – an initial value ``x ∈ \mathcal M``

For more options, especially [`Stepsize`](@ref)s for ``s_k``, see [`gradient_descent`](@ref)
"""
function gradient_descent!(
    M::AbstractManifold,
    F::TF,
    gradF::TDF,
    x;
    stepsize::Stepsize=ConstantStepsize(1.0),
    retraction_method::AbstractRetractionMethod=default_retraction_method(M),
    stopping_criterion::StoppingCriterion=StopAfterIteration(200) |
                                          StopWhenGradientNormLess(10.0^-8),
    direction=IdentityUpdateRule(),
    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
    return_options=false,
    kwargs..., #collect rest
) where {TF,TDF}
    p = GradientProblem(M, F, gradF; evaluation=evaluation)
    o = GradientDescentOptions(
        M,
        x;
        stopping_criterion=stopping_criterion,
        stepsize=stepsize,
        direction=direction,
        retraction_method=retraction_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
#
# Solver functions
#
function initialize_solver!(p::GradientProblem, o::GradientDescentOptions)
    o.gradient = get_gradient(p, o.x)
    return o
end
function step_solver!(p::GradientProblem, o::GradientDescentOptions, iter)
    s, o.gradient = o.direction(p, o, iter)
    retract!(p.M, o.x, o.x, -s * o.gradient, o.retraction_method)
    return o
end
get_solver_result(o::GradientDescentOptions) = o.x


text:
"""
    ngramfirstcombine(x::Strings)

A function that takes a Vector of strings and outputs a vector of string vectors with ngram combinations of the
    first member with all the other elements of the vector. It is the first step for implementing the lexical gravity association measure.
    For more details, cf. (Daudaravičius & Marcinkevičienė 2004).

# Examples
```julia-repl
julia> ngramfirstcombine(["σε", "νόμιμες", "ενέργειες", "κατά"], 4)
3-element Vector{Vector{String}}:
["σε", "νόμιμες"]
["σε", "ενέργειες"]
["σε", "κατά"]
```
"""
@memoize function ngramfirstcombine(str::Strings)
ngramArray::Vector{Vector{String}} = []
for i in 1:(length(str) - 1)
    push!(ngramArray, [str[1], str[1 + i]])
end
return ngramArray
end

@memoize function ngramfirstcombine(str::Strings)
    println("am fine")
    ngramArray = Dictionary{Vector{String}, Float64}()
    str = tokenize(str)
    #println(str)
    for i in 1:(length(str) - 1)
        insert!(ngramArray, [str[1], str[1 + i]], 1/i)
        #println(ngramArray)
    end
    return ngramArray
end

"""
    tokenize(s::StringDocument)

A function that takes a StringDocument and tokenizes it using using the tokenize() function created in the TextAnalysis package.

# Examples
```julia-repl
julia> nicestringdoc = StringDocument("This is a nice string.")
julia> tokenize(nicestringdoc)
6-element Vector{String}:
 "This"
 "is"
 "a"
 "nice"
 "string"
 "."
```
"""
tokenize(s::StringDocument) =  tokenize(text(s))

#ngram_simple(s::Strings_simple1, n::Integer) = [s[i:i+n-1] for i=1:length(s)-n+1]
# s must be first tokenized
ngram(s::Strings, n::Integer) = [s[i:i+n-1] for i=1:length(s)-n+1]

"""
    ngramdirection(x::Strings)

A function that takes a Vector of strings and outputs a vector of string vectors with ngram combinations of the
first member with all the other elements of the vector. It is the first step for implementing the lexical gravity association measure.
For more details, cf. (Daudaravičius & Marcinkevičienė 2004).

# Examples
```julia-repl
julia> ngramfirstcombine(["σε", "νόμιμες", "ενέργειες", "κατά"], 4)
3-element Vector{Vector{String}}:
["σε", "νόμιμες"]
["σε", "ενέργειες"]
["σε", "κατά"]
```
"""
@memoize function ngramdirection(mydata::Strings; n = 4, direction = "forward")
if direction == "forward"
    str = ngram(mydata, n)
    ngramfirstcombine.(str)
elseif  direction == "backward"
    str = reverse.(ngram(mydata, n))
    ngramfirstcombine.(str)
end
end

# tokenizedtext

"""
    bigramsinside(x::Strings)

A function that takes a Vector of strings and outputs a vector of string vectors with ngram combinations of the
first member with all the other elements of the vector. It is the first step for implementing the lexical gravity association measure.
For more details, cf. (Daudaravičius & Marcinkevičienė 2004).

# Examples
```julia-repl
julia> bigramsinside(["σε", "νόμιμες", "ενέργειες", "κατά"], 4)
3-element Vector{Vector{String}}:
["σε", "νόμιμες"]
["σε", "ενέργειες"]
["σε", "κατά"]
```
"""
@memoize function bigramsinside(word, text::Strings, ngram)
typeof(text) == Vector{String} || (text = tokenize(text))
wordbigrams = filter(x -> x[1][1] == word, ngram)
return wordbigrams
end