text:
"""
Compute B-Spline basis
    bs(x, df, norder, nderiv=0)
    :param x: time samples
    :param df: degree of freedom
    :param norder: order of splines
    :param nderiv: derivative number
"""
function bs(x::Vector, df, norder, nderiv=0)
    n = length(x);

    nbreaks = df - norder + 2;
    breaks = LinRange(minimum(x), maximum(x), nbreaks);
    # Check Inputs
    if df <= 0
        error("Degrees of freedom <= 0.")
    end

    if norder <= 0
        error("order of splines <= 0.")
    end

    if minimum(diff(x)) < 0
        isrt = sortperm(x);
        sort!(x);
        sortwrd = true;
    else
        sortwrd = false;
    end

    if x[1] - breaks[1] < -1e-10 || x[n] - breaks[nbreaks] > 1e-10
        println([x[1], x[n]])
        error("Argument values out of range")
    end

    # set some abbreviations
    k = norder;  # order of splines
    km1 = k-1;
    nb = length(breaks);  # number of break points
    nx = length(x);  # number of argument values
    nd = nderiv+1;  # ND is order of derivative plus one
    ns = nb-2+k;  # number of splines to compute
    if ns < 1
        println("There are no B-splines for the given input")
        bsplinemat = [];
        return bsplinemat
    end
    onenx = ones(nx);
    onenb = ones(k);
    onens = ones(ns);

    # augment break sequence to get knots by adding a  K-fold knot at
    # each end.

    knots = [breaks[1].*ones(km1); breaks; breaks[nb].*ones(km1)];
    nbasis = length(knots) - k;

    # For each i, determine left(i) so that K<= left(i) < nbasis+1, and
    # within that restriction, knots(left(i)) <= pts(i) < knot(left(i)+1)

    knotslower = knots[1:nbasis];
    index = sortperm([knotslower; x]);
    point = findall(index .> nbasis) - collect(1:length(x));
    left = maximum(hcat(point, k*onenx),dims=2);

    # compute bspline values and derivatives, if needed:

    # initialize the b array.
    temp = transpose([1; zeros(km1)]);
    b = repeat(temp, nd*nx, 1);
    nxs = nd.*collect(1:nx);

    # run the recurrence simultaneously for all x(i)

    # First, bring it upt to the intended level:

    for j in 1:(k-nd)
        saved = zeros(nx);
        for r in 1:j
            leftpr = round.(Integer, left .+ r);
            tr = knots[leftpr] - x;
            tl = x .- knots[leftpr.-j];
            term = b[nxs,r]./(tr+tl);
            b[nxs,r] = saved + tr.*term;
            saved = tl.*term;
        end
        b[nxs,j+1] = saved;
    end

    # save the B-spline values in successive blocks in b.

    for jj in 1:(nd-1)
        j = k - nd + jj;
        saved = zeros(nx);
        nxn = nxs - 1;
        for r in 1:j
            leftpr = left+r;
            tr = knots[leftpr] - x;
            tl = x-knots[leftpr-j];
            term = b[nxs,r]./(tr+tl);
            b[nxn,r] = saved + tr.*term;
            saved = tl.*term;
        end
        b[nxn, j+1] = saved;
        nxs = nxn;
    end

    # now use the fact that derivative values can be obtained by differencing

    for jj in (nd-1):-1:1
        j = k - jj;
        temp = collect(jj:(nd-1)).*onenx + ones(nd-jj)*nxn;
        nxs = reshape(temp, (nd-1-jj+1)*nx, 1);
        for r in j:-1:1
            leftpr = left + r;
            temp = ones(nd-jj) * (knots[leftpr] - knots[leftpr-j])/j;
            b[nxs,r] = -1*b[nxs,r]./temp;
            b[nxs,r+1] = b[nxs,r+1] - b[nxs,r];
        end
    end

    # Finally, zero out all rows of b corresponding to x outside the basic
    # interval, [breaks[1] ... breaks[nb]]

    index = findall((x.<breaks[1]) .| (x .> breaks[nb]));
    if !isempty(index)
        temp = [(1-nd):0].*ones(length(index))+nd*ones(nd).*index;
        b[temp, :] = zeros(nd*length(index),k);
    end

    # set up output matrix bsplinemat
    width = maximum([ns, nbasis]) + km1 + km1;
    cc = zeros(nx*width);
    index = collect((1-nx):0)*transpose(onenb) + nx * (left*transpose(onenb) + onenx*transpose(collect(-km1:0)));
    index = round.(Integer, index);
    cc[index] = b[nd*collect(1:nx), :];
    # (This uses the fact that for a column vector v and a matrix A,
    #  v(A)(i,j) = v(A(i,j)), all i, j.)
    index2 = round.(Integer, collect((1-nx):0)*transpose(onens) + nx * onenx*transpose(collect(1:ns)));
    bsplinemat = reshape(cc[index2],nx,ns);

    if sortwrd
        temp = copy(bsplinemat);
        bsplinemat[isrt,:] = temp;
    end

    return bsplinemat
end


text:
# NOTE distance in potentials is ångströms when fitting.
# But bohrs when CP2K calculates them. So that is why there are
# 0.52917721090 -terms in few places.

"""
GeneralPowers

General potential with customizable powers.

# Fields
- `constants::Vector{Float64}` : potential constants
- `powers::Vector{Int}`        : powers for radius
"""
mutable struct GeneralPowers <: AbstractPairPotential
    constants::Vector{Float64}
    powers::Vector{Int}
    GeneralPowers() = new([0.0, 0.0], [-6, -12])
    GeneralPowers(i...) = new(zeros(length(i)), collect(i))
    GeneralPowers(ur::UnitRange{Integer}) = new(zeros(length(ur)), collect(ur))
end

function Base.show(io::IO, potential::GeneralPowers; energy_unit="hartree")
    for (c,p) in zip(potential.constants, potential.powers)
        print(io, "r^$(p)  ")
    end
    print(io, " *** ")
    for (c,p) in zip(potential.constants, potential.powers)
        print(io, "$(c)  ")
    end
end

function potential_variables(potential::GeneralPowers,
                            cluster::AbstractCluster, indices)
    # Convert from Å to bohr
    r = distances(cluster, indices[1], indices[2]) ./ 0.52917721090
    return [r.^x for x in potential.powers]'
end

function get_potential!(potential::GeneralPowers, constants...)
    @assert length(constants) == length(potential.constants)  "length of constants with GeneralJones potential do not match"
    return potential.constants .= constants
end

function (p::GeneralPowers)(r::Number)::Float64
    y = 0.0
    rr = r/0.52917721090 # r is in Å and potential is in bohr
    for (a,b) in zip(p.constants, p.powers)
        y += a*rr^b
    end
    return y
end

ext: jl


text:
export InsecureDebug

struct InsecureDebug{P<:SHEShemeParams} <: PassthroughParams{P}
    params::P
end

scheme_name(::Type{InsecureDebug{P}}) where P = string("Insecure ", scheme_name(P))

struct ZeroSampler <: Random.Sampler{Zero}; end
Random.rand(rng::AbstractRNG, zs::ZeroSampler) = Zero()
Base.:+(a::NTT.RingElement, b::Zero) = a
Base.:+(a::Zero, b::NTT.RingElement) = a

𝒩(params::InsecureDebug) = ZeroSampler()


text:
@testset "Functions" begin

    @testset "funcdef" begin
        ex1 = :(f(x) = 1)
        ex2 = :(@inline f(x) = 1)
        ex3 = :(Base.@pure @inline f(x) = 1)
        ex4 = quote
            f(x) = 1
        end
        ex5 = :(function f(x) 1 end)
        ex6 = quote
            function f(x)
                1
            end
        end
        ex7 = quote
            @inline function f(x)
                1
            end
        end
        ex8 = :(immutable NotAFunction end)

        for ex in (ex1, ex5)
            @test is_funcdef_expr(ex)
        end
        for ex in (ex2, ex3, ex4, ex6, ex7, ex8)
            @test !is_funcdef_expr(ex)
        end
        for (ex,contains) in ((ex1,false), (ex2,true), (ex3,true), (ex4,true),
                              (ex5,false), (ex6,true), (ex7,true))
            def, c, container = get_funcdef_expr(ex)
            @test c == contains
            @test is_funcdef_expr(def)
            @test isa(container, Expr)
            if contains
                @test container.args[end] === def
            end
            lf = funcdef_longform(def)
            @test lf.head == :function
        end
        @test_throws ArgumentError get_funcdef_expr(ex8)
        def, c, container = get_funcdef_expr(ex8, false)
        @test isa(def, Expr)
        @test !is_funcdef_expr(def)
    end

    @testset "signature" begin
        # These are all canonicalized, which is OK because we've tested funcdef_longform
        ex1 = :(function f(x) 1 end)
        ex2 = :(function f(x::Int) 1 end)
        ex3 = :(function f{T}(x::T) 1 end)
        ex4 = :(function f{R,S}(x::R, y::S, z) 1 end)
        ex5 = :(function f{T<:Integer}(x::T, y::String) 1 end)
        ex6 = :(function f{T}(x::AbstractArray{T}) 1 end)
        for (ex, params, args, names, types) in ((ex1, Symbol[], [:x], [:x], [:Any]),
                                                 (ex2, Symbol[], [:(x::Int)], [:x], [:Int]),
                                                 (ex3, [:T], [:(x::T)], [:x], [:T]),
                                                 (ex4, [:R,:S], [:(x::R), :(y::S), :z], [:x, :y, :z], [:R, :S, :Any]),
                                                 (ex5, [Expr(:(<:), :T, :Integer)], [:(x::T), :(y::String)], [:x, :y], [:T, :String]),
                                                 (ex6, [:T], [:(x::AbstractArray{T})], [:x], [:(AbstractArray{T})]))
            sig = ex.args[1]
            @test funcdef_name(sig) == :f
            @test funcdef_params(sig) == params
            @test funcdef_args(sig) == args
            @test funcdef_argnames(sig) == names
            @test funcdef_argtypeexprs(sig) == types
        end
    end

end


text:
using IterableTables
using StatsModels
using GLM
using DataValues
using Test

@testset "StatsModels" begin

source_array = [(a=1,b=1.,c="A"), (a=2,b=2.,c="B"), (a=3,b=3.,c="C")]

mf_array = StatsModels.ModelFrame(StatsModels.@formula(a~b), source_array)

@test mf_array isa StatsModels.ModelFrame

x = lm(StatsModels.@formula(a~b), source_array)

end


text:
const PowerMeasure{M,N,R} = ProductMeasure{Kernel{Returns{M}, typeof(identity)}, CartesianIndices{N, R}}

Base.:^(μ::AbstractMeasure, n::Integer) = powermeasure(Returns(μ), (n,))

Base.:^(μ::AbstractMeasure, dim::Tuple) = powermeasure(Returns(μ), dim)

marginals(pm::PowerMeasure) = Fill(pm.κ.value, size(pm))

basemeasure(pm::PowerMeasure) = basemeasure(pm.κ.value) ^ size(pm)

function logdensity_def(μ::PowerMeasure, x::AbstractArray)
    sum(x) do xⱼ
        logdensity_def(μ.κ.value, xⱼ)
    end
end

# function Base.rand(rng::AbstractRNG, ::Type{T}, pm::PowerMeasure) where {T} end
#     x₁ = rand(rng, T, pm.κ.value)

#     # We need an array to pass to `similar`. `Fill` is nice here because it doesn't allocate
#     # TODO: Generalize this so the container type can be chosen at construction time
#     result = similar(Fill(x₁, 1), pm.xs)
#     indices = eachindex(x₁)
#     next = iterate(indices)

#     isnothing(next) || result[first(next)] = x₁

#     while next !== nothing
#         (j, state) = next
#         result[j] = rand(rng, T, pm.parent)
#         next = iterate(indices, state)
#     end

#     return result
# end


text:
using BanditOpt
using Test

macro test_nothrow(ex)
    quote
        try
            $(esc(ex))
            true
        catch e
            print( "ERROR: " )
            if isa( e, KeyError )
                println( "KeyError: key ", e.key, " not found" )
            else
                println( e )
            end
            false
        end
    end
end

tests = [
    # Agent specific tests
    "Agents/EpsGreedy",
    # Arm specific tests
    "Arms/StationaryArms"
    # Experiment specific tests
];

#=
Use
    ARGS = [ "sanity_check" ]; include( joinpath(Pkg.dir("MAB"),"test","runtests.jl") )
or
    ARGS = [ "Algorithms/epsGreedy" ]; include( joinpath(Pkg.dir("MAB"),"test","runtests.jl") )
from REPL.
=#
if length(ARGS) > 0
    tests = ARGS
end

@testset "BanditOpt Tests" begin
for test_script in tests
    fp = joinpath( dirname(@__FILE__), "$(test_script).jl" )
    # println( "Testing : ", test_script )
    # @time include( fp )
    include( fp )
end
end

nothing


text:
import Media
import Juno: @render, Tree, SubTree, Inline, Model,  render
import Atom: span, c, UNDEF
getfield_(x, f) = isdefined(x, f) ? getfield(x, f) : UNDEF

@render i::Inline x::InputData begin
    Tree(x.name, [SubTree(Text("$f → "), getfield(x, f)) for f in fieldnames(typeof(x))[2:end]])
end

