text:
"""
shoot(phi,xspan,lval,lder,rval,rder,init)

Use the shooting method to solve a two-point boundary value problem. The ODE is
u'' = `phi`(x,u,u') for x in `xspan`. Specify a function value or derivative at
the left endpoint using `lval` and `lder`, respectively, and similarly for the
right endpoint  using `rval` and `rder`. (Use an empty array to denote an
unknown quantity.) The value `init` is an initial guess for whichever value is
missing at the left endpoint.

Return vectors for the nodes, the values of u, and the values of u'.
    """
function shoot(phi,xspan,lval,lder,rval,rder,init)

    # Tolerances for IVP solver and rootfinder.
    ivp_opt = 1e-6
    optim_opt = 1e-5

    # Evaluate the difference between computed and target values at x=b.
    function objective(s)
        # Combine s with the known left endpoint value.
        v_init = isempty(lder) ?  [ lval; s[1] ]  :  [ s[1]; lder ]

        # ODE posed as a first-order equation in 2 variables.
        function shootivp(v,p,x)
            [ v[2]; phi(x,v[1],v[2]) ]
        end

        IVP = ODEProblem(shootivp,v_init,xspan)
        sol = solve(IVP,abstol=ivp_opt,reltol=ivp_opt)
        x = sol.t;  v = sol;

        return isempty(rder) ? v[1,end] - rval  :  v[2,end] - rder
    end

    # Find the unknown quantity at x=a by rootfinding.
    x = [];  v = [];   # the values will be overwritten
    s = nlsolve(objective,[init],xtol=optim_opt).zero

    # Don't need to solve the IVP again. It was done within the
    # objective function already.
    u = v[1,:]            # solution
    dudx = v[2,:]         # derivative

    return x,u,dudx
end

"""
    diffmat2(n,xspan)

Compute 2nd-order-accurate differentiation matrices on `n`+1 points in the
interval `xspan`. Return a vector of nodes, and the matrices for the first
and second derivatives.
"""
function diffmat2(n,xspan)
    a,b = xspan
    h = (b-a)/n
    x = [ a + i*h for i=0:n ]   # nodes

    # Define most of Dx by its diagonals.
    dp = fill(0.5/h,n)        # superdiagonal
    dm = fill(-0.5/h,n)       # subdiagonal
    Dx = diagm(-1=>dm,1=>dp)

    # Fix first and last rows.
    Dx[1,1:3] = [-1.5,2,-0.5]/h
    Dx[n+1,n-1:n+1] = [0.5,-2,1.5]/h

    # Define most of Dxx by its diagonals.
    d0 =  fill(-2/h^2,n+1)    # main diagonal
    dp =  ones(n)/h^2         # superdiagonal and subdiagonal
    Dxx = diagm(-1=>dp,0=>d0,1=>dp)

    # Fix first and last rows.
    Dxx[1,1:4] = [2,-5,4,-1]/h^2
    Dxx[n+1,n-2:n+1] = [-1,4,-5,2]/h^2

    return x,Dx,Dxx
end

"""
    diffcheb(n,xspan)

Compute Chebyshev differentiation matrices on `n`+1 points in the
interval `xspan`. Return a vector of nodes, and the matrices for the first
and second derivatives.
"""
function diffcheb(n,xspan)
    x = [ -cos( k*pi/n ) for k=0:n ]    # nodes in [-1,1]
    Dx = zeros(n+1,n+1)
    c = [2; ones(n-1); 2];    # endpoint factors

    # Off-diagonal entries
    Dx = [ (-1)^(i+j)*c[i+1]/(c[j+1]*(x[i+1]-x[j+1])) for i=0:n, j=0:n ]

    # Diagonal entries
    Dx[isinf.(Dx)] .= 0              # fix divisions by zero on diagonal
    s = sum(Dx,dims=2)
    Dx -= diagm(0=>s[:,1])           # "negative sum trick"

    # Transplant to [a,b]
    a,b = xspan
    x = @. a + (b-a)*(x+1)/2
    Dx = 2*Dx/(b-a)

    # Second derivative
    Dxx = Dx^2

    return x,Dx,Dxx
end

"""
     bvplin(p,q,r,xspan,lval,rval,n)

Use finite differences to solve a linear bopundary value problem. The ODE is
u''+`p`(x)u'+`q`(x)u = `r`(x) on the interval `xspan`, with endpoint function
values given as `lval` and `rval`. There will be `n`+1 equally spaced nodes,
including the endpoints.

Return vectors of the nodes and the solution values.
"""
function bvplin(p,q,r,xspan,lval,rval,n)
    x,Dx,Dxx = diffmat2(n,xspan)

    P = diagm(0=>p.(x))
    Q = diagm(0=>q.(x))
    L = Dxx + P*Dx + Q     # ODE expressed at the nodes

    # Replace first and last rows using boundary conditions.
    I = Diagonal(ones(n+1))
    A = [ I[[1],:]; L[2:n,:]; I[[n+1],:] ]
    b = [ lval; r.(x[2:n]); rval ]

    # Solve the system.
    u = A\b

    return x,u
end

"""
    bvp(phi,xspan,lval,lder,rval,rder,init)

Use finite differences to solve a two-point boundary value problem. The ODE is
u'' = `phi`(x,u,u') for x in `xspan`. Specify a function value or derivative at
the left endpoint using `lval` and `lder`, respectively, and similarly for the
right endpoint  using `rval` and `rder`. (Use an empty array to denote an
unknown quantity.) The value `init` is an initial guess for whichever value is
missing at the left endpoint.

Return vectors for the nodes and the values of u.
"""
function bvp(phi,xspan,lval,lder,rval,rder,init)
    n = length(init) - 1
    x,Dx,Dxx = diffmat2(n,xspan)
    h = x[2]-x[1]

    function residual(u)
        # Compute the difference between u'' and phi(x,u,u') at the
        # interior nodes and appends the error at the boundaries.
        dudx = Dx*u                   # discrete u'
        d2udx2 = Dxx*u                # discrete u''
        f = d2udx2 - phi.(x,u,dudx)

        # Replace first and last values by boundary conditions.
        f[1] = isempty(lder) ? (u[1] - lval)/h^2 : (dudx[1] - lder)/h
        f[n+1] = isempty(rder) ? (u[n+1] - rval)/h^2 : (dudx[n+1] - rder)/h
        return f
    end

    u = levenberg(residual,init)
    return x,u[:,end]
end

"""
    fem(c,s,f,a,b,n)

Use a piecewise linear finite element method to solve a two-point boundary
value problem. The ODE is (`c`(x)u')' + `s`(x)u = `f`(x) on the interval
[`a`,`b`], and the boundary values are zero. The discretization uses `n` equal
subintervals.

Return vectors for the nodes and the values of u.
"""
function fem(c,s,f,a,b,n)
    # Define the grid.
    h = (b-a)/n
    x = @. a + h*(0:n)

    # Templates for the subinterval matrix and vector contributions.
    Ke = [1 -1; -1 1]
    Me = (1/6)*[2 1; 1 2]
    fe = (1/2)*[1; 1]

    # Evaluate coefficent functions and find average values.
    cval = c.(x);   cbar = (cval[1:n]+cval[2:n+1]) / 2;
    sval = s.(x);   sbar = (sval[1:n]+sval[2:n+1]) / 2;
    fval = f.(x);   fbar = (fval[1:n]+fval[2:n+1]) / 2;

    # Assemble global system, one interval at a time.
    K = zeros(n-1,n-1);  M = zeros(n-1,n-1);  f = zeros(n-1);
    K[1,1] = cbar[1]/h;  M[1,1] = sbar[1]*h/3;  f[1] = fbar[1]*h/2;
    K[n-1,n-1] = cbar[n]/h;  M[n-1,n-1] = sbar[n]*h/3;  f[n-1] = fbar[n]*h/2;
    for k = 2:n-1
      K[k-1:k,k-1:k] += (cbar[k]/h) * Ke
      M[k-1:k,k-1:k] += (sbar[k]*h) * Me
      f[k-1:k] += (fbar[k]*h) * fe
    end

    # Solve system for the interior values.
    u = (K+M) \ f
    u = [0; u; 0]      # put the boundary values into the result

    return x,u
end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
module DfpWrapper

using JavaCall

const CalculusFieldElement = @jimport org.hipparchus.CalculusFieldElement
const Class = @jimport java.lang.Class
const Dfp = @jimport org.hipparchus.dfp.Dfp
const DfpDec = @jimport org.hipparchus.dfp.DfpDec
const DfpField = @jimport org.hipparchus.dfp.DfpField
const DfpField_RoundingMode = @jimport org.hipparchus.dfp.DfpField$RoundingMode
const DfpMath = @jimport org.hipparchus.dfp.DfpMath
const FieldSinCos = @jimport org.hipparchus.util.FieldSinCos
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const RealFieldElement = @jimport org.hipparchus.RealFieldElement

export Dfp
export DfpDec
export DfpField
export DfpMath
export abs
export acos
export acosh
export add
export asin
export asinh
export atan
export atan2
export atanh
export cbrt
export ceil
export classify
export clear_ieee_flags
export compute_exp
export compute_ln
export copy_sign
export copysign
export cos
export cosh
export divide
export dotrap
export equals
export exp
export expm1
export floor
export get_class
export get_deg_to_rad
export get_e
export get_e_split
export get_exponent
export get_extended_field
export get_field
export get_ieee_flags
export get_ln10
export get_ln2
export get_ln2_split
export get_ln5
export get_ln5_split
export get_one
export get_pi
export get_pi_split
export get_rad_to_deg
export get_radix_digits
export get_real
export get_rounding_mode
export get_runtime_class
export get_sqr2
export get_sqr2_reciprocal
export get_sqr2_split
export get_sqr3
export get_sqr3_reciprocal
export get_two
export get_zero
export greater_than
export hash_code
export hypot
export int_log10
export int_value
export is_infinite
export is_na_n
export is_zero
export less_than
export linear_combination
export log
export log10
export log10k
export log1p
export multiply
export negate
export negative_or_null
export new_dfp
export new_instance
export next_after
export notify
export notify_all
export positive_or_null
export pow
export power10
export power10k
export reciprocal
export remainder
export rint
export root_n
export round
export scalb
export set_ieee_flags
export set_ieee_flags_bits
export set_rounding_mode
export signum
export sin
export sin_cos
export sinh
export sqrt
export strictly_negative
export strictly_positive
export subtract
export tan
export tanh
export to_degrees
export to_double
export to_radians
export to_split_double
export to_string
export unequal
export wait

