text:
import DifferentialEquations;	const DE = DifferentialEquations
using LinearAlgebra
using StaticArrays
using Combinatorics
using FastClosures

"""
	g(t, tf, coeffs, basis_functions)
Correction pulse
# Arguments
- `t::Real or t::Array{<:Real,1}`: time.
- `tf::Real`: total time of the gate.
- `coeffs::Array{<:Real, 1}`: coefficients of the correction pulse.
- `basis_functions::`: .
"""
function g(t::Real, tf, coeffs, basis)
	basis_vector = basis.vector
	basis.function!(basis_vector, t)
	cost, sint = cos(t), sin(t)
	@inbounds gt = (coeffs[1]*cost + coeffs[2]*sint)*basis_vector[1]
	for i = 2: length(basis_vector)
		@inbounds gt += (coeffs[2*i - 1]*cost + coeffs[2*i]*sint)*basis_vector[i]
	end
	return gt
end

function g_envelopes(t, tf, coeffs, basis)
	basis_vector = basis.vector
	basis.function!(basis_vector, t)
	gt = zeros(2)
	for i = 1: length(basis_vector)
		gt[1] += coeffs[2*i - 1]*basis_vector[i]
		gt[2] += coeffs[2*i]*basis_vector[i]
	end
	return gt
end

"""
	Δ(t, tf, coeffs, basis_functions)
Detuning of the pulse frequency with respect to the transition frequency between the zero'th
and first levels.
# Arguments
- `t::Real`: time.
- `tf::Real`: total time of the gate.
- `coeffs::Array{<:Real, 1}`: coefficients of the detuning.
- `basis_functions::`: .
"""
function Δ(t::Real, tf, coeffs, basis)
	basis_vector = basis.vector
	basis.function!(basis_vector, t)
	@inbounds Δt = coeffs[1]*basis_vector[1]
	for i = 2: length(basis_vector)
		@inbounds Δt += coeffs[i]*basis_vector[i]
	end
	return Δt
end

"""
	Wharmonic_I(t, tf, α, η, coeffs, Δ_basis, pulse_basis)
The harmonic correction Hamiltonian in the interaction picture.
# Arguments
- `t::Real`: time.
- `tf::Real`: total time of the gate.
- `α::Real`: the anarmonicity in units of the transition frequency between the zeroth
and first energy levels.
- `η::Real`: parameter that characterizes the strength of the 1<->2 transition relative
to the 0<->1 transition.
- `coeffs::Array{Array{<:Real, 1}, 1}`: The first and second elements are the arrays with the
Δ coefficients and pulse coefficients, respectively.
- `Δ_basis::`: .
- `pulse_basis::`: .
"""
function Wharmonic_I(t, tf, α, η, c, coeffs, Δ_basis, pulse_basis)
	ct = c(t)
	cosc, sinc = cos(ct), sin(ct)
	cos2c, sin2c = cos(2*ct), sin(2*ct)
	φ = (1 + α)*t
	cosφ, sinφ = cos(φ), sin(φ)
	sint = sin(t)

	Δt = Δ(t, tf, coeffs[1], Δ_basis)/2
	gt = g(t, tf, coeffs[2], pulse_basis)
	ηg = η*gt

	WIt = @SVector [gt*cos(t),
				    -gt*sint*cos2c + Δt*sin2c,
				    gt*sint*sin2c + Δt*cos2c,
				    ηg*cosφ*cosc,
				    ηg*sinφ*cosc,
				    ηg*sinφ*sinc,
				    -ηg*cosφ*sinc,
				    -sqrt(3.0)*Δt,
				    Δt]

	return WIt
end
"""
	Wharmonic(t, tf, α, η, coeffs, Δ_basis, pulse_basis)
The harmonic correction Hamiltonian.
# Arguments
- `t::Real`: time.
- `tf::Real`: total time of the gate.
- `α::Real`: the anarmonicity in units of the transition frequency between the zeroth
and first energy levels.
- `η::Real`: parameter that characterizes the strength of the 1<->2 transition relative
to the 0<->1 transition.
- `coeffs::Array{Array{<:Real, 1}, 1}`: The first and second elements are the arrays with the
Δ coefficients and pulse coefficients, respectively.
- `Δ_basis::Function`: Detuning basis functions as a function of t.
- `pulse_basis::Function`: Pulse basis functions as a function of t.
"""
function Wharmonic(t, tf, η, coeffs, Δ_basis, pulse_basis)
	Δt = Δ(t, tf, coeffs[1], Δ_basis)/2
	gt = g(t, tf, coeffs[2], pulse_basis)
	cost, sint = cos(t), sin(t)
	ηg = η*gt

	Wt = @SVector [gt*cost,
				   -gt*sint,
				   Δt,
				   ηg*cost,
				   ηg*sint,
				   0,
				   0,
				   -sqrt(3.0)*Δt,
				   Δt]
	return Wt
end
# correction matrix obtained from the selected envelopes
function du!(du, u, args, t)
	tf = args[1]
	α = args[2]
	η = args[3]
	c = args[4]
	coeffs = args[5]
	Δ_basis = args[6]
	pulse_basis = args[7]

	W0t = Wharmonic_I(t, tf, α, η, c, coeffs, Δ_basis, pulse_basis)
	du .= view(W0t, 1: 7)
	return nothing
end
"""
	correction_matrix(tf, α, η, Δ_basis, pulse_basis; mask=nothing)
Return the correction matrix.
# Arguments
- `tf::Real`: total time of the gate.
- `α::Real`: anarmonicity.
- `η::Real`: parameter that characterizes the strength of the 1<->2 transition relative
to the 0<->1 transition.
- `coeffs::Array{Array{<:Real, 1}, 1}`: The first and second elements are the arrays with the
Δ coefficients and pulse coefficients, respectively.
- `Δ_basis::Function`: detuning basis functions as a function of t.
- `pulse_basis::Function`: pulse basis functions as a function of t.
- `mask::Array{Bool, 1}`: boolean array with a selection of the coefficients.
"""
function correction_matrix(system_params, Δ_basis, pulse_basis)
	α = system_params["alpha"]
    η = system_params["eta"]
    tf = system_params["tf"]
	pulse_integral = system_params["pulse_integral"]

	reltol, abstol = 1e-8, 1e-13
	npulse_basis = length(pulse_basis.vector)
	nΔ_basis = length(Δ_basis.vector)
	M = zeros(7, nΔ_basis + 2*npulse_basis)

	v = zeros(nΔ_basis + 2*npulse_basis)
	for i = 1: nΔ_basis + 2*npulse_basis
		v .= 0.0
		v[i] = 1.0
		coeffs = (v[1: nΔ_basis], v[nΔ_basis + 1: end])

		u0 = zeros(7)
		tspan = (zero(tf), tf)
		args = (tf, α, η, pulse_integral, coeffs, Δ_basis, pulse_basis)
		prob = DE.ODEProblem(du!, u0, tspan, args)
		sol = DE.solve(prob, DE.Vern9(), reltol=reltol, abstol=abstol, saveat=[tf])

		M[:, i] .= sol(tf)
	end
	return M
end
"""
	coefficients_W(Omega, M, Δ_basis, pulse_basis)
Return the coefficients that solve the linear system `M*coeffs = Omega`.
# Arguments
- `Omega::Array{<:Real, 1}`: the error vector, obtained via the Magnus expansion.
- `M::Array{<:Real, 2}`: correction matrix.
- `Δ_basis::Function`: detuning basis functions as a function of t.
- `pulse_basis::Function`: pulse basis functions as a function of t.
"""
function coefficients_W(Omega, M, Δ_basis, pulse_basis)
	ih_vec = Omega[1: 7]
	nΔ_basis = length(Δ_basis.vector)
	coeffs = pinv(M)*ih_vec

	return (coeffs[1: nΔ_basis], coeffs[nΔ_basis + 1: end])
end
"""
	get_coeffs(M, order, system_params, Δ_basis, pulse_basis)
Return the sum of correction coefficients up to order `order`.
# Arguments
- `M::Array{<:Real, 2}`: correction matrix.
- `order::Integer`: the order of the Magnus correction.
- `system_params::Dictionary`: dictionary with the parameters of the system.
- `Δ_basis::Function`: detuning basis functions as a function of t.
- `pulse_basis::Function`: pulse basis functions as a function of t.
"""
function get_coeffs(M, order, system_params, Δ_basis, pulse_basis; max_norm=nothing)
    nΔ_basis = length(Δ_basis.vector)
    npulse_basis = length(pulse_basis.vector)
    coeffs = (zeros(nΔ_basis), zeros(2*npulse_basis))
    for k = 1: order
        coeffs2 = get_coeffs(M, coeffs, k, system_params, Δ_basis, pulse_basis)
		#coeffs2 = get_coeffs(M, coeffs, order, system_params, Δ_basis, pulse_basis)
        coeffs = coeffs .+ coeffs2
		if max_norm != nothing
			if norm(vcat(coeffs...)) > max_norm
				return (zeros(nΔ_basis), zeros(2*npulse_basis))
			end
		end
    end
    return coeffs
end
"""
	get_coeffs(M, order, system_params, Δ_basis, pulse_basis)
Return the correction coefficients of order `order`.
# Arguments
- `M::Array{<:Real, 2}`: correction matrix.
- `coeffs::Tuple{Array{<:Real, 1}, Array{<:Real, 1}}`: correction coefficients.
- `order::Integer`: the order of the Magnus correction.
- `system_params::Dictionary`: dictionary with the parameters of the system.
- `Δ_basis::Function`: detuning basis functions as a function of t.
- `pulse_basis::Function`: pulse basis functions as a function of t.
"""
function get_coeffs(M, coeffs, order, system_params, Δ_basis, pulse_basis)
	α = system_params["alpha"]
    η = system_params["eta"]
    tf = system_params["tf"]
	pulse = system_params["pulse"]
	pulse_integral = system_params["pulse_integral"]

	A = @closure (
		t-> - V_I(t, α, η, pulse, pulse_integral)
			- Wharmonic_I(t, tf, α, η, pulse_integral, coeffs, Δ_basis, pulse_basis))

    Ndims = 9
    Omegaf = solveOmega(tf, order, A, comm!, Ndims)
    vec = sum(Omegaf, dims=2)
    coeffs = coefficients_W(vec, M, Δ_basis, pulse_basis)
    return coeffs
end

function find_best_coeffs(ncoeffs, order, system_params, Δ_basis, pulse_basis, θ, λ)
    nΔ_basis = length(Δ_basis.vector)
    npulse_basis = length(pulse_basis.vector)
    if ncoeffs > nΔ_basis + 2*npulse_basis
        error("ncoeffs cannot be larger than ", nΔ_basis + 2*npulse_basis, ".")
    end

    α = system_params["alpha"]
    η = system_params["eta"]
    tf = system_params["tf"]
    # find the correction matrix
    #M = correction_matrix(tf, α, η, Δ_basis, pulse_basis)
	M = correction_matrix(system_params, Δ_basis, pulse_basis)
    # Ideal evolution operator
    U_ideal = U(system_params, λ; ideal_dynamics=true)
    # mask vector with the selected `ncoeffs` correction coefficients.
    coeffs_mask = trues(nΔ_basis + 2*npulse_basis)
    coeffs_mask[1: ncoeffs] = falses(ncoeffs)
    # all possible permutations of `ncoeffs` coefficients
    permutations = multiset_permutations(coeffs_mask, length(coeffs_mask))
    permutations = collect(permutations)
	println("There are ", length(permutations), " possible combinations of coefficients. "*
			"Starting to calculate ...")
    # select the best subgroup of `ncoeffs` coefficients
    coeffs_opt = (zeros(nΔ_basis), zeros(2*npulse_basis))
    Fopt = 0.0
    for i = 1: length(permutations)
        # Apply mask to the matrix M and obtain coefficients of the ith permutation
        Mp = copy(M)
        Mp[:, permutations[i]] .= 0.0
        coeffs = get_coeffs(Mp, order, system_params, Δ_basis, pulse_basis;
							max_norm= 2*θ/tf)

        Ucorrected = U(system_params, λ, coeffs, Δ_basis, pulse_basis)
        Fcorrected = fidelity(Ucorrected, U_ideal, 2)
        if Fcorrected > Fopt
            Fopt = Fcorrected
            coeffs_opt = coeffs
        end
    end
    return coeffs_opt
end
"""
	U(system_params, coeffs, λ, Δ_basis, pulse_basis)
Corrected evolution operator.
# Arguments
- `system_params::Dictionary`: dictionary with the parameters of the system.
- `λ::Array{Array{<:Complex, 2}}`: Array with the Gell-Mann
- `coeffs::Tuple{Array{<:Real, 1},Array{<:Real, 1}}`: correction coefficients.
- `Δ_basis::Function`: detuning basis functions as a function of t.
- `pulse_basis::Function`: pulse basis functions as a function of t.
"""
function U(system_params, λ, coeffs, Δ_basis, pulse_basis)
    α = system_params["alpha"]
    η = system_params["eta"]
    tf = system_params["tf"]
	pulse = system_params["pulse"]

	H_corrected = @closure t-> (H0(t, α, η, pulse)
						  		+ Wharmonic(t, tf, η, coeffs, Δ_basis, pulse_basis))

    sol = solveU(H_corrected, tf, λ)
    return sol(tf)
end

function solveρ(system_params, λ, coeffs, Δ_basis, pulse_basis, c1, c2, ρ0)
	tf = system_params["tf"]
	α = system_params["alpha"]
	η = system_params["eta"]
	pulse = system_params["pulse"]

	H = @closure t-> (H0(t, α, η, pulse)
					  + Wharmonic(t, tf, η, coeffs, Δ_basis, pulse_basis))

	tspan = (zero(tf), tf)
	args = (H, c1, c2, λ)
	prob = DE.ODEProblem(dρ, ρ0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-14)
	if sol.retcode == :Success
		return sol
	else
		error("Could not solve the master equation.")
	end
end

function average_fidelity(system_params, λ, coeffs, Δ_basis, pulse_basis, c1, c2)
	tf = system_params["tf"]
	Uideal = U(system_params, λ; ideal_dynamics = true)

	dict = Dict("+x"=> (1.0, 0.0, 0.0), "-x"=> (-1.0, 0.0, 0.0),
				"+y"=> (0.0, 1.0, 0.0), "-y"=> (0.0, -1.0, 0.0),
				"+z"=> (0.0, 0.0, 1.0), "-z"=> (0.0, 0.0, -1.0))
	F = 0.0
	for tag in ("+x", "-x", "+y", "-y", "+z", "-z")
		a = dict[tag]
		ρ0 = ρbloch(a)
		ρf = solveρ(system_params, λ, coeffs, Δ_basis, pulse_basis, c1, c2, ρ0)
		ρf = ρf(tf)
		ρfideal = Uideal*ρ0*adjoint(Uideal)
		F += tr(ρfideal*ρf)/6
	end
	if imag(F)>10.0^(-10.0)
		error("Complex fidelity!")
	else
		return real(F)
	end
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
using Documenter, ITensorBenchmarks

#DocMeta.setdocmeta!(ITensors, :DocTestSetup, :(using ITensors); recursive=true)

sitename = "ITensorBenchmarks"

settings = Dict(
  :modules => [ITensorBenchmarks],
  :pages => [
    "Introduction" => "index.md",
    "TeNPy and ITensor Comparisons" => "tenpy_itensor/index.md",
    #"Getting Started with ITensor" => [
    #  "Installing Julia and ITensor" => "getting_started/Installing.md",
    #  "Running ITensor and Julia Codes" => "getting_started/RunningCodes.md",
    #  "Next Steps" => "getting_started/NextSteps.md",
    #],
  ],
  :format => Documenter.HTML(; assets=["assets/favicon.ico"], prettyurls=false),
  :doctest => true,
  :checkdocs => :none,
)

ext: jl


text:
function byrow(y)
	# data by rows
	# if vector, converts to a rowmatrix

	if ndims(y) == 1
		# vector
		nt = length(y)
		ny = 1
		y = reshape(y,1,nt)
	else
		# matrix
		(ny,nt) = size(y)

		if nt < ny
			y=y'
			(ny,nt) = size(y)
		end
	end

	return y,ny,nt
end


text:
using Base.Iterators
using Test

using Arbitrary

const alltypes = [Nothing, Bool, Char,
                  Int8, Int16, Int32, Int64, Int128,
                  UInt8, UInt16, UInt32, UInt64, UInt128,
                  Float16, Float32, Float64,
                  BigInt, BigFloat,
                  Rational{BigInt},
                  Tuple{},
                  Tuple{Int8, Int16},
                  Tuple{Int8, Int16, Int32},
                  Tuple{Int8, Int16, Int32, Int64}]
hasonevalue(T) = T in [Nothing, Tuple{}]

const Array0{T} = Array{T, 0}
const Array1{T} = Array{T, 1}
const Array2{T} = Array{T, 2}
const allcontainers = [nothing, Ref, Tuple, Array0, Array1, Array2]
hasfixedshape(C) = C in [nothing, Ref, Tuple, Array0]

myequal(x::T, y::T) where {T} = isequal(x, y)
myequal(x::Ref{T}, y::Ref{T}) where {T} = isequal(x[], y[])

for C in allcontainers
    for E in alltypes
        T = C === nothing ? E : C{E}
        @testset "Basic functionality for type $T" begin
            # Generate arbitrary values
            arb = arbitrary(T)
            values = collect(take(arb, 100))
            # Generate some other arbitrary values
            arb2 = arbitrary(T)
            values2 = collect(take(arb2, 100))
            # Ensure they are different
            @test all(myequal.(values, values2)) ==
                (hasonevalue(E) && hasfixedshape(C))
            # Generate values from a known RNG
            arb3 = arbitrary(T, UInt(42))
            arb4 = arbitrary(T, UInt(42))
            @test all(myequal.(collect(take(arb3, 100)),
                               collect(take(arb4, 100))))
        end
    end
end

# Floating-point numbers do NOT satisfy the usual arithmetic laws
const arithmetic_types = [Int8, Int16, Int32, Int64, Int128,
                          UInt8, UInt16, UInt32, UInt64, UInt128,
                          # Float16, Float32, Float64,
                          BigInt, BigFloat,
                          Rational{BigInt}]
const division_types = [# Float16, Float32, Float64,
                        # BigFloat,
                        Rational{BigInt}]
for T in arithmetic_types
    @testset "Arithmetic identities for type $T" begin
        # Generate arbitrary values
        xs = collect(take(arbitrary(T), 100))
        ys = collect(take(arbitrary(T), 100))
        zs = collect(take(arbitrary(T), 100))
        ds = collect(take(Iterators.filter(x -> !isequal(x, T(0)),
                                           arbitrary(T)), 100))
        # Addition:
        # Commutativity
        @test all(isequal.(xs .+ ys, ys .+ xs))
        # Associativity
        @test all(isequal.((xs .+ ys) .+ zs, xs .+ (ys .+ zs)))
        # Neutral element
        @test all(isequal.(xs .+ T(0), xs))
        @test all(isequal.(T(0) .+ xs, xs))
        # Inverse
        @test all(isequal.(xs .+ (-xs), T(0)))
        @test all(isequal.(xs .- ys, xs .+ (-ys)))
        # Multiplication:
        # Commutativity
        @test all(isequal.(xs .* ys, ys .* xs))
        # Associativity
        @test all(isequal.((xs .* ys) .* zs, xs .* (ys .* zs)))
        # Neutral element
        @test all(isequal.(xs .* T(1), xs))
        @test all(isequal.(T(1) .* xs, xs))
        # Inverse
        if T in division_types
            @test all(isequal.(ds .* inv.(ds), T(1)))
            @test all(isequal.(xs ./ ds, xs .* inv.(ds)))
        end
    end
end

const function_types = [Int, Float64, Char, UInt]

@testset "Function identities" begin
    T, U, V, R = function_types
    fs = collect(take(arbitrary(Fun{T, U}), 100))
    gs = collect(take(arbitrary(Fun{U, V}), 100))
    hs = collect(take(arbitrary(Fun{V, R}), 100))

    # Neutral element
    for f in fs
        @test fun_isequal(T, f ∘ identity, f)
        @test fun_isequal(T, identity ∘ f, f)
    end

    # Associativity
    for (f, g, h) in zip(fs, gs, hs)
        @test fun_isequal(T, (h ∘ g) ∘ f, h ∘ (g ∘ f))
    end
end


text:
include("imports.jl")
include("initializers.jl")
include("classical_chan_vese.jl")


text:
chanids() = json(get(ch_ids))["ids"]
chinfo(id) = json(get(string(ch_info,id)))
recids() = json(get(rec_ids))["ids"]
recinfo(id) = json(get(string(rec_info,id)))

function episodetable(ids::Vector)
    ans = Dict{AbstractString,Int}()
    for id in ids
        try
            info = recinfo(id)
            tt = string(info["meta"]["recSeries"]["jsonForClient"]["title"])
            ans[tt] = get(ans,tt,0) + 1
        catch
            next
        end
    end
    ans
end


text:
# device capability handling

# select the highest capability that is supported by both the toolchain and device
function supported_capability(dev::CuDevice)
    dev_cap = capability(dev)
    compat_caps = filter(cap -> cap <= dev_cap, target_support[])
    isempty(compat_caps) &&
        error("Device capability v$dev_cap not supported by available toolchain")

    return maximum(compat_caps)
end

# return the capability of the current context's device, or a sane fall-back
function current_capability()
    if CuCurrentContext() !== nothing
        return supported_capability(device())
    else
        # newer devices tend to support cleaner code (higher-level instructions, etc)
        # so target the most recent device as supported by this toolchain
        return maximum(target_support[])
    end
end


text:
####### Payoffs definition

### Spot
include("payoffs/spot.jl")

### European Payoffs
include("payoffs/forward.jl")
include("payoffs/european_option.jl")
include("payoffs/binary_european_option.jl")

### Barrier Payoffs
include("payoffs/barrier_do_option.jl")
include("payoffs/barrier_di_option.jl")
include("payoffs/barrier_uo_option.jl")
include("payoffs/barrier_ui_option.jl")
include("payoffs/double_barrier_option.jl")

### American Payoffs
include("payoffs/american_option.jl")
include("payoffs/binary_american_option.jl")

### Path Dependent Payoffs
include("payoffs/bermudan_option.jl")

### Asian Payoffs
include("payoffs/asian_fixed_strike_option.jl")
include("payoffs/asian_floating_strike_option.jl")

### Basket Payoffs
include("payoffs/basket/n_european_option.jl")


text:
#Tensors in Hierarchical Tucker format + functions
export htensor, randhtensor, cat, display, dropdims, full, hrank, htrunc, innerprod, isequal, minus, mtimes, ndims, plus, reorth, reorth!
export size, trten2mat, trten2ten, ttm, ttv, uminus, norm

"""
    htensor(tree,trten,fmat)
    htensor(tree,trten,U1,U2,...)

Hierarchical Tucker tensor.

## Arguments:
- `tree::dimtree`: Dimension tree.
- `trten::TensorCell`: Transfer tensors. One for each internal node of the tree.
- `fmat::MatrixCell`: Frame matrices. One for each leaf of the tree.
For htensor X, X.isorth=true if frame matrices are othonormal.
"""
mutable struct htensor
	tree::dimtree
	trten::TensorCell
  	fmat::MatrixCell
  	isorth::Bool
	function htensor(tree::dimtree,trten::TensorCell,fmat::MatrixCell,isorth::Bool)
    @assert(length(fmat)==length(tree.leaves),"Dimension mismatch.") #order of a tensor
    @assert(length(trten)==length(tree.internal_nodes),"Dimension mismatch.")
		for U in fmat
			if norm(U'*U-eye(size(U,2)))>(size(U,1)^2)*eps()
				isorth=false
			end
		end
		new(tree,trten,fmat,isorth)
	end
end
#htensor(tree::dimtree{D1,D2},trten::TensorCell{T1},fmat::MatrixCell{T2},isorth::Bool) where {D1<:Integer,D2<:Integer,T1<:Number,T2<:Number}=htensor{D1,D2,T1,T2}(tree,trten,fmat,isorth)
htensor(tree::dimtree,trten::TensorCell,fmat::MatrixCell)=htensor(tree,trten,fmat,true)
htensor(tree::dimtree,trten::TensorCell,mat::AbstractMatrix{<:Number})=htensor(tree,trten,collect(mat),true)
# function htensor(tree::dimtree,trten::TensorCell,mat...)
#   fmat=MatrixCell(undef,0)
#   for M in mat
# 			push!(fmat,M)
# 	end
# 	htensor(tree,trten,fmat,true)
# end
#If array of matrices and tensors aren't defined as TensorCell and MatrixCell, but as M=[M1,M2,...,Mn]:
htensor(tree::dimtree,trten::Array{Array{T1,3}},fmat::Array{Matrix{T2}},isorth::Bool) where {T1<:Number,T2<:Number}=htensor(tree,TensorCell(trten),MatrixCell(fmat),isorth)
htensor(tree::dimtree,trten::Array{Array{T1,3}},fmat::Array{Matrix{T2}}) where {T1<:Number,T2<:Number}=htensor(tree,TensorCell(trten),MatrixCell(fmat),true)

"""
    randhtensor(I::Vector,R::Vector)
    randhtensor(I::Vector,R::Vector,T::dimtree)

Create random htensor of size I, hierarchical rank R and balanced tree T. If T not specified use balanced tree.
"""
function randhtensor(sz::AbstractVector{<:Integer},rnk::AbstractVector{<:Integer},T::dimtree)
  N=length(sz)
  @assert(length(rnk)==non(T),"Size and rank dimension mismatch.")
  In=length(T.internal_nodes) #number of internal nodes
  cmps=TensorCell(undef,N+In)
  [cmps[T.leaves[n]]=rand(sz[n],rnk[T.leaves[n]]) for n=1:N];
  for n in reverse(T.internal_nodes)
    c=children(T,n)
    is_leaf(T,c[1]) ? lsz=size(cmps[c[1]],2) : lsz=size(cmps[c[1]],3)
    is_leaf(T,c[2]) ? rsz=size(cmps[c[2]],2) : rsz=size(cmps[c[2]],3)
    cmps[n]=rand(lsz,rsz,rnk[n])
  end
  htensor(T,cmps[T.internal_nodes],MatrixCell(cmps[T.leaves]))
end
function randhtensor(sz::AbstractVector{<:Integer},rnk::AbstractVector{<:Integer})
    T=dimtree(length(sz))
    randhtensor(sz,rnk,T)
end
function randhtensor(sz::AbstractVector{<:Integer},T::dimtree)
  N=length(sz)
  rnk=ones(Int,non(T))
  for l in T.leaves
    rnk[l]=sz[node2ind(T,l)]
  end
  for n in reverse(T.internal_nodes)[1:end-1]
    c=children(T,n)
    ind=setdiff(1:N,node2ind(T,c))
    rnk[n]=prod(sz[ind])
  end
  randhtensor(sz,rnk)
end
function randhtensor(sz::AbstractVector{<:Integer})
  T=dimtree(length(sz))
  randhtensor(sz,T)
end

"""
    cat(X,Y,n)

Concatente htensors X and Y by mode n.
"""
function cat(X1::htensor,X2::htensor,n::Integer)
    @assert(X1.tree == X2.tree, "htensors must have equal dimension trees.")
    @assert(0<n<=ndims(X1),"Invalid concatenation mode.")
    eqm=setdiff(1:ndims(X1),n) #equal modes
    @assert(size(X1)[eqm] == size(X2)[eqm],"Dimensions mismatch.")
    Xc=deepcopy(X1)
    rnk1=hrank(X1)
    rnk2=hrank(X2)
    rnk=rnk1+rnk2
    rnk[1]=1
    for l=1:ndims(X1)
        if l!=n
            Xc.fmat[l]=[X1.fmat[l] X2.fmat[l]]
        else
            Xc.fmat[l]=blockdiag(X1.fmat[l],X2.fmat[l])
        end
    end
    for i in X1.tree.internal_nodes
        ind=node2ind(X1.tree,i)
        c=children(X1.tree,i)
        Xc.trten[ind]=zeros(rnk[c[1]],rnk[c[2]],rnk[i])
        if i!=1
            Xc.trten[ind][1:rnk1[c[1]],1:rnk1[c[2]],1:rnk1[i]]=X1.trten[ind]
            Xc.trten[ind][rnk1[c[1]]+1:end,rnk1[c[2]]+1:end,rnk1[i]+1:end]=X2.trten[ind]
        else
            Xc.trten[ind][1:rnk1[c[1]],1:rnk1[c[2]]]=X1.trten[ind]
            Xc.trten[ind][rnk1[c[1]]+1:end,rnk1[c[2]]+1:end]=X2.trten[ind]
        end
    end
    Xc
end

#Display a htensor. **Documentation in ttensor.jl
function display(X::htensor,name="htensor")
  println("Hierarchical Tucker tensor of size ",size(X),":\n")
  println("$name.tree: ")
  flush(stdout)
  display(X.tree)
  for n=1:length(X.trten)
    println("\n\n$name.trten[$n]:")
	flush(stdout)
	display(X.trten[n])
  end
  for n=1:length(X.fmat)
    println("\n\n$name.fmat[$n]:")
	flush(stdout)
	display(X.fmat[n])
  end
end

#Makes full tensor out of a ktensor. **Documentation in ttensor.jl.
function full(X::htensor)
  t,tl,tr=structure(X.tree)
  In=length(X.tree.internal_nodes)
  V=MatrixCell(undef,In)
  j=In
  for i=In:-1:1
    if length(tr[i])==1
      Ur=X.fmat[tr[i]...]
    else
      Ur=V[j]
      j-=1;
    end
    if length(tl[i])==1
      Ul=X.fmat[tl[i]...]
    else
      Ul=V[j];
      j-=1;
    end
     V[i]=kron(Ur,Ul)*trten2mat(X.trten[i]);
     #V[i]=krontv(Ur,Ul,B[i]);
  end
  reshape(V[1],size(X))
end

"""
    fmat(X)

Hierarchical ranks of a htensor.
"""
function hrank(X::htensor)
  order=[X.tree.internal_nodes;X.tree.leaves]
  r1=[size(B,3) for B in X.trten]
  r2=[size(U,2) for U in X.fmat]
  [r1;r2][invperm(order)]
end

"""
    htrunc(X[,tree])

Truncate full tensor X into a htensor for a given tree. If tree not specified use balanced tree.
"""
function htrunc(X::Array{<:Number,N},tree::dimtree;method="svd",maxrank=[],atol=1e-8,rtol=0) where N
  @assert(N==length(tree.leaves),"Dimension mismatch.")
  maxrank=check_vector_input(maxrank,2*N-1,10);
  In=length(tree.internal_nodes)
  U=MatrixCell(undef,N)
	for n=1:N
    Xn=float(tenmat(X,n))
    U[n]=colspace(Xn,method=method,maxrank=maxrank[tree.leaves[n]],atol=atol,rtol=rtol) #[U[n]=svdfact(tenmat(X,n))[:U][1:maxrank[tree.leaves[n]]] for n=1:N]
	end
  B=TensorCell(undef,In)
  t,tl,tr=structure(tree)
  for i in reverse(tree.internal_nodes)
    ind=node2ind(tree,i)
    if i!=1
      Xt=tenmat(X,row=t[ind])
      Ut=colspace(Xt,method=method,maxrank=maxrank[ind],atol=atol,rtol=rtol) #Ut=svdfact(Xt)[:U][1:maxrank[ind]]
    else
      Ut=vec(X)
    end
    Xl=tenmat(X,row=tl[ind])
    Xr=tenmat(X,row=tr[ind])
    Ul=colspace(Xl,method=method,maxrank=maxrank[ind],atol=atol,rtol=rtol) #Ul=svdfact(Xl)[:U][1:maxrank[ind]]
    Ur=colspace(Xr,method=method,maxrank=maxrank[ind],atol=atol,rtol=rtol) #Ur=svdfact(Xr)[:U][1:maxrank[ind]]
    #B=krontv(Ur',Ul',Ut)
    #B[ind]=trten2ten(kron(Ur',Ul')*Ut,size(Ul,2),size(Ur,2))

    B[ind]=trten2ten(krontv(copy(Ur'),copy(Ul'),Ut),size(Ul,2),size(Ur,2))
  end
  htensor(tree,B,U)
end
htrunc(X::Array{<:Number,N};method="svd",maxrank=[],atol=1e-8,rtol=0) where N=htrunc(X,dimtree(ndims(X)),method=method,maxrank=maxrank,atol=atol,rtol=rtol)

#Inner product of two htensors. **Documentation in tensor.jl.
function innerprod(X1::htensor,X2::htensor)
	@assert(size(X1) == size(X2),"Dimension mismatch.")
  @assert(X1.tree == X2.tree, "Dimension trees must be equal.")
  N=ndims(X1)
  In=length(X1.tree.internal_nodes) #number of internal nodes
  M=MatrixCell(undef,N+In)
  [M[X1.tree.leaves[n]]=X1.fmat[n]'*X2.fmat[n] for n=1:N]
  for n in reverse(X1.tree.internal_nodes)
    c=children(X1.tree,n)
    M[n]=trten2mat(X1.trten[node2ind(X1.tree,n)])'*kron(M[c[2]],M[c[1]])*trten2mat(X2.trten[node2ind(X2.tree,n)])
  end
  M[1][1]
end

#True if htensors have equal components. **Documentation in ttensor.jl.
function isequal(X1::htensor,X2::htensor)
  if (X1.trten == X2.trten) && (X1.fmat == X2.fmat)
    true
  else
    false
  end
end
==(X1::htensor,X2::htensor)=isequal(X1,X2)

#Subtraction of two htensors. **Documentation in ttensor.jl
function minus(X1::htensor,X2::htensor)
  X1+(-1)*X2
end
-(X1::htensor,X2::htensor)=minus(X1,X2)

#Scalar times htensor. **Documentation in ttensor.jl
function mtimes(α::Number,X::htensor)
  B=deepcopy(X.trten)
  B[1]=α*B[1]
	htensor(X.tree,B,X.fmat)
end
*(α::Number,X::htensor)=mtimes(α,X)
*(X::htensor,α::Number)=*(α,X)

#Number of modes of a htensor. **Documentation in Base.
function ndims(X::htensor)
	length(X.fmat)
end

#Used by dropdims
function next_single_node(T::dimtree,dims::AbstractVector{<:Integer})
    #println("dims = $dims")
    nodes=T.leaves[dims]
    #println("nodes = $nodes")
    if isempty(nodes)
        return -1
    end
    parents=parent(T,nodes)
    #println("parents = $parents")
    if allunique(parents)
       return nodes[1]
    end
    for p in parents
        inds=findall(parents.==p)
        if length(inds)>1
            return nodes[inds[1]]
        end
    end
end
next_single_node(T::dimtree,dims::Integer)=next_single_node(X,[dims])

#Addition of two htensors. **Documentation in ttensor.jl
function plus(X1::htensor,X2::htensor)
  @assert(X1.tree == X2.tree,"Input mismatch.")
  fmat=[[X1.fmat[n] X2.fmat[n]] for n=1:ndims(X1)] #concatenate fmat
  trten=TensorCell(undef,length(X1.tree.internal_nodes))
  tmp=[dropdims(X1.trten[1],dims=3) zeros(size(X1.trten[1],1),size(X2.trten[1],2)); zeros(size(X1.trten[1],2),size(X2.trten[1],1)) dropdims(X2.trten[1],dims=3)]
  trten[1]=reshape(tmp,size(tmp,1),size(tmp,2),1)
  for i=2:length(X1.tree.internal_nodes)
    trten[i]=blockdiag(X1.trten[i],X2.trten[i])
  end
  htensor(X1.tree,trten,MatrixCell(fmat))
end
+(X1::htensor,X2::htensor)=plus(X1,X2)

#Orthogonalize factor matrices of a tensor. **Documentation in ttensor.jl.
function reorth(X::htensor)
	if X.isorth
		return X
  end
  N=ndims(X)
  In=length(X.tree.internal_nodes)
  trten=TensorCell(undef,In)
  fmat=MatrixCell(undef,N)
	R=MatrixCell(undef,N+In)
  n=1;
	for U in X.fmat
		Q,R[X.tree.leaves[n]]=qr(U)
		fmat[n]=Matrix(Q)
    n+=1
	end
  for n in reverse(X.tree.internal_nodes)
    c=children(X.tree,n)
    Rl=R[c[1]]
    Rr=R[c[2]]
    ind=node2ind(X.tree,n)
    trten[ind]=kron(Rr,Rl)*trten2mat(X.trten[ind])
    if n!=1
      Q,R[n]=qr(trten[ind])
	  trten[ind]=Matrix(Q)
    end
    trten[ind]=trten2ten(trten[ind],size(Rl,2),size(Rr,2))
  end
	htensor(X.tree,trten,fmat)
end

#Orthogonalize factor matrices of a tensor. Rewrite tensor. **Documentation in ttensor.jl.
function reorth!(X::htensor)
  if X.isorth==true
    return X
  end
  N=ndims(X)
  In=length(X.tree.internal_nodes)
  R=MatrixCell(undef,N+In)
  n=1;
  for U in X.fmat
    Q,R[X.tree.leaves[n]]=qr(U)
	X.fmat[n]=Matrix(Q)
    n+=1
  end
  for n in reverse(X.tree.internal_nodes)
    c=children(X.tree,n)
    Rl=R[c[1]]
    Rr=R[c[2]]
    ind=node2ind(X.tree,n)
    X.trten[ind]=kron(Rr,Rl)*trten2mat(X.trten[ind])
    if n!=1
      Q,R[n]=qr(X.trten[ind])
	  X.trten[ind]=Matrix(Q)
    end
    X.trten[ind]=trten2ten(X.trten[ind],size(Rl,2),size(Rr,2))
  end
  X.isorth=true
  X
end

function Base.show(io::IO,X::htensor)
    display(X)
end

#Size of a htensor. **Documentation in Base.
function size(X::htensor)
	tuple([size(X.fmat[n],1) for n=1:ndims(X)]...)
end
#Size of n-th mode of a htensor.
function size(X::htensor,n::Integer)
  size(X)[n]
end

"""
TensorToolbox:

    dropdims(X)

Remove singleton dimensions from htensor.
"""
function dropdims(X::htensor)
  B=deepcopy(X.trten)
  U=deepcopy(X.fmat)
  T=deepcopy(X.tree)
  In=copy(X.tree.internal_nodes)
  L=copy(X.tree.leaves)
  sz=[size(X)...]
  sdims=findall(sz.==1) #singleton dimensions
  while length(findall(sz.!=-1)) >2
    node=next_single_node(T,sdims)
    if node == -1
        break
    end
    sibling_node=sibling(T,node)
    parent_node=parent(T,node)
    ind=node2ind(T,node)
    inds=node2ind(T,sibling_node)
    indp=node2ind(T,parent_node)
    if is_left(T,node)
      tmp=ttm(B[indp],U[ind],1)
      tmp=reshape(tmp,size(tmp,2),size(tmp,3))
      lft=1
    else
      tmp=ttm(B[indp],U[ind],2)
      tmp=reshape(tmp,size(tmp,1),size(tmp,3))
      lft=0
    end
    if is_leaf(T,sibling_node)
      v=setdiff(findall(L.>L[inds]),ind)
      L[v].-=2
      L[inds]=parent_node
      #deleteat!(sdims,findall(sdims.==inds))
      sz[ind]=-1
      if lft==1
        U[ind]=U[ind+1]*tmp
        deleteat!(U,ind+1)
      else
        U[ind-1]=U[ind-1]*tmp
        deleteat!(U,ind)
      end
      deleteat!(B,indp)
      deleteat!(In,indp)
    else
      sibling_subnodes=subnodes(T,sibling_node)
      len=length(nodes_on_lvl(T,lvl(T,node)))
      h=height(subtree(T,parent_node))
      hsib=height(subtree(T,sibling(T,parent_node)))
      for n in sibling_subnodes[2:end]
        if lvl(T,n)>hsib+1 && is_left(T,node)
           break
        end
        if is_leaf(T,n)
           L[node2ind(T,n)]-=len;
        else
           In[node2ind(T,n)]-=len;
        end
      end
      if h<=hsib && is_left(T,node)
        len=length(nodes_on_lvl(subtree(T,parent_node),h-1))
        for level=lvl(T,parent_node)+h-1:height(T)-1
          lvl_nodes=nodes_on_lvl(T,level)
           for n in lvl_nodes[findall((in)(subnodes(T,sibling(T,parent_node))),lvl_nodes)]
             if is_leaf(T,n)
                L[node2ind(T,n)]-=len
              else
                In[node2ind(T,n)]-=len
              end
           end
        end
      end
      B[indp]=ttm(B[inds],copy(tmp'),3)
      deleteat!(B,inds)
      deleteat!(In,inds)
      deleteat!(U,ind)
    end
    deleteat!(L,ind)
    deleteat!(sdims,findall(sdims.==ind))
    sdims[findall(sdims.>ind)].-=1
    T=dimtree(copy(L))
    #display(T)
  end
  H=htensor(T,B,U)
end

"""
    trten2mat(B::Array)
    trten2mat(B::TensorCell)

Transfer tensor to matrix. If transfer tensor is given as a tensor of order 3 and size `(r1,r2,r3)`, reshape it into a matrix of size `(r1r2,r3)`.
"""
function trten2mat(B::Array{<:Number,3})
  (r1,r2,r3)=size(B)
  reshape(B,r1*r2,r3)
end
function trten2mat(B::TensorCell)
  In=length(B)
  @assert(any([ndims(B[i])==3 for i=1:In]),"Transfer tensors should be tensors of order 3.")
  Bmat=MatrixCell(undef,In)
  for i=1:In
    (r1,r2,r3)=size(B[i])
    Bmat[i]=reshape(B[i],r1*r2,r3)
  end
  Bmat
end
trten2mat(B::Array{Vector{T}}) where T<:Number= trten2mat(TensorCell(B))

"""
    trten2ten(B::Matrix,r1::Integer,r2::Integer)
    trten2ten(B::MatrixCell,r1::Vector,r2::Vector)

Transfer tensor to tensor. If transfer tensor is given as a matrix, reshape it into a tensor of order 3 and size r1×r2×r3, where `r3=size(B,2)`.
"""
function trten2ten(B::AbstractMatrix{<:Number},r1::Integer,r2::Integer)
  @assert(r1*r2==size(B,1),"Dimension mismatch.")
  r3=size(B,2)
  reshape(B,r1,r2,r3)
end
trten2ten(B::AbstractVector{<:Number},r1::Integer,r2::Integer)=trten2ten(B[:,:],r1,r2)
function trten2ten(B::MatrixCell,r1::AbstractVector{<:Integer},r2::AbstractVector{<:Integer})
  L=length(B)
  @assert([r1[i]*r2[i]==size(B[i],1) for i=1:L],"Dimension mismatch.")
  Bten=TensorCell(undef,L)
  for i=1:L
    r3=size(B[i],2)
    Bten[i]=reshape(B[i],r1[i],r2[i],r3)
  end
  Bten
end
function trten2ten(B::AbstractMatrix{<:Number})
  reshape(B,size(B,1),size(B,2),1)
end

#htensor times matrix (n-mode product). **Documentation in tensor.jl.
#t='t' transposes matrices
function ttm(X::htensor,M::MatrixCell,modes::AbstractVector{<:Integer},t='n')
  if t=='t'
	 [M[n]=M[n]' for n=1:length(M)]
	end
	@assert(length(modes)<=length(M),"Too few matrices")
	@assert(length(M)<=ndims(X),"Too many matrices")
  if length(modes)<length(M)
    M=M[modes]
  end
  U=copy(X.fmat)
  for n=1:length(modes)
    U[modes[n]]=M[n]*X.fmat[modes[n]]
  end
    htensor(X.tree,X.trten,U)
end
ttm(X::htensor,M::MatrixCell,modes::AbstractRange{<:Integer},t='n')=ttm(X,M,collect(modes),t)
ttm(X::htensor,M::AbstractMatrix{<:Number},n::Integer,t='n')=ttm(X,MatrixCell([M]),[n],t)
ttm(X::htensor,M::MatrixCell,t::Char)=ttm(X,M,1:length(M),t)
ttm(X::htensor,M::MatrixCell)=ttm(X,M,1:length(M))
function ttm(X::htensor,M::MatrixCell,n::Integer,t='n')
	if n>0
		ttm(X,M[n],n,t)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,M,modes,t)
	end
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
ttm(X::htensor,M::Array{Matrix{T},1},modes::AbstractVector{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::htensor,M::Array{Matrix{T},1},modes::AbstractRange{<:Integer},t='n') where {T<:Number}=ttm(X,MatrixCell(M),modes,t)
ttm(X::htensor,M::Array{Matrix{T},1},t::Char) where {T<:Number}=ttm(X,MatrixCell(M),t)
ttm(X::htensor,M::Array{Matrix{T},1}) where {T<:Number}=ttm(X,MatrixCell(M))
ttm(X::htensor,M::Array{Matrix{T},1},n::Integer,t='n') where {T<:Number}=ttm(X,MatrixCell(M),n,t)

#htensor times vector (n-mode product). **Documentation in tensor.jl.
#t='t' transposes matrices
function ttv(X::htensor,V::VectorCell,modes::AbstractVector{<:Integer})
  N=ndims(X)
  remmodes=setdiff(1:N,modes)
  U=deepcopy(X.fmat)
  if length(modes) < length(V)
    V=V[modes]
  end
  for n=1:length(modes)
    U[modes[n]]=V[n]'*X.fmat[modes[n]]
  end
  #reshape(htensor(X.tree,X.trten,U),tuple(sz))
  htensor(X.tree,X.trten,U)
end
ttv(X::htensor,v::AbstractVector{<:Number},n::Integer)=ttv(X,VectorCell([v]),[n])
ttv(X::htensor,V::VectorCell,modes::AbstractRange{<:Integer})=ttv(X,V,collect(modes))
ttv(X::htensor,V::VectorCell)=ttv(X,V,1:length(V))
function ttv(X::htensor,V::VectorCell,n::Integer)
	if n>0
		ttm(X,V[n],n)
	else
		modes=setdiff(1:ndims(X),-n)
		ttm(X,A,modes)
	end
end
#If array of vectors isn't defined as VectorCell, but as V=[v1,v2,...,vn]:
ttv(X::htensor,V::Array{Vector{T}},modes::AbstractVector{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::htensor,V::Array{Vector{T}},modes::AbstractRange{<:Integer}) where {T<:Number}=ttv(X,VectorCell(V),modes)
ttv(X::htensor,V::Array{Vector{T}}) where {T<:Number}=ttv(X,VectorCell(V),1:length(V))
ttv(X::htensor,V::Array{Vector{T}},n::Integer) where {T<:Number}=ttv(X,VectorCell(V),n)

#**Documentation in ttensor.jl.
uminus(X::htensor)=mtimes(-1,X)
-(X::htensor)=uminus(X)

#Frobenius norm of a htensor. **Documentation in Base.
function norm(X::htensor;orth=true)
  if orth
    Y=reorth(X)
    norm(Y.trten[1])
  else
    sqrt(innerprod(X,X))
  end
end


text:
# example 5: write all at once with writefasta

module FastaIOExample5
using FastaIO
using GZip

# notes about the data:
#   1. the sequence type could be
#      different, e.g. Vector{Uint8},
#      Vector{AminoAcid} etc.
#   2. the container could be any iterable
#      which returns a tuple (desc,name) when
#      iterated, e.g. a Dict:
#
#      fastadata = Dict("ITEM1"=>"seq1", ...)
#
fastadata = [
    ("ITEM1","""
     -----------STVELTKEN-F--D-Q
     E--F-V--LI--------D-----F--
     P----C-----R----Q---F
     """),
    ("ITEM2","""
     ----------MVKEITDAT-F--E-Q-
     --L-V--LT--------D-----F--W
     -------M-------D-----
     """),
    ("ITEM3","""
     ----------NLESVEQFD--------
     -------G--K-S--VF--------M-
     -D------L------A-----
     """)]

function write_fasta_data(filename::String, data)
    println("Writing $(length(data)) entries to file $filename")

    # note: writefasta can take a writable IO descriptor as
    #       argument instead of a filename; also, it has
    #       an optional mode argument which defaults to "w"
    #       (use "a" for appending to an exisitng file)
    writefasta(filename, data)

    println("This is the content of the file:")
    gzopen(filename) do f
        println(read(f, String))
    end
end

# note: gzip compression is used because the filename ends with .gz
write_fasta_data(joinpath(@__DIR__, "example_out.fasta.gz"), fastadata)

# remove the file
rm(joinpath(@__DIR__, "example_out.fasta.gz"))

end # module


text:
register(DataDep(
    "DS19 presentation materials",
    """
	Dataset: DS19 presentation materials
	Website: https://figshare.com/articles/DS19_presentation_materials/8146106
	Author: Sebastian Micluța-Câmpeanu
	Date of Publication: 2019-05-22T15:36:08Z
	License: MIT (https://opensource.org/licenses/MIT)

	Here are the materials for the SIAM DS19 presentation entitled "Large-Scale Numerical Investigations into the Dynamics of Nonlinear Classical Systems". The dataset used for all the computations is included as the binary file graph.jls and it requires julia 1.1.0 and the two julia packages linked bellow to be correctly read. See the DS19Presentation GitHub repository linked bellow for more details on how to reproduce the presented figures.

	Please cite this dataset: Micluța-Câmpeanu, Sebastian (2019): DS19 presentation materials. figshare. Conference contribution.
	""",
	Any["https://ndownloader.figshare.com/files/15181136", "https://ndownloader.figshare.com/files/15221714"],
	[(md5, "cefed4ae0c85586a8a31305d63da31d4"), (md5, "193df57d9f1c46aed3d8521e0b7d0028")]
))


text:
"""
    MomentsAttribute(N)
    MomentsAttribute()

A constraint attribute for the vector of moments corresponding to the
constraint that a polynomial is zero in the full space in result `N`. If the
polynomial is constrained to be zero in an algebraic set, it is the moments
for the constraint once it is rewritten into an constraint on the full space.
If `N` is omitted, it is 1 by default.

## Examples

Consider the following program:
```julia
@variable(model, α)
@variable(model, β ≤ 1)

using DynamicPolynomials
@polyvar x y
cref = @constraint(model, α * x - β * y == 0, domain = @set x == y)
```
The constraint is equivalent to
```julia
@constraint(model, (α - β) * y == 0)
```
for which the dual is the 1-element vector with the moment of `y` of value `-1`.
This is the result of `moments(cref)`. However, the dual of `cref` obtained by
`dual(cref)` is the 2-elements vector with both the moments of `x` and `y`
of value `-1`.
"""
struct MomentsAttribute <: MOI.AbstractConstraintAttribute
    N::Int
end
MomentsAttribute() = MomentsAttribute(1)

function MOI.is_set_by_optimize(::MomentsAttribute)
    return true
end


text:
module QuantumPropagators

# low-level routines

include("./arnoldi.jl")

include("./specrad.jl")
export specrange

include("./cheby.jl")
export cheby_coeffs, cheby_coeffs!, ChebyWrk, cheby!, cheby

include("./newton.jl")
export NewtonWrk, newton!

include("./expprop.jl")
export ExpPropWrk, expprop!

include("./storage.jl")
export init_storage, map_observables, map_observable, write_to_storage!
export get_from_storage!

# high-level interface
include("./propagate.jl")
export initpropwrk, init_storage, propstep!, propstep, propagate

end


text:
abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end

@noinline function _slice_rmul!(u_temp::AbstractArray{T,N}, A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,N}
    for J in post
        for I in pre
            u_temp[I, :, J] = A*u[I, :, J]
        end
    end
    u_temp
end

"""
slice_rmul lets you multiply each vector like strip of an array `u` with a linear operator `A`, sliced along dimension `dim`
"""
function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
    @assert N != 1
    u_temp = zero(u)

    _slice_rmul!(u_temp, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))

    return u_temp
end

@noinline function _slice_rmul!(lower::AbstractArray, upper::AbstractArray, A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int, pre, post) where {T,B,N,M}
    for J in post
        for I in pre

            tmp = A[I,J]*u[I, :, J]
            lower[I,J], upper[I,J] = tmp.l, tmp.r
        end
    end
    return (lower, upper)
end

function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) where {T, B, N,M}
    @assert N != 1
    @assert M == N-1
    lower = zeros(T,perpsize(u,dim))
    upper = zeros(T,perpsize(u,dim))

    _slice_rmul!(lower, upper, A, u, dim, CartesianIndices(axes(u)[1:dim-1]), CartesianIndices(axes(u)[(dim+1):end]))

    return (lower, upper)
end

struct MultiDimDirectionalBC{T<:Number, B<:AtomicBC{T}, D, N, M} <: MultiDimensionalBC{T, N}
    BCs::Array{B,M} #dimension M=N-1 array of BCs to extend dimension D
end

struct ComposedMultiDimBC{T, B<:AtomicBC{T}, N,M} <: MultiDimensionalBC{T, N}
    BCs::Vector{Array{B, M}} # Why isn't this a vector of MultiDimBCs?
end

"""
A multiple dimensional BC, supporting arbitrary BCs at each boundary point.
To construct an arbitrary BC, pass an Array of BCs with dimension `N-1`, if `N` is the dimensionality of your domain `u`
with a size of `size(u)[setdiff(1:N, dim)]`, where dim is the dimension orthogonal to the boundary that you want to extend.

It is also possible to call
    Q_dim = MultiDimBC(YourBC, size(u), dim)
to use YourBC for the whole boundary orthogonal to that dimension.

Further, it is possible to call
Qx, Qy, Qz... = MultiDimBC(YourBC, size(u))
to use YourBC for the whole boundary for all dimensions. Valid for any number of dimensions greater than 1.
However this is only valid for Robin/General type BCs (including neummann/dirichlet) when the grid steps are equal in each dimension - including uniform grid case.

In the case where you want to extend the same Robin/GeneralBC to the whole boundary with a non unifrom grid, please use
    Qx, Qy, Qz... = RobinBC(l, r, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
or
    Qx, Qy, Qz... = GeneralBC(αl, αr, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))

There are also constructors for NeumannBC, DirichletBC and Dirichlet0BC. Simply replace `dx` in the call with the tuple dxyz... as above, and append `size(u)`` to the argument signature.
The order is a required argument in this case.

where dx, dy, and dz are vectors of grid steps.

For Neumann0BC, please use
    Qx, Qy, Qz... = Neumann0BC(T::Type, (dx::Vector, dy::Vector, dz::Vector ...), approximation_order, size(u))
where T is the element type of the domain to be extended
"""
struct MultiDimBC{N} end

MultiDimBC{dim}(BC::Array{B,N}) where {N, B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC[1]), B, dim, N+1, N}(BC)
#s should be size of the domain
MultiDimBC{dim}(BC::B, s) where  {B<:AtomicBC, dim} = MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)]))

#Extra constructor to make a set of BC operators that extend an atomic BC Operator to the whole domain
#Only valid in the uniform grid case!
MultiDimBC(BC::B, s) where {B<:AtomicBC} = Tuple([MultiDimDirectionalBC{gettype(BC), B, dim, length(s), length(s)-1}(fill(BC, s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])

# Additional constructors for cases when the BC is the same for all boundarties
PeriodicBC{dim}(T,s) where dim = MultiDimBC{dim}(PeriodicBC(T), s)
PeriodicBC(T,s) = MultiDimBC(PeriodicBC(T), s)

NeumannBC{dim}(α::NTuple{2,T}, dx, order, s) where {T,dim} = RobinBC{dim}((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order, s)
NeumannBC(α::NTuple{2,T}, dxyz, order, s) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dxyz, order, s)

DirichletBC{dim}(αl::T, αr::T, s) where {T,dim} = RobinBC{dim}((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2.0, s)
DirichletBC(αl::T, αr::T, s) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), fill(one(T), length(s)), 2.0, s)

Dirichlet0BC{dim}(T::Type, s) where {dim} = DirichletBC{dim}(zero(T), zero(T), s)
Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)

Neumann0BC{dim}(T::Type, dx, order, s) where {dim} = NeumannBC{dim}((zero(T), zero(T)), dx, order, s)
Neumann0BC(T::Type, dxyz, order, s) = NeumannBC((zero(T), zero(T)), dxyz, order, s)

RobinBC{dim}(l::NTuple{3,T}, r::NTuple{3,T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(RobinBC(l, r, dx, order), s)
RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), perpindex(s,dim))) for dim in 1:length(s)])

GeneralBC{dim}(αl::AbstractVector{T}, αr::AbstractVector{T}, dx, order, s) where {T,dim} = MultiDimBC{dim}(GeneralBC(αl, αr, dx, order), s)
GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order),perpindex(s,dim))) for dim in 1:length(s)])

"""
Q = compose(BCs...)

-------------------------------------------------------------------------------------

Example:
Q = compose(Qx, Qy, Qz) # 3D domain
Q = compose(Qx, Qy) # 2D Domain

Creates a ComposedMultiDimBC operator, Q, that extends every boundary when applied to a `u` with compatible size and number of dimensions.

Qx Qy and Qz can be passed in any order, as long as there is exactly one BC operator that extends each dimension.
"""
function compose(BCs::MultiDimDirectionalBC...)
    T = gettype(BCs[1])
    N = ndims(BCs[1])
    Ds = getaxis.(BCs)
    (length(BCs) == N) || throw(ArgumentError("There must be enough BCs to cover every dimension - check that the number of MultiDimBCs == N"))
    for D in Ds
        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple boundary conditions that extend along $D - make sure every dimension has a unique extension"))
    end
    BCs = BCs[sortperm([Ds...])]

    ComposedMultiDimBC{T, Union{eltype.(BC.BCs for BC in BCs)...}, N,N-1}([condition.BCs for condition in BCs])
end

Base.:+(BCs::MultiDimDirectionalBC...) = compose(BCs...)

"""
Qx, Qy,... = decompose(Q::ComposedMultiDimBC{T,N,M})

-------------------------------------------------------------------------------------

Decomposes a ComposedMultiDimBC in to components that extend along each dimension individually
"""
decompose(Q::ComposedMultiDimBC) = Tuple([MultiDimBC(Q.BCs[i], i) for i in 1:ndims(Q)])

getaxis(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = D
getboundarytype(Q::MultiDimDirectionalBC{T, B, D, N, K}) where {T, B, D, N, K} = B

Base.ndims(Q::MultiDimensionalBC{T,N}) where {T,N} = N

Base.:*(BC::AtomicBC, u::AbstractArray) = MultiDimBC{1}(BC, size(u)) * u

function Base.:*(Q::MultiDimDirectionalBC{T, B, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
    @assert perpsize(u, D) == size(Q.BCs) "Size of the BCs array in the MultiDimBC is incorrect, needs to be $(perpsize(u,D)) to extend dimension $D, got $(size(Q.BCs))"
    lower, upper = slice_rmul(Q.BCs, u, D)
    return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
end

function Base.:*(Q::MultiDimDirectionalBC{T, PeriodicBC{T}, D, N, K}, u::AbstractArray{T, N}) where {T, B, D, N, K}
    lower = selectdim(u, D, 1)
    upper = selectdim(u, D, size(u,D))
    return BoundaryPaddedArray{T, D, N, K, typeof(u), Union{typeof(lower), typeof(upper)}}(lower, upper, u)
end

function Base.:*(Q::ComposedMultiDimBC{T, B, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
    for dim in 1:N
        @assert perpsize(u, dim) == size(Q.BCs[dim]) "Size of the BCs array for dimension $dim in the MultiDimBC is incorrect, needs to be $(perpsize(u,dim)), got $(size(Q.BCs[dim]))"
    end
    out = slice_rmul.(Q.BCs, fill(u, N), 1:N)
    lower = [A[1] for A in out]
    upper = [A[2] for A in out]
    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
end

function Base.:*(Q::ComposedMultiDimBC{T, PeriodicBC{T}, N, K}, u::AbstractArray{T, N}) where {T, B, N, K}
    lower = [selectdim(u, d, 1) for d in 1:N]
    upper = [selectdim(u, d, size(u, d)) for d in 1:N]
    return ComposedBoundaryPaddedArray{T, N, K, typeof(u), Union{typeof.(lower)..., typeof.(upper)...}}(lower, upper, u)
end


text:
#=

  Scheduling speakers in Julia ConstraintSolver.jl

  From Rina Dechter, Constraint Processing, page 72
  Scheduling of 6 speakers in 6 slots.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function scheduling_speakers(available_slots,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n = length(available_slots)
    @variable(model, 1 <= xs[1:n] <= n, Int)

    @constraint(model, xs in CS.AllDifferent())

    for (x,a) in zip(xs,available_slots)
        len = length(a)
        b = @variable(model,[1:len],Bin)
        for i in 1:len
            @constraint(model,b[i] := {a[i] == x})
        end
        @constraint(model,sum(b) >= 1)

        # This version is too slow for showin all (2) solutions:
        # ix = @variable(model,[1:len],CS.Integers(1:len))
        # my_element(model,ix[1],a,x)
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                xs_val = convert.(Integer,JuMP.value.(xs; result=sol))
                println("xs:$xs_val")
            end
        end
    else
        println("status:$status")
    end

    return status
end

available_slots = [
                       # Reasoning:
         [3,4,5,6],    # 2) the only one with 6 after speaker F -> 1
         [3,4],        # 5) 3 or 4
         [2,3,4,5],    # 3) only with 5 after F -> 1 and A -> 6
         [2,3,4],      # 4) only with 2 after C -> 5 and F -> 1
         [3,4],        # 5) 3 or 4
         [1,2,3,4,5,6] # 1) the only with 1

]
@time scheduling_speakers(available_slots,true,true)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text: generaterange(min, max, step) = Array(min:step:max)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
function constant_mortality_mean_age_corrected(mx, nx, boundary = 1e-6)
    result = copy(mx)
    upper = mx .>= boundary
    large_mx  = mx[upper]
    expx = exp.(-large_mx .* nx[upper])
    result[upper] = (1 ./ large_mx) .- (nx[upper] .* expx) ./ (1 .- expx)

    lower = .!upper
    small_mx = mx[lower]
    small_nx = nx[lower]
    result[lower] = small_nx .* (0.5 .+ small_nx .* (
        -small_mx / 12 + small_mx.^3 .* small_nx.^2 / 720))
    result
end

function constant_mortality_mean_age(mx, nx)
    expx = exp.(-mx .* nx)
    (1 ./ mx) .- (nx .* expx) ./ (1 .- expx)
end

# The Langevin function is coth(x) - 1/x
function langevin(x)
    x / 3 / (
        1 + x^2 / 15 / (
            1 + x^2 / 35 / (
                1 + x^2 / 63 / (
                    1 + x^2 / 99
                )
            )
        )
    )
end

function langevin_mean_age(mx, nx)
    (nx / 2) .* (1 .- langevin.(mx .* nx / 2))
end

function test_naive_mean_age()
    mx = vcat([0.9 * 10.0^i for i in 0:-1:-17], [0.0])
    nx = ones(size(mx))
    ax = constant_mortality_mean_age(mx, nx)
    println(ax)
    for check in 2:length(ax)
        @assert ax[check] > ax[check - 1]
    end
    @assert all(ax .< .5)
end

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function ARKodeCreate()
    ccall((:ARKodeCreate, libsundials_arkode), ARKODEMemPtr, ())
end

function ARKodeSetDefaults(arkode_mem)
    ccall((:ARKodeSetDefaults, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetOptimalParams(arkode_mem)
    ccall((:ARKodeSetOptimalParams, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetErrHandlerFn(arkode_mem, ehfun::ARKErrHandlerFn, eh_data)
    ccall((:ARKodeSetErrHandlerFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKErrHandlerFn, Ptr{Cvoid}), arkode_mem, ehfun, eh_data)
end

function ARKodeSetErrFile(arkode_mem, errfp)
    ccall((:ARKodeSetErrFile, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{FILE}), arkode_mem, errfp)
end

function ARKodeSetUserData(arkode_mem, user_data)
    ccall((:ARKodeSetUserData, libsundials_arkode), Cint, (ARKODEMemPtr, Any), arkode_mem, user_data)
end

function ARKodeSetDiagnostics(arkode_mem, diagfp)
    ccall((:ARKodeSetDiagnostics, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{FILE}), arkode_mem, diagfp)
end

function __ARKodeSetOrder(arkode_mem, maxord::Cint)
    ccall((:ARKodeSetOrder, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, maxord)
end

function ARKodeSetOrder(arkode_mem, maxord)
    __ARKodeSetOrder(arkode_mem, convert(Cint, maxord))
end

function __ARKodeSetDenseOrder(arkode_mem, dord::Cint)
    ccall((:ARKodeSetDenseOrder, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, dord)
end

function ARKodeSetDenseOrder(arkode_mem, dord)
    __ARKodeSetDenseOrder(arkode_mem, convert(Cint, dord))
end

function __ARKodeSetLinear(arkode_mem, timedepend::Cint)
    ccall((:ARKodeSetLinear, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, timedepend)
end

function ARKodeSetLinear(arkode_mem, timedepend)
    __ARKodeSetLinear(arkode_mem, convert(Cint, timedepend))
end

function ARKodeSetNonlinear(arkode_mem)
    ccall((:ARKodeSetNonlinear, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function __ARKodeSetFixedPoint(arkode_mem, fp_m::Clong)
    ccall((:ARKodeSetFixedPoint, libsundials_arkode), Cint, (ARKODEMemPtr, Clong), arkode_mem, fp_m)
end

function ARKodeSetFixedPoint(arkode_mem, fp_m)
    __ARKodeSetFixedPoint(arkode_mem, convert(Clong, fp_m))
end

function ARKodeSetNewton(arkode_mem)
    ccall((:ARKodeSetNewton, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetExplicit(arkode_mem)
    ccall((:ARKodeSetExplicit, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetImplicit(arkode_mem)
    ccall((:ARKodeSetImplicit, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetImEx(arkode_mem)
    ccall((:ARKodeSetImEx, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function __ARKodeSetERKTable(arkode_mem, s::Cint, q::Cint, p::Cint, c, A, b, bembed)
    ccall((:ARKodeSetERKTable, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, Cint, Cint, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), arkode_mem, s, q, p, c, A, b, bembed)
end

function ARKodeSetERKTable(arkode_mem, s, q, p, c, A, b, bembed)
    __ARKodeSetERKTable(arkode_mem, convert(Cint, s), convert(Cint, q), convert(Cint, p), c, A, b, bembed)
end

function __ARKodeSetIRKTable(arkode_mem, s::Cint, q::Cint, p::Cint, c, A, b, bembed)
    ccall((:ARKodeSetIRKTable, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, Cint, Cint, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), arkode_mem, s, q, p, c, A, b, bembed)
end

function ARKodeSetIRKTable(arkode_mem, s, q, p, c, A, b, bembed)
    __ARKodeSetIRKTable(arkode_mem, convert(Cint, s), convert(Cint, q), convert(Cint, p), c, A, b, bembed)
end

function __ARKodeSetARKTables(arkode_mem, s::Cint, q::Cint, p::Cint, ci, ce, Ai, Ae, bi, be, b2i, b2e)
    ccall((:ARKodeSetARKTables, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, Cint, Cint, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), arkode_mem, s, q, p, ci, ce, Ai, Ae, bi, be, b2i, b2e)
end

function ARKodeSetARKTables(arkode_mem, s, q, p, ci, ce, Ai, Ae, bi, be, b2i, b2e)
    __ARKodeSetARKTables(arkode_mem, convert(Cint, s), convert(Cint, q), convert(Cint, p), ci, ce, Ai, Ae, bi, be, b2i, b2e)
end

function __ARKodeSetERKTableNum(arkode_mem, itable::Cint)
    ccall((:ARKodeSetERKTableNum, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, itable)
end

function ARKodeSetERKTableNum(arkode_mem, itable)
    __ARKodeSetERKTableNum(arkode_mem, convert(Cint, itable))
end

function __ARKodeSetIRKTableNum(arkode_mem, itable::Cint)
    ccall((:ARKodeSetIRKTableNum, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, itable)
end

function ARKodeSetIRKTableNum(arkode_mem, itable)
    __ARKodeSetIRKTableNum(arkode_mem, convert(Cint, itable))
end

function __ARKodeSetARKTableNum(arkode_mem, itable::Cint, etable::Cint)
    ccall((:ARKodeSetARKTableNum, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, Cint), arkode_mem, itable, etable)
end

function ARKodeSetARKTableNum(arkode_mem, itable, etable)
    __ARKodeSetARKTableNum(arkode_mem, convert(Cint, itable), convert(Cint, etable))
end

function __ARKodeSetMaxNumSteps(arkode_mem, mxsteps::Clong)
    ccall((:ARKodeSetMaxNumSteps, libsundials_arkode), Cint, (ARKODEMemPtr, Clong), arkode_mem, mxsteps)
end

function ARKodeSetMaxNumSteps(arkode_mem, mxsteps)
    __ARKodeSetMaxNumSteps(arkode_mem, convert(Clong, mxsteps))
end

function __ARKodeSetMaxHnilWarns(arkode_mem, mxhnil::Cint)
    ccall((:ARKodeSetMaxHnilWarns, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, mxhnil)
end

function ARKodeSetMaxHnilWarns(arkode_mem, mxhnil)
    __ARKodeSetMaxHnilWarns(arkode_mem, convert(Cint, mxhnil))
end

function ARKodeSetInitStep(arkode_mem, hin::realtype)
    ccall((:ARKodeSetInitStep, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, hin)
end

function ARKodeSetMinStep(arkode_mem, hmin::realtype)
    ccall((:ARKodeSetMinStep, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, hmin)
end

function ARKodeSetMaxStep(arkode_mem, hmax::realtype)
    ccall((:ARKodeSetMaxStep, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, hmax)
end

function ARKodeSetStopTime(arkode_mem, tstop::realtype)
    ccall((:ARKodeSetStopTime, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, tstop)
end

function ARKodeSetFixedStep(arkode_mem, hfixed::realtype)
    ccall((:ARKodeSetFixedStep, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, hfixed)
end

function ARKodeSetCFLFraction(arkode_mem, cfl_frac::realtype)
    ccall((:ARKodeSetCFLFraction, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, cfl_frac)
end

function ARKodeSetSafetyFactor(arkode_mem, safety::realtype)
    ccall((:ARKodeSetSafetyFactor, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, safety)
end

function ARKodeSetErrorBias(arkode_mem, bias::realtype)
    ccall((:ARKodeSetErrorBias, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, bias)
end

function ARKodeSetMaxGrowth(arkode_mem, mx_growth::realtype)
    ccall((:ARKodeSetMaxGrowth, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, mx_growth)
end

function ARKodeSetFixedStepBounds(arkode_mem, lb::realtype, ub::realtype)
    ccall((:ARKodeSetFixedStepBounds, libsundials_arkode), Cint, (ARKODEMemPtr, realtype, realtype), arkode_mem, lb, ub)
end

function __ARKodeSetAdaptivityMethod(arkode_mem, imethod::Cint, idefault::Cint, pq::Cint, adapt_params)
    ccall((:ARKodeSetAdaptivityMethod, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, Cint, Cint, Ptr{realtype}), arkode_mem, imethod, idefault, pq, adapt_params)
end

function ARKodeSetAdaptivityMethod(arkode_mem, imethod, idefault, pq, adapt_params)
    __ARKodeSetAdaptivityMethod(arkode_mem, convert(Cint, imethod), convert(Cint, idefault), convert(Cint, pq), adapt_params)
end

function ARKodeSetAdaptivityFn(arkode_mem, hfun::ARKAdaptFn, h_data)
    ccall((:ARKodeSetAdaptivityFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKAdaptFn, Ptr{Cvoid}), arkode_mem, hfun, h_data)
end

function ARKodeSetMaxFirstGrowth(arkode_mem, etamx1::realtype)
    ccall((:ARKodeSetMaxFirstGrowth, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, etamx1)
end

function ARKodeSetMaxEFailGrowth(arkode_mem, etamxf::realtype)
    ccall((:ARKodeSetMaxEFailGrowth, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, etamxf)
end

function __ARKodeSetSmallNumEFails(arkode_mem, small_nef::Cint)
    ccall((:ARKodeSetSmallNumEFails, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, small_nef)
end

function ARKodeSetSmallNumEFails(arkode_mem, small_nef)
    __ARKodeSetSmallNumEFails(arkode_mem, convert(Cint, small_nef))
end

function ARKodeSetMaxCFailGrowth(arkode_mem, etacf::realtype)
    ccall((:ARKodeSetMaxCFailGrowth, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, etacf)
end

function ARKodeSetNonlinCRDown(arkode_mem, crdown::realtype)
    ccall((:ARKodeSetNonlinCRDown, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, crdown)
end

function ARKodeSetNonlinRDiv(arkode_mem, rdiv::realtype)
    ccall((:ARKodeSetNonlinRDiv, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, rdiv)
end

function ARKodeSetDeltaGammaMax(arkode_mem, dgmax::realtype)
    ccall((:ARKodeSetDeltaGammaMax, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, dgmax)
end

function __ARKodeSetMaxStepsBetweenLSet(arkode_mem, msbp::Cint)
    ccall((:ARKodeSetMaxStepsBetweenLSet, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, msbp)
end

function ARKodeSetMaxStepsBetweenLSet(arkode_mem, msbp)
    __ARKodeSetMaxStepsBetweenLSet(arkode_mem, convert(Cint, msbp))
end

function __ARKodeSetPredictorMethod(arkode_mem, method::Cint)
    ccall((:ARKodeSetPredictorMethod, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, method)
end

function ARKodeSetPredictorMethod(arkode_mem, method)
    __ARKodeSetPredictorMethod(arkode_mem, convert(Cint, method))
end

function ARKodeSetStabilityFn(arkode_mem, EStab::ARKExpStabFn, estab_data)
    ccall((:ARKodeSetStabilityFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKExpStabFn, Ptr{Cvoid}), arkode_mem, EStab, estab_data)
end

function __ARKodeSetMaxErrTestFails(arkode_mem, maxnef::Cint)
    ccall((:ARKodeSetMaxErrTestFails, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, maxnef)
end

function ARKodeSetMaxErrTestFails(arkode_mem, maxnef)
    __ARKodeSetMaxErrTestFails(arkode_mem, convert(Cint, maxnef))
end

function __ARKodeSetMaxNonlinIters(arkode_mem, maxcor::Cint)
    ccall((:ARKodeSetMaxNonlinIters, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, maxcor)
end

function ARKodeSetMaxNonlinIters(arkode_mem, maxcor)
    __ARKodeSetMaxNonlinIters(arkode_mem, convert(Cint, maxcor))
end

function __ARKodeSetMaxConvFails(arkode_mem, maxncf::Cint)
    ccall((:ARKodeSetMaxConvFails, libsundials_arkode), Cint, (ARKODEMemPtr, Cint), arkode_mem, maxncf)
end

function ARKodeSetMaxConvFails(arkode_mem, maxncf)
    __ARKodeSetMaxConvFails(arkode_mem, convert(Cint, maxncf))
end

function ARKodeSetNonlinConvCoef(arkode_mem, nlscoef::realtype)
    ccall((:ARKodeSetNonlinConvCoef, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, nlscoef)
end

function ARKodeSetRootDirection(arkode_mem, rootdir)
    ccall((:ARKodeSetRootDirection, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Cint}), arkode_mem, rootdir)
end

function ARKodeSetNoInactiveRootWarn(arkode_mem)
    ccall((:ARKodeSetNoInactiveRootWarn, libsundials_arkode), Cint, (ARKODEMemPtr,), arkode_mem)
end

function ARKodeSetPostprocessStepFn(arkode_mem, ProcessStep::ARKPostProcessStepFn)
    ccall((:ARKodeSetPostprocessStepFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKPostProcessStepFn), arkode_mem, ProcessStep)
end

function __ARKodeInit(arkode_mem, fe::ARKRhsFn, fi::ARKRhsFn, t0::realtype, y0::N_Vector)
    ccall((:ARKodeInit, libsundials_arkode), Cint, (ARKODEMemPtr, ARKRhsFn, ARKRhsFn, realtype, N_Vector), arkode_mem, fe, fi, t0, y0)
end

function ARKodeInit(arkode_mem, fe, fi, t0, y0)
    __y0 = convert(NVector, y0)
    __ARKodeInit(arkode_mem, ARKRhsFn_wrapper(fe), ARKRhsFn_wrapper(fi), t0, convert(N_Vector, __y0))
end

function __ARKodeReInit(arkode_mem, fe::ARKRhsFn, fi::ARKRhsFn, t0::realtype, y0::N_Vector)
    ccall((:ARKodeReInit, libsundials_arkode), Cint, (ARKODEMemPtr, ARKRhsFn, ARKRhsFn, realtype, N_Vector), arkode_mem, fe, fi, t0, y0)
end

function ARKodeReInit(arkode_mem, fe, fi, t0, y0)
    __y0 = convert(NVector, y0)
    __ARKodeReInit(arkode_mem, ARKRhsFn_wrapper(fe), ARKRhsFn_wrapper(fi), t0, convert(N_Vector, __y0))
end

function __ARKodeResize(arkode_mem, ynew::N_Vector, hscale::realtype, t0::realtype, resize::ARKVecResizeFn, resize_data)
    ccall((:ARKodeResize, libsundials_arkode), Cint, (ARKODEMemPtr, N_Vector, realtype, realtype, ARKVecResizeFn, Ptr{Cvoid}), arkode_mem, ynew, hscale, t0, resize, resize_data)
end

function ARKodeResize(arkode_mem, ynew, hscale, t0, resize, resize_data)
    __ynew = convert(NVector, ynew)
    __ARKodeResize(arkode_mem, convert(N_Vector, __ynew), hscale, t0, resize, resize_data)
end

function ARKodeSStolerances(arkode_mem, reltol::realtype, abstol::realtype)
    ccall((:ARKodeSStolerances, libsundials_arkode), Cint, (ARKODEMemPtr, realtype, realtype), arkode_mem, reltol, abstol)
end

function __ARKodeSVtolerances(arkode_mem, reltol::realtype, abstol::N_Vector)
    ccall((:ARKodeSVtolerances, libsundials_arkode), Cint, (ARKODEMemPtr, realtype, N_Vector), arkode_mem, reltol, abstol)
end

function ARKodeSVtolerances(arkode_mem, reltol, abstol)
    __abstol = convert(NVector, abstol)
    __ARKodeSVtolerances(arkode_mem, reltol, convert(N_Vector, __abstol))
end

function ARKodeWFtolerances(arkode_mem, efun::ARKEwtFn)
    ccall((:ARKodeWFtolerances, libsundials_arkode), Cint, (ARKODEMemPtr, ARKEwtFn), arkode_mem, efun)
end

function ARKodeResStolerance(arkode_mem, rabstol::realtype)
    ccall((:ARKodeResStolerance, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, rabstol)
end

function __ARKodeResVtolerance(arkode_mem, rabstol::N_Vector)
    ccall((:ARKodeResVtolerance, libsundials_arkode), Cint, (ARKODEMemPtr, N_Vector), arkode_mem, rabstol)
end

function ARKodeResVtolerance(arkode_mem, rabstol)
    __rabstol = convert(NVector, rabstol)
    __ARKodeResVtolerance(arkode_mem, convert(N_Vector, __rabstol))
end

function ARKodeResFtolerance(arkode_mem, rfun::ARKRwtFn)
    ccall((:ARKodeResFtolerance, libsundials_arkode), Cint, (ARKODEMemPtr, ARKRwtFn), arkode_mem, rfun)
end

function __ARKodeRootInit(arkode_mem, nrtfn::Cint, g::ARKRootFn)
    ccall((:ARKodeRootInit, libsundials_arkode), Cint, (ARKODEMemPtr, Cint, ARKRootFn), arkode_mem, nrtfn, g)
end

function ARKodeRootInit(arkode_mem, nrtfn, g)
    __ARKodeRootInit(arkode_mem, convert(Cint, nrtfn), g)
end

function __ARKode(arkode_mem, tout::realtype, yout::N_Vector, tret, itask::Cint)
    ccall((:ARKode, libsundials_arkode), Cint, (ARKODEMemPtr, realtype, N_Vector, Ptr{realtype}, Cint), arkode_mem, tout, yout, tret, itask)
end

function ARKode(arkode_mem, tout, yout, tret, itask)
    __yout = convert(NVector, yout)
    __ARKode(arkode_mem, tout, convert(N_Vector, __yout), tret, convert(Cint, itask))
end

function __ARKodeGetDky(arkode_mem, t::realtype, k::Cint, dky::N_Vector)
    ccall((:ARKodeGetDky, libsundials_arkode), Cint, (ARKODEMemPtr, realtype, Cint, N_Vector), arkode_mem, t, k, dky)
end

function ARKodeGetDky(arkode_mem, t, k, dky)
    __dky = convert(NVector, dky)
    __ARKodeGetDky(arkode_mem, t, convert(Cint, k), convert(N_Vector, __dky))
end

function ARKodeGetWorkSpace(arkode_mem, lenrw, leniw)
    ccall((:ARKodeGetWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrw, leniw)
end

function ARKodeGetNumSteps(arkode_mem, nsteps)
    ccall((:ARKodeGetNumSteps, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nsteps)
end

function ARKodeGetNumExpSteps(arkode_mem, expsteps)
    ccall((:ARKodeGetNumExpSteps, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, expsteps)
end

function ARKodeGetNumAccSteps(arkode_mem, accsteps)
    ccall((:ARKodeGetNumAccSteps, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, accsteps)
end

function ARKodeGetNumStepAttempts(arkode_mem, step_attempts)
    ccall((:ARKodeGetNumStepAttempts, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, step_attempts)
end

function ARKodeGetNumRhsEvals(arkode_mem, nfe_evals, nfi_evals)
    ccall((:ARKodeGetNumRhsEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, nfe_evals, nfi_evals)
end

function ARKodeGetNumLinSolvSetups(arkode_mem, nlinsetups)
    ccall((:ARKodeGetNumLinSolvSetups, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nlinsetups)
end

function ARKodeGetNumMassSolves(arkode_mem, nMassSolves)
    ccall((:ARKodeGetNumMassSolves, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nMassSolves)
end

function ARKodeGetNumMassMultiplies(arkode_mem, nMassMult)
    ccall((:ARKodeGetNumMassMultiplies, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nMassMult)
end

function ARKodeGetNumErrTestFails(arkode_mem, netfails)
    ccall((:ARKodeGetNumErrTestFails, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, netfails)
end

function ARKodeGetActualInitStep(arkode_mem, hinused)
    ccall((:ARKodeGetActualInitStep, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{realtype}), arkode_mem, hinused)
end

function ARKodeGetLastStep(arkode_mem, hlast)
    ccall((:ARKodeGetLastStep, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{realtype}), arkode_mem, hlast)
end

function ARKodeGetCurrentStep(arkode_mem, hcur)
    ccall((:ARKodeGetCurrentStep, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{realtype}), arkode_mem, hcur)
end

function ARKodeGetCurrentTime(arkode_mem, tcur)
    ccall((:ARKodeGetCurrentTime, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{realtype}), arkode_mem, tcur)
end

function ARKodeGetCurrentButcherTables(arkode_mem, s, q, p, Ai, Ae, ci, ce, bi, be, b2i, b2e)
    ccall((:ARKodeGetCurrentButcherTables, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), arkode_mem, s, q, p, Ai, Ae, ci, ce, bi, be, b2i, b2e)
end

function ARKodeGetTolScaleFactor(arkode_mem, tolsfac)
    ccall((:ARKodeGetTolScaleFactor, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{realtype}), arkode_mem, tolsfac)
end

function __ARKodeGetErrWeights(arkode_mem, eweight::N_Vector)
    ccall((:ARKodeGetErrWeights, libsundials_arkode), Cint, (ARKODEMemPtr, N_Vector), arkode_mem, eweight)
end

function ARKodeGetErrWeights(arkode_mem, eweight)
    __eweight = convert(NVector, eweight)
    __ARKodeGetErrWeights(arkode_mem, convert(N_Vector, __eweight))
end

function __ARKodeGetResWeights(arkode_mem, rweight::N_Vector)
    ccall((:ARKodeGetResWeights, libsundials_arkode), Cint, (ARKODEMemPtr, N_Vector), arkode_mem, rweight)
end

function ARKodeGetResWeights(arkode_mem, rweight)
    __rweight = convert(NVector, rweight)
    __ARKodeGetResWeights(arkode_mem, convert(N_Vector, __rweight))
end

function __ARKodeGetEstLocalErrors(arkode_mem, ele::N_Vector)
    ccall((:ARKodeGetEstLocalErrors, libsundials_arkode), Cint, (ARKODEMemPtr, N_Vector), arkode_mem, ele)
end

function ARKodeGetEstLocalErrors(arkode_mem, ele)
    __ele = convert(NVector, ele)
    __ARKodeGetEstLocalErrors(arkode_mem, convert(N_Vector, __ele))
end

function ARKodeGetNumGEvals(arkode_mem, ngevals)
    ccall((:ARKodeGetNumGEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, ngevals)
end

function ARKodeGetRootInfo(arkode_mem, rootsfound)
    ccall((:ARKodeGetRootInfo, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Cint}), arkode_mem, rootsfound)
end

function ARKodeGetIntegratorStats(arkode_mem, nsteps, expsteps, accsteps, step_attempts, nfe_evals, nfi_evals, nlinsetups, netfails, hinused, hlast, hcur, tcur)
    ccall((:ARKodeGetIntegratorStats, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), arkode_mem, nsteps, expsteps, accsteps, step_attempts, nfe_evals, nfi_evals, nlinsetups, netfails, hinused, hlast, hcur, tcur)
end

function ARKodeGetNumNonlinSolvIters(arkode_mem, nniters)
    ccall((:ARKodeGetNumNonlinSolvIters, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nniters)
end

function ARKodeGetNumNonlinSolvConvFails(arkode_mem, nncfails)
    ccall((:ARKodeGetNumNonlinSolvConvFails, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nncfails)
end

function ARKodeGetNonlinSolvStats(arkode_mem, nniters, nncfails)
    ccall((:ARKodeGetNonlinSolvStats, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, nniters, nncfails)
end

function __ARKodeLoadButcherTable(imethod::Cint, s, q, p, A, b, c, b2)
    ccall((:ARKodeLoadButcherTable, libsundials_arkode), Cint, (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), imethod, s, q, p, A, b, c, b2)
end

function ARKodeLoadButcherTable(imethod, s, q, p, A, b, c, b2)
    __ARKodeLoadButcherTable(convert(Cint, imethod), s, q, p, A, b, c, b2)
end

function __ARKodeGetReturnFlagName(flag::Clong)
    ccall((:ARKodeGetReturnFlagName, libsundials_arkode), Cstring, (Clong,), flag)
end

function ARKodeGetReturnFlagName(flag)
    __ARKodeGetReturnFlagName(convert(Clong, flag))
end

function ARKodeWriteParameters(arkode_mem, fp)
    ccall((:ARKodeWriteParameters, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{FILE}), arkode_mem, fp)
end

function ARKodeWriteButcher(arkode_mem, fp)
    ccall((:ARKodeWriteButcher, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{FILE}), arkode_mem, fp)
end

function ARKodeFree(arkode_mem)
    ccall((:ARKodeFree, libsundials_arkode), Cvoid, (Ref{ARKODEMemPtr},), arkode_mem)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode_spils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function ARKSpilsSetLinearSolver(arkode_mem, LS::SUNLinearSolver)
    ccall((:ARKSpilsSetLinearSolver, libsundials_arkode), Cint, (ARKODEMemPtr, SUNLinearSolver), arkode_mem, LS)
end

function __ARKSpilsSetMassLinearSolver(arkode_mem, LS::SUNLinearSolver, time_dep::Cint)
    ccall((:ARKSpilsSetMassLinearSolver, libsundials_arkode), Cint, (ARKODEMemPtr, SUNLinearSolver, Cint), arkode_mem, LS, time_dep)
end

function ARKSpilsSetMassLinearSolver(arkode_mem, LS, time_dep)
    __ARKSpilsSetMassLinearSolver(arkode_mem, LS, convert(Cint, time_dep))
end

function ARKSpilsSetEpsLin(arkode_mem, eplifac::realtype)
    ccall((:ARKSpilsSetEpsLin, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, eplifac)
end

function ARKSpilsSetMassEpsLin(arkode_mem, eplifac::realtype)
    ccall((:ARKSpilsSetMassEpsLin, libsundials_arkode), Cint, (ARKODEMemPtr, realtype), arkode_mem, eplifac)
end

function ARKSpilsSetPreconditioner(arkode_mem, psetup::ARKSpilsPrecSetupFn, psolve::ARKSpilsPrecSolveFn)
    ccall((:ARKSpilsSetPreconditioner, libsundials_arkode), Cint, (ARKODEMemPtr, ARKSpilsPrecSetupFn, ARKSpilsPrecSolveFn), arkode_mem, psetup, psolve)
end

function ARKSpilsSetMassPreconditioner(arkode_mem, psetup::ARKSpilsMassPrecSetupFn, psolve::ARKSpilsMassPrecSolveFn)
    ccall((:ARKSpilsSetMassPreconditioner, libsundials_arkode), Cint, (ARKODEMemPtr, ARKSpilsMassPrecSetupFn, ARKSpilsMassPrecSolveFn), arkode_mem, psetup, psolve)
end

function ARKSpilsSetJacTimes(arkode_mem, jtsetup::ARKSpilsJacTimesSetupFn, jtimes::ARKSpilsJacTimesVecFn)
    ccall((:ARKSpilsSetJacTimes, libsundials_arkode), Cint, (ARKODEMemPtr, ARKSpilsJacTimesSetupFn, ARKSpilsJacTimesVecFn), arkode_mem, jtsetup, jtimes)
end

function ARKSpilsSetMassTimes(arkode_mem, msetup::ARKSpilsMassTimesSetupFn, mtimes::ARKSpilsMassTimesVecFn, mtimes_data)
    ccall((:ARKSpilsSetMassTimes, libsundials_arkode), Cint, (ARKODEMemPtr, ARKSpilsMassTimesSetupFn, ARKSpilsMassTimesVecFn, Ptr{Cvoid}), arkode_mem, msetup, mtimes, mtimes_data)
end

function ARKSpilsGetWorkSpace(arkode_mem, lenrwLS, leniwLS)
    ccall((:ARKSpilsGetWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwLS, leniwLS)
end

function ARKSpilsGetNumPrecEvals(arkode_mem, npevals)
    ccall((:ARKSpilsGetNumPrecEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, npevals)
end

function ARKSpilsGetNumPrecSolves(arkode_mem, npsolves)
    ccall((:ARKSpilsGetNumPrecSolves, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, npsolves)
end

function ARKSpilsGetNumLinIters(arkode_mem, nliters)
    ccall((:ARKSpilsGetNumLinIters, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nliters)
end

function ARKSpilsGetNumConvFails(arkode_mem, nlcfails)
    ccall((:ARKSpilsGetNumConvFails, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nlcfails)
end

function ARKSpilsGetNumJTSetupEvals(arkode_mem, njtsetups)
    ccall((:ARKSpilsGetNumJTSetupEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, njtsetups)
end

function ARKSpilsGetNumJtimesEvals(arkode_mem, njvevals)
    ccall((:ARKSpilsGetNumJtimesEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, njvevals)
end

function ARKSpilsGetNumRhsEvals(arkode_mem, nfevalsLS)
    ccall((:ARKSpilsGetNumRhsEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nfevalsLS)
end

function ARKSpilsGetLastFlag(arkode_mem, flag)
    ccall((:ARKSpilsGetLastFlag, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, flag)
end

function ARKSpilsGetMassWorkSpace(arkode_mem, lenrwMLS, leniwMLS)
    ccall((:ARKSpilsGetMassWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwMLS, leniwMLS)
end

function ARKSpilsGetNumMassPrecEvals(arkode_mem, nmpevals)
    ccall((:ARKSpilsGetNumMassPrecEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmpevals)
end

function ARKSpilsGetNumMassPrecSolves(arkode_mem, nmpsolves)
    ccall((:ARKSpilsGetNumMassPrecSolves, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmpsolves)
end

function ARKSpilsGetNumMassIters(arkode_mem, nmiters)
    ccall((:ARKSpilsGetNumMassIters, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmiters)
end

function ARKSpilsGetNumMassConvFails(arkode_mem, nmcfails)
    ccall((:ARKSpilsGetNumMassConvFails, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmcfails)
end

function ARKSpilsGetNumMtimesEvals(arkode_mem, nmvevals)
    ccall((:ARKSpilsGetNumMtimesEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmvevals)
end

function ARKSpilsGetLastMassFlag(arkode_mem, flag)
    ccall((:ARKSpilsGetLastMassFlag, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, flag)
end

function __ARKSpilsGetReturnFlagName(flag::Clong)
    ccall((:ARKSpilsGetReturnFlagName, libsundials_arkode), Cstring, (Clong,), flag)
end

function ARKSpilsGetReturnFlagName(flag)
    __ARKSpilsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode_bandpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function ARKBandPrecInit(arkode_mem, N::sunindextype, mu::sunindextype, ml::sunindextype)
    ccall((:ARKBandPrecInit, libsundials_arkode), Cint, (ARKODEMemPtr, sunindextype, sunindextype, sunindextype), arkode_mem, N, mu, ml)
end

function ARKBandPrecGetWorkSpace(arkode_mem, lenrwLS, leniwLS)
    ccall((:ARKBandPrecGetWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwLS, leniwLS)
end

function ARKBandPrecGetNumRhsEvals(arkode_mem, nfevalsBP)
    ccall((:ARKBandPrecGetNumRhsEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nfevalsBP)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode_bbdpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function ARKBBDPrecInit(arkode_mem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dqrely::realtype, gloc::ARKLocalFn, cfn::ARKCommFn)
    ccall((:ARKBBDPrecInit, libsundials_arkode), Cint, (ARKODEMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, ARKLocalFn, ARKCommFn), arkode_mem, Nlocal, mudq, mldq, mukeep, mlkeep, dqrely, gloc, cfn)
end

function ARKBBDPrecReInit(arkode_mem, mudq::sunindextype, mldq::sunindextype, dqrely::realtype)
    ccall((:ARKBBDPrecReInit, libsundials_arkode), Cint, (ARKODEMemPtr, sunindextype, sunindextype, realtype), arkode_mem, mudq, mldq, dqrely)
end

function ARKBBDPrecGetWorkSpace(arkode_mem, lenrwLS, leniwLS)
    ccall((:ARKBBDPrecGetWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwLS, leniwLS)
end

function ARKBBDPrecGetNumGfnEvals(arkode_mem, ngevalsBBDP)
    ccall((:ARKBBDPrecGetNumGfnEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, ngevalsBBDP)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function ARKDlsSetLinearSolver(arkode_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:ARKDlsSetLinearSolver, libsundials_arkode), Cint, (ARKODEMemPtr, SUNLinearSolver, SUNMatrix), arkode_mem, LS, A)
end

function __ARKDlsSetMassLinearSolver(arkode_mem, LS::SUNLinearSolver, M::SUNMatrix, time_dep::Cint)
    ccall((:ARKDlsSetMassLinearSolver, libsundials_arkode), Cint, (ARKODEMemPtr, SUNLinearSolver, SUNMatrix, Cint), arkode_mem, LS, M, time_dep)
end

function ARKDlsSetMassLinearSolver(arkode_mem, LS, M, time_dep)
    __ARKDlsSetMassLinearSolver(arkode_mem, LS, M, convert(Cint, time_dep))
end

function ARKDlsSetJacFn(arkode_mem, jac::ARKDlsJacFn)
    ccall((:ARKDlsSetJacFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKDlsJacFn), arkode_mem, jac)
end

function ARKDlsSetMassFn(arkode_mem, mass::ARKDlsMassFn)
    ccall((:ARKDlsSetMassFn, libsundials_arkode), Cint, (ARKODEMemPtr, ARKDlsMassFn), arkode_mem, mass)
end

function ARKDlsGetWorkSpace(arkode_mem, lenrwLS, leniwLS)
    ccall((:ARKDlsGetWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwLS, leniwLS)
end

function ARKDlsGetNumJacEvals(arkode_mem, njevals)
    ccall((:ARKDlsGetNumJacEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, njevals)
end

function ARKDlsGetNumRhsEvals(arkode_mem, nfevalsLS)
    ccall((:ARKDlsGetNumRhsEvals, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nfevalsLS)
end

function ARKDlsGetLastFlag(arkode_mem, flag)
    ccall((:ARKDlsGetLastFlag, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, flag)
end

function ARKDlsGetMassWorkSpace(arkode_mem, lenrwMLS, leniwMLS)
    ccall((:ARKDlsGetMassWorkSpace, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}, Ptr{Clong}), arkode_mem, lenrwMLS, leniwMLS)
end

function ARKDlsGetNumMassSetups(arkode_mem, nmsetups)
    ccall((:ARKDlsGetNumMassSetups, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmsetups)
end

function ARKDlsGetNumMassSolves(arkode_mem, nmsolves)
    ccall((:ARKDlsGetNumMassSolves, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmsolves)
end

function ARKDlsGetNumMassMult(arkode_mem, nmmults)
    ccall((:ARKDlsGetNumMassMult, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, nmmults)
end

function ARKDlsGetLastMassFlag(arkode_mem, flag)
    ccall((:ARKDlsGetLastMassFlag, libsundials_arkode), Cint, (ARKODEMemPtr, Ptr{Clong}), arkode_mem, flag)
end

function __ARKDlsGetReturnFlagName(flag::Clong)
    ccall((:ARKDlsGetReturnFlagName, libsundials_arkode), Cstring, (Clong,), flag)
end

function ARKDlsGetReturnFlagName(flag)
    __ARKDlsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/arkode/arkode_impl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0


text:
# generic functions
Base.hash(x::Sym, v::UInt64) = hash(getpy(x), v)
Base.hash(x::Sym) = hash(getpy(x))
Base.:(==)(x::Sym, y::Sym) = Bool(getpy(x) == getpy(y))

Base.isless(a::Sym, b::Sym) = (a != sympy.nan && b != sympy.nan) && sympy.Lt(a,b) == sympy.True
Base.isless(a::Sym, b::Number) = isless(promote(a,b)...)
Base.isless(a::Number, b::Sym) = isless(promote(a,b)...)

zero(::Sym) = zero(Sym)
zero(::Type{Sym}) = Sym(0)
one(::Sym) = one(Sym)
one(::Type{Sym}) = Sym(1)

# math functions which are not evaled into instanc
Base.abs2(x::SymbolicObject) = x * conj(x)
Base.cbrt(x::Sym) = x^(1//3)
Base.asech(z::Sym) = log(sqrt(1/z-1)*sqrt(1/z+1) + 1/z)
Base.acsch(z::Sym) = log(sqrt(1+1/z^2) + 1/z) ## http://mathworld.wolfram.com/InverseHyperbolicCosecant.html
Base.atan(y::Sym, x) = sympy.atan2(y,x)

Base.sinc(x::Sym) = iszero(x) ? one(x) : sin(PI*x)/(PI*x)
cosc(x::Sym) = diff(sinc(x))

Base.sincos(x::Sym) = (sin(x), cos(x))
Base.sinpi(x::Sym) = sin(PI * x)
Base.cospi(x::Sym) = cos(PI * x)
# XXX do cosd, ...?

Base.rad2deg(x::Sym) = (x * 180) / PI
Base.deg2rad(x::Sym) = (x * PI) / 180

Base.hypot(x::Sym, y::Number) = hypot(promote(x,y)...)
Base.hypot(xs::Sym...) = sqrt(sum(abs(xᵢ)^2 for xᵢ ∈ xs))

## exponential/Log
Base.log1p(x::Sym) = log(1 + x)
Base.log(b::Number, x::Sym) = sympy.log(x, b)
Base.log2(x::SymbolicObject) = log(2, x)
Base.log10(x::SymbolicObject) = log(10, x)

## --------------------------------------------------
## floating point stuff

# Floating point bits
Base.eps(::Type{Sym}) = zero(Sym)
Base.eps(::Sym) = zero(Sym)
Base.signbit(x::Sym) = x < 0
Base.copysign(x::Sym, y::Number) = abs(x) * sign(y)
Base.flipsign(x::Sym, y) = signbit(y) ? -x : x
Base.typemax(::Type{Sym}) = oo
Base.typemin(::Type{Sym}) = -oo

Base.fld(x::SymbolicObject, y) = floor(x/y)
Base.cld(x::SymbolicObject, y) = ceil(x/y)

#

Base.isfinite(x::Sym) = !is_(:infinite, x)
Base.isinf(x::Sym) = is_(:infinite, x)
Base.isnan(x::Sym) = x == sympy.nan
Base.isinteger(x::Sym) = is_(:integer, x)
Base.iseven(x::Sym) = is_(:even, x)
Base.isreal(x::Sym) = is_(:real, x)
Base.isodd(x::Sym) = is_(:odd, x)

# XXX right semantics?
function Base.float(x::Sym)
    u = N(x)
    u === x && return x
    float(u)
end

Base.Integer(x::Sym) = is_(:integer, x) ? N(x) : throw(DomainError("x can not be converted to an integer"))

Base.complex(::Type{Sym}) = Sym
Base.complex(r::Sym) = real(r) + imag(r) * im
function Base.complex(r::Sym, i)
    isreal(r) || throw(ArgumentError("r and i must not be complex"))
    isreal(i) || throw(ArgumentError("r and i must not be complex"))
    N(r) + N(i) * im
end
Base.complex(xs::AbstractArray{Sym}) = complex.(xs) # why is this in base?

Base.real(::Type{Sym}) = Sym
Base.real(x::Sym) = sympy.re(x)
Base.imag(x::Sym) =  sympy.im(x)

Base.angle(z::SymbolicObject) = atan(sympy.im(z), sympy.re(z))

Base.divrem(x::Sym, y) = sympy.div(x, y)
# needed for #390; but odd
Base.div(x::Sym, y::Union{Sym,Number}) = convert(Sym, sympy.floor(x/convert(Sym,y)))
Base.rem(x::Sym, y::Union{Sym,Number}) = x-Sym(y)*Sym(sympy.floor.(x/y))

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
function _precompile_()
  ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
  @assert Base.precompile(Tuple{typeof(print_timer), typeof(stdout), TimerOutput})
  @assert Base.precompile(Tuple{typeof(print_timer), TimerOutput})
  @assert Base.precompile(Tuple{typeof(push!), TimerOutput, String})
  @assert Base.precompile(Tuple{typeof(reset_timer!), TimerOutput})
  @assert Base.precompile(Tuple{typeof(disable_timer!), TimerOutput})
  @assert Base.precompile(Tuple{typeof(enable_timer!), TimerOutput})
  @assert Base.precompile(Tuple{typeof(complement!), TimerOutput})
  @assert Base.precompile(Tuple{typeof(do_accumulate!), TimeData, UInt64, Int64})
  @assert Base.precompile(Tuple{Type{TimerOutput}, String})
end


text:
module OpenTelemetryAPI

using TOML

const PKG_VERSION =
    VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])

include("common/common.jl")
include("trace/trace.jl")
include("propagator/propagator.jl")
include("metric/metric.jl")
include("log.jl")

function __init__()
    push!(GLOBAL_PROPAGATOR, TraceContextTextMapPropagator())
end

end # module


text:
using Pkg

function install(name)
    installed_packages = keys(Pkg.installed())
    if name ∉ installed_packages
        Pkg.add(name)
    end
end

function includedir(dir)
	for (root, dirs, files) in walkdir(dir)
		for f in files
			include(root*"/"*f)
		end
	end
end


text:
#
# implementation of Spline-based pair potentials
# and other functionality
#

import Interpolations
using Interpolations: interpolate, BSpline, Cubic, Line, OnGrid,
                      ScaledInterpolation, BSplineInterpolation
using OffsetArrays: OffsetVector
using ForwardDiff

const CubicBSpline{T} = ScaledInterpolation{T, 1, BSplineInterpolation{T, 1, OffsetVector{T, Vector{T}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{T, Base.TwicePrecision{T}, Base.TwicePrecision{T}}}}

"""
`type SplinePairPotential`

Pair potential defined via splines. In the construction of these splines
atomic units (Angstrom) are assumed. If this is used in different units,
then something will probably break.

### Constructors:
```
SplinePairPotential(xdat, ydat; kwargs...)   # fit spline  from data points
SplinePairPotential(fname; kwargs...)        # load data points from file
```
Keyword arguments:
* `fixcutoff = true`: if true, then the fit will be artificially modified at the cut-off to ensure that the transition to zero is smooth
"""
mutable struct SplinePairPotential{T} <: SimplePairPotential
   spl::CubicBSpline{T}      # The actual spline object
   rcut::T                   # cutoff radius (??? could just use spl.t[end] ???)
end

@pot SplinePairPotential

function SplinePairPotential(xdat, ydat, rcut = xdat[end]; fixcutoff=true)
   # check that the data is equi-spaced. If not, then we need to generalize the
   # code to allow splines on arbitrary grids (griddedinterpolations?)
   if maximum(abs, diff(diff(xdat))) > 1e-12
      @show xdat[1:10]
      @show xdat[end-10:end]
      @show maximum(abs, diff(diff(xdat)))
      error("spline knots arent equispaced")
   end
   h = xdat[2] - xdat[1]
   # add xdat points to the left to make sure zero is included
   addn = ceil(Int, xdat[1]/h)
   xrg = range(xdat[1] - addn * h, xdat[end]+1.5*h, step=h)
   ydat = [ fill(ydat[1], addn); ydat; [0.0] ]
   if fixcutoff
      # set the last few spline data-points to 0 to get a guaranteed
      # smooth transition to 0 at the cutoff.
      ydat[end-2:end] .= 0.0   # end-order=1:end
   end
   # construct the spline
   spl_pre = interpolate(ydat, BSpline(Cubic(Line(OnGrid()))))
   # ... and rescale it to the correct grid
   spl = Interpolations.scale(spl_pre, xrg)

   return SplinePairPotential(spl, rcut)
end

SplinePairPotential(fname::AbstractString; kwargs...) =
      SplinePairPotential(load_ppfile(fname)...; kwargs...)

cutoff(V::SplinePairPotential) = V.rcut

evaluate(V::SplinePairPotential{T}, r::Number) where {T} =
      r < V.rcut ? V.spl(r) : zero(T)

evaluate_d(V::SplinePairPotential{T}, r::Number) where {T} =
      r < V.rcut ? Interpolations.gradient(V.spl, r)[1] : zero(T)

evaluate_dd(V::SplinePairPotential, r::Number) =
      r < V.rcut ? Interpolations.hessian(V.spl, r)[1] : zero(T)

function load_ppfile(fname::AbstractString)
   if fname[end-3:end] == ".plt"
      return load_plt(fname)
   else
      error("unknown file type: $(fname[end-3:end])")
   end
   nothing
end

function load_plt(fname::AbstractString)
   data = readdlm(fname, Float64; comments=true)
   return data[:, 1], data[:, 2]
end