@render i::Inline x::DFInput begin
    runstring = x.run ? span(".syntax--constant.syntax--other.syntax--symbol", name(x)) : span(".syntax--comment", name(x))
    fields = filter(x->x != run, [fieldnames(typeof(x))[2:end]...])
    Tree(runstring , [SubTree(Text("$f → "), getfield(x, f)) for f in fields])
end

@render i::Inline x::AbstractAtom begin
    pos = position_cart(x)
    Tree(Text("$(name(x))  ang: $(round(pos[1], digits=5)) $(round(pos[2], digits=5)) $(round(pos[3], digits=5))"), [SubTree(Text("$f → "), getfield_(x, f)) for f in fieldnames(typeof(x))[2:end]])
end

@render i::Inline x::DFJob begin
    Tree(x.name, [SubTree(Text("$f → "), getfield(x, f)) for f in fieldnames(typeof(x))[3:end]])
end

@render i::Inline x::AbstractStructure begin
    Tree(name(x) , [SubTree(Text("$f → "), getfield(x, f)) for f in fieldnames(typeof(x))])
end


text:
export Chebyshev

function Chebyshev(nlp_stop,
				   PData :: PDataFact,
				   dₙ :: Vector,
				   gt :: Vector)
	# printstyled("on est dans Chebyshev  \n", color = :yellow)
	T = eltype(gt)
	# @show T
	nlp_at_x = nlp_stop.current_state
	x = nlp_at_x.x
	# xdemi = NaN * rand(length(dₙ))
	# @show x
	# @show ∇f³xuv
	# @show typeof(nlp_stop.pb)
	Tdndn = ∇f³xuv(nlp_stop.pb, x, dₙ, dₙ)
	# printstyled("on a Tdndn  \n", color = :yellow)

	ϵ2 =  sqrt(eps(eltype(gt)))
   	Γ = max.(abs.(PData.Δ), ϵ2)

	d̃ = PData.L \ Tdndn[PData.pp]
    d̂ = PData.L' \ (PData.Q * (PData.Q' * d̃ ./ Γ))
    d = -d̂[invperm(PData.pp)]

	dC = dₙ + T(0.5) * d
	return dC
end


text:
using InplaceLinalg
using Test, LinearAlgebra

@test @macroexpand(@inplace C = A * B)             == :(InplaceLinalg.C_AB!(C, 0, +, 1, A, B))
@test @macroexpand(@inplace C = α * A * B)         == :(InplaceLinalg.C_AB!(C, 0, +, α, A, B))
#@test @macroexpand(@inplace C = A * 7B)            == :(InplaceLinalg.C_AB!(C, 0, +, 7, A, B))
#@test @macroexpand(@inplace C = 3A * B)            == :(InplaceLinalg.C_AB!(C, 0, +, 3, A, B))
#@test @macroexpand(@inplace C = 3A * 7B)           == :(InplaceLinalg.C_AB!(C, 0, +, 3 * 7, A, B))

@test @macroexpand(@inplace C = C + α * A * B)     == :(InplaceLinalg.C_AB!(C, 1, +, α, A, B))
@test @macroexpand(@inplace C = β * C + α * A * B) == :(InplaceLinalg.C_AB!(C, β, +, α, A, B))
@test @macroexpand(@inplace C = β * C + 2π * A * B) == :(InplaceLinalg.C_AB!(C, β, +, 2π, A, B))
@test @macroexpand(@inplace C = β * C + 2 * π * A * B) == :(InplaceLinalg.C_AB!(C, β, +, 2π, A, B))
@test @macroexpand(@inplace C = β * C + A * B)     == :(InplaceLinalg.C_AB!(C, β, +, 1, A, B))
@test @macroexpand(@inplace C = 2π * C + A * B)    == :(InplaceLinalg.C_AB!(C, 2π, +, 1, A, B))
@test @macroexpand(@inplace C = 2*π * C + A * B)   == :(InplaceLinalg.C_AB!(C, 2π, +, 1, A, B))

@test @macroexpand(@inplace C += A * B)            == :(InplaceLinalg.C_AB!(C, 1, +, 1, A, B))
@test @macroexpand(@inplace C += α * A * B)        == :(InplaceLinalg.C_AB!(C, 1, +, α, A, B))
@test @macroexpand(@inplace C -= A * B)            == :(InplaceLinalg.C_AB!(C, 1, -, 1, A, B))
@test @macroexpand(@inplace C -= α * A * B)        == :(InplaceLinalg.C_AB!(C, 1, -, α, A, B))

@test @macroexpand(@inplace C += 2C + A * B)       == :(InplaceLinalg.C_AB!(C, 3, +, 1, A, B))
@test @macroexpand(@inplace C -= 2C + A * B)       == :(InplaceLinalg.C_AB!(C, -1, -, 1, A, B))
@test @macroexpand(@inplace C += 2.5C + π * A * B) == :(InplaceLinalg.C_AB!(C, 3.5, +, π, A, B))
@test @macroexpand(@inplace C -= 2.5C + 2exp(1) * A * B) == :(InplaceLinalg.C_AB!(C, -1.5, -, 2 * exp(1), A, B))

@test @macroexpand(@inplace C += A)                == :(InplaceLinalg.C_AB!(C, 1, +, 1, 1, A))
@test @macroexpand(@inplace C = C + A)             == :(InplaceLinalg.C_AB!(C, 1, +, 1, 1, A))
@test @macroexpand(@inplace C = 0.1C + 0.2A)       == :(InplaceLinalg.C_AB!(C, 0.1, +, 1, 0.2, A))

@test @macroexpand(@inplace B *= A)                == :(InplaceLinalg.mult_update!(B, A, Val(1)))

@test @macroexpand(@inplace B /= A)                == :(InplaceLinalg.div_update!(B, /, A))

@test @macroexpand(@inplace B = B / A)             == :(InplaceLinalg.div_update!(B, /, A))
@test @macroexpand(@inplace B = α * B / A)         == :(InplaceLinalg.div_update!(B, α, /, A))
@test @macroexpand(@inplace B = 2π * B / A)        == :(InplaceLinalg.div_update!(B, 2π, /, A))
@test @macroexpand(@inplace B = (2*π) * B / A)     == :(InplaceLinalg.div_update!(B, 2π, /, A))
@test @macroexpand(@inplace B = A \ B)             == :(InplaceLinalg.div_update!(B, \, A))
@test @macroexpand(@inplace B = A \ 2B)            == :(InplaceLinalg.div_update!(B, 2, \, A))
@test @macroexpand(@inplace B = A \ (α * B))       == :(InplaceLinalg.div_update!(B, α, \, A))
@test @macroexpand(@inplace B = A \ (2π * B))      == :(InplaceLinalg.div_update!(B, 2π, \, A))

@test @macroexpand(@inplace B = B * a * b)         == :(InplaceLinalg.mult_update!(B, a, b, Val(1)))
@test @macroexpand(@inplace B = a * B * b)         == :(InplaceLinalg.mult_update!(B, a, b, Val(2)))
@test @macroexpand(@inplace B = a * b * B)         == :(InplaceLinalg.mult_update!(B, a, b, Val(3)))
@test @macroexpand(@inplace B = B * a)             == :(InplaceLinalg.mult_update!(B, a, Val(1)))
@test @macroexpand(@inplace B = a * B)             == :(InplaceLinalg.mult_update!(B, a, Val(2)))

## esoteric cases
@test @macroexpand(@inplace B = A*-B)              == :(InplaceLinalg.mult_update!(B, -1, A, Val(3)))
@test @macroexpand(@inplace B = -A*B)              == :(InplaceLinalg.mult_update!(B, -1, A, Val(3)))
@test @macroexpand(@inplace B = -2A*B)             == :(InplaceLinalg.mult_update!(B, -2, A, Val(3)))
@test @macroexpand(@inplace B = 2A*B)              == :(InplaceLinalg.mult_update!(B, 2, A, Val(3)))

@test @macroexpand(@inplace B = -B*A)              == :(InplaceLinalg.mult_update!(B, -1, A, Val(2)))
@test @macroexpand(@inplace B = B*-A)              == :(InplaceLinalg.mult_update!(B, -1, A, Val(2)))
@test @macroexpand(@inplace B = B*-2A)             == :(InplaceLinalg.mult_update!(B, -2, A, Val(1)))

#@test_throws InplaceException try @eval @macroexpand(@inplace C += A \ B) catch err; throw(err.error) end
#@test_throws InplaceException try @eval @inplace(C += B / A) catch err; throw(err.error) end

#@test @macroexpand(@inplace C = B / A)             == :((InplaceLinalg.ip_error)("LHS must be equal to numerator in updating divide"))
#@test @macroexpand(@inplace C = 2B / A)            == :(InplaceLinalg.C_div!(C, 2, B, $(/), A))

A = randn(100, 200)
B = randn(200, 100)

C = zeros(100, 100)
C1 = A * B

## Basic matrix multiplication, gemm!()
@time @inplace C = A * B
@test C ≈ C1

@time @inplace C += A * B
@test C ≈ 2C1

## With scaling factor
@time @inplace C = 3.0 * A * B
@test C ≈ 3C1

@time @inplace C += 5 * A * B
@test C ≈ 8C1

@time @inplace C += 2 * 4.0f0 * A * B
@test C ≈ 16C1

## Transpose etc.