include("dfp.jl")
include("dfp_dec.jl")
include("dfp_field.jl")
include("dfp_math.jl")

end

ext: jl


text:
using StatsPlots, DelimitedFiles

# generate the dats
#=
include("../SBEM/SVlib.jl")
σe = exp(-0.736/2.0)
ρ = 0.9
σu = 0.363
θtrue = [σe, ρ, σu] # true param values, on param space
n = 1000 # sample size
burnin = 100
shocks_u = randn(n+burnin,1)
shocks_e = randn(n+burnin,1)
m = SVmodel(σu, ρ, σe, n, shocks_u, shocks_e, true)
=#
y = readdlm("svdata.txt")
p1 = plot(y)
p2 = density(y)
plot(p1, p2, layout=(2,1))
gui()
savefig("svdata.svg")
dstats(y)

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
module TestWoodburyFactorizations
using WoodburyFactorizations
using LazyInverses
using LinearAlgebra
using Test

function getW(n, m; diagonal = true, symmetric = false)
    A = diagonal ? Diagonal(exp.(.1randn(n))) + I(n) : randn(n, n)
    A = symmetric && !diagonal ? A'A : A
    U = randn(n, m)
    C = diagonal ? Diagonal(exp.(.1randn(m))) + I(m) : randn(m, m)
    C = symmetric && !diagonal ? C'C : C
    V = symmetric ? U' : randn(m, n)
    W = Woodbury(A, U, C, V)
end

Base.AbstractMatrix(F::BunchKaufman) = Symmetric(Matrix(F))

