text:
abstract type AbstractExpansion{T} end

struct GradientExpansion{T,N,M} <: AbstractExpansion{T}
	x::SizedVector{N,T,Vector{T}}
	u::SizedVector{M,T,Vector{T}}
	function GradientExpansion{T}(n::Int,m::Int) where T
		new{T,n,m}(SizedVector{n}(zeros(T,n)), SizedVector{m}(zeros(T,m)))
	end
end

# TODO: Move to ALTRO
"""
	DynamicsExpansion{T,N,N̄,M}

Stores the dynamics expansion for a single time instance.
For a `LieGroupModel`, it will provide access to both the state and state
error Jacobians.

# Constructors
```julia
DynamicsExpansion{T}(n0, n, m)
DynamicsExpansion{T}(n, m)
```
where `n0` is the size of the full state, and `n` is the size of the error state.

# Methods
To evaluate the dynamics Jacobians, use

	dynamics_expansion!(::Type{Q}, D::DynamicsExpansion, model, Z)

To compute the Jacobians for the error state, use

	error_expansion!(D::DynamicsExpansion, model, G)

where `G` is a vector of error-state Jacobians. These can be computed using
`RobotDynamics.state_diff_jacobian(G, model, Z)`.

# Extracting Jacobians
The Jacobians should be extracted using

	fdx, fdu = error_expansion(D::DynamicsExpansion, model)

This method will provide the error state Jacobians for `LieGroupModel`s, and
	the normal Jacobian otherwise. Both `fdx` and `fdu` are a `SizedMatrix`.
"""
struct DynamicsExpansion{T,N,N̄,M}
	∇f::Matrix{T} # n × (n+m)
	∇²f::Matrix{T}  # (n+m) × (n+m)
	A_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	B_::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	A::SizedMatrix{N̄,N̄,T,2,Matrix{T}}
	B::SizedMatrix{N̄,M,T,2,Matrix{T}}
	fxx::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	fuu::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	fux::SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}
	tmpA::SizedMatrix{N,N,T,2,Matrix{T}}
	tmpB::SizedMatrix{N,M,T,2,Matrix{T}}
	tmp::SizedMatrix{N,N̄,T,2,Matrix{T}}
	function DynamicsExpansion{T}(n0::Int, n::Int, m::Int) where T
		∇f = zeros(n0,n0+m)
		∇²f = zeros(n0+m,n0+m)
		ix = 1:n0
		iu = n0 .+ (1:m)
		A_ = view(∇f, ix, ix)
		B_ = view(∇f, ix, iu)
		A = SizedMatrix{n,n}(zeros(n,n))
		B = SizedMatrix{n,m}(zeros(n,m))
		fxx = view(∇²f, ix, ix)
		fuu = view(∇²f, iu, iu)
		fux = view(∇²f, iu, ix)
		tmpA = SizedMatrix{n0,n0}(zeros(n0,n0))
		tmpB = SizedMatrix{n0,m}(zeros(n0,m))
		tmp = SizedMatrix{n0,n}(zeros(n0,n))
		new{T,n0,n,m}(∇f,∇²f,A_,B_,A,B,fxx,fuu,fux,tmpA,tmpB,tmp)
	end
	function DynamicsExpansion{T}(n::Int, m::Int) where T
		∇f = zeros(n,n+m)
		∇²f = zeros(n+m,n+m)
		ix = 1:n
		iu = n .+ (1:m)
		A_ = view(∇f, ix, ix)
		B_ = view(∇f, ix, iu)
		A = SizedMatrix{n,n}(zeros(n,n))
		B = SizedMatrix{n,m}(zeros(n,m))
		fxx = view(∇²f, ix, ix)
		fuu = view(∇²f, iu, iu)
		fux = view(∇²f, iu, ix)
		tmpA = A
		tmpB = B
		tmp = zeros(n,n)
		new{T,n,n,m}(∇f,∇²f,A_,B_,A,B,fxx,fuu,fux,tmpA,tmpB,tmp)
	end
end

function save_tmp!(D::DynamicsExpansion)
	D.tmpA .= D.A_
	D.tmpB .= D.B_
end

function dynamics_expansion!(Q, D::Vector{<:DynamicsExpansion}, model::AbstractModel,
		Z::Traj, args...)
	for k in eachindex(D)
		RobotDynamics.discrete_jacobian!(Q, D[k].∇f, model, Z[k], args...)
		# save_tmp!(D[k])
		# D[k].tmpA .= D[k].A_  # avoids allocations later
		# D[k].tmpB .= D[k].B_
	end
end

function error_expansion!(D::DynamicsExpansion,G1,G2)
    mul!(D.tmp, D.tmpA, G1)
    mul!(D.A, Transpose(G2), D.tmp)
    mul!(D.B, Transpose(G2), D.tmpB)
end

@inline error_expansion(D::DynamicsExpansion, model::LieGroupModel) = D.A, D.B
@inline error_expansion(D::DynamicsExpansion, model::AbstractModel) = D.tmpA, D.tmpB

@inline DynamicsExpansion(model::AbstractModel) = DynamicsExpansion{Float64}(model)
@inline function DynamicsExpansion{T}(model::AbstractModel) where T
	n,m = size(model)
	n̄ = state_diff_size(model)
	DynamicsExpansion{T}(n,n̄,m)
end

function error_expansion!(D::Vector{<:DynamicsExpansion}, model::AbstractModel, G)
	for d in D
		save_tmp!(d)
	end
end

function error_expansion!(D::Vector{<:DynamicsExpansion}, model::LieGroupModel, G)
	for k in eachindex(D)
		save_tmp!(D[k])
		error_expansion!(D[k], G[k], G[k+1])
	end
end

struct StaticExpansion{T,N,M,NN,MM,NM}
	x::SVector{N,T}
	xx::SMatrix{N,N,T,NN}
	u::SVector{M,T}
	uu::SMatrix{M,M,T,MM}
	ux::SMatrix{M,N,T,NM}
end

function StaticExpansion(E::AbstractExpansion)
	StaticExpansion(SVector(E.x), SMatrix(E.xx),
		SVector(E.u), SMatrix(E.uu), SMatrix(E.ux))
end

function StaticExpansion(x,xx,u,uu,ux)
	StaticExpansion(SVector(x), SMatrix(xx), SVector(u), SMatrix(uu), SMatrix(ux))
end

"""
	Expansion{n,m,T}

Stores a full second-order expansion of a scalar-valued function (e.g. the cost function)
with respect to both the state and control. The Hessian and gradient with respect to the
concatenated state and control, available via the `hess` and `grad` fields.

The pieces with respect to the state and control separately are available via
`x`, `xx`, `u`, `uu`, and `ux`, which are aliased (for backward compatibility) with
`q`, `Q`, `r`, `R`, and `H`.
"""
struct Expansion{n,m,T}
    # not sure why calling it with
    res::DiffResults.MutableDiffResult{2,T,Tuple{Vector{T},Matrix{T}}}
    hess::Matrix{T}
    grad::Vector{T}
    x::SizedVector{n,T,SubArray{T,1,Vector{T},Tuple{UnitRange{Int}},true}}
    u::SizedVector{m,T,SubArray{T,1,Vector{T},Tuple{UnitRange{Int}},true}}
    xx::SizedMatrix{n,n,T,2,SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
    uu::SizedMatrix{m,m,T,2,SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
    ux::SizedMatrix{m,n,T,2,SubArray{T,2,Matrix{T},Tuple{UnitRange{Int},UnitRange{Int}},false}}
    function Expansion{T}(n::Int, m::Int) where T
        ix,iu = 1:n, n .+ (1:m)
        res = DiffResults.HessianResult(zeros(T,n+m))
        grad = DiffResults.gradient(res)
        x = SizedVector{n}(view(grad,ix))
        u = SizedVector{m}(view(grad,iu))
        hess = DiffResults.hessian(res)
        xx = SizedMatrix{n,n}(view(hess,ix,ix))
        uu = SizedMatrix{m,m}(view(hess,iu,iu))
		ux = SizedMatrix{m,n}(view(hess,iu,ix))
		hess .= I(n+m)
		grad .= 0
        new{n,m,T}(res, hess, grad, x, u, xx, uu, ux)
    end
end
function Base.getproperty(E::Expansion, field::Symbol)
    if field == :q
        getfield(E, :x)
    elseif field == :r
        getfield(E, :u)
    elseif field == :Q
        getfield(E, :xx)
    elseif field == :R
        getfield(E, :uu)
    elseif field == :H
        getfield(E, :ux)
    else
        getfield(E, field)
    end
end

static_expansion(E::Expansion) = StaticExpansion(E.x, E.xx, E.u, E.uu, E.ux)


text:
using Mimi
using MimiDICE2016

DICE = MimiDICE2016.get_model()
run(DICE)

explore(DICE)

ext: jl


text:
using Optim
using ODEOpt
using ODE
using PyPlot
using Base.Test
using NLsolve
using ForwardDiff
PLOT = true
let
    ystart = [0.0;1.0]
    trange = linspace(0,2π,100)
    a, b, c = 1.0, 2.0, 3.0
    foo = secondorder(a,b,c)
    @test foo(0,[0;1]) == [1.0;-2.0]
end

let
    # Test based on the example at
    # http://www.sosmath.com/diffeq/second/constantcof/constantcof.html
    #
    #     x'' + 2x' + 2x = 0
    #     x(π/4) = 2; x'(π/4) = -2
    #
    #     convert to first order using the system
    #     y[1] = x; y[2] = x'
    #
    #     descretize at linspace(π/4,2π,100)
    #
    # solve the ode numerically then plug into the formulas to
    # check that the formulas match the ode45 solution.
    ystart = [2.0;-2.0]
    trange = linspace(π/4,2π,100)
    a, b, c = 1.0, 2.0, 2.0
    foo  = secondorder(a,b,c)
    sol = solve45(foo, ystart, collect(trange))
    result, formula, deriv = ivp(a,b,c, ystart, π/4)
    @show params = result.zero
    @show result.trace
    @show result
    # figure()
    # test that the computed solution matches the formula
    form = [formula(t,params...) for t in sol.time]
    deri = [deriv(t,params...) for t in sol.time]
    @test norm(form-[y[1] for y in sol.y]) <= 1e-3
    @test norm(deri-[y[2] for y in sol.y]) <= 1e-3
    # plot general solution.
    if PLOT
        plot(sol.time, form, label="formula")
        plot(sol.time, deri, label="f'")
        title("parameter fitting \$ $a \\ddot{x} + $b \\dot{x} + $c x = 0\$")
        solutionplot(sol)
    end
end

let
    # parameter fitting for 2nd order linear ode
    # the equation is ax'' + bx' + cx = 0
    # we want to recover the coefficients a,b,c as well as the initial values x₀,v₀
    # generate "exact" solution
    info("Testing Parameter Fitting")
    ystart = [2.0;-2.0]
    trange = linspace(π/4,2π,100)
    a, b, c = 1.0, 2.0, 2.0
    params = [a,b,c]
    # initial parameter guess must be fairly close to avoid local optima
    param₀ = 1.5*ones(3)
    #param₀ = params .+ randn(3)/5

    foo  = secondorder(a,b,c)
    sol = solve45(foo, ystart, collect(trange))
    result, formula, deriv = ivp(a,b,c, ystart, π/4)
    consts = result.zero
    form = Float64[formula(t,consts...) for t in sol.time]
    deri = Float64[deriv(t,consts...) for t in sol.time]

    # add noise to get data as if measured
    srand(1)
    noisev = 0.10 * randn(length(form))
    data = form .+ noisev
    # x = [y[1] for y in sol.y]
    # make sure the soserror matches the norm of the measurement noise.
    # @show sqrt(soserr(sol.time,form, sol.time, data)) norm(noisev) 1e-2
    tsolve_opt, x⁺_opt, v_opt, d²_opt = obj(foo, ystart, sol.time, data)
    @test_approx_eq_eps sqrt(d²_opt) norm(noisev) 1e-2

    # fit the parameters to the data
    @show param₀
    γ(x) = sqrt(obj(secondorder(x...), ystart, sol.time, data)[end])
    @show res = optimize(γ,
                         param₀,
                         xtol=1e-12,
                         ftol=1e-32,
                         method=ConjugateGradient(),
                         #method=GradientDescent(),
                         #method=LBFGS(),
                         autodiff=false,
                         show_trace=true,
                         iterations=1000,
                         )
    # report the quality of the fit
    @show param⁼ = res.minimum
    @show res.f_minimum d²_opt
    tsolve, x⁺, v⁺, d² = obj(secondorder(param⁼...), ystart, sol.time, data)
    if PLOT
        figure()
        plot(tsolve, x⁺, label="\$x^+_o\$")
        plot(sol.time,form,label="\$x_f\$")
        plot(tsolve, x⁺, label="\$x^+\$")
        plot(sol.time, data, ".",  label="data")
        title("parameter fitting \$ $a \\ddot{x} + $b \\dot{x} + $c x = 0\$")
        legend()
    end
    # check that residual is small
    @show sqrt(soserr(tsolve, x⁺, sol.time, form))
    @test_approx_eq_eps γ(res.minimum) sqrt(d²_opt) 2e-1
    #check that error is small
    @show err = norm(param⁼ - params)
    @test err < 3e-1
end

max_stars_repo_name: jpfairbanks/ODEOpt

max_issues_repo_name: jpfairbanks/ODEOpt

max_forks_repo_name: jpfairbanks/ODEOpt


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

using Test

Test.@testset "densities" begin
    include("test_logdval.jl")
    include("test_abstract_density.jl")
    include("test_generic_density.jl")
    include("test_const_density.jl")
    include("test_density_product.jl")
    include("test_distribution_density.jl")
    include("test_parameter_mapped_density.jl")
    include("test_renormalize_density.jl")
    include("test_truncate_density.jl")
    include("test_transformed_density.jl")
    include("test_external_density.jl")
end


text:
module PlotHelper
using Plots
using Statistics

export plotBSTsample, saveBSTSplot, plotBSTsummary, plotPolarSDM
export format_wavename

function plotBSTsample(ind::Int, SDMs, SDM, SDM_RD, SDM_RD_err)
    bts = [real(s[ind,ind]) for s in SDMs]
    stephist(bts, bins=20, lab="bootstrap", title="SDM[$(ind),$(ind)]")
    vline!([real(SDM[ind,ind])], l=(:red, 2), lab="main fit")
    vline!(real.([SDM_RD[ind,ind]] .+ SDM_RD_err[ind,ind] .* [-1,1]), l=(:orange, 1.5), lab="official fit")
    vline!(quantile(bts,[0.16,0.84]), lab="quantile")
    plot!()
end

function plotBSTsample(func, SDMs, SDM, SDM_RD, SDM_RD_err)
    bts = [func(s) for s in SDMs]
    stephist(bts, bins=20, lab="bootstrap", title="")
    vline!([func(SDM)], l=(:red, 2), lab="main fit")
    vline!([func(SDM_RD)], l=(:orange, 1.5), lab="official fit")
    vline!(quantile(bts,[0.16,0.84]), lab="quantile")
    plot!()
end

function saveBSTSplot(fout, SDMs, SDM, SDM_RD, SDM_RD_err)
    Nw = size(SDM,2)
    onames = Array{String}(Nw)
    for i in 1:Nw
        plotBSTsample(i, SDMs, SDM, SDM_RD, SDM_RD_err)
        onames[i] = fout * "-tmp-$(i).pdf"
        savefig(onames[i])
    end
    run(`pdfunite $onames $fout`)
    run(`rm $onames`)
end

function construct_values(ind, SDMs, SDM, SDM_RD, SDM_RD_err)
    bts = [real(s[ind,ind]) for s in SDMs]
    v_main = real(SDM[ind,ind])
    v_off = real(SDM_RD[ind,ind])
    v_off_err = real(SDM_RD_err[ind,ind])
    v_qnt = quantile(bts,[0.16,0.84]) # sigma to every direction
    [ind v_main v_off v_off_err v_qnt...]
end

function plotBSTsummary(SDMs, SDM, SDM_RD, SDM_RD_err; tosort=false, toannotate=false)
    combres = vcat([construct_values(i, SDMs, SDM, SDM_RD, SDM_RD_err) for i in 1:size(SDM,1)]...);
    plotBSTsummary(combres; tosort=tosort, toannotate=toannotate)
end

function plotBSTsummary(combres; tosort=false, toannotate=false)
    Nw = size(combres,1)
    mat = tosort ? sortrows(combres, by=x->x[2], rev=true) : combres
    bar(mat[:,5], fillrange=mat[:,6], c=:green, l=nothing, lab="bstrap quantiles",
        xlab="# wave", ylab="Wave Intensity, SDM[#,#]")
    tosort && plot!(xaxis=nothing)
    toannotate && annotate!([(i+0.4,mat[i,6]+50,text(Int64(mat[i,1]),4)) for i in 1:Nw])
    scatter!(mat[:,3], yerr=mat[:,4], xerr=[0.5 for i in 1:size(mat,1)],
        m=(1.5,:orange, :d, stroke(0.0)),lab="")
    scatter!(mat[:,3],m=(1.5,:orange, :d, stroke(0.0)), lab="official fit")
    scatter!(mat[:,2], m=(1.5, :red, :d, stroke(0.0)), lab="main fit")
    #scatter!(mat[:,3],m=(1.5,:orange, :d, stroke(0.0)), lab="official fit")
end

function format_wavename(name)
    name == "flat" && return "\$\\mathrm{FLAT}\$"
    tmps = "\$$(name[4])^{$(name[5:6])}$(name[8])^{$(name[9])}\\,$(name[10:end-1])\\,$(name[end])\$"
    for (k,v) in Dict("rho" => "\\rho",
                       "f2" => "f_2",
                       "pi" => "\\pi",
                       "(\\pi\\pi)_S" => "\\sigma")
        tmps = replace(tmps,k,v)
    end
    tmps
end

function plotPolarSDM(SDM; cutoff_scale=0.4)
    plot(proj=:polar)
    arg(z) = atan2(imag(z), real(z))
    SDM_with_sign = sqrt.(diag(SDM)).*cis.([arg(SDM[2,i]) for i in 1:size(SDM,1)]);
    cutoff = cutoff_scale*max(abs.(SDM_with_sign)...)
    for (i,v) in enumerate(SDM_with_sign)
        plot!([0, arg(v)], [0, abs(v)], lab="", l=(1))
        (abs(v) > cutoff) && annotate!(
            [(arg(v), 1.1*abs(v), text(format_wavename(wavenames[i]),10))])
    end
    plot!()
end

end


text:
module TestInfo

using BenchmarkConfigSweeps
using BenchmarkConfigSweeps.Internal: asconfig, asinfo, julia_info
using Test

function test_external_julia_info()
    julia = BenchmarkConfigSweeps.julia(Base.julia_cmd())
    @test julia_info(julia) == julia_info()
end

function test_asinfo()
    jl = Base.julia_cmd().exec[1]
    specs = [
        BenchmarkConfigSweeps.julia(`$jl --startup-file=no`),
        BenchmarkConfigSweeps.julia(`$jl --startup-file=no --check-bounds=yes`),
        BenchmarkConfigSweeps.julia(`$jl --startup-file=no`),
        BenchmarkConfigSweeps.julia(`$jl --startup-file=no --check-bounds=yes`),
    ]
    configs = asconfig.(specs)
    info = asinfo(configs)
    @test length(info.julia) == 2
end

end  # module


text:
export iplpAffineScaling, iplpPathFollowing , iplpwrapper

"""
function F(x,y,s,Ae,ae,c;tau)

Set up KKT conditions for solving

min c'x
s.t. Ax = b, x >= 0
(primal problem)

or equivalently

max b'y
s.t. A'y + s = c, s>= 0
(dual problem)

F(x,y,s) = | Ae*x - ae    | = |0|
           | Ae'y + s - c |   |0|
           | x.*s - tau   |   |0|

Consult Nocedal & Wright "Numerical Optimization" Chapter 14.1 for more details

    Input:
        1. x  - starting point, nx1
        2. y  - Lagrange multipliers, nx1
        3. s  - slack variables, nx1
        4. Ae - matrix of weights encoding equality constraints,mxn
        5. ae - equality constraints to satisfy,nx1
        6. c  - starting vector, nx1

    Keyword arguments:
        1. tau - non-negativity constraint such that (x,s) >= 0

    Output:
        1. r1 - Ae*x - ae
        2. r2 - Ae'y + s - c
        3. r3 - XS - tau
"""
function F(x,y,s,Ae,ae,c;tau=0.0)
    r1 = Ae*x - ae
    r2 = Ae'*y + s - c
    r3 = x.*s - tau

    return (r1,r2,r3)
end

"""
Split output into components correspond to x variables, lagrange multipliers, slack variables
"""
function splitVec(d,nx,ny,ns)
    return d[1:nx], d[nx+1:nx+ny], d[nx+ny+1:end]
end

"""
function iplpAffineScaling(A,b,c,x,y,s;kwargs...)

Wrapper function for iplp()
"""

iplpAffineScaling(A,b,c,x,y,s;kwargs...) = iplp(A,b,c,x,y,s;kwargs...)

"""
function iplpPathFollowing(A,b,c,x,y,s;kwargs...)

Wrapper function for iplp()
"""

iplpPathFollowing(A,b,c,x,y,s;kwargs...) = iplp(A,b,c,x,y,s;sig=iter->max(min(0.5^iter,0.75),1e-2), kwargs...)

function stepLengthIPLP(x,s,dx,ds;eta=0.999)
    alphax = -1./min(minimum(dx./x),-1)
    alphax = min(1,eta*alphax)

    alphas = -1./min(minimum(ds./s),-1)
    alphas = min(1,eta*alphas)

    return min(alphax,alphas), min(alphax,alphas)
end

function solveBigSystem(x,y,s,A,b,c,r1,r2,r3)
    (m,n) = size(A)
    KKT = [A zeros(m,m) zeros(m,n); zeros(n,n) A' eye(n); diagm(s) zeros(n,m) diagm(x)]
    step = KKT\[r1;r2;r3]

    return splitVec(step,n,m,n)
end

function solveNormEq(x,y,s,A,b,c,r1,r2,r3)
    NormEq = A*sdiag(x./s)*A'
    b      = r1 + A*(x.*(r2-r3./x)./s)

    dy3 =  NormEq\b
    dx3 =  - x.*(r2-r3./x-A'*dy3)./s
    ds3 = - s.*(dx3./x) + r3./x
   return [dx3;dy3;ds3]
end

function solveSymInd(x,y,s,A,b,c,r1,r2,r3)
    m,n = size(A)
    SymInd = [-sdiag(s./x) A'; A zeros(m,m)]
    b = [r2-r3./x;r1]

    xt =  SymInd\b
    dx2 = xt[1:n]
    dy2 = xt[n+1:n+m]
    ds2 = - s.*(dx2./x) + r3./x
    return [dx2;dy2;ds2]
end

"""
function iplp(A,b,c,x,y,s;maxIter,out,eta,ftol,mutol,storeInterm,sig,Jinv)

Interior Point Linear Programming (IPLP) solver

    Input:
        1. A - matrix of weights for encoding equality constraints, mxn
        2. b - equality constraints to satisfy, nx1
        3. c - starting vector, nx1
        4. x - starting point, nx1
        5. y - Lagrange multipliers, nx1
        6. s - slack variables, nx1

    Keyword arguments:
        7. maxIter      - number of iterations for Interior Point method
        8. out          - print intermediate outputs for Interior Point method for linear programming if out = 2, final iteration results if out >= 1 or = 0
        9. eta          - step size for updating search direction
        10. ftol        - tolerance for breaking out of iterations if change in sum of search directions (outputs of F()) < ftol AND duality measure < mutol
        11. mutol       - tolerance for breaking out of iterations if change in sum of norm of search directions (outputs of F()) < ftol AND duality measure < mutol
        12. storeInterm - return intermediate history
        13. sig         - function handle for setting centering parameter; if sigma = 0 then all steps are standard Newton steps
        14. Jinv        - function handle for solving

                          |0 A' I| |dx| = |- (A'y + s - c)  |
                          |A 0  0| |dy| = |- (Ax - b)       |
                          |S 0  X| |ds| = |-XSe + sigma*mu*e|

                          S = diag(s), X = diag(x), sigma = centering parameter in (0,1), mu = centrality measure = (x.*s)/n (n = length(s) = length(x)), e = unit vector of ones

    Output:
        1. x   - solution vector
        2. y   - Lagrange multipliers
        3. s   - slack variables
        4. his - history of iterations storing sum of norm of search directions, duality measure, etc

"""
function iplp(A,b,c,x::Vector,y::Vector,s::Vector; maxIter::Int=10,out::Int=1,eta::Real=0.999,
                    ftol::Real=1e-4,mutol::Real=1e-4,storeInterm=0,sig::Function=iter->0.0,Jinv::Function=solveBigSystem)

    nx = length(x)
    ny = length(y)
    ns = length(s)

    (f1,f2,f3) = F(x,y,s,A,b,c,tau=sig(1)*dot(x,s)/nx)

    his = zeros(maxIter,8)

    if storeInterm == 1
        X = zeros(nx,maxIter)
        Y = zeros(ny,maxIter)
        S = zeros(ns,maxIter)
    end

    if out==2
        @printf "=== Interior Point for LP ===\n"
        @printf "%4s\t%3s\t\t%2s\t\t%4s\t\t%4s\t\t%4s\t\t%5s\t\t%3s\n" "iter" "|F|" "mu" "|dx|" "|dy|" "|ds|" "alpha" "sig"
    end
    flag = -1
    iter = 1
    for iter=1:maxIter
        (dx,dy,ds) = Jinv(x,y,s,A,b,c,-f1,-f2,-f3)
        aMaxx,aMaxs = stepLengthIPLP(x,s,dx,ds,eta=eta)

        # update
        x += aMaxx*dx
        y += aMaxs*dy
        s += aMaxs*ds

        if storeInterm == 1
            X[:,iter] = x
            Y[:,iter] = y
            S[:,iter] = s
        end

        (f1,f2,f3) = F(x,y,s,A,b,c,tau=sig(iter)*dot(x,s)/nx)
        # his[iter,:] = [norm(f1)+norm(f2)+norm(f3) dot(s,x)/nx norm(dx) norm(dy) norm(ds) aMaxx aMaxs sig(iter)]
        his[iter,1:7] = [norm(f1)+norm(f2)+norm(f3) dot(s,x)/nx norm(dx) norm(dy) norm(ds) aMaxx aMaxs ]

        if out==2
            @printf "%4d\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\t%1.2e\n" iter his[iter,1] his[iter,2] his[iter,3] his[iter,4] his[iter,5]  his[iter,6] his[iter,8]
        end
        if his[iter,1]<ftol && (his[iter,2]<mutol)
            flag = 0
            break
        end
    end

    if out>=0
        if flag==-1
            @printf "iplp iterated maxIter(=%d) times but reached only |F(x,y,s)|=%1.2e (ftol=%1.2e) and mu=%1.2e (mutol=%1.2e)\n" maxIter his[end,1] ftol his[end,2] mutol
        elseif flag==0 && out>=1
            @printf "iplp achieved ftol=%1.2e and mutol=%1.2e at iteration %d. Returned result satisfies |F(x,y,s)|=%1.2e and mu=%1.2e.\n" ftol  mutol iter his[iter,1] his[iter,2]
        end
    end

    (x,y,s) = (storeInterm==1) ? (X[:,1:iter],Y[:,1:iter],S[:,1:iter]) : (x,y,s)
    return x,y,s,his[1:iter,:]
end


text:
using MacroTools: postwalk, @capture

export @dice, compile, rundice

struct DiceCode
    interpret::Function
end

"""
process dice code before running it
currently it makes if-then-else, &&, and || polymorphic
"""
function to_dice(code)

    # manual hygiene
    mgr = gensym(:mgr)
    flip = gensym(:flip)
    ite = gensym(:ite)

    transformed_code = postwalk(esc(code)) do x
        # TODO: replace by custom desugaring that is still lazy for boolean guards
        # this will not work in general
        # for example when control flow has `return`` inside
        if x isa Expr && (x.head == :if || x.head == :elseif)
            @assert length(x.args) == 3 "@dice macro only supports purely functional if-then-else"
            return :($ite($(x.args...)))
        end
        @capture(x, flip(P_)) && return :($flip($P))
        @capture(x, A_ || B_) && return :($A | $B)
        @capture(x, A_ && B_) && return :($A & $B)
        return x
    end

    return quote

        DiceCode( $(esc(mgr)) -> begin

            $(esc(flip))(prob::Number) =
                DistBool($(esc(mgr)), prob)

            $(esc(ite))(args...) =
                ifelse(args...)

            # transformed user code
            $transformed_code
        end)
    end
end

macro dice(code)
    to_dice(code)
end

function compile(code::DiceCode)
    code.interpret(CuddMgr())
end

function to_dice_ir(code::DiceCode)
    interpretation = code.interpret(IrMgr())
    treeify(interpretation.bit)
end

function rundice(code::DiceCode)
    rundice(to_dice_ir(code))
end

function infer(code::DiceCode, algo)
    x = if algo == :ocaml
        to_dice_ir(code)
    elseif algo == :bdd
        compile(code)
    end
    infer(x)
end


text:
using GAP, Hecke, Printf, ArgParse

include(joinpath(Hecke.pkgdir,"examples/FieldEnumeration.jl"))

###############################################################################
#
#  Read-Write
#
###############################################################################

function _write_fields(list::Array{Tuple{AnticNumberField, fmpz},1}, filename::String)
  f=open(filename, "a")
  for L in list
    x=([coeff(L[1].pol, i) for i=0:degree(L[1].pol)], L[2])
    Base.write(f, "$x\n")
  end
  close(f)
end

function _read_fields(filename::String)
  f=open(filename, "r")
  Qx,x=PolynomialRing(FlintQQ,"x")
  pols=Tuple{fmpq_poly, fmpz}[]
  for s in eachline(f)
    a=eval(parse(s))
    push!(pols,(Qx(a[1]), a[2]))
  end
  close(f)
  return pols
end

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--order", "-o"
      help = "Order of the group"
      arg_type = Int
      default = -1
    "--id", "-i"
      help = "Id of the group"
      arg_type = Int
      default = -1
    "--number", "-n"
      help = "Number of the group"
      arg_type = Int
      default = -1
    "--tamely-ramified"
      help = "Only tamely ramified fields"
      action = :store_true
    "--only-real"
      help = "Only totally real fields"
      action = :store_true
    "--disc-bound"
      help = "Discriminant bound"
      arg_type = fmpz
      default = fmpz(-1)
  end

  return parse_args(s)
end

function main()
  parsed_args = parse_commandline()

  local grp_order::Int
  local grp_id::Int
  local dbound::fmpz
  local only_real::Bool
  local only_tame::Bool
  local grp_no::Int

  for (arg, val) in parsed_args
    println("$arg => $val")
    if arg == "order"
      grp_order = val
    elseif arg == "id"
      grp_id = val
    elseif arg == "disc-bound"
      dbound = val
    elseif arg == "only-real"
      only_real = val
    elseif arg == "tamely-ramified"
      only_tame = val
    elseif arg == "number"
      grp_no = val
    end
  end

  sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

  if grp_id != -1
    @assert grp_order != -1
    n = grp_order
    i = grp_id
    grp_no = findfirst(isequal((n, i)), small_solvable_groups)
  else
    @assert grp_no != -1
    n, i = small_solvable_groups[grp_no]
  end

  file = sprint_formatted("%0$(length(string(length(small_solvable_groups))))d", grp_no) * "-$n-$i."* "log"

  @show grp_order
  @show grp_id
  @show grp_no
  @show dbound
  @show only_real
  @show only_tame
  @show file

  if isfile(file)
    throw(error("File $file does already exist"))
  end

  println("========================================")
  println("Group: $grp_no: $n $i")
  println("========================================")

  if dbound == -1
    dbound = Hecke.lower_discriminant_bound(n, n)
  end

  dbound = 10 * dbound

  println("========================================")
  println("Discriminant bound: $dbound")
  println("========================================")

  set_verbose_level(:FieldsNonFancy, 1)

  l = []

  complex_fields = []
  real_fields = []
  all = []

  found_complex = false
  found_real = false

  if mod(n, 2) != 0
    found_complex = true
  end

  while true

    @show only_real

    l = fields(n, i, dbound, only_real = only_real)

    for v in l
      x = v.field
      if !found_complex
        if signature(x)[2] != 0
          push!(complex_fields, v)
        end
      end

      if !found_real
        if signature(x)[2] == 0
          push!(real_fields, v)
        end
      end
    end

    if !found_complex && length(complex_fields) > 0
      append!(all, complex_fields)
      found_complex = true
      println("== FOUND SOMETHING COMPLEX ==")
      empty!(complex_fields)
    end

    if !found_real && length(real_fields) > 0
      append!(all, real_fields)
      found_real = true
      println("== FOUND SOMETHING REAL ==")
      empty!(real_fields)
    end

    if found_complex
      if !found_real
        only_real = true
        println("== I AM IN ONLY_REAL = TRUE MODE ==")
      else
        println("== EXITING ==")
        break
      end
    end

    dbound = 100 * dbound
    println("========================================")
    println("Increasing discriminant bound: $dbound")
    println("========================================")
    println("only_real = $only_real")
    println("found_complex = $found_complex")
    println("found_real = $found_real")
  end

  l = all

  ll = map(v -> v.field, l)
  ffields = [ (x, discriminant(maximal_order(x))) for x in ll ]
  sort!(ffields, lt = (x, y) -> abs(x[2]) <= abs(y[2]))

  @show length(ffields)

  _write_fields(ffields, file)
end

main()


text:
using LinearAlgebra:
    Symmetric,
    Hermitian,
    UpperTriangular,
    UnitUpperTriangular,
    LowerTriangular,
    UnitLowerTriangular,
    UpperHessenberg,
    Tridiagonal,
    SymTridiagonal,
    Bidiagonal,
    Diagonal,
    UniformScaling

## Types for which a `generate` method is implemented. For parametric
## types, the parameters are arbitrary.
types_with_generate = [Bool,
                       UInt8, UInt16, UInt32, UInt64, UInt128,
                       Int128, Int16, Int32, Int64, Int8,
                       Float16, Float32, Float64,
                       Complex{Float64},
                       String, Char,
                       Vector{Float64}, Matrix{Int}, Array{Float32, 3},
                       BitVector, BitMatrix, BitArray{5},
                       SquareMatrix{Float64},
                       Symmetric{Float64, SquareMatrix{Float64}},
                       Hermitian{Float64, SquareMatrix{Float64}},
                       UpperTriangular{Float64},
                       UpperTriangular{Float64, SquareMatrix{Float64}},
                       UnitUpperTriangular{Float64},
                       UnitUpperTriangular{Float64, SquareMatrix{Float64}},
                       LowerTriangular{Float64},
                       LowerTriangular{Float64, SquareMatrix{Float64}},
                       UnitLowerTriangular{Float64},
                       UnitLowerTriangular{Float64, SquareMatrix{Float64}},
                       UpperHessenberg{Float64},
                       UpperHessenberg{Float64, Matrix{Float64}},
                       Tridiagonal{Float64},
                       Tridiagonal{Float64, Vector{Float64}},
                       SymTridiagonal{Float64},
                       SymTridiagonal{Float64, Vector{Float64}},
                       Bidiagonal{Float64},
                       Bidiagonal{Float64, Vector{Float64}},
                       Diagonal{Float64},
                       Diagonal{Float64, Vector{Float64}},
                       UniformScaling{Float64},
                       Union{Complex{Float64}, String, Char,
                             Vector{Float64}, Matrix{Int},
                             Array{Float32, 3}, BitVector, BitMatrix,
                             BitArray{5}, SquareMatrix{Float64},
                             Symmetric{Float64, SquareMatrix{Float64}}}]


text:
module LogClustering

    include("KATE.jl")

end # module


text:
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# minellipse.jl
#
# This example is from the Boyd & Vandenberghe book "Convex Optimization".
# If PyPlot is installed, passing the argument `plot` will plot
# the solution, e.g. julia minellipse.jl plot
#
# Given a set of ellipses centered on the origin
# E(A) = { u | u^T inv(A) u <= 1 }
# find a "minimal" ellipse that contains the provided ellipses
#
# We can formulate this as an SDP:
#     minimize  trace(WX)
#   subject to  X >= A_i,    i = 1,...,m
#               X PSD
# where W is a PD matrix of weights to choose between different solutions
#############################################################################

using JuMP
using Mosek
using PyPlot  # Comment out if not installed

# We will use three ellipses
m = 3
# Two "simple" ones
As = Any[ [2.0  0.0;
           0.0  1.0],
          [1.0  0.0;
           0.0  3.0]]
# and a random one
randA = rand(2,2)
push!(As, (randA' * randA) * (rand()*2+1))

# We change the weights to see different solutions, if they exist
W = [1.0 0.0
     0.0 1.0];

mod = Model()
@variable(mod, X[1:2,1:2], SDP)
@objective(mod, Min, trace(W*X))
for i = 1:m
    @SDconstraint(mod, X >= As[i])
end
solve(mod)

X_val = getvalue(X)
println(X_val)

# Plot it, if desired (e.g. julia minellipse.jl plot)
if length(ARGS) > 0 && Pkg.installed("PyPlot") !== nothing
    # Setup the figure
    fig = figure()
    ax = fig[:gca]()
    ax[:set_xticks]([-4:+4])
    ax[:set_yticks]([-4:+4])

    # Draw provided ellipses
    for i = 1:m
        xs = Float64[]
        ys = Float64[]
        for angle in linspace(0, 2*pi, 100)
            u = [cos(angle),sin(angle)]
            x = As[i] * u
            push!(xs, x[1])
            push!(ys, x[2])
        end
        plot(xs, ys, "b", linewidth=2.0)
    end

    # Draw bounding ellipse
    xs = Float64[]
    ys = Float64[]
    for angle in linspace(0, 2*pi, 100)
        u = [cos(angle),sin(angle)]
        x = X_val * u
        push!(xs, x[1])
        push!(ys, x[2])
    end
    plot(xs, ys, "r", linewidth=2.0)

    fig[:canvas][:draw]()
    readline()
end


text:
"""
    calc_sun_position_hor(datetime, lat, long)

Compute the Sun position at given time and observer coordinates in horizontal coordinates.

# Arguments:
- `datetime`: time: Either a `ZonedDateTime`, or `DateTime` assumed in UTC
- `lat`, `long`: latitude and longitude in degree

# Value
`NamedTuple`: sun position with entries of the same type as `lat` and `long`
- `altitude`: angle above the horizon [rad].
- `azimuth`: angle ange the horizon plane eastwards of north [rad]
- `hourangle`: [rad] as output by AstroLib.eq2hor
   Seems to represent time [day/2pi] after solar noon.
   Value at local timezone noon provdes (local time - solar time).
"""
function calc_sun_position_hor(datetime::ZonedDateTime, lat, long)
  datetimeUTC = DateTime(datetime,UTC)
  calc_sun_position_hor(datetimeUTC, lat, long)
end
function calc_sun_position_hor(datetime::DateTime, lat::FT, long::FT) where FT
    deg2rad = FT(π/180)
    jd = FT(datetime2julian(datetime))
    pos_eq = calc_sun_position_MOD(jd)
    # precession is already account for in MOD
    pos_hor = eq2hor(pos_eq.α/deg2rad, pos_eq.δ/deg2rad, jd, lat, long; precession = false) .* deg2rad
    (altitude = FT(pos_hor[1]), azimuth = FT(pos_hor[2]), hourangle = FT(pos_hor[3]))
end

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the sun position.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] The Astronomical Almanac for the year 2000 (p. C24).
#
#   [2] http://aa.usno.navy.mil/faq/docs/SunApprox.php
#
#   [3] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorne, CA.
#
#   [4] The Astronomical Almanac for the year 2006.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# modified from https://github.com/JuliaSpace/SatelliteToolbox.jl/blob/c5244c50e55e85f6b3a1702de16f00a9c726ee24/src/sun/sun_position.jl#L24-L30
#
# https://en.wikipedia.org/wiki/Ecliptic_coordinate_system

"""
    calc_sun_position_MOD(JD::Number)

Compute the Sun position at the Julian Day `JD`.

Results are represented in the Mean Equinox of Date (MOD),
i.e. accounting for precession but not for nutation and smaller pertubation
of the polar axes, in spherical ecliptic and equatorial coordinates.
The algorithm was adapted from [Vallado 2013, p. 277-279].

# Arguments:
- `JD`: time given as Julian Day .

# Value
`NamedTuple`: sun position where
- Ecliptic coordinates (1:3)
  - λ: Ecliptic longitude of the Sun [rad].
  - β: Ecliptic latitude of the Sun [rad] is assumed 0.
  - r: Distance of the Sun from Earth [m].
- Equatorial coordinate (4:5)
  - α: ascention [rad]
  - δ: declination [rad]
- ϵ: Obliquity of the ecliptic [rad].

# References
- Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
       Microcosm Press, Hawthorne, CA.
"""
function calc_sun_position_MOD(JD::Number)
    # Constants (apdated from SatelliteToolbox)
    deg2rad = π/180
    au2m = 149597870700.0
    JD_J2000 = 2.451545e6

    # Number of Julian centuries from J2000 epoch.
    T_UT1 = (JD-JD_J2000)/36525.0

    # Here, we will assume that T_TBD ≈ T_UT1. This can be done because we are
    # using a low-precision algorithm [3].
    T_TBD = T_UT1

    # Mean anomaly of the Sun [deg].
    Ms = 357.529_109_2 + 35_999.050_34T_TBD

    # Convert Ms to [rad] and limit to the interval [0,2π].
    Ms = mod2pi(Ms*deg2rad)

    # Compute auxiliary variables.
    sinMs,  cosMs  = sincos(Ms)
    sin2Ms, cos2Ms = sincos(2Ms)

    # Mean longitude of the Sun [deg].
    λ_m = 280.460 + 36_000.771T_UT1

    # Ecliptic latitude of the Sun [deg]. # twutz 2110 longitude?
    λ_e = λ_m + 1.914_666_471sinMs + 0.019_994_643sin2Ms

    # Obliquity of the ecliptic [deg].
    ϵ = 23.439_291 - 0.013_004_2T_TBD

    # Convert λ_e and ϵ to [rad] and limit to the interval [0,2π].
    λ_e = mod2pi(λ_e*deg2rad)
    ϵ   = mod2pi(  ϵ*deg2rad)

    # Auxiliary variables.
    sinϵ   , cosϵ   = sincos(ϵ)
    sinλ_e , cosλ_e = sincos(λ_e)

    # Distance of the Sun from Earth [m].
    r = (1.000_140_612 - 0.016_708_617cosMs - 0.000_139_589cos2Ms )*au2m

    # ascension
    α = atan(cosϵ * sinλ_e, cosλ_e)

    # declination
    δ = asin(sinϵ * sinλ_e)

    S_MOD_rad = (
        λ = λ_e, β = 0.0, r = r,
        α = α, δ = δ,
        ϵ = ϵ
        )
end


text:
@testset "Frobenius extensions" begin
  # S3 on 3 points
  # Data taken from Cohen & Diaz y Diaz & Olivier, Counting discriminants of number fields
  flds = Hecke.primitive_frobenius_extensions(QQ, (6, 1), fmpz(10)^2)
  @test length(flds) == 7
  flds = Hecke.primitive_frobenius_extensions(QQ, (6, 1), fmpz(10)^3, only_real = true)
  @test length(flds) == 22
  @test all(count(L -> isisomorphic(L, K)[1], flds) == 1 for K in flds)
  flds = Hecke.primitive_frobenius_extensions(QQ, (6, 1), fmpz(10)^3, only_non_real = true)
  @test length(flds) == 127
  flds = Hecke.primitive_frobenius_extensions(QQ, (6, 1), fmpz(10)^3)
  @test length(flds) == 149

  QQx, x = QQ["x"]

  # D5 on 5 points, according to Malle-Klüners
  flds = Hecke.primitive_frobenius_extensions(QQ, (10, 1), fmpz(2208), only_non_real = true)
  @test length(flds) == 0
  flds = Hecke.primitive_frobenius_extensions(QQ, (10, 1), fmpz(2209), only_non_real = true)
  @test length(flds) == 1
  @test isisomorphic(flds[1], number_field(x^5 - 2*x^4 + 2*x^3 - x^2 + 1)[1])[1]

  flds = Hecke.primitive_frobenius_extensions(QQ, (10, 1), fmpz(106800), only_real = true)
  @test length(flds) == 0
  flds = Hecke.primitive_frobenius_extensions(QQ, (10, 1), fmpz(160801), only_real = true)
  @test length(flds) == 1
  @test isisomorphic(flds[1], number_field(x^5 - x^4 - 5*x^3 + 4*x^2 + 3*x - 1)[1])[1]

  @test_throws ArgumentError Hecke.primitive_frobenius_extensions(QQ, (1, 1), fmpz(10))
  @test_throws ArgumentError Hecke.primitive_frobenius_extensions(QQ, (10, 1), fmpz(160801), only_real = true, only_non_real = true)
end


text:
#=
Alignment methods. Each method needs
1) Type with fields necessary for algorithm
2) function "align" to operate on sort with type field defined above
3) any other necessary functions for alignment algorithm
=#

export AlignMax, AlignMin, AlignMinMax, AlignCOM, AlignProm

#=
Maximum signal
=#
mutable struct AlignMax <: Align
end

function align(a::AlignMax, sort::Sorting)
    mymax=sort.p_temp[sort.s.win2]
    sort.cent=sort.s.win2
    @inbounds for i=(sort.s.win2+1):sort.s.s_end
        if sort.p_temp[i]>mymax
            sort.cent=i+1
            mymax=sort.p_temp[i]
        end
    end
    nothing
end

mysize(align::AlignMax,win)=win

#=
Minimum signal
=#
mutable struct AlignMin <: Align
    shift::Int64
end

AlignMin()=AlignMin(10)

function align(a::AlignMin,sort::Sorting)
    mymin=sort.p_temp[sort.s.win2]
    sort.cent=sort.s.win2
    @inbounds for i=(sort.s.win2+1):sort.s.s_end
        if sort.p_temp[i]<mymin
            sort.cent=i+1
            mymin=sort.p_temp[i]
        end
    end
    sort.cent += a.shift
    nothing
end

mysize(align::AlignMin,win)=win

#=
Middle of Minimum and Max
=#
mutable struct AlignMinMax <: Align
    shift::Int64
end

AlignMinMax()=AlignMinMax(5)

function align(a::AlignMinMax,sort::Sorting)
    mymin=sort.p_temp[sort.s.win2]
    mymax=sort.p_temp[sort.s.win2]
    minind=sort.s.win2
    maxind=sort.s.win2
    @inbounds for i=(sort.s.win2+1):sort.s.s_end
        if sort.p_temp[i]<mymin
            minind=i+1
            mymin=sort.p_temp[i]
        end
        if sort.p_temp[i]>mymax
            maxind=i+1
            mymax=sort.p_temp[i]
        end
    end
    sort.cent = a.shift+round(Int64,(minind+maxind)/2)
    nothing
end

mysize(align::AlignMinMax,win)=win

#=
Center of Mass Alignment
=#

mutable struct AlignCOM <: Align
    shift::Int64
end

AlignCOM()=AlignCOM(5)

function align(a::AlignCOM,sort::Sorting)

    com=0.0
    mysum=0.0

    for i=(sort.s.win2+1):sort.s.s_end
        com += i*abs(sort.p_temp[i])
        mysum += abs(sort.p_temp[i])
    end

    sort.cent = a.shift + round(Int64, com/mysum)

end

mysize(align::AlignCOM,win)=win

#=
Prominence Alignment
=#

mutable struct AlignProm <: Align
    shift::Int64
end

AlignProm()=AlignProm(5)

function align(a::AlignProm,sort::Sorting)

    indprom=sort.s.win2+1
    prom=sort.p_temp[indprom]

    for i=(sort.s.win2+1):sort.s.s_end

        test = (sort.p_temp[i-1]-sort.p_temp[i-2])
        test += (sort.p_temp[i]-sort.p_temp[i-1])
        test += (sort.p_temp[i]-sort.p_temp[i+1])
        test += (sort.p_temp[i+1]-sort.p_temp[i+2])

        if abs(test) > prom
            prom=abs(test)
            indprom=i
        end
    end
    sort.cent = a.shift + indprom
end

mysize(align::AlignProm,win)=win

#=
Maximum Magnitude
=#

#=
FFT upsampling
=#

#=
type AlignFFT <: Align
    M::Int64
    x_int::Array{Complex{Float64},1}
    fout::Array{Complex{Float64},1}
    upsamp::Array{Float64,1}
    align_range::UnitRange{Int64}
end

function AlignFFT(M::Int64)

    AlignFFT(M,zeros(Complex{Float64},M*2*window),
             zeros(Complex{Float64},2*window),zeros(Float64,M*2*window),
             (window_half*M+1):((window+window_half)*M))

end

function align(a::AlignFFT, sort::Sorting)

    sort.a.fout[:]=fft(sort.p_temp)

    sort.a.x_int[1:window]=sort.a.fout[1:window]
    sort.a.x_int[window+1]=sort.a.fout[window+1]/2
    sort.a.x_int[(window+2):(sort.a.M*2*window-window)]=zeros(Complex{Float64},2*sort.a.M*window-2*window-1)
    sort.a.x_int[sort.a.M*2*window-window+1]=sort.a.fout[window+1]/2
    sort.a.x_int[(sort.a.M*2*window-window+2):end]=sort.a.fout[(window+2):end]

    ifft!(sort.a.x_int)
    sort.a.upsamp[:]=sort.a.M.*real(sort.a.x_int)

    j=indmax(sort.a.upsamp[sort.a.align_range])+sort.a.M*window_half
    sort.waveform=view(sort.a.upsamp,j-sort.a.M*window_half:j+sort.a.M*window_half-1)

    return j-window_half:j+window_half-1
end

mysize(align::AlignFFT)=window*align.M
=#
#=
FFT upsampling + temporal order of peaks

Rutishauser 2006
=#
#=
type AlignOsort <: Align

end

function align(a::AlignOsort, sort::Sorting)

end

function mysize(align::AlignOsort)
    window*align.M
end
=#


text:
@enum coordstype auto griddata pointdata

# Insertion point for defining my own functions. They are mainly abstract
# over the full type and it's view, so another possibility would be to
# implement with unions
abstract type SESpatialData{T<:EcoBase.AbstractLocationData} <: EcoBase.AbstractPlaces{T} end
abstract type SELocations{T<:EcoBase.AbstractLocationData} <: EcoBase.AbstractPlaces{T} end
abstract type SEThings{D <: Real} <: EcoBase.AbstractThings end
abstract type SEGrid <: EcoBase.AbstractGrid end
abstract type SEPoints <: EcoBase.AbstractPoints end

abstract type AbstractComMatrix{D<:Real} end

# I could implement sitestats as a Dict with several DataFrames to make space for big data sets, but I prefer to not do this now. Example below.

# I could do a lot more with immutable types if I had a clearer view/copy implementation
mutable struct GridTopology <: EcoBase.AbstractGrid
    xs::StepRangeLen{Float64}
    ys::StepRangeLen{Float64}
end

mutable struct Bbox
    xmin::Number
    xmax::Number
    ymin::Number
    ymax::Number
end

# Do I need sitenames here? I think so, they should match those in sitestats, and be separate
mutable struct PointData <: SEPoints
    coords::Matrix{Float64}
end

mutable struct GridData <: SEGrid
    indices::Matrix{Int}
    grid::GridTopology
end

mutable struct Locations{T<:Union{GridData, PointData}} <: SELocations{T}
    coords::T
    sitestats::DataFrames.DataFrame
    function Locations{T}(coords, sitestats = DataFrames.DataFrame(sites = string.(1:nsites(coords)))) where T
        DataFrames.nrow(sitestats) == nsites(coords) || throw(DimensionMismatch("Wrong number of rows in sitestat")) # a little check for the right number
        new(coords, sitestats)
    end
end

Locations(coords, sitestats) = Locations{typeof(coords)}(coords, sitestats)

#it's a question how many of these structs need to be mutable, as opposed to when you want to allocate a new object
mutable struct ComMatrix{D} <: AbstractComMatrix{D}
    occurrences::SparseMatrixCSC{D}
    speciesnames::Vector{String}
    sitenames::Vector{String}
    ComMatrix{D}(occ::SparseMatrixCSC{D}, spn::Vector{String}, sin::Vector{String}) where {D} = new(dropzeros!(occ), spn, sin)
end

# likewise, do I need a speciesnames here? Should traits have a :series field (like now) or all matching be done on the speciesnames?
mutable struct SpeciesData{D} <: SEThings{D}
    commatrix::ComMatrix{D}
    traits::DataFrames.DataFrame
    function SpeciesData{D}(commatrix::ComMatrix{D}, traits::DataFrames.DataFrame) where D <: Real
        DataFrames.nrow(traits) ==  nspecies(commatrix) || throw(DimensionMismatch("Wrong number of species in traits"))
        new(commatrix, traits)
    end
end
SpeciesData(commatrix::ComMatrix{D}, traits) where D<:Real = SpeciesData{D}(commatrix, traits)

# Not really sure what this type is for
mutable struct SiteData{T<:Union{GridData, PointData}} <: SESpatialData{T}
    site::T
end

abstract type SEAssemblage{D<:Real, T<:SEThings, P<:SELocations} <: EcoBase.AbstractAssemblage{D, T, P} end

mutable struct Assemblage{D<:Real, P<:Locations} <: SEAssemblage{D, SpeciesData{D}, P} # A type to keep subtypes together, ensuring that they are all aligned at all times
    site::P
    occ::SpeciesData{D}

    # inner constructor
    function Assemblage{D, P}(site::P, occ::SpeciesData{D}) where {P <: SELocations, D <: Real}
        size(occurrences(occ), 2) == size(coordinates(site), 1) || error("Length mismatch between occurrence matrix and coordinates")
        #TODO activate this # sitenames(occ) == sitenames(site) || error("sitenames do not match") #I need a constructor that matches them up actively
        new(site, occ)
    end
end


text:
#    Problem : GROUP A
#    *********
#    The extended Dennis and Schnabel problem, as defined by Li.
#
#    Original SIF Source:
#    G. Li,
#    "The secant/finite difference algorithm for solving sparse
#    nonlinear systems of equations",
#    SIAM Journal on Optimization, (to appear), 1990.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    EDENSCH.SIF classification  OUR2-AN-V-0
#
#    Number of variables (at least 2)
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        item1 = x[i] - 2
        item2 = x[i]x[i+1] - 2x[i+1]
        item3 = x[i+1] + 1
        fx += 16 + item1^4 + item2^2 + item3^2
    end
    return fx
end

g! = (g, x) -> begin
    for i in firstindex(x):lastindex(x)-1
        item1 = x[i] - 2
        item2 = x[i]*x[i+1] - 2x[i+1]
        item3 = x[i+1] + 1
        g[i] += 4item1^3 + 2x[i+1]item2
        g[i+1] += 2item2*(x[i] - 2) + 2item3
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)-1
        item1 = x[i] - 2
        item2 = x[i]*x[i+1] - 2x[i+1]
        item3 = x[i+1] + 1
        fx += 16 + item1^4 + item2^2 + item3^2
        g[i] += 4item1^3 + 2x[i+1]item2
        g[i+1] += 2item2*(x[i] - 2) + 2item3
    end
    return fx, g
end

init = (n::Int=2000) -> begin

    return n, 8.0*ones(n)
end

@warn "TODO: Debug EDENSCH"
# TestSet["EDENSCH"] = UncProgram("EDENSCH", f, g!, fg!, init)


text:
##########################################################################
# Simulation results
abstract Results # super class for PrecodingResults and AssignmentResults

##########################################################################
# Simulation types
typealias SimulationParams Dict{ASCIIString, Any}
typealias PlotParams Dict{ASCIIString, Any}

type SingleSimulationResults{T <: Results}
    precoding_results::Dict{ASCIIString, T}
end
SingleSimulationResults{T <: Results}(::Type{T}) =
    SingleSimulationResults(Dict{ASCIIString, T}())
Base.getindex(s::SingleSimulationResults, k) =
    getindex(s.precoding_results, k)
Base.setindex!(s::SingleSimulationResults, v, k) =
    setindex!(s.precoding_results, v, k)

type MultipleSimulationResults{N}
    simulation_results::Array{SingleSimulationResults, N}
end
MultipleSimulationResults(dims...) =
    MultipleSimulationResults(Array(SingleSimulationResults, dims...))
Base.getindex(m::MultipleSimulationResults, inds...) =
    getindex(m.simulation_results, inds...)
Base.setindex!(m::MultipleSimulationResults, v, inds...) =
    setindex!(m.simulation_results, v, inds...)

##########################################################################
# Simulate the entire system, including assignment and precoding.
# If loop_over == :precoding_methods, we run the first assignment method
# provided and then loop over all precoding methods.
# If loop_over == :assignment_methods, we loop over all assignment methods
# and run the first provided precoding method.
# If there are several methods of the type not looped over specified,
# a warning is given.
function simulate(network, simulation_params;
    loop_over::Symbol=:precoding_methods)

    # Number of drops and small scale fading realizations
    Ndrops = get(simulation_params, "Ndrops", 1)::Int
    Nsim = get(simulation_params, "Nsim", 1)::Int

    # Main independent variable
    idp_func = simulation_params["independent_variable"][1]::Function
    idp_vals = simulation_params["independent_variable"][2]
    idp_vals_length = length(idp_vals)

    # Auxiliary independent variables
    Naux, aux_idp_funcs, aux_idp_vals, aux_idp_vals_length = get_aux_idp(simulation_params)

    # Methods that should not be swept as a function of independent and auxiliary variables
    precoding_methods_nosweep = get(simulation_params, "precoding_methods_nosweep", Function[])::Vector{Function}

    # Get assignment/precoding functions depending on looping mode
    precoding_method, assignment_method = get_other_method(simulation_params, loop_over)

    println("-- simulate on $network.")
    println("--- Ndrops: $Ndrops, Nsim: $Nsim.")
    Lumberjack.info("Starting simulation.", @compat Dict(:network => network, :simulation_params => simulation_params))

    # Set initial aux params
    set_initial_aux_params!(simulation_params, network)

    # Ensure that we are not storing all intermediate iterations.
    set_aux_precoding_param!(network, :final_iteration, "output_protocol")

    # Storage container for results
    raw_precoding_results = MultipleSimulationResults(Ndrops, Nsim, idp_vals_length, aux_idp_vals_length)
    raw_assignment_results = MultipleSimulationResults(Ndrops, 1, idp_vals_length, aux_idp_vals_length)

    # Outer simulation loop
    if loop_over == :precoding_methods
        progress = ProgressMeter.Progress(Ndrops*idp_vals_length*aux_idp_vals_length*length(simulation_params["precoding_methods"])*Nsim)
    else
        progress = ProgressMeter.Progress(Ndrops*idp_vals_length*aux_idp_vals_length*length(simulation_params["assignment_methods"])*Nsim)
    end
    for Ndrops_idx = 1:Ndrops
        draw_user_drop!(network)
        channels = [ draw_channel(network) for n = 1:Nsim ]

        # Inner simulation loop depends on loop mode
        if loop_over == :precoding_methods
            # Loop over main independent variable
            for idp_vals_idx = 1:idp_vals_length
                # Set main independent variable
                idp_func(network, idp_vals[idp_vals_idx])

                # Loop over auxiliary variables
                for aux_idp_vals_idx = 1:aux_idp_vals_length
                    # Set all auxiliary independent variables
                    for Naux_idx = 1:Naux
                        aux_idp_funcs[Naux_idx](network, aux_idp_vals[Naux_idx][aux_idp_vals_idx])
                    end

                    # Allocate memory for assignment method, and run it
                    raw_assignment_results[Ndrops_idx, 1, idp_vals_idx, aux_idp_vals_idx] = SingleSimulationResults(AssignmentResults)
                    raw_assignment_results[Ndrops_idx, 1, idp_vals_idx, aux_idp_vals_idx][stringify_method(assignment_method)] = assignment_method(channels[1], network)

                    for precoding_method in simulation_params["precoding_methods"]
                        for Nsim_idx = 1:Nsim
                            # Allocate memory if this is the first method to be run
                            if precoding_method == simulation_params["precoding_methods"][1]
                                raw_precoding_results[Ndrops_idx, Nsim_idx, idp_vals_idx, aux_idp_vals_idx] = SingleSimulationResults(PrecodingResults)
                            end
                            if in(precoding_method, precoding_methods_nosweep) && idp_vals_idx != 1 && aux_idp_vals_idx != 1
                                # We should not sweep this precoding method, so re-store old result.
                                raw_precoding_results[Ndrops_idx, Nsim_idx, idp_vals_idx, aux_idp_vals_idx][stringify_method(precoding_method)] = raw_precoding_results[Ndrops_idx, Nsim_idx, 1, 1][stringify_method(precoding_method)]
                            else
                                raw_precoding_results[Ndrops_idx, Nsim_idx, idp_vals_idx, aux_idp_vals_idx][stringify_method(precoding_method)] = precoding_method(channels[Nsim_idx], network)
                            end

                            ProgressMeter.next!(progress)
                        end
                    end
                end
            end
        elseif loop_over == :assignment_methods
            # Loop over main independent variable
            for idp_vals_idx = 1:idp_vals_length
                # Set main independent variable
                idp_func(network, idp_vals[idp_vals_idx])

                # Loop over auxiliary variables
                for aux_idp_vals_idx = 1:aux_idp_vals_length
                    # Set all auxiliary independent variables
                    for Naux_idx = 1:Naux
                        aux_idp_funcs[Naux_idx](network, aux_idp_vals[Naux_idx][aux_idp_vals_idx])
                    end

                    for assignment_method in simulation_params["assignment_methods"]
                        # Allocate memory if this is the first method to be run
                        if assignment_method == simulation_params["assignment_methods"][1]
                            raw_assignment_results[Ndrops_idx, 1, idp_vals_idx, aux_idp_vals_idx] = SingleSimulationResults(AssignmentResults)
                        end
                        raw_assignment_results[Ndrops_idx, 1, idp_vals_idx, aux_idp_vals_idx][stringify_method(assignment_method)] = assignment_method(channels[1], network)

                        for Nsim_idx = 1:Nsim
                            # Allocate memory if this is the first method to be run
                            if assignment_method == simulation_params["assignment_methods"][1]
                                raw_precoding_results[Ndrops_idx, Nsim_idx, idp_vals_idx, aux_idp_vals_idx] = SingleSimulationResults(PrecodingResults)
                            end
                            raw_precoding_results[Ndrops_idx, Nsim_idx, idp_vals_idx, aux_idp_vals_idx][stringify_method(assignment_method)] = precoding_method(channels[Nsim_idx], network)

                            ProgressMeter.next!(progress)
                        end
                    end
                end
            end
        end
    end

    return raw_precoding_results, raw_assignment_results
end

##########################################################################
# Simulate precoding methods and store their convergence behaviour.
# The loop_over parameter has the same behaviour as in the simulate function.
function simulate_precoding_convergence(network, simulation_params;
    loop_over::Symbol=:precoding_methods)

    # Number of drops and small scale fading realizations
    Ndrops = get(simulation_params, "Ndrops", 1)::Int
    Nsim = get(simulation_params, "Nsim", 1)::Int

    # Auxiliary independent variables
    Naux, aux_idp_funcs, aux_idp_vals, aux_idp_vals_length = get_aux_idp(simulation_params)

    # Get assignment/precoding functions depending on looping mode
    precoding_method, assignment_method = get_other_method(simulation_params, loop_over)

    println("-- simulate_precoding_convergence on $network.")
    println("--- Ndrops: $Ndrops, Nsim: $Nsim.")
    Lumberjack.info("Starting convergence simulation.",
        @compat Dict(:network => network, :simulation_params => simulation_params))

    # Set initial aux params
    set_initial_aux_params!(simulation_params, network)

    # We want to store all intermediate iterations.
    set_aux_precoding_param!(network, :all_iterations, "output_protocol")

    # Warn if we are not using a fixed number of iterations.
    if get_aux_precoding_param(network, "stop_crit") > 0.
        Lumberjack.warn("postprocess_convergence will not be able to run since stop_crit is non-zero, and the algorithms may therefore use different numbers of iterations.")
    end

    # Storage container for results
    raw_results = MultipleSimulationResults(Ndrops, Nsim, aux_idp_vals_length)

    # Outer simulation loop
    if loop_over == :precoding_methods
        progress = ProgressMeter.Progress(Ndrops*aux_idp_vals_length*length(simulation_params["precoding_methods"])*Nsim)
    else
        progress = ProgressMeter.Progress(Ndrops*aux_idp_vals_length*length(simulation_params["assignment_methods"])*Nsim)
    end
    for Ndrops_idx = 1:Ndrops
        draw_user_drop!(network)
        channels = [ draw_channel(network) for n = 1:Nsim ]

        # Inner simulation loop depends on loop mode
        if loop_over == :precoding_methods
            # Loop over auxiliary variables
            for aux_idp_vals_idx = 1:aux_idp_vals_length
                # Set all auxiliary independent variables
                for Naux_idx = 1:Naux
                    aux_idp_funcs[Naux_idx](network, aux_idp_vals[Naux_idx][aux_idp_vals_idx])
                end

                assignment_method(channels[1], network)

                for precoding_method in simulation_params["precoding_methods"]
                    for Nsim_idx = 1:Nsim
                        # Allocate memory if this is the first method to be run
                        if precoding_method == simulation_params["precoding_methods"][1]
                            raw_results[Ndrops_idx, Nsim_idx, aux_idp_vals_idx] = SingleSimulationResults(PrecodingResults)
                        end
                        raw_results[Ndrops_idx, Nsim_idx, aux_idp_vals_idx][stringify_method(precoding_method)] = precoding_method(channels[Nsim_idx], network)

                        ProgressMeter.next!(progress)
                    end
                end
            end
        elseif loop_over == :assignment_methods
            # Loop over auxiliary variables
            for aux_idp_vals_idx = 1:aux_idp_vals_length
                # Set all auxiliary independent variables
                for Naux_idx = 1:Naux
                    aux_idp_funcs[Naux_idx](network, aux_idp_vals[Naux_idx][aux_idp_vals_idx])
                end

                for assignment_method in simulation_params["assignment_methods"]
                    assignment_method(channels[1], network)

                    for Nsim_idx = 1:Nsim
                        # Allocate memory if this is the first method to be run
                        if assignment_method == simulation_params["assignment_methods"][1]
                            raw_results[Ndrops_idx, Nsim_idx, aux_idp_vals_idx] = SingleSimulationResults(PrecodingResults)
                        end
                        raw_results[Ndrops_idx, Nsim_idx, aux_idp_vals_idx][stringify_method(assignment_method)] = precoding_method(channels[Nsim_idx], network)

                        ProgressMeter.next!(progress)
                    end
                end
            end
        end
    end

    return raw_results
end

##########################################################################
# Simulate assignment methods and store their convergence behaviour.
function simulate_assignment_convergence(network, simulation_params)
    # Number of drops and small scale fading realizations
    Ndrops = simulation_params["Ndrops"]::Int

    # Auxiliary independent variables
    Naux, aux_idp_funcs, aux_idp_vals, aux_idp_vals_length = get_aux_idp(simulation_params)

    println("-- simulate_assignment_convergence on $network.")
    println("--- Ndrops: $Ndrops.")
    Lumberjack.info("Starting convergence simulation.",
        @compat Dict(:network => network, :simulation_params => simulation_params))

    # Set initial aux params
    set_initial_aux_params!(simulation_params, network)

    # Storage container for results
    raw_results = MultipleSimulationResults(Ndrops, 1, aux_idp_vals_length)

    # Simulation loop
    progress = ProgressMeter.Progress(Ndrops*aux_idp_vals_length*length(simulation_params["assignment_methods"]))
    for Ndrops_idx = 1:Ndrops
        draw_user_drop!(network)
        channel = draw_channel(network)

        # Loop over auxiliary variables
        for aux_idp_vals_idx = 1:aux_idp_vals_length
            # Set all auxiliary independent variables
            for Naux_idx = 1:Naux
                aux_idp_funcs[Naux_idx](network, aux_idp_vals[Naux_idx][aux_idp_vals_idx])
            end

            for assignment_method in simulation_params["assignment_methods"]
                # Allocate memory if this is the first method to be run
                if assignment_method == simulation_params["assignment_methods"][1]
                    raw_results[Ndrops_idx, 1, aux_idp_vals_idx] = SingleSimulationResults(AssignmentResults)
                end
                raw_results[Ndrops_idx, 1, aux_idp_vals_idx][stringify_method(assignment_method)] = assignment_method(channel, network)

                ProgressMeter.next!(progress)
            end
        end
    end

    return raw_results
end

##########################################################################
# Helper methods

# Given the type of methods that we are looping over, provide the other
# method to be run before/after. If there are several of those methods
# provided in the simulation_params, present warning.
function get_other_method(simulation_params, loop_over)
    # Dummies
    precoding_method() = nothing
    assignment_method() = nothing

    if loop_over == :precoding_methods
        if haskey(simulation_params, "assignment_methods")
            if length(simulation_params["assignment_methods"]) > 1
                Lumberjack.warn("Looping over precoding methods: will only use first assignment method provided.")
            end
            assignment_method(channel, network) = simulation_params["assignment_methods"][1](channel, network)
        else
            assignment_method(channel, network) = IDCellAssignment!(channel, network)
        end
    elseif loop_over == :assignment_methods
        if haskey(simulation_params, "precoding_methods")
            if length(simulation_params["precoding_methods"]) > 1
                Lumberjack.warn("Looping over assignment methods: will only use first precoding method provided.")
            end
            precoding_method(channel, network) = simulation_params["precoding_methods"][1](channel, network)
        else
            Lumberjack.error("No precoding method specified.")
        end
    end

    return precoding_method, assignment_method
end

# Get all auxiliary independent variable values
function get_aux_idp(simulation_params)
    if haskey(simulation_params, "aux_independent_variables")
        Naux = length(simulation_params["aux_independent_variables"])
        aux_idp_funcs = [ simulation_params["aux_independent_variables"][n][1]::Function for n = 1:Naux ]
        aux_idp_vals = [ simulation_params["aux_independent_variables"][n][2] for n = 1:Naux ]

        # Check that all independent variable vectors are the same length
        aux_idp_vals_length = length(simulation_params["aux_independent_variables"][1][2])
        for n = 2:Naux
            aux_idp_vals_length == length(simulation_params["aux_independent_variables"][n][2]) ? nothing : Lumberjack.error("Auxiliary independent variable vectors must have equal length.")
        end
    else
        Naux = 0
        aux_idp_funcs = Function[]
        aux_idp_vals = Vector{Int}[]
        aux_idp_vals_length = 1
    end

    return Naux, aux_idp_funcs, aux_idp_vals, aux_idp_vals_length
end

# Set all initial auxiliary parameters based on provided simulation_params
function set_initial_aux_params!(simulation_params, network)
    haskey(simulation_params, "aux_network_params") && set_aux_network_params!(network, simulation_params["aux_network_params"])
    haskey(simulation_params, "aux_precoding_params") && set_aux_precoding_params!(network, simulation_params["aux_precoding_params"])
    haskey(simulation_params, "aux_assignment_params") && set_aux_assignment_params!(network, simulation_params["aux_assignment_params"])
end


text:
module test_bukdu_router

using Test
using Bukdu
using HTTP.Messages: header

struct ExController <: ApplicationController
    conn::Conn
end

f1(::ExController) = nothing
f2(::ExController) = "a"
f3(::ExController) = SubString("ab")
f4(::ExController) = 0
f5(::ExController) = 3.14
f6(::ExController) = (1, 2)
f7(::ExController) = (a=1,)

g1(::ExController) = render(JSON, nothing)
g2(::ExController) = render(JSON, "a")
g3(::ExController) = render(JSON, SubString("ab"))
g4(::ExController) = render(JSON, 0)
g5(::ExController) = render(JSON, 3.14)
g6(::ExController) = render(JSON, (1, 2))
g7(::ExController) = render(JSON, (a=1,))

routes() do
    for func in [f1 f2 f3 f4 f5 f6 f7
                 g1 g2 g3 g4 g5 g6 g7]
        get(string("/", nameof(func)), ExController, func)
    end
end

result = Router.call(get, "/f2")
@test result.got == "a"
@test header(result.resp, "Content-Type") == "application/julia; charset=utf-8"
@test result.resp.body == Vector{UInt8}(repr("a"))

for (f, g) in zip([f1 f2 f3 f4 f5 f6 f7],
                  [g1 g2 g3 g4 g5 g6 g7])
    result1 = Router.call(get, string("/", nameof(f)))
    result2 = Router.call(get, string("/", nameof(g)))
    @test result1.got == result2.got
end

Routing.empty!()

end # module test_bukdu_router


text:
module CreepLaw

# This implements viscous creep laws and routines to compute with them
#
# Note that various simple creep laws are defined in this file;
# more complex ones (such as DislocationCreep) are in separate files
# in the same directory
#
# In case you want to add new creep laws, have a look at how the ones
# here are implemented. Please add tests as well!

using Base: Float64
using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam
using BibTeX
using ..MaterialParameters: MaterialParamsInfo
import Base.show, GeoParams.param_info

abstract type AbstractCreepLaw{T} <: AbstractMaterialParam end

export  computeCreepLaw_EpsII, computeCreepLaw_TauII,       # calculation routines
        CreepLawVariables,                                     # holds additional parameters required for calculations
        LinearViscous,
        PowerlawViscous,
        param_info

# NOTE: we will likely have to remove this, in favor of multiple dispatch options
"""
 Struct that holds parameters required for creep law calculations (such as P,T, grain size etc.)

You can assign the struct as
```julia-repl
 p = CreepLawVariables(P=0.0, T=0.0, d=1.0)
```
where you can also pass vectors or arrays as values.

Note that, if needed, this can be extended, w/out interfering with existing calculation
"""

@with_kw struct CreepLawVariables{_T,U1,U2,U3,U4}
    P::GeoUnit{_T,U1}  =   100.0MPa   # pressure
    T::GeoUnit{_T,U2}  =   500.0C     # temperature
    d::GeoUnit{_T,U3}  =   1.0cm      # grainsize
    f::GeoUnit{_T,U4}  =   0.0MPa     # water-fugacity
end
CreepLawVariables(args...) = CreepLawVariables(convert.(GeoUnit,args)...)

include("DislocationCreep.jl")
include("DiffusionCreep.jl")
include("Viscosity.jl")

# Linear viscous rheology ------------------------------------------------
"""
    LinearViscous(η=1e20Pa*s)

Defines a linear viscous creeplaw

The (isotopic) linear viscous rheology is given by
```math
    \\tau_{ij} = 2 \\eta \\dot{\\varepsilon}_{ij}
```
or
```math
    \\dot{\\varepsilon}_{ij}  = {\\tau_{ij}  \\over 2 \\eta }
```

where ``\\eta_0`` is the reference viscosity [Pa*s] at reference strain rate ``\\dot{\\varepsilon}_0``[1/s], and ``n`` the power law exponent [].
"""
@with_kw_noshow struct LinearViscous{T,U} <: AbstractCreepLaw{T}
    η::GeoUnit{T,U}     =   1e20Pa*s                # viscosity
    η_val::T            =   1.0
end
LinearViscous(args...) = LinearViscous(convert(GeoUnit,args[1]), args[2])

function param_info(s::LinearViscous) # info about the struct
    return MaterialParamsInfo(Equation = L"\tau_{ij} = 2 \eta  \dot{\varepsilon}_{ij}")
end

# Calculation routines for linear viscous rheologies
function computeCreepLaw_EpsII(TauII, s::LinearViscous, p::CreepLawVariables)
    @unpack η   = s

    return EpsII = (TauII/η)*0.5;
end

function computeCreepLaw_TauII(EpsII, s::LinearViscous, p::CreepLawVariables)
    @unpack η   = s

    return TauII = 2.0*(η*EpsII);
end

# Print info
function show(io::IO, g::LinearViscous)
    print(io, "Linear viscosity: η=$(g.η.val)")
end
#-------------------------------------------------------------------------

# Powerlaw viscous rheology ----------------------------------------------
"""
    PowerlawViscous(η0=1e18Pa*s, n=2.0NoUnits, ε0=1e-15/s)

Defines a power law viscous creeplaw as:

```math
    \\tau_{ij}^n  = 2 \\eta_0 \\left( \\dot{\\varepsilon}_{ij} \\over \\dot{\\varepsilon}_0 \\right)
```
where ``\\eta`` is the effective viscosity [Pa*s].
"""
@with_kw_noshow struct PowerlawViscous{T,U1,U2,U3} <: AbstractCreepLaw{T}
    η0::GeoUnit{T,U1}     =  1e18Pa*s            # reference viscosity
    n::GeoUnit{T,U2}      =  2.0*NoUnits         # powerlaw exponent
    ε0::GeoUnit{T,U3}     =  1e-15*1/s;          # reference strainrate
end
PowerlawViscous(a...) = PowerlawViscous(convert.(GeoUnit,a)...)

# Print info
function show(io::IO, g::PowerlawViscous)
    print(io, "Powerlaw viscosity: η0=$(g.η0.val), n=$(g.n.val), ε0=$(g.ε0.val) ")
end
#-------------------------------------------------------------------------

for myType in (:DiffusionCreep, :DislocationCreep)
    @eval begin
        computeCreepLaw_EpsII(TauII, a::$(myType), args) = computeCreepLaw_EpsII(TauII, a; args...)
        computeCreepLaw_TauII(EpsII, a::$(myType), args) = computeCreepLaw_TauII(EpsII, a; args...)
    end
end

# Help info for the calculation routines
"""
    computeCreepLaw_EpsII(TauII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the strainrate invariant ``\\dot{\\varepsilon}_{II}`` for a given deviatoric stress
invariant ``\\tau_{II}`` for any of the viscous creep laws implemented.
``p`` is a struct that can hold various parameters (such as ``P``, ``T``) that the creep law
may need for the calculations

```math
    \\dot{\\varepsilon}_{II}   = f( \\tau_{II} )
```

"""
computeCreepLaw_EpsII

"""
    computeCreepLaw_TauII(EpsII, s:<AbstractCreepLaw, p::CreepLawVariables)

Returns the deviatoric stress invariant ``\\tau_{II}`` for a given strain rate
invariant ``\\dot{\\varepsilon}_{II}`` for any of the viscous creep laws implemented.
``p`` is a struct that can hold various parameters (such as ``P``, ``T``) that the creep law
may need for the calculations

```math
    \\tau_{II}  = f( \\dot{\\varepsilon}_{II} )
```

"""
computeCreepLaw_TauII

end


text:
abstract type Spec end

specialize(::Type{ECGroup{P}}, spec::Spec; name = nothing) where P <: ECPoint = ECGroup{specialize(P, spec; name)}
specialize(::Type{ECGroup}, spec::Spec; name = nothing) = ECGroup{specialize(ECPoint, spec; name)}

function specialize(::Type{ECPoint{P}}, curve::Spec; name=nothing) where P <: AbstractPoint

    Q = specialize(P, curve)

    _order = order(curve)
    _cofactor = 1 # Need to update this one

    R = specialize(ECPoint{Q}, _order, _cofactor, name)

    return R
end

specialize(::Type{ECPoint}, spec::Spec; name = nothing) = specialize(ECPoint{AffinePoint}, spec; name)

##################### Macro for defining curve as also a group ###################

# macro def(constant_name, type, group_spec)

#     name_str = string(constant_name)

#     M = @__MODULE__

#     return esc(quote
#         const $constant_name = $M.specialize($type, $group_spec)

#         $M.name(::Type{$constant_name}) = $name_str

#         local ORDER = $group_spec.n
#         $M.order(::Type{$constant_name}) = ORDER

#         local GENERATOR = $constant_name($M.generator($group_spec)...)
#         $M.generator(::Type{$constant_name}) = GENERATOR
#     end)
# end

macro def(constant_name, type, group_spec)

    M = @__MODULE__

    return esc(quote
        const $constant_name = $M.specialize($type, $group_spec; name = $(QuoteNode(constant_name)))
    end)
end

include("curve_specs.jl")
include("modp_specs.jl")

### For now keeping things simple
# Need to think about a more sustainable strategy
function spec(x::Symbol)
    if x == :P_192
        return Curve_P_192
    elseif x == :P_244
        return Curve_P_244
    elseif x == :P_256
        return Curve_P_256
    elseif x == :P_384
        return Curve_P_384
    elseif x == :P_521
        return Curve_P_521
    else
        error("$x not implemented")
    end
end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
module IntervalLinearEquations
using ValidatedNumerics;

include("Types.jl")
include("sti.jl")
include("F.jl")
include("G.jl")

function solve{T}(name :: AbstractString, A :: Array{Interval{T}, 2}, B :: Array{Interval{T}, 1}, precision :: AbstractFloat, max_iteration :: Int, scale :: AbstractFloat)
    @assert precision > 0.0
    @assert scale <= 1.0 && scale > 0

    case_module = get_module(name)
    solver = initialize(case_module, A, B)
    condition = Types.StopConditions(precision, max_iteration)

    is_initial = true
    while is_initial || is_stop_iteration(solver, condition)
        solver.previous = solver.current
        calculatedSubDifferential = case_module.subDifferential(solver)
        equationValue = case_module.equation(solver)
        if det(calculatedSubDifferential) == 0
            break
        end
        solver.current = iterate(solver.previous, scale, calculatedSubDifferential, equationValue)
        solver.roots = vcat(solver.roots, [solver.current])
        solver.iteration += 1
        is_initial = false
    end
    solver.current.intervals
end

function initialize{T}(case_module, A :: Array{Interval{T}, 2}, B :: Array{Interval{T}, 1})
    @assert size(A, 1) == size(A, 2)
    @assert size(A, 1) == size(B, 1)

    system = Types.Configuration{T}(
        A, Types.IntervalVector(B, IntervalLinearEquations.sti.STI(B)),
        size(A, 1), size(A, 1)*2
    )
    initialInNumbers = case_module.initialConditions(system)
    initial = Types.IntervalVector(sti.reverseSTI(initialInNumbers), initialInNumbers)
    solver = Types.Solver{T}(
        initial,
        initial,
        initial,
        system,
        [initial],
        1
    )
    solver
end

function is_stop_iteration{T}(solver :: Types.Solver{T}, condition :: Types.StopConditions)
    get_iter_diff(solver) > condition.precision && solver.iteration > condition.max_iterations
end

function get_iter_diff{T}(solver :: Types.Solver{T})
    diff_for_dimensions = map(interval_diff, solver.current.intervals, solver.previous.intervals)
    diff = sqrt(sum(diff_for_dimensions))
    diff
end

function interval_diff{T}(current :: Interval{T}, previous :: Interval{T})
    hi_diff = current.hi - previous.hi
    lo_diff = current.lo - previous.lo
    sqrt(hi_diff^2 + lo_diff^2)
end

function iterate(previous, scale, subdiff, equation_value)
    new_sti = previous.sti - scale * inv(subdiff) * equation_value
    Types.IntervalVector(IntervalLinearEquations.sti.reverseSTI(new_sti), new_sti)
end

function get_module(name :: AbstractString)
    @assert name == "F" || name == "G"

    if name == "F"
        F
    else
        G
    end
end

end


text:
using ThreeBodyDecay

ThreeBodySpins(1//2, 1//2, 3//2, 3//2) |> x2

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module FastTransforms

using ArrayLayouts, FastGaussQuadrature, FillArrays, LinearAlgebra,
      Reexport, SpecialFunctions, ToeplitzMatrices

import DSP

@reexport using AbstractFFTs
@reexport using FFTW

import Base: convert, unsafe_convert, eltype, ndims, adjoint, transpose, show,
             *, \, inv, length, size, view, getindex

import Base.GMP: Limb

import AbstractFFTs: Plan, ScaledPlan,
                     fft, ifft, bfft, fft!, ifft!, bfft!, rfft, irfft, brfft,
                     plan_fft, plan_ifft, plan_bfft, plan_fft!, plan_ifft!,
                     plan_bfft!, plan_rfft, plan_irfft, plan_brfft,
                     fftshift, ifftshift, rfft_output_size, brfft_output_size,
                     plan_inv, normalization

import DSP: conv

import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!,
             plan_dct, plan_idct, fftwNumber

import FastGaussQuadrature: unweightedgausshermite

import FillArrays: AbstractFill, getindex_value

import LinearAlgebra: mul!, lmul!, ldiv!

export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
       lag2lag, jac2ultra, ultra2jac, jac2cheb,
       cheb2jac, ultra2cheb, cheb2ultra, associatedjac2jac,
       sph2fourier, sphv2fourier, disk2cxf, rectdisk2cheb, tri2cheb, tet2cheb,
       fourier2sph, fourier2sphv, cxf2disk, cheb2rectdisk, cheb2tri, cheb2tet

export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
       plan_lag2lag, plan_jac2ultra, plan_ultra2jac, plan_jac2cheb,
       plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra, plan_associatedjac2jac,
       plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis,
       plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
       plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,
       plan_rectdisk2cheb, plan_rectdisk_synthesis, plan_rectdisk_analysis,
       plan_tri2cheb, plan_tri_synthesis, plan_tri_analysis,
       plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis,
       plan_spinsph2fourier, plan_spinsph_synthesis, plan_spinsph_analysis

include("clenshaw.jl")

include("libfasttransforms.jl")

export nufft, nufft1, nufft2, nufft3, inufft1, inufft2

export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3,
       plan_inufft1, plan_inufft2

include("nufft.jl")
include("inufft.jl")

export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!,
       paduapoints

export plan_paduatransform!, plan_ipaduatransform!

include("PaduaTransform.jl")

export chebyshevtransform, ichebyshevtransform,
       chebyshevtransform!, ichebyshevtransform!,
       chebyshevutransform, ichebyshevutransform,
       chebyshevutransform!, ichebyshevutransform!, chebyshevpoints

export plan_chebyshevtransform, plan_ichebyshevtransform,
       plan_chebyshevtransform!, plan_ichebyshevtransform!,
       plan_chebyshevutransform, plan_ichebyshevutransform,
       plan_chebyshevutransform!, plan_ichebyshevutransform!

include("chebyshevtransform.jl")

export clenshawcurtisnodes, clenshawcurtisweights, fejernodes1, fejerweights1,
       fejernodes2, fejerweights2

export plan_clenshawcurtis, plan_fejer1, plan_fejer2

include("clenshawcurtis.jl")
include("fejer.jl")

export weightedhermitetransform, iweightedhermitetransform

include("hermite.jl")

include("fftBigFloat.jl")

export gaunt

include("gaunt.jl")

export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
       sphvones, sphvzeros, sphvrand, sphvrandn,
       diskones, diskzeros, diskrand, diskrandn,
       rectdiskones, rectdiskzeros, rectdiskrand, rectdiskrandn,
       triones, trizeros, trirand, trirandn, trievaluate,
       tetones, tetzeros, tetrand, tetrandn,
       spinsphones, spinsphzeros, spinsphrand, spinsphrandn

include("specialfunctions.jl")

end # module


text:
module TestDescentUpdates
using Test
using LinearAlgebra
using OptimizationAlgorithms
using OptimizationAlgorithms: Gradient, ScaledGradient, StoppingCriterion, fixedpoint!
using OptimizationAlgorithms: PolyakMomentum, Nesterov, ADAM, AMSGRAD
using ForwardDiff: gradient

ε = 1e-4
@testset "PolyakMomentum" begin
    n = 2
    A = randn(n, n)
    A = A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    x = randn(n)
    G = Gradient(f, x)
    α = .1
    β = .5
    G = PolyakMomentum(G, copy(x), α, β)
    fixedpoint!(G, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε
end

@testset "Nesterov" begin
    # TODO: better test for nesterov
    n = 2
    A = randn(n, n)
    A = A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    x = randn(n)
    α = 5e-2
    β = .5
    G = Gradient(f, x)
    G = Nesterov(G, x, α, β)
    fixedpoint!(G, x, StoppingCriterion(x, dx = 1e-2ε))
    @test norm(gradient(f, x)) < ε
end

@testset "ADAM" begin
    n = 2
    A = randn(n, n)
    A = A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    x = randn(n)
    G = Gradient(f, x)
    α = .1
    β₁ = 0.5
    β₂ = 0.999 # this should be large
    G = ADAM(G, x, α, β₁, β₂)
    fixedpoint!(G, x, StoppingCriterion(x, dx = 1e-2ε, maxiter = 1024))
    @test norm(gradient(f, x)) < ε
end

@testset "AMSGRAD" begin
    n = 2
    A = randn(n, n)
    A = A'A + I
    μ = randn(n)
    f = x -> (x-μ)'A*(x-μ)
    x = randn(n)
    G = Gradient(f, x)
    α = .1
    β₁ = 0.5
    β₂ = 0.999 # this should be large
    G = AMSGRAD(G, x, α, β₁, β₂)
    fixedpoint!(G, x, StoppingCriterion(x, dx = 1e-2ε, maxiter = 1024))
    @test norm(gradient(f, x)) < ε
end

end # TestDescentUpdates


text:
struct MutualInformationContainer{H}
    hist::H
    pxy::Array{Float32,2}
    px::Array{Float32,2}
    py::Array{Float32,2}
    px_py::Array{Float32,2}
    nzs::BitArray{2}

    function MutualInformationContainer(hist::H) where {H}
        pxy = counts(hist) ./ sum(counts(hist))
        px = sum(pxy, dims = 2)
        py = sum(pxy, dims = 1)
        px_py = px * py
        nzs = pxy .> 0

        new{H}(hist, pxy, px, py, px_py, nzs)
    end
end

function _mutual_information!(mi::MutualInformationContainer)
    mi.pxy .= counts(mi.hist) ./ sum(counts(mi.hist))
    sum!(mi.px, mi.pxy)
    sum!(mi.py, mi.pxy)
    mi.px_py .= mi.px * mi.py
    mi.nzs .= mi.pxy .> 0
    pxys = mi.pxy[mi.nzs]
    sum(pxys .* log.(pxys ./ mi.px_py[mi.nzs]))
end

"""
    mutual_information!(mi::MutualInformationContainer, x, y)

Computes the mutual information between the two variables `x` and `y`. The histogram within `mi` must be of the correct
type to handle the formats of `x` and `y`.
"""
function mutual_information!(mi::MutualInformationContainer, x, y)
    zero!(mi.hist)
    increment_bins!(mi.hist, x, y)
    _mutual_information!(mi)
end

"""
    mutual_information!(
        mi::MutualInformationContainer,
        fixed,
        buffer,
        full_image,
        moving_bbox,
        range_x,
        range_y,
        ::Missing;
        set_buffer!,
        get_buffer_crop,
        prefilter_frame_crop! = x -> nothing,
    )

Calculates the mutual information of two images at all shifts within the `range_x` and `range_y`. The `fixed` image
must already be filtered. This will set the `buffer` and filter its contents using `prefilter_frame_crop!`.
"""
function mutual_information!(
    mi::MutualInformationContainer,
    fixed,
    buffer,
    full_image,
    moving_bbox,
    range_x,
    range_y,
    ::Missing;
    set_buffer!,
    get_buffer_crop,
    prefilter_frame_crop! = x -> nothing,
)
    w, h = size(fixed)

    mis = OffsetArray(
        Array{Float32}(undef, length(range_x), length(range_y)),
        range_x,
        range_y,
    )
    fixed_vec = vec(fixed)

    # Crop and prefilter a section of `current_frame` big enough to handle the shift extents.
    set_buffer!(buffer, full_image, moving_bbox)
    prefilter_frame_crop!(buffer)

    @inbounds for shift_x in range_x
        @inbounds for shift_y in range_y
            moving_vec = vec(get_buffer_crop(buffer, moving_bbox, shift_x, shift_y))
            mis[shift_x, shift_y] = mutual_information!(mi, fixed_vec, moving_vec)
        end
    end

    return mis
end

"""
    mutual_information!(
        mi::MutualInformationContainer,
        fixed,
        buffer,
        ::Any,
        moving_bbox,
        range_x,
        range_y,
        prev_mis::AbstractArray{Float32,2};
        get_buffer_crop,
        kwargs...
    )

Calculates the mutual information of two images at all shifts within the `range_x` and `range_y`. Warm-starts the
evaluation using previous results (`prev_mis`; the return value from a previous call of this function) and using the
previously set and filtered contents of the `buffer`.
"""
function mutual_information!(
    mi::MutualInformationContainer,
    fixed,
    buffer,
    ::Any,
    moving_bbox,
    range_x,
    range_y,
    prev_mis::AbstractArray{Float32,2};
    get_buffer_crop,
    kwargs...,
)
    w, h = size(fixed)

    mis = OffsetArray(
        Array{Float32}(undef, length(range_x), length(range_y)),
        range_x,
        range_y,
    )

    prev_range_x = axes(prev_mis, 1)
    prev_range_y = axes(prev_mis, 2)

    fixed_vec = vec(fixed)

    # No need to extract and prefilter a crop to fill the buffer here because it is done in `mutual_information!`
    # where `prev_mis` is of type `Missing`. We just reuse the prefiltered crop here.

    @inbounds for shift_x in range_x
        @inbounds for shift_y in range_y
            if shift_x ∈ prev_range_x && shift_y ∈ prev_range_y
                mis[shift_x, shift_y] = prev_mis[shift_x, shift_y]
            else
                moving_vec = vec(get_buffer_crop(buffer, moving_bbox, shift_x, shift_y))
                mis[shift_x, shift_y] = mutual_information!(mi, fixed_vec, moving_vec)
            end
        end
    end

    return mis
end


text:
type receiver
    nr::Int64
    loc::Array{<:Real}
    nloc::Array{Int64}
    BDnloc::Array{Int64}
end

#=== receiver ===#
function initrec(
    nr::Int64,
    loc::Array{<:Real},
    model::Union{elmod2d,acmod2d})

    ext = model.medium.ext
    iflag = model.medium.iflag
    if length(loc) == 2 #=== receiver 2d ===#
        dz = model.medium.dz
        dx = model.medium.dx
        nloc = zeros(Int64,nr,2)
        for i = 1:nr
            nloc[i,1] = Int64(round(loc[i,1]./dz))
            nloc[i,2] = Int64(round(loc[i,2]./dx))
        end
        if iflag == 1 # free surface
            BDnloc = [nloc[:,1] nloc[:,2].+ext]
        elseif iflag == 2 # unlimited space
            BDnloc = [nloc[:,1].+ext nloc[:,2].+ext]
        end
    elseif length(loc) == 3 #=== receiver 3d ===#
        dz = model.medium.dz
        dx = model.medium.dx
        dy = model.medium.dy
        nloc = zeros(Int64,nr,3)
        for i = 1:nr
            nloc[i,1] = Int64(round(loc[i,1]./dz))
            nloc[i,2] = Int64(round(loc[i,2]./dx))
            nloc[i,3] = Int64(round(loc[i,3]./dy))
        end
        if iflag == 1 # free surface
            BDnloc = [nloc[:,1] nloc[:,2].+ext nloc[:,3].+ext]
        elseif iflag == 2 # unlimited space
            BDnloc = [nloc[:,1].+ext nloc[:,2].+ext nloc[:,3].+ext]
        end
    end
    return receiver(nr,loc,nloc,BDnloc)
end


text:
module DataReaders
using FileIO
using CSVFiles
using FeatherFiles
using Feather
using DataFrames
using HDF5
using JLD
using Parquet
using Dates

export DataReader, fit!, transform!
export Transformer
export datareaderrun

using TSML
using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

const FILEFMT = Dict("csv"=>Val(:csv),"feather"=>Val(:feather),"hdf5"=>Val(:hdf5),
          "h5"=>Val(:hdf5),"parquet"=>Val(:parquet),"jld"=>Val(:jld))

const DATEFMT = "d/m/yyyy HH:MM"

mutable struct DataReader <: Transformer
    model
    args

    function DataReader(args=Dict())
        default_args=Dict(
            :filename => "",
            :dateformat => DATEFMT,
            :impl_args => Dict()
        )
        new(nothing,mergedict(default_args,args))
    end
end

function fit!(dtr::DataReader,x::T=[],y::Vector=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fname = dtr.args[:filename]
    fmt = dtr.args[:dateformat]
    (fname != "" && isfile(fname) && fmt != "") || error("missing filename or date format: (",fname,"), (",fmt,")")
    dtr.model = dtr.args
end

function transform!(dtr::DataReader,x::T=[]) where {T<:Union{DataFrame,Vector,Matrix}}
    fullname = dtr.args[:filename]
    fmt = dtr.args[:dateformat]
    (fullname != "" && isfile(fullname) && fmt != "") || error("missing filename or date format: (",fullname,"), (",fmt,")")
    fname = basename(fullname)
    fname != "" || error("filename is empty")
    fn,ext=split(fname,".")
    ext != "" || error("no filename extension format")
    ext in keys(FILEFMT)  || error("extension not recognized "*ext)
    # dispatch based on extension
    readfmt(FILEFMT[ext],fullname,fmt)
end

function readfmt(::T,fname::String,datefmt::String) where {T<:Union{Val{:csv},Val{:feather},Val{:parquet}}}
    df = load(fname) |> DataFrame
    ncol(df) == 2 || error("dataframe should have only two columns: Date,Value")
    rename!(df,names(df)[1]=>:Date,names(df)[2]=>:Value)
    if eltype(df[:Date]) <: DateTime
        return df
    else
        df[:Date] = DateTime.(df[:Date],datefmt)
        return df
    end
end

function readfmt(atype::Union{Val{:hdf5},Val{:jld}},fname::String,datefmt::String)
    # R library(hdf5r); dat=H5File$new(fname);
    # dat[["dateval/date"]];dat[["dateval/value"]]
    if atype == Val{:hdf5}
        fileopen = h5open
    else
        fileopen = jldopen
    end
    ldate = fileopen(fname, "r") do file
        read(file, "dateval/date")
    end
    lvalue = fileopen(fname, "r") do file
        read(file, "dateval/value")
    end
    DataFrame(date=DateTime.(ldate,datefmt),Value=lvalue)
end

end


text:
#=
general contract for search functions f in SearchBenchmarks:

    f(pattern::SearchSequence)::Tuple

initializes a search structure as Tuple, which includes pattern and
derived data for higher order search functions.
It is put in a separate method to allow for simpler timing
measurement of initialization code.

It is used in the search "kernel", which consists of the main
search loop:

    f(ss::Tuple,text::SearchSequence,start::Int)

=#

"""
This is the original implementation in julia Base,
method _searchindex(s::SearchSequence, t::SearchSequence, i::Integer).

Some technical changes are:
 * split into initialization and kernel
 * benchmark measurements, compiled away if last parameter is nothing
 * code is valid for String and ByteArray sequences
"""
function bloom_v0 end

function bloom_v0(t::SearchSequence)
    n = sizeof(t)
    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j - 1
        end
    end
    return t,bloom_mask,skip,tlast
end

function bloom_v0(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_v0(s,bloom_v0(t),i)
end

function bloom_v0(s::SearchSequence, p::Tuple,i::Integer,sv::MaybeVector=nothing)
    ( t,bloom_mask,skip,tlast) = p
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)

    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end

    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

#=
    bloom_mask = UInt64(0)
    skip = n - 1
    tlast = _nthbyte(t,n)
    for j in 1:n
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast && j < n
            skip = n - j - 1
        end
    end
=#
    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i -= 1
    while i <= w
        if DOSTATS loops += 1 end
        if _nthbyte(s,i+n) == tlast
            # check candidate
            j = 0
            while j < n - 1
                if _nthbyte(s,i+j+1) != _nthbyte(t,j+1)
                    break
                end
                j += 1
            end

            # match found
            if j == n - 1
                if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                return i+1
            end

            # no match, try to rule out the next character
            if DOSTATS bloomtests += 1 end
            if i < w && bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                if DOSTATS bloomskips += 1 end
                i += n
            else
                i += skip
            end
        elseif i < w
            if DOSTATS bloomtests += 1 end
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i+n+1)) == 0
                if DOSTATS bloomskips += 1 end
                i += n
            end
        end
        i += 1
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
    0
end


text:
using OED
using Base.Test
using jInv.Mesh
using jInv.Utils

domain = [-1 1 -1 1.]; n = 4*[8 8]
Mr     = getRegularMesh(domain,n)
x,y    = getCellCenteredAxes(Mr)
xc     = getCellCenteredGrid(Mr)

dataT     = exp.(-xc[:,1].^2- xc[:,2].^2)
dataT     = reshape(xc*rand(2),tuple(n...))

Tc,dT = linearInter(dataT,domain,xc,doDerivative=true)
A0 = getLinearInterMatrix(domain,n,xc)
Tt = A0*vec(dataT)
@test norm(vec(Tc)-vec(dataT))/norm(vec(dataT)) < 1e-8
@test norm(vec(Tt)-vec(dataT))/norm(vec(dataT)) < 1e-8

function linearInterTestFun(x,v=[])
	Tc,dT = linearInter(dataT,domain,x,doDerivative=true)

	if isempty(v)
		return Tc
	else
		return Tc,dT*v
	end
end
x0 = vec(xc+1e-3*randn(size(xc)))
chkDer, = checkDerivative(linearInterTestFun,x0)

@test chkDer
println("check getLinearInterMatrix with points in domain")
xt = reshape(copy(x0),Mr.nc,2);
 xt[:,1]=xc[:,1];  xt[:,2]*=2
A = getLinearInterMatrix(copy(domain),copy(n),copy(xt))
T = linearInter(copy(dataT),copy(domain),copy(xt))
Tt = A*vec(dataT)
@test norm(Tt-T)/norm(T) < 1e-10
println("check getLinearInterMatrix with points out of domain")
A = getLinearInterMatrix(.5*domain,n,x0)
T = linearInter(dataT*0+1,.5*domain,x0)
Tt = A*vec(0*dataT+1)
@test norm(Tt-T)/norm(T) < 1e-10


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using GlobalSensitivity
using Distributions
using Plots

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

tmax = 10000.0
tspan = (0.0,tmax);

cb_ss = TerminateSteadyState();

N = 1000.0;
u0 = [990.0,10.0,0.0];
p = [0.05,10.0,0.25]; # β,c,γ

n_samples = 1000 # Number of samples
# Parameters are β, c, γ, I₀
lb = [0.01, 5.0, 0.1, 1.0]
ub = [0.1, 20.0, 1.0, 50.0]
n_params = 4;

prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

f1 = function(pu0)
  p = pu0[1:3]
  I0 = pu0[4]
  u0 = [N-I0,I0,0.0]
  prob = remake(prob_ode;p=p,u=u0)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  [maximum(sol[2,:]), sol.t[argmax(sol[2,:])], sol[end][3]]
end;

m_morris = gsa(f1, Morris(num_trajectory=n_samples), [[lb[i],ub[i]] for i in 1:n_params]);

m_morris.means

m_morris.variances

m_sobol = gsa(f1, Sobol(), [[lb[i],ub[i]] for i in 1:n_params],N=n_samples);

m_sobol.ST

m_sobol.S1

m_regression = gsa(f1, RegressionGSA(rank=true), [[lb[i],ub[i]] for i in 1:n_params]; samples = n_samples);

m_regression.partial_correlation

m_regression.partial_rank_correlation

m_efast = gsa(f1, eFAST(), [[lb[i],ub[i]] for i in 1:n_params]; n = n_samples);

m_efast.ST

m_efast.S1

pf1 = function (pu0)
  p = pu0[1:3,:]
  I0 = pu0[4,:]
  prob_func(prob,i,repeat) = remake(prob;p=p[:,i],u=[N-I0[i],I0[i],0.0])
  ensemble_prob = EnsembleProblem(prob_ode,prob_func=prob_func)
  sol = solve(ensemble_prob,ROS34PW3(),EnsembleThreads();trajectories=size(p,2))
  out = zeros(3,size(p,2))
  for i in 1:size(p,2)
    out[1,i] = maximum(sol[i][2,:])
    out[2,i] = sol[i].t[argmax(sol[i][2,:])]
    out[3,i] = sol[i][end][3]
  end
  out
end;

m_efast_parallel = gsa(pf1, eFAST(), [[lb[i],ub[i]] for i in 1:n_params]; n = n_samples, batch = true);