Bt = collect(B')
At = collect(A')

@time @inplace C = At' * B
@test C ≈ C1

@time @inplace C += At' * B
@test C ≈ 2C1

@time @inplace C = A * Bt'
@test C ≈ C1

@time @inplace C += A * Bt'
@test C ≈ 2C1

@time @inplace C = At' * Bt'
@test C ≈ C1

@time @inplace C += At' * Bt'
@test C ≈ 2C1

## scaling scal!()
#@time @inplace C *= 2
#@test C ≈ 4C1

#@time @inplace C *= π
#@test C ≈ 4π * C1

## copy blascopy!()
#@inplace C *= 0
#@time @inplace C = C1
#@test C == C1

## simple in-place add sxpy!()
#@time @inplace C += C1
#@test C == 2C1

# TRSM =========================================================================
m,n = 10,20
B0 = randn(m,n)
AL = LowerTriangular(randn(m,m))
AR = UpperTriangular(randn(n,n))
AL1 = UnitLowerTriangular(randn(m,m))
AR1 = UnitUpperTriangular(randn(n,n))
DL = Diagonal(randn(m))
DR = Diagonal(randn(n))
rI = UniformScaling(randn())
α = randn()

# test basics: \,/; Upper and Lower triangles

B = copy(B0); @inplace B = AL \ B
@test B ≈ AL \ B0

B = copy(B0); @inplace B = B / AR
@test B ≈ B0 / AR

#test scaling
B = copy(B0); @inplace B = AL \ (α*B)
@test B ≈ AL \ (α*B0)

B = copy(B0); @inplace B = AL \ (B*α)
@test B ≈ AL \ (α*B0)

B = copy(B0); @inplace B = AL \ 2B
@test B ≈ AL \ 2B0

B = copy(B0); @inplace B = AL \ (2*B)
@test B ≈ AL \ 2B0

B = copy(B0); @inplace B = α*B / AR
@test B ≈ α*B0 / AR

B = copy(B0); @inplace B = B*α / AR
@test B ≈ α*B0 / AR

B = copy(B0); @inplace B = 2B / AR
@test B ≈ 2B0 / AR

B = copy(B0); @inplace B = 2*B / AR
@test B ≈ 2B0 / AR

#test unit-diagonal triangles and UniformScaling
B = copy(B0); @inplace B = AL1 \ B
@test B ≈ AL1 \ B0

B = copy(B0); @inplace B = rI \ B
@test B ≈ rI \ B0

B = copy(B0); @inplace B = B / AR1
@test B ≈ B0 / AR1

B = copy(B0); @inplace B = B / rI
@test B ≈ B0 / rI

#test Diagonal
B = copy(B0); @inplace B = DL \ B
@test B ≈ DL \ B0

B = copy(B0); @inplace B = B / DR
@test B ≈ B0 / DR

#Diagonal solve with prescaling not allowed
B = copy(B0);
@test_throws InplaceException @inplace B = DL \ 2B
@test_throws InplaceException @inplace B = 2B / DR

#test /=
B = copy(B0); @inplace B /= AR
@test B ≈ B0 / AR

#test some disallowed stuff
B = copy(B0)

@test_throws InplaceException @inplace C += AL \ B
@test_throws InplaceException @inplace C = C + AL \ B
@test_throws InplaceException @inplace C += B / AR
@test_throws InplaceException @inplace C = C + B / AR

Am = randn(m,m)
An = randn(n,n)
@test_throws InplaceException @inplace B = Am*B
@test_throws InplaceException @inplace B = B*An
@test_throws InplaceException @inplace C = AL / B

@test_throws InplaceException @inplace C = A*B + C*D

# TRSV =========================================================================
m = 10
B0 = randn(m)
AL = LowerTriangular(randn(m,m))
AU = UpperTriangular(randn(m,m))
AL1 = UnitLowerTriangular(randn(m,m))
AU1 = UnitUpperTriangular(randn(m,m))
rI = UniformScaling(randn())
α = randn()

# test basics: triangle variants

B = copy(B0); @inplace B = AL \ B
@test B ≈ AL \ B0

B = copy(B0); @inplace B = AU \ B
@test B ≈ AU \ B0

B = copy(B0); @inplace B = AL1 \ B
@test B ≈ AL1 \ B0

B = copy(B0); @inplace B = AU1 \ B
@test B ≈ AU1 \ B0

B = copy(B0); @inplace B = rI \ B
@test B ≈ rI \ B0

#scaling is disallowed
@test_throws Exception @inplace B = AL \ α*B
#@test_throws Exception @inplace B = AL \ B*α
@test_throws InplaceException @inplace B = AL \ 2B
@test_throws Exception @inplace B = AL \ 2*B

# this too ...
@test_throws InplaceException @inplace B = B \ AL

# AXPY!, AXPBY! and copyto!
α = randn()
β = randn()
m = 3
x = randn(m)
y = randn(m)
y0 = copy(y)

y = copy(y0)
@inplace y += x
@test y ≈ x + y0

y = copy(y0)
@inplace y += 2x
@test y ≈ 2x + y0

y = copy(y0)
@inplace y += α*x
@test y ≈ α*x + y0

y = copy(y0)
@inplace y -= α*x
@test y ≈ y0 - α*x

y = copy(y0)
@inplace y = α*x
@test y ≈ α*x

y = copy(y0)
@inplace y = β*y + α*x
@test y ≈ β*y0 + α*x

y = copy(y0)
@inplace y += β*y + α*x
@test y ≈ (1+β)*y0 + α*x

y = copy(y0)
@inplace y -= β*y + α*x
@test y ≈ (1-β)*y0 - α*x

y = copy(y0)
@inplace y = β*y + x
@test y ≈ β*y0 + x

y = copy(y0)
@inplace y = β*y - α*x
@test y ≈ β*y0 - α*x

y = copy(y0)
@inplace y = β*y - x
@test y ≈ β*y0 - x

y = copy(y0)
@inplace y = 2x
@test y ≈ 2x

y = copy(y0)
@inplace y = x
@test y == x
@test !(y === x)

# symmetric updates
α = randn()
β = randn()
m = 3
x = Symmetric(randn(m,m),:L)
y = Symmetric(randn(m,m),:L)
# x = LowerTriangular(randn(m,m))
# y = LowerTriangular(randn(m,m))
y0 = copy(y)

y = copy(y0)
@inplace y += x
@test y ≈ x + y0

y = copy(y0)
@inplace y += 2x
@test y ≈ 2x + y0

y = copy(y0)
@inplace y += α*x
@test y ≈ α*x + y0

y = copy(y0)
@inplace y -= α*x
@test y ≈ y0 - α*x

y = copy(y0)
@inplace y = α*x
@test y ≈ α*x

y = copy(y0)
@inplace y = β*y + α*x
@test y ≈ β*y0 + α*x

y = copy(y0)
@inplace y += β*y + α*x
@test y ≈ (1+β)*y0 + α*x

y = copy(y0)
@inplace y -= β*y + α*x
@test y ≈ (1-β)*y0 - α*x

y = copy(y0)
@inplace y = β*y + x
@test y ≈ β*y0 + x

y = copy(y0)
@inplace y = β*y - α*x
@test y ≈ β*y0 - α*x

y = copy(y0)
@inplace y = β*y - x
@test y ≈ β*y0 - x

y = copy(y0)
@inplace y = 2x
@test y ≈ 2x

y = copy(y0)
@inplace y = x
@test y == x
@test !(y === x)

#test scalar multiplication and division
B = copy(B0); @inplace B *= 2
@test B ≈ 2*B0

B = copy(B0); @inplace B = B*2
@test B ≈ 2*B0

B = copy(B0); @inplace B = 2*B
@test B ≈ 2*B0

B = copy(B0); @inplace B /= 2
@test B ≈ B0/2

B = copy(B0); @inplace B = B/2
@test B ≈ B0/2

B = copy(B0); @inplace B = 2\B
@test B ≈ B0/2

include("triangle_tests.jl")


text:
module Immortal
    # Подключаем библиотеку, связанную с временем
    include("./utils/Utils.jl")
    import .Utils

    # Подключаем библиотеку для работы с RabbitMQ и каналами
    include("./rabbitmq/Rabbit.jl")
    import .Rabbit

    export init, declares, get_adapter

    """ Функция получения адаптера для сервера очередей

    value: словарь настроек брокера или имя адаптера сервера очередей
    """
    function get_adapter(value)
        @info Utils.Chronometer.message_with_time("Вызвана функция получения адаптера брокера сообщений/очередей...")
        @debug "Параметры вызова функции:"
        @debug "> параметр для получения адаптера брокера" value
        adapter = nothing
        @debug "Инициализирован параметр adapter с пустым значением"
        adapter_name = nothing
        @debug "Инициализирован параметр adapter_name с пустым значением"
        value_type = typeof(value)
        @debug "Получили тип переданного в функцию параметра" value_type
        @debug "Сопоставляем тип переданного параметра. В зависимости от него продоим вычисление параматра adapter_name"
        if value_type === String
            @debug "Тип переданного значения определился как строка. Считаем переданное значение именем адаптера" value
            adapter_name = value
        elseif value_type == Dict{String, Any}
            @debug "Тип переданного значения опредилился как словарь (String, Any). Получаем имя адаптера из словаря"
            adapter_name = get(value, "adapter", "rabbit")
        else
            @debug "Нет корретного значеня типа. Считаем имя адаптера - пустым"
            adapter_name = nothing
        end

        @debug "По имени адаптера присваиваем модуль"
        if adapter_name === "rabbit"
            @debug "Адаптеру соответствует модуль для работы с RabbitMQ"
            adapter = Rabbit
        else
            @debug "Нет совпадений модулей адаптеров"
            adapter = nothing
        end

        @debug "Возвразаем адаптер" adapter
        return adapter
    end

    function init(config_path::String)
        try
            connection = nothing
            broker_config, channels_config = Utils.ConfigReader.get_configs(config_path)
            adapter, mechanism, login, password, host, port, virtual_host = Utils.Getters.connection_params(broker_config)

            Adapter = get_adapter(adapter)

            if Adapter !== nothing
                auth_params = Adapter.get_auth_params(login, password, mechanism)
                connection = Adapter.get_connection(virtual_host, host, port, auth_params)
            end

            return connection, broker_config, channels_config
        catch error
            @error error
            throw(error)
        end
    end

    function declares(connection, broker_config, channels_config)
        try
            Adapter = get_adapter(broker_config)

            results = nothing

            names = keys(channels_config)
            for name in names
                channel_config = get(channels_config, name, nothing)

                channel = Utils.Getters.channel(connection, name, true, Adapter)
                exchangers, queues = Utils.Declares.in_channel(channel, channel_config, Adapter)
                channel_to_add = Utils.Getters.channel_to_add(channel, exchangers, queues)
                if channel_to_add !== nothing
                    if results === nothing
                       results = Dict{String, Dict{String, Any}}(name => channel_to_add)
                    else
                        push!(resultsg, name => channel_to_add)
                    end
                end
            end

            return results
        catch error
            @error error
            throw(error)
        end
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions related Data and Time conversion.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] https://support.microsoft.com/en-us/help/214019/method-to-determine-whether-a-year-is-a-leap-year
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export is_leap_year
export JD_UT1toUTC, JD_UTCtoUT1
export JD_UTCtoTT,  JD_TTtoUTC
export get_ΔAT

################################################################################
#                                  Constants
################################################################################

# Table containing the leap seconds between UTC and the International Atomic
# Time (TAI). Notice that only the dates in which an increment in the leap
# seconds occurred should be added to this table.
#
#           Julian Day         ΔAT [s]
const ΔAT_Data = [2441499.500000     11.0;
                  2441683.500000     12.0;
                  2442048.500000     13.0;
                  2442413.500000     14.0;
                  2442778.500000     15.0;
                  2443144.500000     16.0;
                  2443509.500000     17.0;
                  2443874.500000     18.0;
                  2444239.500000     19.0;
                  2444786.500000     20.0;
                  2445151.500000     21.0;
                  2445516.500000     22.0;
                  2446247.500000     23.0;
                  2447161.500000     24.0;
                  2447892.500000     25.0;
                  2448257.500000     26.0;
                  2448804.500000     27.0;
                  2449169.500000     28.0;
                  2449534.500000     29.0;
                  2450083.500000     30.0;
                  2450630.500000     31.0;
                  2451179.500000     32.0;
                  2453736.500000     33.0;
                  2454832.500000     34.0;
                  2456109.500000     35.0;
                  2457204.500000     36.0;
                  2457754.500000     37.0;
                 ]

"""
    function get_ΔAT(JD::Number)

Get the accumulated leap seconds (ΔAT) [s] between UTC and International Atomic
Time (TAI) in the given `JD`. This function search for ΔAT in the array
`ΔAT_Data`.

# Remarks

If `JD` is before `ΔAT_Data[1,1]`, then 10 will be returned. **Notice that this
can lead to errors.**

If `JD` is after `ΔAT_Data[end,1]`, then `ΔAT_Data[end,2]` will be returned,
because it is not possible yet to predict when leap seconds will be added.

"""
function get_ΔAT(JD::Number)
    # If `JD` is before `ΔAT_Data[1,1]`, then return 10.0.
    @inbounds if JD < ΔAT_Data[1,1]
        return 10.0
    else
        for i = 2:size(ΔAT_Data,1)
            (JD < ΔAT_Data[i,1]) && return ΔAT_Data[i-1,2]
        end
    end

    # In this case, `JD` is after `ΔAT_Data[end,1]`.
    return @inbounds ΔAT_Data[end,2]
end

################################################################################
#                                  Functions
################################################################################

"""
    function JD_UTCtoUT1(JD_UTC::Number, ΔUT1::Number)

Convert the Julian Day in UTC `JD_UTC` to the Julian Day in UT1 using the
accumulated difference `ΔUT1`, which is provided by IERS EOP Data.

"""
JD_UTCtoUT1(JD_UTC::Number, ΔUT1::Number) = JD_UTC + ΔUT1/86400

"""
    function JD_UT1toUTC(JD_UT1::Number, ΔUT1::Number)

Convert the Julian Day in UT1 `JD_UT1` to the Julian Day in UTC using the
accumulated difference `ΔUT1`, which is provided by IERS EOP Data.

"""
JD_UT1toUTC(JD_UT1::Number, ΔUT1::Number) = JD_UT1 - ΔUT1/86400

"""
    function JD_UTCtoUT1(JD_UTC::Number, eop::Union{EOPData_IAU1980,EOPData_IAU2000A})

Convert the Julian Day in UTC `JD_UTC` to the Julian Day in UT1 using the
accumulated difference given by the EOP Data `eop` (see `get_iers_eop`). Notice
that the accumulated difference will be interpolated.

"""
function JD_UTCtoUT1(JD_UTC::Number, eop::Union{EOPData_IAU1980,EOPData_IAU2000A})
	JD_UTCtoUT1(JD_UTC, eop.UT1_UTC(JD_UTC))
end

"""
    function JD_UTCtoUT1(JD_UTC::Number, eop::Union{EOPData_IAU1980,EOPData_IAU2000A})

Convert the Julian Day in UT1 `JD_UT1` to the Julian Day in UTC using the
accumulated difference given by the EOP Data `eop` (see `get_iers_eop`). Notice
that the accumulated difference will be interpolated.

"""
function JD_UT1toUTC(JD_UT1::Number, eop::Union{EOPData_IAU1980,EOPData_IAU2000A})
	JD_UT1toUTC(JD_UT1, eop.UT1_UTC(JD_UT1))
end

"""
    function JD_UTCtoTT(JD_UTC::Number [, ΔAT::Number])

Convert the Julian Day in UTC `JD_UTC` to the Julian Day in TT (Terrestrial
Time) using the accumulated difference `ΔAT` between UTC and the International
Atomic Time (TAI). If no value is provided, then the leap seconds will be
obtained from the table `ΔAT_Data`. **Notice that, in this case, if a date
previous to 1973 is provided, then a fixed value of 10 will be used, leading to
wrong computations.**

"""
JD_UTCtoTT(JD_UTC::Number, ΔAT::Number) = JD_UTC + (ΔAT + 32.184)/86400

function JD_UTCtoTT(JD_UTC::Number)
    ΔAT = get_ΔAT(JD_UTC)
    JD_UTCtoTT(JD_UTC,ΔAT)
end

"""
    function JD_TTtoUTC(JD_TT::Number, ΔAT::Number = 37)

Convert the Julian Day in TT `JD_TT` (Terrestrial Time) to the Julian Day in UTC
(Terrestrial Time) using the accumulated difference `ΔAT` between UTC and the
International Atomic Time (TAI). If no value is provided, then the leap seconds
will be obtained from the table `ΔAT_Data`. **Notice that, in this case, if a
date previous to 1973 is provided, then a fixed value of 10 will be used,
leading to wrong computations.**

"""
JD_TTtoUTC(JD_TT::Number, ΔAT::Number) = JD_TT - (ΔAT + 32.184)/86400

function JD_TTtoUTC(JD_UTC::Number)
    ΔAT = get_ΔAT(JD_UTC)
    JD_TTtoUTC(JD_UTC,ΔAT)
end

"""
    function is_leap_year(year::Int)

Check if the year `year` is a leap year. It returns `true` if `year` is a leap
year, or `false` otherwise.

# Remarks

This algorithm was based on [3].

"""
function is_leap_year(year::Int)
    # Check if `year` is positive. This algorithm does not handle negative
    # years.
    (year < 0) && throw(ArgumentError("The year must be positive."))

    if (year % 4) == 0
        if (year % 100) == 0
            if (year % 400) == 0
                return true
            else
                return false
            end
        else
            return true
        end
    else
        return false
    end
end


text:
using Gadfly

set_default_plot_size(21cm, 8cm)

y = [0.46, 0.13, 0.4, 0.73, 0.43]

p1 = plot(x=[0,10], y=[0,1], Geom.blank)
p2 = plot()
p3 = push!(plot(), layer(x=1:5, y=y, Geom.point))

hstack(p1, p2, p3)


text:
include("Graph.jl")
include("Compress.jl")
include("PlanarPolyDraw.jl")
include("MakePicture.jl")
include("AddLine.jl")

using Printf

function DrawNetwork(ags1, ags2)
	lg = open("../logs/DrawLog.txt", "a")

	# read network data
	filePath = "../data/"
	fileName = string(filePath, ags1, ".txt")
	G0 = readGraph(fileName, "unweighted")
	println(lg, "Network ", ags1, " has ", G0.n, " nodes and ", G0.m, " edges.")

	# get compressed network
	buf = split(ags2, ',')
	selectedNodes = (size(buf, 1) > 0) ? parse(Int, buf[1]) : 20
	selectedEdges = (size(buf, 1) > 1) ? parse(Int, buf[2]) : 20
	selectMode = (size(buf, 1) > 2) ? buf[3] : "ne"
	addMode = (size(buf, 1) > 3) ? buf[4] : "KeepOne"
	G, NodeC, maxNodeC = Compress(G0, selectedNodes, selectedEdges, selectMode, addMode)
	nci, eci, ci = CentralityInformation(G, NodeC, maxNodeC)
	println(lg, "Compress Mode : ", selectedNodes, " ", selectedEdges, " ", selectMode, " ", addMode)
	println(lg, "Compressed Graph has ", G.n, " nodes and ", G.m, " edges.")
	println(lg, "Node Centrality Information : ", @sprintf("%.2f", nci*100), " %")
	println(lg, "Edge Centrality Information : ", @sprintf("%.2f", eci*100), " %")
	println(lg, "Centrality Information : ", @sprintf("%.2f", ci*100), " %")

	# draw a planar subgraph of G by using polyline drawing
	X, Y, PATH, deletedEdgeList, L, H, preservedEdgeC, deletedEdgeC = PlanarPolyDraw(G)
	println(lg, "The planar subgraph has ", size(PATH, 1), " edges.")
	println(lg, "Planar Subgraph ECI : ", @sprintf("%.2f", (preservedEdgeC/eci)*100), " %")

	# try to add some deleted edges
	X, Y, PATH, PATH2, L, H, dropedEdges = AddLine(X, Y, PATH, deletedEdgeList, L, H)
	println(lg, "Total selected edges : ", size(PATH, 1)+size(PATH2, 1)-dropedEdges, "(+", size(deletedEdgeList, 1)-dropedEdges, ")")
	println(lg, "Dropped edges : ", dropedEdges)

	nodeInfo, edgeInfo, totalInfo = CalculatePreservedCentralityInformation(G, NodeC, maxNodeC, preservedEdgeC, deletedEdgeC, PATH2)
	println(lg, "Preserved ECI : ", @sprintf("%.2f", (edgeInfo/eci)*100), " %")
	println(lg)

	# generate picture
	generatePicture(X, Y, G.V, PATH, PATH2, NodeC, ags1, max(L, H))
end


text:
function AngularAzEl(arg0::GroundStation, arg1::AbsoluteDate, arg2::Vector{jdouble}, arg3::Vector{jdouble}, arg4::Vector{jdouble}, arg5::ObservableSatellite)
    return AngularAzEl((GroundStation, AbsoluteDate, Vector{jdouble}, Vector{jdouble}, Vector{jdouble}, ObservableSatellite), arg0, arg1, arg2, arg3, arg4, arg5)
end

function get_station(obj::AngularAzEl)
    return jcall(obj, "getStation", GroundStation, ())
end


text:
# export get_pha!

function irisws(cha::String, d0::String, d1::String;
                fmt::String   = KW.fmt,
                to::Int       = KW.to,
                opts::String  = KW.opts,
                v::Int        = KW.v,
                w::Bool       = KW.w)

  # init
  parse_err = false
  S = SeisData()
  Ch = SeisChannel()
  parsable = false
  if fmt == "mseed"
    fmt = "miniseed"
  end

  # parse channel string cha
  c = (parse_chstr(cha)[1,:])[1:min(end,4)]
  if isempty(c[3])
    c[3] = "--"
  end
  ID = join([c[1], c[2], strip(c[3],'-'), c[4]], '.')
  setfield!(Ch, :id, ID)

  # Build query URL
  url = "http://service.iris.edu/irisws/timeseries/1/query?" *
          build_stream_query(c,d0,d1) * "&scale=AUTO&output=" * fmt
  v > 0 && println(url)
  Ch.src = url
  req_info_str = datareq_summ("IRISWS data", ID, d0, d1)

  # Do request
  (R, parsable) = get_http_req(url, req_info_str, to)
  if parsable
    if w
      savereq(R, fmt, Ch.id, d0)
    end
    if fmt == "sacbl"
      Ch = read_sac_stream(IOBuffer(R), BUF.sac_fv, BUF.sac_iv, BUF.sac_cv, false, false)
      if isempty(Ch.name)
        Ch.name = deepcopy(Ch.id)
      end
    elseif fmt == "miniseed"
      parsemseed!(S, IOBuffer(R), v, KW.nx_add, KW.nx_add)
      Ch = S[1]
      if isempty(Ch.loc)
        Ch.loc = GeoLoc()
      end
    elseif fmt == "geocsv"
      read_geocsv_tspair!(S, IOBuffer(R))
      Ch = S[1]
    else
      # other parsers not yet written
      @warn(string("Unsupported data format", req_info_str, "\nFORMAT = ", fmt,
            "\n\nUnparsed request data in .misc[\"data\"]"))
      Ch.misc["data"] = R
    end
  else
    parse_err = true
    Ch.misc["data"] = String(R)
  end
  return parse_err, Ch
end

function IRISget(C::Array{String,1}, d0::String, d1::String;
                  fmt::String   = KW.fmt,
                  to::Int       = KW.to,
                  opts::String  = KW.opts,
                  v::Int        = KW.v,
                  w::Bool       = KW.w)

  parse_err = false
  S = SeisData()
  K = size(C,1)
  v > 0 && println("IRISWS data request begins...")
  for k = 1:K
    (p, Ch) = irisws(C[k], d0, d1, fmt = fmt, opts = opts, to = to, v = v, w = w)
    S += Ch
    parse_err = max(parse_err, p)
  end
  return parse_err, S
end

# Programming note: if this method is the default, and S is only modified
# within the for loop, then S is copied to a local scope and the newly-added
# data are deleted upon return
function IRISget!(S::SeisData, C::Array{String,1}, d0::String, d1::String;
                  fmt::String   = KW.fmt,
                  to::Int       = KW.to,
                  opts::String  = KW.opts,
                  v::Int        = KW.v,
                  w::Bool       = KW.w)

  (parse_err, U) = IRISget(C, d0, d1, fmt = fmt, to = to, opts = opts, v = v, w = w)
  append!(S,U)
  return parse_err
end


text:
table = dbc_table(
    # using the same table as in the above example
    table_header * table_body,
    bordered = true,
    dark = true,
    hover = true,
    responsive = true,
    striped = true,
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
using LinearAlgebra, StaticArrays
using OrdinaryDiffEq
using OrdinaryDiffEq: ODEIntegrator
using DynamicalSystemsBase: MinimalDiscreteIntegrator

export lyapunovs, lyapunov

#####################################################################################
#                               Lyapunov Spectum                                    #
#####################################################################################
"""
    lyapunovs(ds::DynamicalSystem, N [, k::Int | Q0]; kwargs...) -> λs

Calculate the spectrum of Lyapunov exponents [1] of `ds` by applying
a QR-decomposition on the parallelepiped matrix `N` times. Return the
spectrum sorted from maximum to minimum.

The third argument `k` is optional, and dictates how many lyapunov exponents
to calculate (defaults to `dimension(ds)`).
Instead of passing an integer `k` you can pass
a pre-initialized matrix `Q0` whose columns are initial deviation vectors (then
`k = size(Q0)[2]`).

## Keyword Arguments
* `u0 = get_state(ds)` : State to start from.
* `Ttr = 0` : Extra "transient" time to evolve the system before application of the
  algorithm. Should be `Int` for discrete systems. Both the system and the
  deviation vectors are evolved for this time.
* `dt` : Time of individual evolutions
  between successive orthonormalization steps. Defaults to `1`. For continuous
  systems this is approximate.
* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

## Description
The method we employ is "H2" of [2], originally stated in [3]. The deviation vectors
defining a `D`-dimensional parallepiped in tangent space
are evolved using the tangent dynamics of the system.
A QR-decomposition at each step yields the local growth rate for each dimension
of the parallepiped. The growth rates are
then averaged over `N` successive steps, yielding the lyapunov exponent spectrum
(at each step the parallepiped is re-normalized).

## Performance Notes
This function uses a [`tangent_integrator`](@ref).
For loops over initial conditions and/or
parameter values one should use the low level method that accepts
an integrator, and `reinit!` it to new initial conditions.
See the "advanced documentation" for info on the integrator object.
The low level method is
```
lyapunovs(tinteg, N, dt::Real, Ttr::Real)
```

If you want to obtain the convergence timeseries of the Lyapunov spectrum,
use the method
```
ChaosTools.lyapunovs_convergence(tinteg, N, dt, Ttr)
```
(not exported).

## References

[1] : A. M. Lyapunov, *The General Problem of the Stability of Motion*,
Taylor & Francis (1992)

[2] : K. Geist *et al.*, Progr. Theor. Phys. **83**, pp 875 (1990)

[3] : G. Benettin *et al.*, Meccanica **15**, pp 9-20 & 21-30 (1980)
"""
lyapunovs(ds::DS, N, k::Int = dimension(ds); kwargs...) =
lyapunovs(ds, N, orthonormal(dimension(ds), k); kwargs...)

function lyapunovs(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix; Ttr::Real = 0,
    dt::Real = 1, u0 = get_state(ds), diffeq...) where {IIP, S, D}

    T = stateeltype(ds)
    # Create tangent integrator:
    if typeof(ds) <: DDS
        @assert typeof(Ttr) == Int
        integ = tangent_integrator(ds, Q0; u0 = u0)
    else
        integ = tangent_integrator(ds, Q0; u0 = u0, diffeq...)
    end

    λ::Vector{T} = lyapunovs(integ, N, dt, Ttr)
    return λ
end

function lyapunovs(integ, N, dt::Real, Ttr::Real = 0.0)

    T = stateeltype(integ)
    t0 = integ.t
    if Ttr > 0
        while integ.t < t0 + Ttr
            step!(integ, dt)
            qrdec = LinearAlgebra.qr(get_deviations(integ))
            set_deviations!(integ, _get_Q(qrdec))
        end
    end
    k = size(get_deviations(integ))[2]
    λ::Vector{T} = zeros(T, k)
    t0 = integ.t

    for i in 2:N
        step!(integ, dt)
        qrdec = LinearAlgebra.qr(get_deviations(integ))
        for j in 1:k
            @inbounds λ[j] += log(abs(qrdec.R[j,j]))
        end
        set_deviations!(integ, _get_Q(qrdec))
    end
    λ ./= (integ.t - t0)
    return λ
end

_get_Q(qrdec::StaticArrays.QR) = qrdec.Q
_get_Q(qrdec::LinearAlgebra.QRCompactWY) = Matrix(qrdec.Q)

lyapunovs(ds::DynamicalSystem{IIP, T, 1}, a...; kw...) where {IIP, T} = error(
"For 1D systems, only discrete & out-of-place method is implemented.")

function lyapunovs(ds::DDS{false, T, 1}, N; Ttr = 0) where {T}

    x = get_state(ds); f = ds.f
    p = ds.p; t0 = ds.t0
    t = 0
    if Ttr > 0
        for i in t0:(Ttr+t0)
            x = f(x, p, i)
        end
    end
    λ = zero(T)
    for i in (t0+Ttr):(t0+Ttr+N)
        x = f(x, p, i)
        λ += log(abs(ds.jacobian(x, p, i)))
    end
    return λ/N
end

lyapunov(ds::DDS{false, T, 1}, N; Ttr = 0) where {T} = lyapunovs(ds, N; Ttr = Ttr)

#####################################################################################
#                           Maximum Lyapunov Exponent                               #
#####################################################################################
inittest_default(D) = (state1, d0) -> state1 .+ d0/sqrt(D)

"""
    lyapunov(ds::DynamicalSystem, Τ; kwargs...) -> λ

Calculate the maximum Lyapunov exponent `λ` using a method due to Benettin [1],
which simply
evolves two neighboring trajectories (one called "given" and one called "test")
while constantly rescaling the test one.
`T`  denotes the total time of evolution (should be `Int` for discrete systems).

## Keyword Arguments

* `Ttr = 0` : Extra "transient" time to evolve the trajectories before
  starting to measure the expontent. Should be `Int` for discrete systems.
* `d0 = 1e-9` : Initial & rescaling distance between the two neighboring trajectories.
* `upper_threshold = 1e-6` : Upper distance threshold for rescaling.
* `lower_threshold = 1e-12` : Lower distance threshold for rescaling (in order to
   be able to detect negative exponents).
* `dt = 1` : Time of evolution between each check of
  distance exceeding the thresholds. For continuous
  systems this is approximate.
* `inittest = (u1, d0) -> u1 .+ d0/sqrt(D)` :
  A function that given `(u1, d0)`
  initializes the test state with distance
  `d0` from the given state `u1` (`D` is the dimension
  of the system). This function can be used when you want to avoid
  the test state appearing in a region of the phase-space where it would have
  e.g. different energy or escape to infinity.
* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
  See [`trajectory`](@ref) for examples. Only valid for continuous systems.

## Description
Two neighboring trajectories with initial distance `d0` are evolved in time.
At time ``t_i`` their distance ``d(t_i)`` either exceeds the `upper_threshold`,
or is lower than `lower_threshold`, which initializes
a rescaling of the test trajectory back to having distance `d0` from
the given one, while the rescaling keeps the difference vector along the maximal
expansion/contraction direction: `` u_2 \\to u_1+(u_2−u_1)/(d(t_i)/d_0)``.

The maximum
Lyapunov exponent is the average of the time-local Lyapunov exponents
```math
\\lambda = \\frac{1}{t_{n}}\\sum_{i=1}^{n}
\\ln\\left( a_i \\right),\\quad a_i = \\frac{d(t_{i})}{d_0}.
```

## Performance Notes
This function uses a [`parallel_integrator`](@ref).
For loops over initial conditions and/or
parameter values one should use the low level method that accepts
an integrator, and `reinit!` it to new initial conditions.
See the "advanced documentation" for info on the integrator object.
The low level method is
```
lyapunov(pinteg, T, Ttr, dt, d0, ut, lt)
```

## References
[1] : G. Benettin *et al.*, Phys. Rev. A **14**, pp 2338 (1976)
"""
function lyapunov(ds::DS, T;
                  Ttr = 0,
                  d0=1e-9,
                  upper_threshold = 1e-6,
                  lower_threshold = 1e-12,
                  inittest = inittest_default(dimension(ds)),
                  dt = 1,
                  diffeq...
                  )

    ST = stateeltype(ds)
    lower_threshold ≤ d0 ≤ upper_threshold || throw(ArgumentError(
    "d0 must be between thresholds!"))
    D = dimension(ds)
    if typeof(ds) <: DDS
        pinteg = parallel_integrator(ds,
            [deepcopy(get_state(ds)), inittest(get_state(ds), d0)])
    else
        pinteg = parallel_integrator(ds,
            [deepcopy(get_state(ds)), inittest(get_state(ds), d0)];
            diffeq...)
    end
    λ::ST = lyapunov(pinteg, T, Ttr, dt, d0, upper_threshold, lower_threshold)
    return λ
end

function lyapunov(pinteg, T, Ttr, dt, d0, ut, lt)
    # transient
    t0 = pinteg.t
    while pinteg.t < t0 + Ttr
        step!(pinteg, dt)
        d = λdist(pinteg)
        lt ≤ d ≤ ut || rescale!(pinteg, d/d0)
    end

    t0 = pinteg.t
    d = λdist(pinteg)
    d == 0 && error("Initial distance between states is zero!!!")
    rescale!(pinteg, d/d0)
    λ = zero(d)
    while pinteg.t < t0 + T
        d = λdist(pinteg)
        #evolve until rescaling:
        while lt ≤ d ≤ ut
            step!(pinteg, dt)
            d = λdist(pinteg)
            pinteg.t ≥ t0 + T && break
        end
        # local lyapunov exponent is simply the relative distance of the trajectories
        a = d/d0
        λ += log(a)
        rescale!(pinteg, a)
    end
    # Do final rescale, in case no other happened
    d = λdist(pinteg)
    a = d/d0
    λ += log(a)
    return λ/(pinteg.t - t0)
end

################ Helper functions that allow a single definition ######################
function λdist(integ::ODEIntegrator{Alg, M}) where {Alg, M<:Matrix}
    d = 0.0
    for i in 1:size(integ.u)[1]
        d += (integ.u[i, 1] - integ.u[i, 2])^2
    end
    return sqrt(d)
    # return norm(view(integ.u, :, 1) .- view(integ.u, :, 2))
end
# No-annotation case is with Vectors
function λdist(integ)
    @inbounds s = zero(eltype(integ.u[1]))
    @inbounds for k in 1:length(integ.u[1])
        x = (integ.u[1][k] - integ.u[2][k])
        s += x*x
    end
    return sqrt(s)
end
function λdist(integ::MinimalDiscreteIntegrator{true, Vector{S}}) where {S<:SVector}
    return norm(integ.u[1] - integ.u[2])
end
function λdist(integ::ODEIntegrator{Alg, Vector{S}}) where {Alg, S<:SVector}
    return norm(integ.u[1] - integ.u[2])
end

# Rescales:
function rescale!(integ::MinimalDiscreteIntegrator{true, Vector{S}}, a) where {S<:SVector}
    integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::MinimalDiscreteIntegrator{true, Vector{S}}, a) where {S<:Vector}
    @. integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::ODEIntegrator{Alg, M}, a) where {Alg, M<:Matrix}
    for i in 1:size(integ.u)[1]
        integ.u[i, 2] = integ.u[i,1] + (integ.u[i,2] - integ.u[i,1])/a
    end
    u_modified!(integ, true)
end
function rescale!(integ::ODEIntegrator{Alg, Vector{S}}, a) where {Alg, S<:Vector}
    @. integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end
function rescale!(integ::ODEIntegrator{Alg, Vector{S}}, a) where {Alg, S<:SVector}
    integ.u[2] = integ.u[1] + (integ.u[2] - integ.u[1])/a
    u_modified!(integ, true)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Printf

using Base.Ryu

export @printf, @sprintf

# format specifier categories
const Ints = Union{Val{'d'}, Val{'i'}, Val{'u'}, Val{'x'}, Val{'X'}, Val{'o'}}
const Floats = Union{Val{'e'}, Val{'E'}, Val{'f'}, Val{'F'}, Val{'g'}, Val{'G'}, Val{'a'}, Val{'A'}}
const Chars = Union{Val{'c'}, Val{'C'}}
const Strings = Union{Val{'s'}, Val{'S'}}
const Pointer = Val{'p'}
const HexBases = Union{Val{'x'}, Val{'X'}, Val{'a'}, Val{'A'}}

"""
Typed representation of a format specifier.

`T` is a `Val{'_'}`, where `_` is a valid format specifier character.

Fields are the various modifiers allowed for various format specifiers.
"""
struct Spec{T} # T => %type => Val{'type'}
    leftalign::Bool
    plus::Bool
    space::Bool
    zero::Bool
    hash::Bool
    width::Int
    precision::Int
end

# recreate the format specifier string from a typed Spec
Base.string(f::Spec{T}; modifier::String="") where {T} =
    string("%", f.leftalign ? "-" : "", f.plus ? "+" : "", f.space ? " " : "",
        f.zero ? "0" : "", f.hash ? "#" : "", f.width > 0 ? f.width : "",
        f.precision == 0 ? ".0" : f.precision > 0 ? ".$(f.precision)" : "", modifier, char(T))
Base.show(io::IO, f::Spec) = print(io, string(f))

floatfmt(s::Spec{T}) where {T} =
    Spec{Val{'f'}}(s.leftalign, s.plus, s.space, s.zero, s.hash, s.width, 0)
ptrfmt(s::Spec{T}, x) where {T} =
    Spec{Val{'x'}}(s.leftalign, s.plus, s.space, s.zero, true, s.width, sizeof(x) == 8 ? 16 : 8)

"""
    Printf.Format(format_str)

Create a C printf-compatible format object that can be used for formatting values.

The input `format_str` can include any valid format specifier character and modifiers.

A `Format` object can be passed to `Printf.format(f::Format, args...)` to produce a
formatted string, or `Printf.format(io::IO, f::Format, args...)` to print the
formatted string directly to `io`.

For convenience, the `Printf.format"..."` string macro form can be used for building
a `Printf.Format` object at macro-expansion-time.
"""
struct Format{S, T}
    str::S # original full format string as CodeUnits
    # keep track of non-format specifier strings to print
    # length(substringranges) == length(formats) + 1
    # so when printing, we start with printing
      # str[substringranges[1]], then formats[1] + args[1]
      # then str[substringranges[2]], then formats[2] + args[2]
      # and so on, then at the end, str[substringranges[end]]
    substringranges::Vector{UnitRange{Int}}
    formats::T # Tuple of Specs
end

# what number base should be used for a given format specifier?
base(T) = T <: HexBases ? 16 : T <: Val{'o'} ? 8 : 10
char(::Type{Val{c}}) where {c} = c

# parse format string
function Format(f::AbstractString)
    isempty(f) && throw(ArgumentError("empty format string"))
    bytes = codeunits(f)
    len = length(bytes)
    pos = 1
    b = 0x00
    while true
        b = bytes[pos]
        pos += 1
        (pos > len || (b == UInt8('%') && pos <= len && bytes[pos] != UInt8('%'))) && break
    end
    strs = [1:pos - 1 - (b == UInt8('%'))]
    fmts = []
    while pos <= len
        b = bytes[pos]
        pos += 1
        # positioned at start of first format str %
        # parse flags
        leftalign = plus = space = zero = hash = false
        while true
            if b == UInt8('-')
                leftalign = true
            elseif b == UInt8('+')
                plus = true
            elseif b == UInt8(' ')
                space = true
            elseif b == UInt8('0')
                zero = true
            elseif b == UInt8('#')
                hash = true
            else
                break
            end
            pos > len && throw(ArgumentError("incomplete format string: '$f'"))
            b = bytes[pos]
            pos += 1
        end
        if leftalign
            zero = false
        end
        # parse width
        width = 0
        while b - UInt8('0') < 0x0a
            width = 10 * width + (b - UInt8('0'))
            b = bytes[pos]
            pos += 1
            pos > len && break
        end
        # parse precision
        precision = 0
        parsedprecdigits = false
        if b == UInt8('.')
            pos > len && throw(ArgumentError("incomplete format string: '$f'"))
            parsedprecdigits = true
            b = bytes[pos]
            pos += 1
            if pos <= len
                while b - UInt8('0') < 0x0a
                    precision = 10precision + (b - UInt8('0'))
                    b = bytes[pos]
                    pos += 1
                    pos > len && break
                end
            end
        end
        # parse length modifier (ignored)
        if b == UInt8('h') || b == UInt8('l')
            prev = b
            b = bytes[pos]
            pos += 1
            if b == prev
                pos > len && throw(ArgumentError("invalid format string: '$f'"))
                b = bytes[pos]
                pos += 1
            end
        elseif b in b"Ljqtz"
            b = bytes[pos]
            pos += 1
        end
        # parse type
        !(b in b"diouxXDOUeEfFgGaAcCsSpn") && throw(ArgumentError("invalid format string: '$f', invalid type specifier: '$(Char(b))'"))
        type = Val{Char(b)}
        if type <: Ints && precision > 0
            zero = false
        elseif (type <: Strings || type <: Chars) && !parsedprecdigits
            precision = -1
        elseif type <: Union{Val{'a'}, Val{'A'}} && !parsedprecdigits
            precision = -1
        elseif type <: Floats && !parsedprecdigits
            precision = 6
        end
        push!(fmts, Spec{type}(leftalign, plus, space, zero, hash, width, precision))
        start = pos
        while pos <= len
            b = bytes[pos]
            pos += 1
            if b == UInt8('%')
                pos > len && throw(ArgumentError("invalid format string: '$f'"))
                if bytes[pos] == UInt8('%')
                    # escaped '%'
                    b = bytes[pos]
                    pos += 1
                else
                    break
                end
            end
        end
        push!(strs, start:pos - 1 - (b == UInt8('%')))
    end
    return Format(bytes, strs, Tuple(fmts))
end

macro format_str(str)
    Format(str)
end

const hex = b"0123456789abcdef"
const HEX = b"0123456789ABCDEF"

# write out a single arg according to format options
# char
@inline function writechar(buf, pos, c)
    u = bswap(reinterpret(UInt32, c))
    while true
        buf[pos] = u % UInt8
        pos += 1
        (u >>= 8) == 0 && break
    end
    return pos
end

@inline function fmt(buf, pos, arg, spec::Spec{T}) where {T <: Chars}
    leftalign, width = spec.leftalign, spec.width
    if !leftalign && width > 1
        for _ = 1:(width - 1)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    pos = writechar(buf, pos, arg isa String ? arg[1] : Char(arg))
    if leftalign && width > 1
        for _ = 1:(width - 1)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    return pos
end

# strings
@inline function fmt(buf, pos, arg, spec::Spec{T}) where {T <: Strings}
    leftalign, hash, width, prec = spec.leftalign, spec.hash, spec.width, spec.precision
    str = string(arg)
    op = p = prec == -1 ? (length(str) + (hash ? arg isa AbstractString ? 2 : 1 : 0)) : prec
    if !leftalign && width > p
        for _ = 1:(width - p)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    if hash
        if arg isa Symbol
            buf[pos] = UInt8(':')
            pos += 1
            p -= 1
        elseif arg isa AbstractString
            buf[pos] = UInt8('"')
            pos += 1
            p -= 1
        end
    end
    for c in str
        p == 0 && break
        pos = writechar(buf, pos, c)
        p -= 1
    end
    if hash && arg isa AbstractString && p > 0
        buf[pos] = UInt8('"')
        pos += 1
    end
    if leftalign && width > op
        for _ = 1:(width - op)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    return pos
end

# integers
toint(x) = x
toint(x::Rational) = Integer(x)

fmt(buf, pos, arg::AbstractFloat, spec::Spec{T}) where {T <: Ints} =
    fmt(buf, pos, arg, floatfmt(spec))

@inline function fmt(buf, pos, arg, spec::Spec{T}) where {T <: Ints}
    leftalign, plus, space, zero, hash, width, prec =
        spec.leftalign, spec.plus, spec.space, spec.zero, spec.hash, spec.width, spec.precision
    bs = base(T)
    arg2 = toint(arg)
    n = i = ndigits(arg2, base=bs, pad=1)
    x, neg = arg2 < 0 ? (-arg2, true) : (arg2, false)
    arglen = n + (neg || (plus | space)) +
        (T == Val{'o'} && hash ? 1 : 0) +
        (T == Val{'x'} && hash ? 2 : 0) + (T == Val{'X'} && hash ? 2 : 0)
    arglen2 = arglen < width && prec > 0 ? arglen + min(max(0, prec - n), width - arglen) : arglen
    if !leftalign && !zero && arglen2 < width
        # pad left w/ spaces
        for _ = 1:(width - arglen2)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    if neg
        buf[pos] = UInt8('-'); pos += 1
    elseif plus # plus overrides space
        buf[pos] = UInt8('+'); pos += 1
    elseif space
        buf[pos] = UInt8(' '); pos += 1
    end
    if T == Val{'o'} && hash
        buf[pos] = UInt8('0')
        pos += 1
    elseif T == Val{'x'} && hash
        buf[pos] = UInt8('0')
        buf[pos + 1] = UInt8('x')
        pos += 2
    elseif T == Val{'X'} && hash
        buf[pos] = UInt8('0')
        buf[pos + 1] = UInt8('X')
        pos += 2
    end
    if zero && arglen2 < width
        for _ = 1:(width - arglen2)
            buf[pos] = UInt8('0')
            pos += 1
        end
    elseif n < prec
        for _ = 1:(prec - n)
            buf[pos] = UInt8('0')
            pos += 1
        end
    elseif arglen < arglen2
        for _ = 1:(arglen2 - arglen)
            buf[pos] = UInt8('0')
            pos += 1
        end
    end
    while i > 0
        @inbounds buf[pos + i - 1] = bs == 16 ?
            (T == Val{'x'} ? hex[(x & 0x0f) + 1] : HEX[(x & 0x0f) + 1]) :
            (48 + (bs == 8 ? (x & 0x07) : rem(x, 10)))
        if bs == 8
            x >>= 3
        elseif bs == 16
            x >>= 4
        else
            x = oftype(x, div(x, 10))
        end
        i -= 1
    end
    pos += n
    if leftalign && arglen2 < width
        # pad right
        for _ = 1:(width - arglen2)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    return pos
end

# floats
"""
    Printf.tofloat(x)

Convert an argument to a Base float type for printf formatting.
By default, arguments are converted to `Float64` via `Float64(x)`.
Custom numeric types that have a conversion to a Base float type
that wish to hook into printf formatting can extend this method like:

```julia
Printf.tofloat(x::MyCustomType) = convert_my_custom_type_to_float(x)
```

For arbitrary precision numerics, you might extend the method like:

```julia
Printf.tofloat(x::MyArbitraryPrecisionType) = BigFloat(x)
```
"""
tofloat(x) = Float64(x)
tofloat(x::Base.IEEEFloat) = x
tofloat(x::BigFloat) = x

@inline function fmt(buf, pos, arg, spec::Spec{T}) where {T <: Floats}
    leftalign, plus, space, zero, hash, width, prec =
        spec.leftalign, spec.plus, spec.space, spec.zero, spec.hash, spec.width, spec.precision
    x = tofloat(arg)
    if x isa BigFloat && isfinite(x)
        ptr = pointer(buf, pos)
        newpos = @ccall "libmpfr".mpfr_snprintf(ptr::Ptr{UInt8}, (length(buf) - pos + 1)::Csize_t, string(spec; modifier="R")::Ptr{UInt8}; arg::Ref{BigFloat})::Cint
        newpos > 0 || error("invalid printf formatting for BigFloat")
        return pos + newpos
    elseif x isa BigFloat
        x = Float64(x)
    end
    if T == Val{'e'} || T == Val{'E'}
        newpos = Ryu.writeexp(buf, pos, x, prec, plus, space, hash, char(T), UInt8('.'))
    elseif T == Val{'f'} || T == Val{'F'}
        newpos = Ryu.writefixed(buf, pos, x, prec, plus, space, hash, UInt8('.'))
    elseif T == Val{'g'} || T == Val{'G'}
        prec = prec == 0 ? 1 : prec
        x = round(x, sigdigits=prec)
        newpos = Ryu.writeshortest(buf, pos, x, plus, space, hash, prec, T == Val{'g'} ? UInt8('e') : UInt8('E'), true, UInt8('.'))
    elseif T == Val{'a'} || T == Val{'A'}
        x, neg = x < 0 ? (-x, true) : (x, false)
        newpos = pos
        if neg
            buf[newpos] = UInt8('-')
            newpos += 1
        elseif plus
            buf[newpos] = UInt8('+')
            newpos += 1
        elseif space
            buf[newpos] = UInt8(' ')
            newpos += 1
        end
        if isnan(x)
            buf[newpos] = UInt8('N')
            buf[newpos + 1] = UInt8('a')
            buf[newpos + 2] = UInt8('N')
            newpos += 3
        elseif !isfinite(x)
            buf[newpos] = UInt8('I')
            buf[newpos + 1] = UInt8('n')
            buf[newpos + 2] = UInt8('f')
            newpos += 3
        else
            buf[newpos] = UInt8('0')
            newpos += 1
            buf[newpos] = T <: Val{'a'} ? UInt8('x') : UInt8('X')
            newpos += 1
            if arg == 0
                buf[newpos] = UInt8('0')
                newpos += 1
                if prec > 0
                    while prec > 0
                        buf[newpos] = UInt8('0')
                        newpos += 1
                        prec -= 1
                    end
                end
                buf[newpos] = T <: Val{'a'} ? UInt8('p') : UInt8('P')
                buf[newpos + 1] = UInt8('+')
                buf[newpos + 2] = UInt8('0')
            else
                if prec > -1
                    s, p = frexp(x)
                    sigbits = 4 * min(prec, 13)
                    s = 0.25 * round(ldexp(s, 1 + sigbits))
                    # ensure last 2 exponent bits either 01 or 10
                    u = (reinterpret(UInt64, s) & 0x003f_ffff_ffff_ffff) >> (52 - sigbits)
                    i = n = (sizeof(u) << 1) - (leading_zeros(u) >> 2)
                else
                    s, p = frexp(x)
                    s *= 2.0
                    u = (reinterpret(UInt64, s) & 0x001f_ffff_ffff_ffff)
                    t = (trailing_zeros(u) >> 2)
                    u >>= (t << 2)
                    i = n = 14 - t
                end
                frac = u > 9 || hash || prec > 0
                while i > 1
                    buf[newpos + i] = T == Val{'a'} ? hex[(u & 0x0f) + 1] : HEX[(u & 0x0f) + 1]
                    u >>= 4
                    i -= 1
                    prec -= 1
                end
                if frac
                    buf[newpos + 1] = UInt8('.')
                end
                buf[newpos] = T == Val{'a'} ? hex[(u & 0x0f) + 1] : HEX[(u & 0x0f) + 1]
                newpos += n + frac
                while prec > 0
                    buf[newpos] = UInt8('0')
                    newpos += 1
                    prec -= 1
                end
                buf[newpos] = T <: Val{'a'} ? UInt8('p') : UInt8('P')
                newpos += 1
                p -= 1
                buf[newpos] = p < 0 ? UInt8('-') : UInt8('+')
                p = p < 0 ? -p : p
                newpos += 1
                n = i = ndigits(p, base=10, pad=1)
                while i > 0
                    buf[newpos + i - 1] = 48 + rem(p, 10)
                    p = oftype(p, div(p, 10))
                    i -= 1
                end
                newpos += n
            end
        end
    end
    if newpos - pos < width
        # need to pad
        if leftalign
            # easy case, just pad spaces after number
            for _ = 1:(width - (newpos - pos))
                buf[newpos] = UInt8(' ')
                newpos += 1
            end
        else
            # right aligned
            n = width - (newpos - pos)
            if zero
                ex = (arg < 0 || (plus | space)) + (T <: Union{Val{'a'}, Val{'A'}} ? 2 : 0)
                so = pos + ex
                len = (newpos - pos) - ex
                copyto!(buf, so + n, buf, so, len)
                for i = so:(so + n - 1)
                    buf[i] = UInt8('0')
                end
                newpos += n
            else
                copyto!(buf, pos + n, buf, pos, newpos - pos)
                for i = pos:(pos + n - 1)
                    buf[i] = UInt8(' ')
                end
                newpos += n
            end
        end
    end
    return newpos
end

# pointers
fmt(buf, pos, arg, spec::Spec{Pointer}) = fmt(buf, pos, Int(arg), ptrfmt(spec, arg))

# old Printf compat
function fix_dec end
function ini_dec end

# generic fallback
function fmtfallback(buf, pos, arg, spec::Spec{T}) where {T}
    leftalign, plus, space, zero, hash, width, prec =
        spec.leftalign, spec.plus, spec.space, spec.zero, spec.hash, spec.width, spec.precision
    buf2 = Base.StringVector(309 + 17 + 5)
    ise = T <: Union{Val{'e'}, Val{'E'}}
    isg = T <: Union{Val{'g'}, Val{'G'}}
    isf = T <: Val{'f'}
    if isg
        prec = prec == 0 ? 1 : prec
        arg = round(arg, sigdigits=prec)
    end
    n, pt, neg = isf ? fix_dec(arg, prec, buf2) : ini_dec(arg, min(prec + ise, length(buf2) - 1), buf2)
    if isg && !hash
        while buf2[n] == UInt8('0')
            n -= 1
        end
    end
    expform = ise || (isg && !(-4 < pt <= prec))
    n2 = n + (expform ? 4 : 0) + (prec > 0 || hash) + (neg || (plus | space)) +
        (isf && pt >= n ? prec + 1 : 0)
    if !leftalign && !zero && n2 < width
        # pad left w/ spaces
        for _ = 1:(width - n2)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    if neg
        buf[pos] = UInt8('-'); pos += 1
    elseif plus # plus overrides space
        buf[pos] = UInt8('+'); pos += 1
    elseif space
        buf[pos] = UInt8(' '); pos += 1
    end
    if zero && n2 < width
        for _ = 1:(width - n2)
            buf[pos] = UInt8('0')
            pos += 1
        end
    end
    if expform
        buf[pos] = buf2[1]
        pos += 1
        if n > 1 || hash
            buf[pos] = UInt8('.')
            pos += 1
            for i = 2:n
                buf[pos] = buf2[i]
                pos += 1
            end
        end
        buf[pos] = T <: Val{'e'} || T <: Val{'g'} ? UInt8('e') : UInt8('E')
        pos += 1
        exp = pt - 1
        buf[pos] = exp < 0 ? UInt8('-') : UInt8('+')
        pos += 1
        exp = abs(exp)
        if exp < 10
            buf[pos] = UInt8('0')
            buf[pos + 1] = 48 + exp
            pos += 2
        else
            buf[pos] = 48 + div(exp, 10)
            buf[pos + 1] = 48 + rem(exp, 10)
            pos += 2
        end
    elseif pt <= 0
        buf[pos] = UInt8('0')
        buf[pos + 1] = UInt8('.')
        pos += 2
        while pt < 0
            buf[pos] = UInt8('0')
            pos += 1
            pt += 1
        end
        for i = 1:n
            buf[pos] = buf2[i]
            pos += 1
        end
    elseif pt >= n
        for i = 1:n
            buf[pos] = buf2[i]
            pos += 1
        end
        while pt > n
            buf[pos] = UInt8('0')
            pos += 1
            n += 1
        end
        if hash || (isf && prec > 0)
            buf[pos] = UInt8('.')
            pos += 1
            while prec > 0
                buf[pos] = UInt8('0')
                pos += 1
                prec -= 1
            end
        end
    else
        for i = 1:pt
            buf[pos] = buf2[i]
            pos += 1
        end
        buf[pos] = UInt8('.')
        pos += 1
        for i = pt+1:n
            buf[pos] = buf2[i]
            pos += 1
        end
    end

    if leftalign && n2 < width
        # pad right
        for _ = 1:(width - n2)
            buf[pos] = UInt8(' ')
            pos += 1
        end
    end
    return pos
end

const UNROLL_UPTO = 16
# if you have your own buffer + pos, write formatted args directly to it
@inline function format(buf::Vector{UInt8}, pos::Integer, f::Format, args...)
    # write out first substring
    for i in f.substringranges[1]
        buf[pos] = f.str[i]
        pos += 1
    end
    # for each format, write out arg and next substring
    # unroll up to 16 formats
    N = length(f.formats)
    Base.@nexprs 16 i -> begin
        if N >= i
            pos = fmt(buf, pos, args[i], f.formats[i])
            for j in f.substringranges[i + 1]
                buf[pos] = f.str[j]
                pos += 1
            end
        end
    end
    if N > 16
        for i = 17:length(f.formats)
            pos = fmt(buf, pos, args[i], f.formats[i])
            for j in f.substringranges[i + 1]
                buf[pos] = f.str[j]
                pos += 1
            end
        end
    end
    return pos
end

plength(f::Spec{T}, x) where {T <: Chars} = max(f.width, 1) + (ncodeunits(x isa AbstractString ? x[1] : Char(x)) - 1)
plength(f::Spec{Pointer}, x) = max(f.width, 2 * sizeof(x) + 2)

function plength(f::Spec{T}, x) where {T <: Strings}
    str = string(x)
    p = f.precision == -1 ? (length(str) + (f.hash ? (x isa Symbol ? 1 : 2) : 0)) : f.precision
    return max(f.width, p) + (sizeof(str) - length(str))
end

function plength(f::Spec{T}, x) where {T <: Ints}
    x2 = toint(x)
    return max(f.width, f.precision + ndigits(x2, base=base(T), pad=1) + 5)
end

plength(f::Spec{T}, x::AbstractFloat) where {T <: Ints} =
    max(f.width, 0 + 309 + 17 + f.hash + 5)
plength(f::Spec{T}, x) where {T <: Floats} =
    max(f.width, f.precision + 309 + 17 + f.hash + 5)

@inline function computelen(substringranges, formats, args)
    len = sum(length, substringranges)
    N = length(formats)
    # unroll up to 16 formats
    Base.@nexprs 16 i -> begin
        if N >= i
            len += plength(formats[i], args[i])
        end
    end
    if N > 16
        for i = 17:length(formats)
            len += plength(formats[i], args[i])
        end
    end
    return len
end

@noinline argmismatch(a, b) =
    throw(ArgumentError("mismatch between # of format specifiers and provided args: $a != $b"))

"""
    Printf.format(f::Printf.Format, args...) => String
    Printf.format(io::IO, f::Printf.Format, args...)

Apply a printf format object `f` to provided `args` and return the formatted string
(1st method), or print directly to an `io` object (2nd method). See [`@printf`](@ref)
for more details on C `printf` support.
"""
function format end

function format(io::IO, f::Format, args...) # => Nothing
    length(f.formats) == length(args) || argmismatch(length(f.formats), length(args))
    buf = Base.StringVector(computelen(f.substringranges, f.formats, args))
    pos = format(buf, 1, f, args...)
    write(io, resize!(buf, pos - 1))
    return
end

function format(f::Format, args...) # => String
    length(f.formats) == length(args) || argmismatch(length(f.formats), length(args))
    buf = Base.StringVector(computelen(f.substringranges, f.formats, args))
    pos = format(buf, 1, f, args...)
    return String(resize!(buf, pos - 1))
end

"""
    @printf([io::IO], "%Fmt", args...)

Print `args` using C `printf` style format specification string, with some caveats:
`Inf` and `NaN` are printed consistently as `Inf` and `NaN` for flags `%a`, `%A`,
`%e`, `%E`, `%f`, `%F`, `%g`, and `%G`. Furthermore, if a floating point number is
equally close to the numeric values of two possible output strings, the output
string further away from zero is chosen.
Optionally, an `IO`
may be passed as the first argument to redirect output.
See also: [`@sprintf`](@ref)
# Examples
```jldoctest
julia> @printf("%f %F %f %F\\n", Inf, Inf, NaN, NaN)
Inf Inf NaN NaN\n
julia> @printf "%.0f %.1f %f\\n" 0.5 0.025 -0.0078125
0 0.0 -0.007812
```
"""
macro printf(io_or_fmt, args...)
    if io_or_fmt isa String
        fmt = Format(io_or_fmt)
        return esc(:($Printf.format(stdout, $fmt, $(args...))))
    else
        io = io_or_fmt
        isempty(args) && throw(ArgumentError("must provide required format string"))
        fmt = Format(args[1])
        return esc(:($Printf.format($io, $fmt, $(Base.tail(args)...))))
    end
end

"""
    @sprintf("%Fmt", args...)

Return `@printf` formatted output as string.
# Examples
```jldoctest
julia> @sprintf "this is a %s %15.1f" "test" 34.567
"this is a test            34.6"
```
"""
macro sprintf(fmt, args...)
    f = Format(fmt)
    return esc(:($Printf.format($f, $(args...))))
end

end # module

max_forks_repo_name: jishnub/julia


text:
"""
    LU solver
"""
function lu_solver(A::SparseMatrixCSC{T,Int}) where T
    A = Array(A)
    lu_solver(A)
end

function linear_solve!(solver::LUSolver{T}, x::Vector{T}, A::SparseMatrixCSC{T,Int},
        b::Vector{T}; reg::T = 0.0, fact::Bool = true) where T
    linear_solve!(solver, x, Array(A), b, reg=reg, fact=fact)
end

"""
    LU (sparse) solver
"""
mutable struct LUSparseSolver{T} <: LinearSolver
    A::SparseMatrixCSC{T,Int}
    F::SuiteSparse.UMFPACK.UmfpackLU{T, Int}
    umf_ctrl::Vector{T}
    umf_info::Vector{T}
    tmp::Vector{Ptr{Cvoid}}
end

function lu_sparse_solver(A::SparseMatrixCSC{T,Int}) where T
    F = lu(copy(A))
    UMFPACK_INFO = 90
    UMFPACK_CONTROL = 20
    umf_ctrl = Vector{Float64}(undef, UMFPACK_CONTROL)
    umf_info = Vector{Float64}(undef, UMFPACK_INFO)
    tmp = Vector{Ptr{Cvoid}}(undef, 1)
    LUSparseSolver(copy(A), F, umf_ctrl, umf_info, tmp)
end

function factorize!(s::LUSparseSolver{T}, A::SparseMatrixCSC{T,Int}) where T
    _lu!(s.F, A, s.tmp, s.umf_ctrl, s.umf_info)
    # s.F = lu(A)
end

function linear_solve!(s::LUSparseSolver{T}, x::Vector{T}, A::SparseMatrixCSC{T,Int},
        b::Vector{T}; reg::T = 0.0, fact::Bool = true) where T
    fact && factorize!(s, A)
    _solve!(x, s.F, b, s.umf_ctrl, s.umf_info)
    # x .= A \ b
end

function linear_solve!(s::LUSparseSolver{T}, x::Matrix{T}, A::Matrix{T},
    b::Matrix{T}; reg::T = 0.0, fact::Bool = true) where T
    fill!(x, 0.0)
    n, m = size(x)
    r_idx = 1:n
    fact && factorize!(s, A)
    for j = 1:m
        xv = @views x[r_idx, j]
        bv = @views b[r_idx, j]
        _solve!(xv, s.F, bv, s.umf_ctrl, s.umf_info)
    end
end

function _umfpack_numeric!(U::SuiteSparse.UMFPACK.UmfpackLU{Float64,Int64}, tmp, umf_ctrl, umf_info; reuse_numeric = true)
    status = ccall((:umfpack_di_numeric, :libumfpack), Int64,
                   (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
                    Ptr{Float64}, Ptr{Float64}),
                   U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
                   umf_ctrl, umf_info)
    return nothing
end

function _lu!(F::SuiteSparse.UMFPACK.UmfpackLU, S::SparseMatrixCSC{<:SuiteSparse.UMFPACK.UMFVTypes,<:SuiteSparse.UMFPACK.UMFITypes}, tmp, umf_ctrl, umf_info; check::Bool=true)
    zerobased = SparseArrays.getcolptr(S)[1] == 0
    F.m = size(S, 1)
    F.n = size(S, 2)

    F.colptr .= SparseArrays.getcolptr(S)
    F.rowval .= SparseArrays.rowvals(S)
    if !zerobased
        F.colptr .-= 1
        F.rowval .-= 1
    end

    F.nzval .= nonzeros(S)

    _umfpack_numeric!(F, tmp, umf_ctrl, umf_info, reuse_numeric = false)
end

function _solve!(x, lu::SuiteSparse.UMFPACK.UmfpackLU{T, Int}, b, umf_ctrl, umf_info) where T
    ccall((:umfpack_di_solve, :libumfpack), Int64,
        (Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Float64},
         Ptr{Float64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Float64},
         Ptr{Float64}),
        0, lu.colptr, lu.rowval, lu.nzval,
        x, b, lu.numeric, umf_ctrl,
        umf_info)
    return nothing
end

# # test
# n = length(p.newton.res.r)
# S = copy(p.newton.jac.R)#sprand(n, n, 0.1)
# b = rand(n)
# b
# xsol1 = S \ b
# solver = lu_sparse_solver(S)
# xsol2 = zeros(n)
# linear_solve!(solver, xsol2, S, b)

# for i = 1:10
#     @show err = b - S * xsol1
#     Δ = S \ err
#     xsol1 += Δ
# end
# # for i = 1:10
# #     @show err = b - S * xsol2
# #     Δ = S \ err
# #     xsol2 += Δ
# # end
# norm(xsol1 - xsol2, Inf)

# @benchmark linear_solve!($solver, $xsol2, $S, $b)
# @benchmark factorize!($solver, $S)

# # UMFPACK_INFO = 90
# # UMFPACK_CONTROL = 20
# # umf_ctrl = Vector{Float64}(undef, UMFPACK_CONTROL)
# # umf_info = Vector{Float64}(undef, UMFPACK_INFO)
# # tmp = Vector{Ptr{Cvoid}}(undef, 1)

# # # function umfpack_numeric!(U, tmp, umf_ctrl, umf_info)
# # #     status = ccall((:umfpack_di_numeric, :libumfpack), Int64,
# # #                    (Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Cvoid}, Ptr{Cvoid},
# # #                     Ptr{Float64}, Ptr{Float64}),
# # #                    U.colptr, U.rowval, U.nzval, U.symbolic, tmp,
# # #                    umf_ctrl, umf_info)
# # #     return nothing
# # # end

# # # @benchmark umfpack_numeric!($F, $tmp, $umf_ctrl, $umf_info)
# # # F.colptr
# # # typeof(F)
# # S
# # function sparse_lu!(F::SuiteSparse.UMFPACK.UmfpackLU{T,I}, S::SparseMatrixCSC{T,I}, tmp, umf_ctrl, umf_info) where {T,I}
# #     F.m = size(S, 1)
# #     F.n = size(S, 2)
# #     F.colptr = SparseArrays.getcolptr(S)
# #     F.rowval = SparseArrays.rowvals(S)
# #     F.nzval .= S.nzval
# #     umfpack_numeric!(F, tmp, umf_ctrl, umf_info)
# # end

# my_lu!(F, S, tmp, umf_ctrl, umf_info)
# xsol2 = F \ b
# @benchmark my_lu!($F, $S, $tmp, $umf_ctrl, $umf_info)

# # F.m = size(S, 1)
# # F.n = size(S, 2)
# # F.colptr .= SparseArrays.getcolptr(S)# zerobased ? copy(getcolptr(S)) : decrement(getcolptr(S))
# # F.rowval .= SparseArrays.rowvals(S)#zerobased ? copy(rowvals(S)) : decrement(rowvals(S))
# # F.nzval .= S.nzval#copy(nonzeros(S))
# # umfpack_numeric!(F, reuse_numeric = false)

# # umfpack_numeric!(F, reuse_numeric = false)
# # check && (issuccess(F) || throw(LinearAlgebra.SingularException(0)))

# # function lu!(F::UmfpackLU, S::SparseMatrixCSC{<:UMFVTypes,<:UMFITypes}; check::Bool=true)
# #     zerobased = getcolptr(S)[1] == 0
# #     F.m = size(S, 1)
# #     F.n = size(S, 2)
# #     F.colptr = zerobased ? copy(getcolptr(S)) : decrement(getcolptr(S))
# #     F.rowval = zerobased ? copy(rowvals(S)) : decrement(rowvals(S))
# #     F.nzval = copy(nonzeros(S))

# #     umfpack_numeric!(F, reuse_numeric = false)
# #     check && (issuccess(F) || throw(LinearAlgebra.SingularException(0)))
# #     return F
# # end

# using SparseArrays, SuiteSparse, LinearAlgebra, BenchmarkTools

# const UMFPACK_INFO = 90
# const UMFPACK_CONTROL = 20
# const umf_ctrl = Vector{Float64}(undef, UMFPACK_CONTROL)
# const umf_info = Vector{Float64}(undef, UMFPACK_INFO)
# tmp = Vector{Ptr{Cvoid}}(undef, 1)

# n = 350
# S = sprand(n, n, 0.4)
# b = rand(n)

# xsol1 = S \ b

# F = lu(copy(S))

# # @benchmark SuiteSparse.UMFPACK.lu!($F, $S)

# _umfpack_numeric!(F, tmp, reuse_numeric=false)

# _lu!(F, S, tmp)
# @benchmark _lu!($F, $S, $tmp)

# xsol2 = F \ b
# norm(xsol1 - xsol2)

# xsol3 = zeros(n)
# # ldiv!(xsol3, F, b)
# # @benchmark ldiv!($xsol3, $F, $b)
# # norm(xsol3 - xsol1)

# sparse_lu_solve!(xsol3, F, b)

# norm(xsol3 - xsol1, Inf)


text:
using Base.Test
using CUDArt
using Knet
include(Pkg.dir("Knet/test/mnist.jl"))
sparse32{T}(a::Array{T})=convert(SparseMatrixCSC{T,Int32}, a)

xtrn = MNIST.xtrn
xtst = MNIST.xtst
ytrn = MNIST.ytrn
ytst = MNIST.ytst
w0 = similar(ytst, size(ytst,1), 0)
d0 = 6.0
c0 = 1.0
g0 = 0.1
niter = 100
nbatch = 128
net = nothing
nc = size(ytrn,1)

for ker in (
            (:kpoly, [c0,d0]),
            (:kgauss, [g0]),
            (:klinear, nothing),
            (:perceptron, nothing),
            )
    aa = ss = nothing
    for prc in (
                :double,
                :single,
                )
        for fmt in (
                    :sparse,
                    :dense,
                    )
            for loc in (
                        :cpu,
                        :gpu,
                        )
                # loc == :gpu && ker[1] == :perceptron && continue
                println("\n$ker, $prc, $fmt, $loc")
                Knet.gpu(loc == :gpu)
                for p in (:xtrn, :xtst, :ytrn, :ytst); @eval $p=copy(MNIST.$p); end
                prc == :double && (for p in (:xtrn, :xtst, :ytrn, :ytst); @eval $p=convert(Array{Float64},$p); end)
                fmt == :sparse && (for p in (:xtrn, :xtst); @eval $p=sparse32($p); end)

                xtrn,ytrn=xtst,ytst # DBG: For quick results

                net = (ker[1] == :perceptron ?
                       Op[Mmul(nc;average=true,init=initzero), PercLoss()] :
                       Op[KPerceptron(nc, Knet.(ker[1]), ker[2])])
                gc(); @date train(net, xtrn, ytrn; iters=niter,batch=nbatch)
                @date a = accuracy(ytst, predict(net, xtst))
                s = isdefined(net[1],:s) ? size(net[1].s) : 0
                @show (s,a)
                aa == nothing && (aa=a)
                ss == nothing && (ss=s)
                @test aa == a
                @test ss == s
            end
        end
    end
end

# with naive hcat!:

# (:klinear0,nothing), cpu, dense
# 2015-06-23T03:19:20 train(net,xtrn,ytrn; iters=niter,batch=nbatch)
# elapsed time: 0.780812858 seconds (789688232 bytes allocated, 42.90% gc time)
# ((:klinear0,nothing),:cpu,:dense,(784,2926),0.877)
# elapsed time: 0.366831263 seconds (120308196 bytes allocated, 11.45% gc time)

# (:klinear0,nothing), cpu, sparse
# 2015-06-23T03:19:22 train(net,xtrn,ytrn; iters=niter,batch=nbatch)
# elapsed time: 4.665061582 seconds (1148766840 bytes allocated, 10.10% gc time)
# ((:klinear0,nothing),:cpu,:sparse,(784,2926),0.877)
# elapsed time: 5.582400282 seconds (864266820 bytes allocated, 6.32% gc time)

# with sparse hcat!:

# (:klinear0,nothing), cpu, sparse
# 2015-06-23T03:22:54 train(net,xtrn,ytrn; iters=niter,batch=nbatch)
# elapsed time: 4.136324843 seconds (834224992 bytes allocated, 8.03% gc time)
# ((:klinear0,nothing),:cpu,:sparse,(784,2926),0.877)
# elapsed time: 5.153346267 seconds (864266820 bytes allocated, 6.83% gc time)

# if false; info("KPerceptron+klinear2")
# for i=1:2
# ftrn = sparse(MNIST.xtrn)
# ftst = sparse(MNIST.xtst)
# @show fnet = Op[KPerceptron(10, Knet.klinear0, [0f0])]
#     gc()
#     @date train(fnet, ftrn, ytrn; iters=niter,batch=nbatch)
#     gc()
#     @time println((i, size(fnet[1].s),
#                    accuracy(ytst, predict(fnet, ftst)),
#                    # accuracy(ytrn, predict(lnet, ltrn)),
#                    ))
# end; end

# if false; info("KPerceptron+klinear")
# for i=1:2
# ltrn = sparse(MNIST.xtrn)
# ltst = sparse(MNIST.xtst)
# @show lnet = Op[KPerceptron(10, Knet.klinear, [0f0])]
# gc()
#     @date train(lnet, ltrn, ytrn; iters=niter,batch=nbatch)
# gc()
#     @time println((i, size(lnet[1].s),
#                    accuracy(ytst, predict(lnet, ltst)),
#                    # accuracy(ytrn, predict(lnet, ltrn)),
#                    ))
# end; end

# if false; info("KPerceptron+kgauss")
# for i=1:2
# qtrn = sparse(MNIST.xtrn)
# qtst = sparse(MNIST.xtst)
# @show qnet = Op[KPerceptron(10, Knet.kgauss, [g0])]
#     @date train(qnet, qtrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(qnet[1].s),
#                    accuracy(ytst, predict(qnet, qtst)),
#                    # accuracy(ytrn, predict(qnet, qtrn)),
#              ))
# end; end

# if false; info("KPerceptron+kgauss0")
# for i=1:2
# qtrn = sparse(MNIST.xtrn)
# qtst = sparse(MNIST.xtst)
# @show qnet = Op[KPerceptron(10, Knet.kgauss0, [g0])]
#     @date train(qnet, qtrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(qnet[1].s),
#                    accuracy(ytst, predict(qnet, qtst)),
#                    # accuracy(ytrn, predict(qnet, qtrn)),
#              ))
# end; end

# if false; info("KPerceptron+kpoly")
# for i=1:2
# ktrn = sparse(MNIST.xtrn)
# ktst = sparse(MNIST.xtst)
# @show knet = Op[KPerceptron(10, Knet.kpoly, [c0,d0])]
#     @date train(knet, ktrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(knet[1].s),
#                    accuracy(ytst, predict(knet, ktst)),
#                    # accuracy(ytrn, predict(knet, ktrn)),
#                    ))
# end; end

# if false; info("KPerceptron+kpoly0")
# for i=1:2
# ktrn = sparse(MNIST.xtrn)
# ktst = sparse(MNIST.xtst)
# @show knet = Op[KPerceptron(10, Knet.kpoly, [c0,d0])]
#     @date train(knet, ktrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(knet[1].s),
#                    accuracy(ytst, predict(knet, ktst)),
#                    # accuracy(ytrn, predict(knet, ktrn)),
#                    ))
# end; end

# if false; info("Perceptron (sparse)")
# ptrn = sparse(MNIST.xtrn)
# ptst = sparse(MNIST.xtst)
# @show pnet = Op[Perceptron(10)]
# for i=1:1
#     @date train(pnet, ptrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i,
#                    accuracy(ytst, predict(pnet, ptst)),
#                    # accuracy(ytrn, predict(pnet, ptrn)),
#                    ))
# end; end

# if false; info("Perceptron (dense)")
# dtrn = copy(MNIST.xtrn)
# dtst = copy(MNIST.xtst)
# @show dnet = Op[Perceptron(10)]
# for i=1:1
#     @date train(dnet, dtrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i,
#                    accuracy(ytst, predict(dnet, dtst)),
#                    # accuracy(ytrn, predict(dnet, dtrn)),
#                    ))
# end; end

# if false; info("Poly+PercLoss")
# strn = sparse(MNIST.xtrn)
# stst = sparse(MNIST.xtst)
# @show snet = Op[Poly(10;c=c0,d=d0,w=w0), PercLoss()]
# for i=1:1
#     @date train(snet, strn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(snet[1].s),
#                    accuracy(ytst, predict(snet, stst)),
#                    # accuracy(ytrn, predict(snet, strn)),
#                    ))
# end; end

# if false
# @show knet[1].s == snet[1].s
# @show knet[1].w0 == snet[1].v
# @show knet[1].w2 == snet[1].w
# end # if false

# if false; info("Rbfk+PercLoss")
# rtrn = sparse(MNIST.xtrn)
# rtst = sparse(MNIST.xtst)
# @show rnet = Op[Rbfk(10;gamma=0.1), PercLoss()]
# for i=1:1
#     @date train(rnet, rtrn, ytrn; iters=niter,batch=nbatch)
#     println((i, size(rnet[1].s),
#              accuracy(ytst, predict(rnet, rtst)),
#              # accuracy(ytrn, predict(rnet, rtrn)),
#              ))
# end; end

# if false; info("KPerceptron+kgauss")
# qtrn = sparse(MNIST.xtrn)
# qtst = sparse(MNIST.xtst)
# @show qnet = Op[KPerceptron(10, Knet.kgauss, [g0])]
# for i=1:1
#     @date train(qnet, qtrn, ytrn; iters=niter,batch=nbatch)
#     @time println((i, size(qnet[1].s),
#                    accuracy(ytst, predict(qnet, qtst)),
#                    # accuracy(ytrn, predict(qnet, qtrn)),
#              ))
# end; end

# if false
# @show qnet[1].s == rnet[1].s
# @show qnet[1].w0 == rnet[1].v
# @show qnet[1].w2 == rnet[1].w
# end # if false

# if false; info("Poly+PercLoss+dense")
# xnet = Op[Poly(c=c0,d=d0,w=w0), PercLoss()]
# for i=1:1
#     @time train(xnet, xtrn, ytrn; iters=niter,batch=nbatch)
#     println((i, size(xnet[1].s),
#              accuracy(ytst, predict(xnet, xtst)),
#              )) # accuracy(ytrn, predict(xnet, xtrn))))
# end; end

# if false; info("Poly+PercLoss+GPU+dense") # mmul, hcat, ctranspose do not work
# Knet.gpu(true)
# cnet = Op[Poly(c=c0,d=d0,w=CudaArray(w0)), PercLoss()]
# for i=1:1
#     @time train(cnet, xtrn, ytrn; iters=niter,batch=nbatch)
#     println((i, size(cnet[1].s),
#              accuracy(ytst, predict(cnet, xtst)),
#              )) # accuracy(ytrn, predict(cnet, xtrn))))
# end; end


text:
# TODO: review and rewrite this

using NBodyIPs, JuLIP, Test, IPFitting
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB, observation, vec_obs
using NBodyIPs: BondLengthDesc, nbpolys
const Lsq = IPFitting.Lsq
using LinearAlgebra: qr, norm

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Cu)
calc = let r0=r0
   LennardJones(r0=r0) * C2Shift(2.5*r0)
end
data = generate_data(:Cu, 3, 0.33*r0, 50, calc)
rcut2 = cutoff(calc)
D2 = BondLengthDesc("exp( - 2 * (r/$r0 - 1) )", (:cos, rcut2-1, rcut2))
B2 = nbpolys(2, D2, 6)

# generate a database so that we generate the row indices
db = LsqDB("", B2, data)

## generate a LSQ system by hand
atlen = length(data[1])
len = 1 + length(observation(data[1], "F"))
Y = zeros(length(data) * len)
Ψ = zeros(length(data) * len, length(B2))
for (okey, evalb, w) in zip(["E", "F"], (energy, forces), (1.345/atlen, 2.987))
   for d in data
      o = observation(d, okey)
      irows = IPFitting.DB.matrows(d, okey)
      leno = length(o)
      Y[irows] = w * o
      for colidx = 1:length(B2)
         Ψ[irows, colidx] = w * vec_obs(okey, evalb(B2[colidx], d.at))
      end
   end
end

## test that Ψ_man, F_man are equivalent to a system generated from a db
Ψ_man, Y_man = Ψ, Y
Ψ, Y = IPFitting.Lsq.get_lsq_system( db,
               E0 = 0.0,
               weights = Dict("default" => Dict("E" => 1.345, "F" => 2.987)),
            )

##
println(@test Y ≈ Y_man)
println(@test Ψ ≈ Ψ_man)


text:
module TestExceptions

using PDFHighlights
using Test

# Print the header
println("\e[1;32mRUNNING\e[0m: TestExceptions.jl")

PDFHighlights.Internal.Exceptions.@exception(
    exception_with_file,
    file::String,
)

@testset "@exception_with_file" begin

    @exception_with_file(
        HelloThere,
        "Exception thrown when Obi-Wan Kenobi joins the server.",
        "Hello there, ", e.file, "!"
    )

    try
        throw(HelloThere("Obi-Wan"))
    catch e
        @test sprint(showerror, e) == """

        Main.TestExceptions.HelloThere:
        Hello there, Obi-Wan!
        """
    end

end

PDFHighlights.Internal.Exceptions.@exception(
    exception_with_line,
    file::String,
    line::Int,
)

@testset "@exception_with_line" begin

    @exception_with_line(
        ObjectsDestroyed,
        "Exception thrown when another god damn Death Star was destroyed.",
        e.file, " destroyed: ", e.line, ".",
    )

    try
        throw(ObjectsDestroyed("Death Stars", 2))
    catch e
        @test sprint(showerror, e) == """

        Main.TestExceptions.ObjectsDestroyed:
        Death Stars destroyed: 2.
        """
    end

end

end


text:
module Example
  # split the model, data, inference, and plotting into separate scripts
  include("model.jl")
  include("data.jl")
  include("inference.jl")
  include("plot.jl")
end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
sets=(
    largest_clique_number,largest_cliques_set
)
for f in sets
    suite[string(f)] = BenchmarkGroup(["set"])
    suite[string(f)]["graphs"] = BenchmarkGroup()
    for (name,g) in GRAPHS
        suite[string(f)]["graphs"][name] = @benchmarkable $(f)($g)
    end
end