@testset "Woodbury" begin
    n = 5
    m = 2
    W = getW(n, m, symmetric = true)
    MW = Matrix(W)
    @test size(W) == (n, n)
    x = randn(size(W, 2))

    # multiplication
    @test W*x ≈ MW*x
    @test x'W ≈ x'MW
    @test dot(x, W, x) ≈ dot(x, MW*x)
    X = randn(size(W, 2), 3)
    @test W*X ≈ MW*X
    @test X'*W ≈ X'*MW
    @test W == copy(W) # testing isequal and copy
    @test W == deepcopy(W) # testing isequal and deepcopy

    # in-place multiplication
    y = randn(size(W, 1)) # with vector
    α, β = randn(2)
    b = α*(W*x) + β*y
    mul!(y, W, x, α, β)
    @test y ≈ b

    Y = randn(size(W, 1), size(X, 2)) # with matrix
    α, β = randn(2)
    B = α*(W*X) + β*Y
    mul!(Y, W, X, α, β)
    @test Y ≈ B

    @test α*W isa Woodbury
    @test Matrix(α * W) ≈ α * Matrix(W) # with scalar
    @test Matrix(W * α) ≈ Matrix(W) * α # with scalar
    @test Matrix(W / α) ≈ Matrix(W) / α # with scalar
    @test Matrix(α \ W) ≈ α \ Matrix(W) # with scalar
    @test logdet(abs(α)*W) ≈ logdet(abs(α)*Matrix(W)) # checking that logdet works correctly

    @test eltype(W) == Float64
    @test issymmetric(W)
    @test ishermitian(W)
    # @test ishermitian(MW)
    # @test ishermitian(Woodbury(W.A, W.U, W.C, copy(W.U)')) # testing edge case
    @test !issymmetric(getW(n, m))
    @test !ishermitian(getW(n, m))

    # test solves
    @test W\(W*x) ≈ x
    @test (x'*W)/W ≈ x'

    # factorization
    n = 1024
    m = 3
    W = getW(n, m, symmetric = true)
    MW = Matrix(W)
    F = factorize(W)
    @test F isa Inverse
    x = randn(n)
    @test W \ x ≈ F \ x
    @test MW \ x ≈ F \ x

    n = 4
    m = 3
    W = getW(n, m, symmetric = true)
    MW = Matrix(W)

    ## determinant
    @test det(W) ≈ det(MW)
    @test logdet(W) ≈ logdet(MW)
    @test logabsdet(W)[1] ≈ logabsdet(MW)[1]
    @test logabsdet(W)[2] ≈ logabsdet(MW)[2]
    @test det(inverse(W)) ≈ det(Inverse(W))
    @test det(inverse(W)) ≈ det(inv(MW))
    @test logdet(inverse(W)) ≈ logdet(Inverse(W))
    @test logdet(inverse(W)) ≈ logdet(inv(MW))
    @test all(logabsdet(inverse(W)) .≈ logabsdet(Inverse(W)))
    @test all(logabsdet(inverse(W)) .≈ logabsdet(inv(MW)))

    # indexing
    for i in 1:n, j in 1:n
        @test W[i, j] ≈ MW[i, j]
    end
    i = 1:2
    j = 3:3
    Wij = W[i, j]
    @test Wij isa Woodbury
    @test Wij[1, 1] == W[1, 3]
    Wij = @view W[i, j]
    @test Wij isa Woodbury
    @test Wij.U isa SubArray

    # trace
    @test tr(W) ≈ tr(MW)

    # factorize
    F = factorize(W)
    @test Matrix(inverse(F)) ≈ inv(MW)
    @test inv(W) ≈ inv(MW)
    @test logdet(F) ≈ logdet(MW)
    @test isposdef(F)

    # factorize and logdet after factorization
    F = factorize(W)
    @test logabsdet(F)[1] ≈ logabsdet(MW)[1]
    @test logabsdet(F)[2] ≈ logabsdet(MW)[2]
    @test logdet(F) ≈ logdet(MW)
    @test Matrix(inverse(F)) ≈ inv(MW)

    # test recursive stacking of woodbury objects
    b = randn(n, 1)
    W2 = Woodbury(W, b, ones(1, 1), b')
    F2 = factorize(W2)
    M2 = Matrix(W2)
    @test logdet(W2) ≈ logdet(M2)
    @test logdet(F2) ≈ logdet(M2)

    # test factorize_D behavior on non-hermitian 1x1 matrix
    using WoodburyFactorizations: factorize_D
    FD = factorize_D(W2, fill(-rand(), (1, 1)))
    @test FD isa Matrix
    @test size(FD) == (1, 1)

    # test rank 1 constructor
    x = randn(n)
    A = Diagonal(2 .+ exp.(randn(n)))
    u = randn(n)
    W = Woodbury(A, u)
    @test Matrix(W) ≈ A + u*u'
    b = W*x
    @test W \ b ≈ x

    v = randn(n)
    W = Woodbury(A, u, v')
    @test Matrix(W) ≈ A + u*v'
    b = W*x
    @test W \ b ≈ x

    # LowRank constructor
    U = randn(n, 1)
    W = Woodbury(I(n), U, U')

    # CholeskyPivoted constructor
    A = randn(1, n)
    A = A'A
    C = cholesky(A, RowMaximum(), check = false)
    W = Woodbury(A, C)
    @test W isa Woodbury
    @test Matrix(W) ≈ 2A
    # tests with α = -
    W = Woodbury(I(n), C, -1)
    # display(W)
    @test Matrix(W) ≈ I(n) - A
    @test inv(W) ≈ inv(I(n) - A)

    FW = factorize(W)
    @test FW isa Inverse
    @test Matrix(FW) ≈ I(n) - A
    @test Matrix(inverse(FW)) ≈ inv(I(n) - A)
    MW = I(n) - A
    # indexing
    for i in 1:n, j in 1:n
        @test W[i, j] ≈ MW[i, j]
    end
    @test tr(W) ≈ tr(MW)
    @test diag(W) ≈ diag(MW)

    # Woodbury structure with 1d outer dimension
    a, b, c, d = randn(4)
    W = Woodbury(a, b, c, d)
    MW = Matrix(W)
    @test size(MW) == (1, 1)
    @test MW[1, 1] ≈ a + *(b, c, d)

    a = randn()
    B = randn(1, n)
    C = randn(n, n)
    D = randn(n, 1)
    W = Woodbury(a, B, C, D)
    MW = Matrix(W)
    @test size(MW) == (1, 1)
    @test MW[1, 1] ≈ a + dot(B, C, D)

    @test inv(W) isa Real
    @test factorize(W) isa Real
end

end # TestWoodburyFactorizations


text:
## Exports
module Expect
export ExpectProc, expect!, sendline
export ExpectTimeout, ExpectEOF

## Support lib
const libttymakeraw = Pkg.dir("Expect", "deps", "libttymakeraw.so")

## Imports
import Base: AsyncStream, Process, TTY, wait_readnb
import Base: eof, read, readbytes!, readuntil, write, close

## Types
type ExpectTimeout <: Exception end
type ExpectEOF <: Exception end

type ExpectProc <: IO
    proc::Process
    timeout::Real
    codec::Function
    in_stream::AsyncStream
    out_stream::AsyncStream
    before
    match
    buffer::Vector{Uint8}

    function ExpectProc(cmd::Cmd, timeout::Real; env::Base.EnvHash=ENV, codec::Function=utf8)
        in_stream, out_stream, proc = _spawn(cmd, env)
        new(proc, timeout, codec, in_stream, out_stream, "", nothing, [])
    end
end

## Support functions
function raw!(tty::TTY, raw::Bool)
    # TODO: Base.Terminals.raw! does not currently set the correct line discipline
    #       See https://github.com/JuliaLang/libuv/pull/27
    #       We use our custom little hack in order to avoid waiting for libuv.
    ccall((:tty_makeraw, libttymakeraw), Int32, (Ptr{Void}, Int32), tty.handle, (raw? 1: 0))
end

raw!(proc::ExpectProc, raw::Bool) = raw!(proc.in_stream, raw)

function _spawn(cmd::Cmd, env::Base.EnvHash=ENV)
    env = copy(ENV)
    env["TERM"] = "dumb"
    setenv(cmd, env)
    detach(cmd)

    @unix? begin
        const O_RDWR = Base.FS.JL_O_RDWR
        const O_NOCTTY = Base.FS.JL_O_NOCTTY

        fdm = RawFD(ccall(:posix_openpt, Cint, (Cint,), O_RDWR|O_NOCTTY))
        fdm == -1 && error("openpt failed: $(strerror())")

        rc = ccall(:grantpt, Cint, (Cint,), fdm)
        rc != 0 && error("grantpt failed: $(strerror())")

        rc = ccall(:unlockpt, Cint, (Cint,), fdm)
        rc != 0 && error("unlockpt failed: $(strerror())")

        pts = ccall(:ptsname, Ptr{Uint8}, (Cint,), fdm)
        fds = RawFD(ccall(:open, Cint, (Ptr{Uint8}, Cint), pts, O_RDWR|O_NOCTTY))
        fds == -1 && error("open failed: $(strerror())")

        ttym = TTY(fdm; readable=true)
        in_stream = out_stream = ttym
        raw!(ttym, true) != 0 && error("raw! failed: $(strerror())")

        close_fds = (_...)->ccall(:close, Cint, (Cint,), fds)
        proc = spawn(true, cmd, (fds, fds, fds), close_fds)
    end : begin
        in_stream, out_stream, proc = readandwrite(cmd)
    end

    Base.start_reading(in_stream)
    return (in_stream, out_stream, proc)
end

# Base IO functions
eof(proc::ExpectProc) = eof(proc.in_stream)
close(proc::ExpectProc) = close(proc.out_stream)

write(proc::ExpectProc, buf::AbstractArray{Uint8}) = write(proc.out_stream, buf)
write(proc::ExpectProc, buf::Uint8) = write(proc.out_stream, buf)

function read(proc::ExpectProc, ::Type{Uint8})
    proc.buffer = []
    proc.before = nothing
    read(proc.in_stream, Uint8)
end

function readbytes!(proc::ExpectProc, b::AbstractArray{Uint8}, nb=length(b))
    proc.buffer = []
    proc.before = nothing
    readbytes!(proc.in_stream, b, nb)
end

function _readuntil(proc::ExpectProc, delim)
    proc.buffer = []
    proc.before = nothing
    readuntil(proc.in_stream, delim)
end

readuntil(proc::ExpectProc, delim::AbstractArray{Uint8}) = _readuntil(proc, delim)
readuntil(proc::ExpectProc, delim::String) = _readuntil(proc, delim)
readuntil(proc::ExpectProc, delim::Uint8) = _readuntil(proc, delim)

# Some helpers
sendline(proc::ExpectProc, line::String) = write(proc, string(line, "\n"))

# Expect
function _expect_search(buf::String, str::String)
    pos = search(buf, str)
    return pos == 0:-1? nothing: (buf[pos], pos)
end

function _expect_search(buf::String, regex::Regex)
    m = match(regex, buf)
    return m == nothing? nothing: (m.match, m.offset:(m.offset+length(m.match)-1))
end

function _expect_search(buf::String, vec::Vector)
    for idx=1:length(vec)
        ret = _expect_search(buf, vec[idx])
        if ret != nothing
            return idx, ret[1], ret[2]
        end
    end
    return nothing
end

function expect!(proc::ExpectProc, vec; timeout::Real=proc.timeout)
    pos = 0:-1
    idx = 0
    while true
        if nb_available(proc.in_stream) > 0
            proc.buffer = vcat(proc.buffer, takebuf_array(proc.in_stream.buffer))
        end
        if length(proc.buffer) > 0
            buffer = try proc.codec(proc.buffer) end
            if buffer != nothing
                ret = _expect_search(buffer, vec)
                if ret != nothing
                    idx, proc.match, pos = ret
                    break
                end
            end
        end
        if nb_available(proc.in_stream) == 0 && (!isopen(proc.in_stream) || process_exited(proc.proc))
            throw(ExpectEOF())
        end
        cond = Condition()
        @schedule begin
            wait_readnb(proc.in_stream, 1)
            notify(cond, true)
        end
        @schedule begin
            sleep(timeout)
            notify(cond, false)
        end
        if wait(cond) == false
            throw(ExpectTimeout())
        end
    end
    proc.before = proc.codec(proc.buffer[1:pos[1]-1])
    proc.buffer = proc.buffer[pos[end]+1:end]
    return idx
end

function expect!(proc::ExpectProc, regex::Regex; timeout::Real=proc.timeout)
    expect!(proc, [regex]; timeout=timeout)
    proc.before
end

function expect!(proc::ExpectProc, str::String; timeout::Real=proc.timeout)
    # TODO: this is worth implementing more efficiently
    expect!(proc, [str]; timeout=timeout)
    proc.before
end

end


text:
"""
    struct TargetInfo <: Any
Hold a pointer to a `clang::TargetInfo` object.
"""
struct TargetInfo
    ptr::CXTargetInfo_
end

Base.unsafe_convert(::Type{CXTargetInfo_}, x::TargetInfo) = x.ptr
Base.cconvert(::Type{CXTargetInfo_}, x::TargetInfo) = x


text:
using Random

@testset "Clgp" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  @testset "class numbers" begin
    @testset "quadratic fields" begin
      classnumbersofquadraticfields = Tuple{Int, Int}[(-50,1),(-49,1),(-48,1),(-47,5),(-46,4)
        ,(-45,2),(-44,1),(-43,1),(-42,4),(-41,8),(-40,2),(-39,4),(-38,6)
        ,(-37,2),(-36,1),(-35,2),(-34,4),(-33,4),(-32,1),(-31,3),(-30,4)
        ,(-29,6),(-28,1),(-27,1),(-26,6),(-25,1),(-24,2),(-23,3),(-22,2)
        ,(-21,4),(-20,2),(-19,1),(-18,1),(-17,4),(-16,1),(-15,2),(-14,4)
        ,(-13,2),(-12,1),(-11,1),(-10,2),(-9,1),(-8,1),(-7,1),(-6,2),(-5,2)
        ,(-4,1),(-3,1),(-2,1),(-1,1),(2,1),(3,1),(5,1),(6,1),(7,1),(8,1)
        ,(10,2),(11,1),(12,1),(13,1),(14,1),(15,2),(17,1),(18,1),(19,1),(20,1)
        ,(21,1),(22,1),(23,1),(24,1),(26,2),(27,1),(28,1),(29,1),(30,2),(31,1)
        ,(32,1),(33,1),(34,2),(35,2),(37,1),(38,1),(39,2),(40,2),(41,1),(42,2)
        ,(43,1),(44,1),(45,1),(46,1),(47,1),(48,1),(50,1)]

      @testset "K = Q[√$d]" for (d,h) in classnumbersofquadraticfields
        K, a = NumberField(x^2-d, "a")
        O = maximal_order(K)
        Cl, mCl = Hecke.class_group(O, redo = true)
        U, mU = Hecke.unit_group(O)
        @test order(Cl) == h

        Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
        U, mU = Hecke.unit_group(O)
        @test order(Cl) == h

        O = Order(K, shuffle(basis(O)), isbasis = true)
        O.ismaximal = 1
        Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
        U, mU = Hecke.unit_group(O)
        @test order(Cl) == h
      end
    end

    @testset "x^3-3*x-1" begin
      f = x^3 - 3*x - 1
      K, a = NumberField(f, "a")
      O = maximal_order(K)
      Cl, mCl = @inferred Hecke.class_group(O, redo = true)
      U, mU = @inferred Hecke.unit_group(O)
      @test order(Cl) == 1

      Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
      U, mU = Hecke.unit_group(O)
      @test order(Cl) == 1

      O = Order(K, shuffle(basis(O)), isbasis = true)
      O.ismaximal = 1
      Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
      U, mU = Hecke.unit_group(O)
      @test order(Cl) == 1
    end

    @testset "29th cyclotomic polynomial" begin
      phi29 = x^28+x^27+x^26+x^25+x^24+x^23+x^22+x^21+x^20+x^19+x^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1
      K, a = NumberField(phi29, "a")
      O = maximal_order(K)
      Cl, mCl = Hecke.class_group(O)
      @test order(Cl) == 8

      O = Order(K, shuffle(basis(O)), isbasis = true)
      O.ismaximal = 1
      Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
      U, mU = Hecke.unit_group(O)
      @test order(Cl) == 8
    end

    @testset "53th cyclotomic polynomial" begin
      #= phi53 = x^52+x^51+x^50+x^49+x^48+x^47+x^46+x^45+x^44+x^43+x^42+x^41+x^40+x^39+x^38+x^37+x^36+x^35+x^34+x^33+x^32+x^31+x^30+x^29+x^28+x^27+x^26+x^25+x^24+x^23+x^22+x^21+x^20+x^19+x^18+x^17+x^16+x^15+x^14+x^13+x^12+x^11+x^10+x^9+x^8+x^7+x^6+x^5+x^4+x^3+x^2+x+1 =#
      #= K, a = NumberField(phi53, "a") =#
      #= O = maximal_order(K) =#
      #= Cl = Hecke.class_group(O) =#
      #= @test Cl.h == 4889 =#
    end

    @testset "121th cyclotomic polynomial" begin
      #= phi121 = x^110+x^99+x^88+x^77+x^66+x^55+x^44+x^33+x^22+x^11+1 =#
      #= K, a = NumberField(phi121, "a") =#
      #= O = maximal_order(K) =#
      #= Cl = Hecke.class_group(O) =#
      #= @test Cl.h == 12188792628211 =#
    end

  end

end

@testset "_class_unit_group" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  AF = ArbField(20, cached = false)

  @testset "K = Q" begin
    K, a = NumberField(x, "a", cached = false)
    O = maximal_order(K)

    Cl, mCl = Hecke.class_group(O, redo = true)
    U, mU = Hecke.unit_group(O)
    @test order(Cl) == 1
  end

  @testset "K = Q[√2]" begin
    K, a = NumberField(x^2-2, "a")
    O = maximal_order(K)

    Cl, mCl = Hecke.class_group(O, redo = true)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 2
    @test mG(G[1]) == O(-1)
    @test order(U) == O
    @test U.torsion_units_order == 2
    @test contains(AF(0.88137358701),U.tentative_regulator)
    @test order(Cl) == 1

    O = Order(K, shuffle(basis(O)), isbasis = true)
    O.ismaximal = 1
    Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 2
    @test mG(G[1]) == O(-1)
    @test order(U) == O
    @test U.torsion_units_order == 2
    @test contains(AF(0.88137358701),U.tentative_regulator)
    @test order(Cl) == 1
  end

  @testset "K = Q[x]/(f), f = x^3 - 2" begin
    K, a = NumberField(x^2 - 3, "a")
    O = maximal_order(K)
    Cl, mCl = Hecke.class_group(O, redo = true)
    UU, mU = Hecke.unit_group(O)

    @test order(Cl) == 1
  end

  @testset "f = Q[x]/(f), f = x^5 - 11^2 * 7" begin
    K, a = NumberField(x^5 - 11^2 * 7, "a")
    O = maximal_order(K)
    Cl, mCl = Hecke.class_group(O, redo = true)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 2
    @test mG(G[1]) == O(-1)
    @test order(U) == O
    @test U.torsion_units_order == 2
    @test contains(AF(2027.9289425180057),U.tentative_regulator)
    @test order(Cl) == 5

    O = Order(K, shuffle(basis(O)), isbasis = true)
    O.ismaximal = 1
    Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 2
    @test mG(G[1]) == O(-1)
    @test order(U) == O
    @test U.torsion_units_order == 2
    @test contains(AF(2027.9289425180057),U.tentative_regulator)
    @test order(Cl) == 5
  end

  @testset "Cyclotomic Field 13" begin
    K, a = cyclotomic_field(13, cached = false)
    O = maximal_order(K)
    Cl, mCl = Hecke.class_group(O, redo = true)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 26
    @test order(U) == O
    @test U.torsion_units_order == 26
    @test order(Cl) == 1

    O = Order(K, shuffle(basis(O)), isbasis = true, cached = false)
    O.ismaximal = 1
    Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
    UU, mU = Hecke.unit_group(O)

    U = Hecke._get_UnitGrpCtx_of_order(O)

    G, mG = torsion_unit_group(O)
    @test order(G) == 26
    @test order(U) == O
    @test U.torsion_units_order == 26
    @test order(Cl) == 1
  end

  @testset "f = Q[x]/(f), f = x^18 + 18*x^16 + 135*x^14 + 192*x^12 - 2961*x^10 - 17334*x^8+ 20361*x^6 +  315108*x^4 + 514944*x^2 + 123904" begin
    K, a = NumberField(x^18 + 18*x^16 + 135*x^14 + 192*x^12 - 2961*x^10 - 17334*x^8+ 20361*x^6 +  315108*x^4 + 514944*x^2 + 123904, "a")
    O = maximal_order(K)

    Cl, mCl = Hecke.class_group(O, redo = true)
    UU, mU = Hecke.unit_group(O)

    @test order(Cl)== 36

    O = Order(K, shuffle(basis(O)), isbasis = true)
    O.ismaximal = 1

    Cl, mCl = Hecke.class_group(O, redo = true, do_lll = false)
    UU, mU = Hecke.unit_group(O)

    @test order(Cl)== 36
  end

  @testset "S3 field" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    f = x^6-24*x^4+157*x^2-162
    K, a = NumberField(f)
    OK = maximal_order(K)
    OK = lll(OK)
    C, mC = class_group(OK, redo = true, use_aut = true)
    @test order(C) == 3
  end

  @testset "Proof" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    f = x^2 - 3 * 5 * 7 * 11
    K, a = NumberField(f)
    OK = maximal_order(K)
    C, mC = class_group(OK, GRH = false)
    @test order(C) == 8
  end

  @testset "Auto but no autos" begin
    K, a = number_field(x^3 - 3*x - 4)
    OK = (maximal_order(K))
    c, mc = class_group(OK, use_aut = true)
    @test order(c) == 1
  end

  @testset "A class group which was hanging" begin
    K, a = number_field(x^4 - 4*x^3 - 1039*x^2 + 2086*x + 241020)
    OK = lll(maximal_order(K))
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
    c, mc = class_group(OK, redo = true)
  end

  @testset "A class group that was failing" begin
    K, a = number_field(x^24 + x^23 - x^19 - x^18 - x^17 - x^16 + x^14 + x^13 + x^12 + x^11 + x^10 - x^8 - x^7 - x^6 - x^5 + x + 1, cached = false, check = false)
    automorphisms(K)
    OK = maximal_order(K)
    @test order(class_group(OK)[1]) == 1
  end

  @testset "Class group proof" begin
    K, a = number_field(x^2 - 2)
    OK = maximal_order(K)
    c, mc = class_group(OK, GRH = false)
    @test isone(order(c))
  end

  @testset "Class group with small generating set" begin
    K, a = number_field(x^6 - 2*x^5 + 756*x^4 - 1006*x^3 + 192535*x^2 - 128528*x + 16515841)
    OK = maximal_order(K)
    for i in 1:20
      c, mc = class_group(OK, redo = true)
      @test order(c) == 892
    end
  end
end


text:
function update_and_get_max_abs_diff!(  fssrb::PotentialSimulationSetupRB{T, N1, N2},
                                        depletion_handling::Val{depletion_handling_enabled},
                                        only2d::Val{only_2d} = Val{false}(),
                                        is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                        use_nthreads::Int = Base.Threads.nthreads()
                                        )::T where {T, N1, N2, depletion_handling_enabled, only_2d, _is_weighting_potential}
    tmp_potential::Array{T, N2} = copy(fssrb.potential)
    for i in 1:10
        update!(fssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
    end
    max_diff::T = maximum(abs.(tmp_potential - fssrb.potential))
    return max_diff
end

function _update_till_convergence!( fssrb::PotentialSimulationSetupRB{T, N1, N2},
                                    convergence_limit::T;
                                    n_iterations_between_checks = 500,
                                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
                                    only2d::Val{only_2d} = Val{false}(),
                                    is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                    use_nthreads::Int = Base.Threads.nthreads(),
                                    max_n_iterations::Int = -1
                                    )::T where {T, N1, N2, depletion_handling_enabled, only_2d, _is_weighting_potential}
    n_iterations::Int = 0
    cf::T = Inf
    cfs::Vector{T} = fill(cf, 10)
    cl::T = _is_weighting_potential ? convergence_limit : abs(convergence_limit * fssrb.bias_voltage) # to get relative change in respect to bias voltage
    prog = ProgressThresh(cl, 0.1, "Convergence: ")
    while cf > cl
        for i in 1:n_iterations_between_checks
            update!(fssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
        end
        cf = update_and_get_max_abs_diff!(fssrb, depletion_handling, only2d, is_weighting_potential, use_nthreads)
        @inbounds cfs[1:end-1] = cfs[2:end]
        @inbounds cfs[end] = cf
        slope::T = abs(mean(diff(cfs)))
        ProgressMeter.update!(prog, cf)
        n_iterations += n_iterations_between_checks
        if slope < cl
            # @info "Slope is basically 0 -> Converged: $slope"
            cf = slope
        end
        if max_n_iterations > 0 && n_iterations > max_n_iterations
            # @show n_iterations_between_checks
            @info "Maximum number of iterations reached. (`n_iterations = $(n_iterations)`)"
            break
        end
    end
    if depletion_handling_enabled
        tmp_pointtypes::Array{PointType, N2} = fssrb.pointtypes .& undepleted_bit
        @showprogress "Checking undepleted regions " for i in 1:10
            update!(fssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
            @inbounds for i in eachindex(fssrb.pointtypes)
                if (fssrb.pointtypes[i] & undepleted_bit == 0) && (tmp_pointtypes[i] > 0)
                    fssrb.pointtypes[i] += undepleted_bit
                elseif (fssrb.pointtypes[i] & undepleted_bit > 0) && (tmp_pointtypes[i] == 0)
                    tmp_pointtypes[i] += undepleted_bit
                end
            end
        end
        @inbounds for i in eachindex(fssrb.pointtypes)
            if (fssrb.pointtypes[i] & update_bit == 0)
                fssrb.pointtypes[i] = PointType(0)
            else
                if (fssrb.pointtypes[i] & pn_junction_bit == 0)
                    if fssrb.pointtypes[i] & undepleted_bit > 0 fssrb.pointtypes[i] -= undepleted_bit end
                end
            end
        end
    end

    ProgressMeter.finish!(prog)
    return cf
end

function _get_refinement_inds( potential::Array{T, 3}, grid::Grid{T, 3, :cylindrical},
                max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T})::NTuple{3, Vector{Int}} where {T <: SSDFloat}
    inds_r::Vector{Int} = Int[]
    inds_φ::Vector{Int} = Int[]
    inds_z::Vector{Int} = Int[]

    ax1::DiscreteAxis{T} = grid[1]
    ax2::DiscreteAxis{T} = grid[2]
    ax3::DiscreteAxis{T} = grid[3]
    int1::Interval = ax1.interval
    int2::Interval = ax2.interval
    int3::Interval = ax3.interval

    refinement_value_r::T = max_diffs[1]
    refinement_value_φ::T = max_diffs[2]
    refinement_value_z::T = max_diffs[3]

    minimum_distance_r::T = minimum_distances[1]
    minimum_distance_φ::T = minimum_distances[2]
    minimum_distance_z::T = minimum_distances[3]

    minimum_distance_φ_deg::T = rad2deg(minimum_distance_φ)

    axr::Vector{T} = collect(grid.axes[1])
    axφ::Vector{T} = collect(grid.axes[2])
    axz::Vector{T} = collect(grid.axes[3])
    cyclic::T = grid.axes[2].interval.right - grid.axes[2].interval.left

    for iz in 1:size(potential, 3)
        z = axz[iz]
        for iφ in 1:size(potential, 2)
            φ = axφ[iφ]
            isum = iz + iφ
            for ir in 1:size(potential, 1)
                r = axr[ir]

                # r
                if ir != size(potential, 1)
                    if abs(potential[ir + 1, iφ, iz] - potential[ir, iφ, iz]) >= refinement_value_r
                        if axr[ir + 1] - axr[ir] > minimum_distance_r
                            if !in(ir, inds_r) push!(inds_r, ir) end
                        end
                    end
                end

                # φ
                if iφ == size(potential, 2)
                    if typeof(int2).parameters[2] == :open # ugly solution for now. should be done over multiple dispatch..
                        if abs(potential[ir, 1, iz] - potential[ir, iφ, iz]) >= refinement_value_φ
                            if cyclic - axφ[iφ] > minimum_distance_φ
                                if !in(iφ, inds_φ) push!(inds_φ, iφ) end
                            end
                        end
                    end
                else
                    if abs(potential[ir, iφ + 1, iz] - potential[ir, iφ, iz]) >= refinement_value_φ
                        if axφ[iφ + 1] - axφ[iφ] > minimum_distance_φ
                            if !in(iφ, inds_φ) push!(inds_φ, iφ) end
                        end
                    end
                end

                # z
                if iz != size(potential, 3)
                    if abs(potential[ir, iφ, iz + 1] - potential[ir, iφ, iz]) >= refinement_value_z
                        if axz[iz + 1] - axz[iz] > minimum_distance_z
                            if !in(iz, inds_z) push!(inds_z, iz) end
                        end
                    end
                end

            end
        end
    end

    if typeof(int2).parameters[2] == :open
        if isodd(length(inds_φ))
            for iφ in 1:size(potential, 2)
                if !in(iφ, inds_φ)
                    if iφ != size(potential, 2)
                        if axφ[iφ + 1] - axφ[iφ] > minimum_distance_φ
                            append!(inds_φ, iφ)
                            break
                        end
                    else
                        if cyclic - axφ[iφ] > minimum_distance_φ
                            append!(inds_φ, iφ)
                            break
                        end
                    end
                end
            end
        end
        if isodd(length(inds_φ))
            b::Bool = true
            diffs = diff(axφ)
            while b
                if length(diffs) > 0
                    idx = findmax(diffs)[2]
                    if !in(idx, inds_φ)
                        append!(inds_φ, idx)
                        b = false
                    else
                        deleteat!(diffs, idx)
                    end
                else
                    for i in eachindex(axφ)
                        if !in(i, inds_φ) push!(inds_φ, i) end
                    end
                    b = false
                end
            end
        end
        @assert iseven(length(inds_φ)) "Refinement would result in uneven grid in φ."
    end

    if isodd(length(inds_z))
        for iz in 1:size(potential, 3)
            if !in(iz, inds_z)
                if iz != size(potential, 3)
                    if axz[iz + 1] - axz[iz] > minimum_distance_z
                        append!(inds_z, iz)
                        break
                    end
                end
            end
        end
    end
    if isodd(length(inds_z)) && ((length(inds_z) + 1) < length(axz))
        b = true
        diffs = diff(axz)
        while b
            if length(diffs) > 0
                idx = findmax(diffs)[2]
                if !in(idx, inds_z)
                    append!(inds_z, idx)
                    b = false
                else
                    deleteat!(diffs, idx)
                end
            else
                for i in eachindex(axz[1:end-1])
                    if !in(i, inds_z) push!(inds_z, i) end
                end
                b = false
            end
        end
    end
    if isodd(length(inds_z)) && (length(inds_z) > 1)
        b = true
        diffs = diff(axz)
        while b && length(diffs) > 1
            idx = findmin(diffs)[2]
            if in(idx, inds_z)
                deleteat!(inds_z, findfirst(i -> i == idx, inds_z))
                b = false
            else
                deleteat!(diffs, idx)
            end
        end
    end
    if isodd(length(inds_z)) inds_z = inds_z[1:end-1] end
    @assert iseven(length(inds_z)) "Refinement would result in uneven grid in z."

    return sort(inds_r), sort(inds_φ), sort(inds_z)
end

function _get_refinement_inds(  potential::Array{T, 3}, grid::Grid{T, 3, :cartesian},
                                max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T})::NTuple{3, Vector{Int}} where {T}
    inds_x::Vector{Int} = Int[]
    inds_y::Vector{Int} = Int[]
    inds_z::Vector{Int} = Int[]

    ax1::DiscreteAxis{T} = grid[1]
    ax2::DiscreteAxis{T} = grid[2]
    ax3::DiscreteAxis{T} = grid[3]
    int1::Interval = ax1.interval
    int2::Interval = ax2.interval
    int3::Interval = ax3.interval

    refinement_value_x::T = max_diffs[1]
    refinement_value_y::T = max_diffs[2]
    refinement_value_z::T = max_diffs[3]

    minimum_distance_x::T = minimum_distances[1]
    minimum_distance_y::T = minimum_distances[2]
    minimum_distance_z::T = minimum_distances[3]

    ax_x::Vector{T} = collect(grid.axes[1])
    ax_y::Vector{T} = collect(grid.axes[2])
    ax_z::Vector{T} = collect(grid.axes[3])

    for iz in 1:size(potential, 3)
        z = ax_z[iz]
        for iy in 1:size(potential, 2)
            y = ax_y[iy]
            isum = iz + iy
            for ix in 1:size(potential, 1)
                x = ax_x[ix]

                # x
                if ix != size(potential, 1)
                    if abs(potential[ix + 1, iy, iz] - potential[ix, iy, iz]) >= refinement_value_x
                        if ax_x[ix + 1] - ax_x[ix] > minimum_distance_x
                            if !in(ix, inds_x) push!(inds_x, ix) end
                        end
                    end
                end

                # y
                if iy != size(potential, 2)
                    if abs(potential[ix, iy + 1, iz] - potential[ix, iy, iz]) >= refinement_value_y
                        if ax_y[iy + 1] - ax_y[iy] > minimum_distance_y
                            if !in(iy, inds_y) push!(inds_y, iy) end
                        end
                    end
                end

                # z
                if iz != size(potential, 3)
                    if abs(potential[ix, iy, iz + 1] - potential[ix, iy, iz]) >= refinement_value_z
                        if ax_z[iz + 1] - ax_z[iz] > minimum_distance_z
                            if !in(iz, inds_z) push!(inds_z, iz) end
                        end
                    end
                end

            end
        end
    end

    if isodd(length(inds_x))
        for ix in 1:size(potential, 1)
            if !in(ix, inds_x)
                if ix != size(potential, 1)
                    if ax_x[ix + 1] - ax_x[ix] > minimum_distance_x
                        append!(inds_x, ix)
                        break
                    end
                end
            end
        end
    end
    if isodd(length(inds_x))
        b = true
        diffs = diff(ax_x)
        while b
            if length(diffs) > 0
                idx = findmax(diffs)[2]
                if !in(idx, inds_x)
                    append!(inds_x, idx)
                    b = false
                else
                    deleteat!(diffs, idx)
                end
            else
                for i in eachindex(ax_x[1:end-1])
                    if !in(i, inds_x) push!(inds_x, i) end
                end
                b = false
            end
        end
    end
    if isodd(length(inds_x)) inds_x = inds_x[1:end-1] end
    @assert iseven(length(inds_x)) "Refinement would result in uneven grid in x. This is not allowed since this is the red black dimension."

    return sort(inds_x), sort(inds_y), sort(inds_z)
end

function refine(p::ScalarPotential, max_diffs::Tuple{<:Real,<:Real,<:Real}, minimum_distances::Tuple{<:Real,<:Real,<:Real})::typeof(p)
    T = eltype(p.grid)
    refine_at_inds = _get_refinement_inds(p.data, p.grid, T.(max_diffs), T.(minimum_distances))
    potential, grid = add_points_and_interpolate(p.data, p.grid, refine_at_inds...)
    return typeof(p)(potential, grid)
end


text:
### Operator Implementation
mutable struct JacVecOperator{T,F,T1,T2,uType,P,tType,O} <:
               DiffEqBase.AbstractDiffEqLinearOperator{T}
    f::F
    cache1::T1
    cache2::T2
    u::uType
    p::P
    t::tType
    autodiff::Bool
    ishermitian::Bool
    opnorm::O

    function JacVecOperator{T}(
        f,
        p = nothing,
        t::Union{Nothing,Number} = nothing;
        autodiff = true,
        ishermitian = false,
        opnorm = true,
    ) where {T}
        p === nothing ? P = Any : P = typeof(p)
        t === nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),Nothing,Nothing,Any,P,tType,typeof(opnorm)}(
            f,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            autodiff,
            ishermitian,
        )
    end
    function JacVecOperator{T}(
        f,
        u::AbstractArray,
        p = nothing,
        t::Union{Nothing,Number} = nothing;
        autodiff = true,
        ishermitian = false,
        opnorm = true,
    ) where {T}
        if autodiff
            cache1 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(u, u)
            cache2 = ForwardDiff.Dual{SparseDiffTools.DeivVecTag}.(u, u)
        else
            cache1 = similar(u)
            cache2 = similar(u)
        end
        p === nothing ? P = Any : P = typeof(p)
        t === nothing ? tType = Any : tType = typeof(t)
        new{
            T,
            typeof(f),
            typeof(cache1),
            typeof(cache2),
            typeof(u),
            P,
            tType,
            typeof(opnorm),
        }(
            f,
            cache1,
            cache2,
            u,
            p,
            t,
            autodiff,
            ishermitian,
            opnorm,
        )
    end
    function JacVecOperator(f, u, args...; kwargs...)
        JacVecOperator{eltype(u)}(f, u, args...; kwargs...)
    end
end

LinearAlgebra.opnorm(L::JacVecOperator, p::Real = 2) = L.opnorm
LinearAlgebra.ishermitian(L::JacVecOperator) = L.ishermitian

Base.size(L::JacVecOperator) = (length(L.cache1), length(L.cache1))
Base.size(L::JacVecOperator, i::Int) = length(L.cache1)
function update_coefficients!(L::JacVecOperator, u, p, t)
    L.u = u
    L.p = p
    L.t = t
    !L.autodiff && L.cache1 !== nothing && L.f(L.cache1, L.u, L.p, L.t)
end

# Interpret the call as df/du * u

function (L::JacVecOperator)(u, p, t::Number)
    update_coefficients!(L, u, p, t)
    L * u
end

function (L::JacVecOperator)(du, u, p, t::Number)
    update_coefficients!(L, u, p, t)
    mul!(du, L, u)
end

function Base.:*(L::JacVecOperator, x::AbstractVector)
    if hasmethod(L.f, typeof.((L.u, L.p, L.t)))
        return vec(
            L.autodiff ? auto_jacvec(_u -> L.f(_u, L.p, L.t), L.u, x) :
            num_jacvec(_u -> L.f(_u, L.p, L.t), L.u, x),
        )
    end
    return mul!(similar(vec(L.u)), L, x)
end

function LinearAlgebra.mul!(
    du::AbstractVector,
    L::JacVecOperator,
    x::AbstractVector,
)
    du = reshape(du, size(L.u))
    let p = L.p, t = L.t
        if L.cache1 === nothing
            if L.autodiff
                auto_jacvec!(du, (_du, _u) -> L.f(_du, _u, p, t), L.u, x)
            else
                num_jacvec!(
                    du,
                    (_du, _u) -> L.f(_du, _u, p, t),
                    L.u,
                    x;
                    compute_f0 = true,
                )
            end
        else
            if L.autodiff
                auto_jacvec!(
                    du,
                    (_du, _u) -> L.f(_du, _u, p, t),
                    L.u,
                    x,
                    L.cache1,
                    L.cache2,
                )
            else
                num_jacvec!(
                    du,
                    (_du, _u) -> L.f(_du, _u, p, t),
                    L.u,
                    x,
                    L.cache1,
                    L.cache2;
                    compute_f0 = true,
                )
            end
        end
    end
    return vec(du)
end

function Base.resize!(J::JacVecOperator, i)
    resize!(J.cache1, i)
    resize!(J.cache2, i)
end

### AnalyticalOperator Implementation

mutable struct AnalyticalJacVecOperator{T,F,uType,P,tType,O} <:
               DiffEqBase.AbstractDiffEqLinearOperator{T}
    f::F
    u::uType
    p::P
    t::tType
    ishermitian::Bool
    opnorm::O
    function AnalyticalJacVecOperator{T}(
        f,
        u = nothing,
        p = nothing,
        t::Union{Nothing,Number} = nothing;
        ishermitian = false,
        opnorm = true,
    ) where {T}
        u === nothing ? uType = Any : uType = typeof(u)
        p === nothing ? P = Any : P = typeof(p)
        t === nothing ? tType = Any : tType = typeof(t)
        new{T,typeof(f),uType,P,tType,typeof(opnorm)}(
            f,
            u,
            p,
            t,
            ishermitian,
            opnorm,
        )
    end
    function AnalyticalJacVecOperator(f, u, args...; kwargs...)
        AnalyticalJacVecOperator{eltype(u)}(f, u, args...; kwargs...)
    end
end

LinearAlgebra.opnorm(L::AnalyticalJacVecOperator, p::Real = 2) = L.opnorm
LinearAlgebra.ishermitian(L::AnalyticalJacVecOperator) = L.ishermitian

Base.size(L::AnalyticalJacVecOperator) = (length(L.u), length(L.u))
Base.size(L::AnalyticalJacVecOperator, i::Int) = length(L.u)
function update_coefficients!(L::AnalyticalJacVecOperator, u, p, t)
    L.u = u
    L.p = p
    L.t = t
    return L
end

# Interpret the call as df/du * u

function (L::AnalyticalJacVecOperator)(u, p, t::Number)
    update_coefficients!(L, u, p, t)
    L * u
end

function (L::AnalyticalJacVecOperator)(du, u, p, t::Number)
    update_coefficients!(L, u, p, t)
    mul!(du, L, u)
end

Base.:*(L::AnalyticalJacVecOperator, x::AbstractVector) = L.f(x, L.u, L.p, L.t)

function LinearAlgebra.mul!(
    du::AbstractVector,
    L::AnalyticalJacVecOperator,
    x::AbstractVector,
)
    L.f(du, x, L.u, L.p, L.t)
end


text:
using DataStructures: OrderedDict
using MLGoodies

@testset "toposort" begin
    nodes = OrderedDict(2 => [1], 5 => [1, 3, 4], 4 => [1, 2], 6 => [5], 1 => [], 3 => [1], 20 => [5, 10], 10 => [1])
    res = toposort(nodes, v -> last(v), v -> first(v))
    @test res == [1, 2, 3, 4, 5, 6, 10, 20]
end

max_stars_repo_name: dhonza/MLGoodies


text:
using StatsBase
using ClusterValidityIndices
using Logging
using Plots
using MetaICVI
using Random

# Plotting options
dpi = 300       # Plotting dots-per-inch
theme(:dark)    # Plotting style
# gr()            # GR backend (default for Plots.jl)
unicodeplots()

# Include the library definitions
# include(projectdir("julia/lib_sim.jl"))
# include("rocket.jl")

# Set the logging level to Info and standardize the random seed
LogLevel(Logging.Info)
Random.seed!(0)
# data_dir(args...) = projectdir("work/data/meta_icvi", args...)
# results_dir(args...) = projectdir("work/results/meta_icvi", args...)
data_dir(args...) = joinpath("../data/training", args...)
results_dir(args...) = joinpath("../data/results", args...)

# Setup the data
data_paths = [
    data_dir("correct_partition.csv"),
    data_dir("over_partition.csv"),
    data_dir("under_partition.csv")
]
n_samples = 2000
n_data = 3

# Setup ICVIs
cvis = [
    cSIL(),
    PS()
    # PS(),
    # GD43()
]
n_cvis = length(cvis)

# Setup the correlation
n = 100                         # Window size
n_corr = n_samples - n + 1      # Number of correlations
corrs = zeros(n_corr, n_data)   # Preallocation of correlations
iter_corr = n:n_samples         # Iterator for the correlation generation

plot_labels = zeros(n_samples, n_data)

# Iterate over all partition types
for dx = 1:n_data
    data_path = data_paths[dx]

    # Load the training data
    data, labels = get_cvi_data(data_path)

    plot_labels[:, dx] = labels

    labels = relabel_cvi_data(labels)
    data_name = splitext(basename(data_path))[1]    # Get the data_name

    # Instantiate the icvis
    local_cvis = deepcopy(cvis)
    criterion_values = zeros(n_samples, n_cvis)

    # Run the ICVIs iteratively
    for i = 1:n_samples
        sample = data[:, i]
        label = labels[i]
        for j = 1:n_cvis
            criterion_values[i, j] = get_icvi!(local_cvis[j], sample, label)
        end
    end

    # Calculate the windowed correlations
    for i = iter_corr
        corrs[i-n+1, dx] = corspearman(criterion_values[i-n+1:i, 1], criterion_values[i-n+1:i, 2])
    end
end

n_window = 100
n_kernels = 2
n_rocket = n_corr - n_window + 1
# corrs = zeros(n_corr, n_data)   # Preallocation of correlations
rocket_features = zeros(n_rocket, n_kernels, 2, n_data)
# rocket_features = []
# rocket_maxes = zeros(n_rocket, n_data)
iter_rocket = n + n_window:n_samples
@info n_rocket
@info iter_rocket

rocket = RocketModule(n_window, n_kernels)
for dx = 1:n_data
    # for i = iter_rocket
    for i = 1:size(corrs)[1]
        if i >= n_window
            # sample_window = corrs[i - n_window + 1: i, dx]
            sample_window = corrs[i - n_window + 1: i, dx]
            # push!(rocket_features, apply_kernels(rocket, sample_window))
            # rocket_features[i-n_window-n+1, :, :, dx] = apply_kernels(rocket, sample_window)
            # index, ppv, max, data
            rocket_features[i - n_window + 1, :, :, dx] = apply_kernels(rocket, sample_window)
        end
    end
end

# Plot
p = plot(dpi=dpi, legend=:bottomleft)
for i = 1:n_data
    plot!(p, iter_corr, corrs[:, i], label=basename(data_paths[i]))
end
# plot!(p, 1:n_samples, labels, label="Labels")
# plot!(p, 1:n_samples, data.train_y, label="True Labels")
# plot!(p, 1:n_corr, criterion_values_p, label="Porcelain")
# title!("Spearman, Data: " * basename(data_path))
title!("Partition ICVI Correlations")
xlabel!("Sample Index")
ylabel!("Spearman")
# xlims!(1, n_corr)
# ylims!(0, Inf)

g = plot(dpi=dpi, legend=:bottomleft)
for i = 1:n_data
    plot!(g, plot_labels[:, i], label=basename(data_paths[i]))
end

h = plot(dpi=dpi, legend=:bottomleft)
for i = 1:n_data
    # for j = 1:size(rocket_features[i])[1]
    for j = 1:n_kernels
        # plot!(h, rocket_features[i][j, 1], label=basename(data_paths[i]))
        plot!(h, 199:2000, rocket_features[:, j, 2, i], label=basename(data_paths[i]))
    end
end

l = @layout [a; b; c]
# pt = plot(p, g, layout = l, size=(500,500))
pt = plot(p, g, h, layout = l, size=(800,500))

# Save and show the plot
# png(pt, results_dir("4_corr_rocket"))
display(pt)


text:
using Makie
include("database.jl")

include("tutorials.jl")
include("examples2d.jl")
include("examples3d.jl")
include("interactive.jl")
include("documentation.jl")
include("implicits.jl")
include("short_tests.jl")
include("recipes.jl")
include("bigdata.jl")
include("layouting.jl")

database


text:
#InspectDR: Basic math operations
#-------------------------------------------------------------------------------
#=NOTE:
These tools should eventually be moved to a separate unit.
=#

#==Useful tests
===============================================================================#

#Verifies that v is strictly increasing (no repeating values):
#TODO: support repeating values (non-strict)
function isincreasing(v::Vector)
	if length(v) < 1; return true; end
	prev = v[1]

	for x in v[2:end]
		if !(x > prev) #Make sure works for NaN
			return false
		end
	end
	return true
end

isincreasing(r::AbstractRange) = (step(r) > 0)

#==Basic operations
===============================================================================#

#Safe version of extrema (returns DNaN on error):
#Also ignores NaNs.  TODO: optimize?
function extrema_nan(v::Array)
	vmin = DInf; vmax = -DInf
	for elem in v
		if !isnan(elem)
			if elem > vmax
				vmax = elem
			end
			if elem < vmin
				vmin = elem
			end
		end
	end
	if vmin>vmax
		vmin = vmax = DNaN
	end
	return (vmin, vmax)
end

#Finds the finite extents of a vector.
#Returns infinite extents if none are found.
#TODO: optimize?
function extents_finite(v::Array)
	vmin = DInf; vmax = -DInf
	vmin_inf = DInf; vmax_inf = -DInf
	for elem in v
		if isfinite(elem)
			if elem > vmax
				vmax = elem
			end
			if elem < vmin
				vmin = elem
			end
		elseif isinf(elem)
			if elem > vmax_inf
				vmax_inf = elem
			end
			if elem < vmin_inf
				vmin_inf = elem
			end
		end
	end
	if vmin>vmax #Did not find finite values
		vmin = vmin_inf; vmax = vmax_inf

		if vmin>vmax #Did not even find infinite values
			vmin = vmax = DNaN
		end
	end
	return (vmin, vmax)
end

#Last line


text:
import ArgParse, RIOPA
import Test: @testset, @test, @test_throws

@testset "args" begin
    ArgError = ArgParse.ArgParseError

    inputs = RIOPA.parse_inputs(["hello"])
    @test inputs["%COMMAND%"] == "hello"

    @test_throws(
        ArgError("unknown command: hello-fail"),
        RIOPA.parse_inputs(
            ["hello-fail"],
            error_handler = ArgParse.debug_handler,
        )
    )

    @test_throws(
        ArgError("unrecognized option --hey"),
        RIOPA.parse_inputs(["--hey"], error_handler = ArgParse.debug_handler)
    )
end


text:
# # [Implicit Nonlinear Dynamics : Autoregulation](@id autoregulation)
#
# The following is another example on how to use the [`ImplicitOptimizer`](@ref) describing a biological autoregulation process using
# two coupled implicit equations.

using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra
#md using Plots
#md gr()

@parameters begin
    t
    α = 1.0
    β = 1.3
    γ = 2.0
    δ = 0.5
end

@variables begin
    x[1:2](t) = [20.0; 12.0]
end

x = collect(x)
D = Differential(t)

eqs = [
    D(x[1]) ~ α/(1+x[2])-β*x[1];
    D(x[2]) ~ γ/(1+x[1])-δ*x[2];
]

sys = ODESystem(eqs, t, x, [α, β, γ, δ], name = :Autoregulation)

x0 = [x[1] => 20.0; x[2] => 12.0]

tspan = (0.0, 5.0)

de_problem = ODEProblem(sys, x0, tspan)
de_solution = solve(de_problem, Tsit5(), saveat = 0.005);
#md plot(de_solution)

# As always, we start by defining a [DataDrivenProblem](@ref problem) and a sufficient basis for sparse regression.

dd_prob = ContinuousDataDrivenProblem(de_solution)

eqs = [
    polynomial_basis(x, 4); D.(x); x .* D(x[1]); x .* D(x[2])
    ]

basis = Basis(eqs, x, independent_variable = t, implicits = D.(x))

#md plot(dd_prob)

# Next to varying over different sparsity penalties, we want also to batch our data using the [DataSampler](@ref).
# We define a train-test split of 80-20 for our data and batch the resulting training data into 10 minibatches, allowing
# shuffled and repeated values. Our goal is to find the model with the minimal error.

sampler = DataSampler(
    Split(ratio = 0.8), Batcher(n = 10, shuffle = true, repeated = true, batchsize_min = 30)
)

res = solve(dd_prob, basis, ImplicitOptimizer(STLSQ(1e-1:1e-1:9e-1)), by = :min, sampler = sampler, digits = 1)
print(res) #hide
print(result(res)) #hide

# We have recovered the correct equations of motion!
# Another visual check using the problem and the result yields

system = result(res)
@named ode = ODESystem(equations(system), t, x, parameters(system));
ode_prob = ODEProblem(ode, x0, tspan, parameter_map(res));

prediction = solve(ode_prob, Tsit5(), saveat = 0.2);
#md plot(de_solution, label = ["Groundtruth" nothing])
#md scatter!(prediction, label = ["Prediction" nothing])

#md # ## [Copy-Pasteable Code](@id autoregulation_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

for r_ in [res] #src
    @test all(l2error(r_) .< 0.5) #src
    @test all(aic(r_) .> 1e3) #src
    @test all(determination(r_) .>= 0.9) #src
end #src


text:
@inline function initialize!(integrator,cache::ExplicitRKConstantCache,f=integrator.f)
  integrator.kshortsize = 2
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.fsalfirst = f(integrator.t,integrator.uprev)
end

@inline function perform_step!(integrator,cache::ExplicitRKConstantCache,f=integrator.f)
  @unpack t,dt,uprev,u = integrator
  @unpack A,c,α,αEEst,stages = cache
  @unpack kk = cache
  # Calc First
  kk[1] = integrator.fsalfirst
  # Calc Middle
  for i = 2:stages-1
    utilde = zero(kk[1])
    for j = 1:i-1
      utilde = @muladd utilde + A[j,i]*kk[j]
    end
    kk[i] = f(@muladd(t+c[i]*dt),@muladd(uprev+dt*utilde));
  end
  #Calc Last
  utilde = zero(kk[1])
  for j = 1:stages-1
    utilde = @muladd utilde + A[j,end]*kk[j]
  end
  kk[end] = f(@muladd(t+c[end]*dt),@muladd(uprev+dt*utilde)); integrator.fsallast = kk[end] # Uses fsallast as temp even if not fsal
  # Accumulate Result
  utilde = α[1]*kk[1]
  for i = 2:stages
    utilde = @muladd utilde + α[i]*kk[i]
  end
  u = @muladd uprev + dt*utilde
  if integrator.opts.adaptive
    uEEst = αEEst[1]*kk[1]
    for i = 2:stages
      uEEst = @muladd uEEst + αEEst[i]*kk[i]
    end
    integrator.EEst = integrator.opts.internalnorm( dt*(utilde-uEEst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
  end
  if isfsal(integrator.alg.tableau)
    integrator.fsallast = kk[end]
  else
    integrator.fsallast = f(t+dt,u)
  end
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  @pack integrator = t,dt,u
end

@inline function initialize!(integrator,cache::ExplicitRKCache,f=integrator.f)
  integrator.kshortsize = 2
  integrator.fsallast = cache.fsallast
  integrator.fsalfirst = cache.kk[1]
  integrator.k = eltype(integrator.sol.k)(integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::ExplicitRKCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack A,c,α,αEEst,stages = cache.tab
  @unpack kk,utilde,tmp,atmp,uEEst = cache
  # Middle
  for i = 2:stages-1
    for l in uidx
      utilde[l] = zero(kk[1][1])
    end
    for j = 1:i-1
      for l in uidx
        utilde[l] = @muladd utilde[l] + A[j,i]*kk[j][l]
      end
    end
    for l in uidx
      tmp[l] = @muladd uprev[l]+dt*utilde[l]
    end
    f(@muladd(t+c[i]*dt),tmp,kk[i])
  end
  #Last
  for l in uidx
    utilde[l] = zero(kk[1][1])
  end
  for j = 1:stages-1
    for l in uidx
      utilde[l] = @muladd utilde[l] + A[j,end]*kk[j][l]
    end
  end
  for l in uidx
    u[l] = @muladd uprev[l]+dt*utilde[l]
  end
  f(@muladd(t+c[end]*dt),u,kk[end]) #fsallast is tmp even if not fsal
  #Accumulate
  if !isfsal(integrator.alg.tableau)
    for i in uidx
      utilde[i] = α[1]*kk[1][i]
    end
    for i = 2:stages
      for l in uidx
        utilde[l] = @muladd utilde[l] + α[i]*kk[i][l]
      end
    end
    for i in uidx
      u[i] = @muladd uprev[i] + dt*utilde[i]
    end
  end
  if integrator.opts.adaptive
    for i in uidx
      uEEst[i] = αEEst[1]*kk[1][i]
    end
    for i = 2:stages
      for j in uidx
        uEEst[j] = @muladd uEEst[j] + αEEst[i]*kk[i][j]
      end
    end
    for i in uidx
      atmp[i] = (dt*(utilde[i]-uEEst[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
    end
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  if !isfsal(integrator.alg.tableau)
    f(t+dt,u,integrator.fsallast)
  end
  @pack integrator = t,dt,u
end


text:
"""
`songlength{T}(syllables::Array{Syllable{T}, 1}, windowStep::Integer=windowstep, samplerate::Integer=samplerate)`
Return the length (in seconds) of a song, given the isolated syllables
"""
function songlength{T}(syllables::Array{Syllable{T}, 1}, windowStep::Integer=windowstep, sampleRate::Integer=samplerate)
   min_t_idx = minimum(syllables[1].times)
   max_t_idx = maximum(syllables[end].times)

   song_length_idx = max_t_idx-min_t_idx
   song_length = timeslice2sec(song_length_idx, windowStep, sampleRate)

   return song_length
end

"""
`minpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate) -> minimum pitch`
Return the minimum pitch attained by any syllable in `syllables`
"""
function minpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    min_harmonic = minimum(syll->minimum(syll.harmonics), syllables)
    min_pitch = harmonic2hz(min_harmonic, lowHarmonic, windowLength, sampleRate)

    return min_pitch
end

"""
`maxpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate) -> maximum pitch`
Return the maximum pitch attained by any syllable in `syllables`
"""
function maxpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    max_harmonic = maximum(syll->maximum(syll.harmonics), syllables)
    max_pitch = harmonic2hz(max_harmonic, lowHarmonic, windowLength, sampleRate)

    return max_pitch
end

"""
`pitch_distribution{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)`
Return the distribution of frequencies in a song
"""
function pitch_distribution{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    pool = pitch_pool(syllables, lowHarmonic, windowLength, sampleRate)
    return fit(Histogram, pool)
end

"""
`SDpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate) -> stdDev`
Compute the standard deviation of pitch in `syllables`
"""
function SDpitch{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    return std(pitch_pool(syllables, lowHarmonic, windowLength, sampleRate))
end

"""
`psong{T}(syllables::AbstractArray{Syllable{T}, 1}) -> proportion`
Return the proportion of time that a bird spent making noise relative to the total song time
"""
function psong{T}(syllables::AbstractArray{Syllable{T}, 1}, windowStep::Integer=windowstep, sampleRate::Integer=samplerate)
    singing_timeslices = sum(syll->maximum(syll.times)-minimum(syll.times)+1, syllables)
    singing_time = timeslice2sec(singing_timeslices, windowStep, sampleRate)

    return singing_time/songlength(syllables, windowStep, sampleRate)
end

"""
`syllrate{T}(syllables::Array{Syllable{T}, 1}, windowStep::Integer=windowstep, sampleRate::Integer=samplerate)`
Return the average number of syllables/second in `syllables`
"""
function syllrate{T}(syllables::Array{Syllable{T}, 1}, windowStep::Integer=windowstep, sampleRate::Integer=samplerate)
    return length(syllables)/songlength(syllables, windowStep, sampleRate)
end

"""
`bandlimits{T}(syllables::AbstractArray{Syllable{T}, band_proportion::Float64=1, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate) -> (min_freq, max_freq)`
return a range of frequencies that contain `band_proportion` of the histogram of frequencies of `syllables`
"""
function bandlimits{T}(syllables::AbstractArray{Syllable{T}, 1}, band_proportion::Float64=1.0, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    @argcheck (band_proportion > 0) && (band_proportion <=1)

    fpool = pitch_pool(syllables, lowHarmonic, windowLength, sampleRate)
    max_freq = quantile(fpool, .5+(band_proportion/2))
    min_freq = quantile(fpool, .5-(band_proportion/2))

    return (min_freq, max_freq)
end

"""
`bandwidth{T}(syllables::AbstractArray{Syllable{T}}, band_proportion::Float64=1, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)`
Return the size of the frequency range containg `band_proportion` of the histogram of frequencies in `syllables`
"""
function bandwidth{T}(syllables::AbstractArray{Syllable{T}, 1}, band_proportion::Float64=1, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    @argcheck cutoff >= 0 && cutoff <= 1

    blim = bandlimits(syllables, band_proportion, lowHarmonic, windowLength, sampleRate)

    return blim[2]-blim[1]
end

"""Utility function to convert harmonics to pitch (in HZ)"""
function harmonic2hz(harmonic::Integer, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    return (harmonic+lowHarmonic) * sampleRate/windowLength
end

"""Utility function to convert time slices (in number of windows) to seconds"""
#### NOTE: the psong function above impliticly assumes that this function distributes over addition (because right now it is just a scalar multiplication)
#### Don't change this without making sure that every other function is updated as necessary.
function timeslice2sec(timeslice::Integer, windowStep::Integer=windowstep, sampleRate::Integer=samplerate)
    return timeslice * windowStep/sampleRate
end

"""Utility function to extract all the frequencies in a song"""
function pitch_pool{T}(syllables::AbstractArray{Syllable{T}, 1}, lowHarmonic::Integer=lowharmonic, windowLength::Integer=windowlength, sampleRate::Integer=samplerate)
    pool = vcat(map(x->x.harmonics, syllables)...)
    pitchPool = map(x->harmonic2hz(x, lowHarmonic, windowLength, sampleRate), pool)

    return pitchPool
end

export songlength
export minpitch
export maxpitch
export SDpitch
export pitch_distribution
export psong
export syllrate
export bandlimits
export bandwidth


text:
#=
choose the frequency of experiments to minimize a given convex spectral function
of the information matrix and satisfy an experiment budget constraint and other
side constraints

adapted from Boyd and Vandenberghe, "Convex Optimization", section 7.5

minimize    f(V × Diagonal(x) × V')
subject to  x ≥ 0
            e'x = k
            A x = b
where k = 2d, variable x ∈ ℝᵏ is the frequency of each experiment, k is the
number of experiments to run, the columns of V ∈ ℝ^(d × k) correspond to each
experiment, f is a convex spectral function, and A x = b are side constraints
=#

struct ExperimentDesignJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    d::Int
    ext::MatSpecExt # formulation specifier
end

function build(inst::ExperimentDesignJuMP{T}) where {T <: Float64}
    d = inst.d
    @assert d >= 1
    @assert is_domain_pos(inst.ext)
    k = 2 * d

    V = randn(T, d, k)
    V .*= d / sum(svdvals(V))
    A = randn(T, round(Int, sqrt(k - 1)), k)
    b = sum(A, dims = 2)

    model = JuMP.Model()
    JuMP.@variable(model, x[1:k] >= 0)
    JuMP.@constraint(model, sum(x) == k)
    JuMP.@constraint(model, A * x .== b)

    vec_dim = Cones.svec_length(d)
    Q = V * diagm(x) * V' # information matrix
    Q_vec = zeros(JuMP.AffExpr, vec_dim)
    Cones.smat_to_svec!(Q_vec, Q, sqrt(T(2)))

    # convex objective
    JuMP.@variable(model, epi)
    JuMP.@objective(model, Min, epi)
    add_homog_spectral(inst.ext, d, vcat(1.0 * epi, Q_vec), model)

    # save for use in tests
    model.ext[:Q_var] = Q

    return model
end

function test_extra(inst::ExperimentDesignJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective
    tol = eps(T)^0.2
    Q_opt = JuMP.value.(model.ext[:Q_var])
    λ = eigvals(Symmetric(Q_opt, :U))
    @test minimum(λ) >= -tol
    obj_result = get_val(pos_only(λ), inst.ext)
    @test JuMP.objective_value(model) ≈ obj_result atol=tol rtol=tol
    return
end


text:
mutable struct ModelSelectionData
    equation::Array{Symbol}
    depvar::Symbol
    expvars::Array{Symbol}
    panel::Union{Symbol, Nothing}
    time::Union{Symbol, Nothing}
    depvar_data::Union{Array{Float64}, Array{Float32}, Array{Union{Float64, Missing}}, Array{Union{Float32, Missing}}}
    expvars_data::Union{Array{Float64}, Array{Float32}, Array{Union{Float64, Missing}}, Array{Union{Float32, Missing}}}
    panel_data::Union{Nothing, Array{Int64}, Array{Int32}, Array{Union{Int64, Missing}}, Array{Union{Int32, Missing}}}
    time_data::Union{Nothing, Array{Float64}, Array{Float32}, Array{Union{Float64, Missing}}, Array{Union{Float32, Missing}}}
    intercept::Bool
    datatype::DataType
    removemissings::Bool
    nobs::Int64
    options::Array{Any}
    extras::Dict
    previous_data::Array{Any}
    results::Array{Any}

    function ModelSelectionData(
            equation,
            depvar::Symbol,
            expvars::Array{Symbol},
            panel::Union{Symbol, Nothing},
            time::Union{Symbol, Nothing},
            depvar_data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}},
            expvars_data::Union{Array{Float64}, Array{Float32}, Array{Union{Float32, Missing}}, Array{Union{Float64, Missing}}},
            panel_data::Union{Nothing, Array{Int64}, Array{Int32}, Array{Union{Int64, Missing}}, Array{Union{Int32, Missing}}},
            time_data::Union{Nothing, Array{Float64}, Array{Float32}, Array{Union{Float64, Missing}}, Array{Union{Float32, Missing}}},
            intercept::Bool,
            datatype::DataType,
            removemissings::Bool,
            nobs::Int64
        )
        extras = Dict()
        options = Array{Any}(undef, 0)
        previous_data = Array{Any}(undef, 0)
        results = Array{Any}(undef, 0)
        new(equation, depvar, expvars, panel, time, depvar_data, expvars_data, panel_data, time_data, intercept, datatype, removemissings, nobs, options, extras, previous_data, results)
    end
end


text:
function SubLine(arg0::Line, arg1::IntervalsSet)
    return SubLine((Line, IntervalsSet), arg0, arg1)
end

function SubLine(arg0::Segment)
    return SubLine((Segment,), arg0)
end

function SubLine(arg0::Vector3D, arg1::Vector3D, arg2::jdouble)
    return SubLine((Vector3D, Vector3D, jdouble), arg0, arg1, arg2)
end

function get_segments(obj::SubLine)
    return jcall(obj, "getSegments", List, ())
end

function intersection(obj::SubLine, arg0::SubLine, arg1::jboolean)
    return jcall(obj, "intersection", Vector3D, (SubLine, jboolean), arg0, arg1)
end


text:
module ProgrammableAPI

export PAPIManager, PAPIInput, PAPIOutput, CreateEngine, EngineCreated, WSAPIService

using Circo
using Plugins

PAPIReqId = UInt64

struct PAPIInput
    respondto::Addr
    program::Vector{UInt8}
    PAPIInput(respondto::Addr, program) = new(respondto, program)
end

struct PAPIOutput
    token::UInt64
    body::Vector{UInt8}
end

struct CreateEngine
    token::UInt64
    respondto::Addr
end

struct EngineCreated
    token::UInt64
    addr::Addr
end

struct DestroyEngine
    token::UInt64
end

mutable struct PAPIManager <: Actor{Any}
    engine_factory
    living::Set{Addr}
    got_ping_from::Set{Addr} # TODO use
    core
    PAPIManager(engine_factory) = new(engine_factory, Set(), Set())
end

function Circo.onmessage(me::PAPIManager, msg::CreateEngine, service)
    engine = me.engine_factory(msg.token, emptycore(service))
    addr = spawn(service, engine)
    push!(me.living, addr)
    send(service, me, msg.respondto, EngineCreated(msg.token, addr))
end

include("wsapi.jl")

end # module


text:
function del_void_rxns_mets!(metnet::MetNet)
    # find mets with no rxns
    for (meti, met) in enumerate(metnet.mets)
        rxns = met_rxns(metnet, meti)
        isempty(rxns) && del_met!(metnet, meti)
    end

    # find rxns with no mets
    for (rxni, rxn) in enumerate(metnet.rxns)
        mets = rxn_mets(metnet, rxni)
        isempty(mets) && del_rxn!(metnet, rxni)
    end

    return metnet
end


text:
module UserSessionsController

using Genie, SearchLight, Genie.Authentication, Genie.Flax, Nullables, Genie.Loggers

function show_login()
  Genie.Flax.include_template("/Users/adrian/Dropbox/Projects/todo_mvc/app/layouts/login.flax.html", partial = false)
end

function login()
  success = do_login(@params(:email), @params(password), session(@params()))
  if isnull(success)
    flash("Unknown username and password combination", @params())
    redirect_to("/login")
  else
    flash("You are logged in", @params())
    redirect_to("/")
  end
end

function do_login(email::String, password::String, session::Sessions.Session) :: Nullable
  users = SearchLight.find(User, SQLQuery(where = [SQLWhere(:email, email), SQLWhere(:password, sha256(password) |> bytes2hex)]))

  if isempty(users)
    log("Failed login: Can't find user")
    return Nullable()
  end
  user = users[1]

  Authentication.login(user, session)
end

function logout(session::Sessions.Session)
  Authentication.logout(session(@params()))
end

end


text:
# This file was generated, do not modify it. # hide
using MLJ, CategoricalArrays, PrettyPrinting
import DataFrames: DataFrame, nrow
using UrlDownload

LinearRegressor = @load LinearRegressor pkg=GLM
LinearBinaryClassifier = @load LinearBinaryClassifier pkg=GLM

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials