text:
# Pseudocode: http://www.stat.rice.edu/%7Escottdw/stat550/HW/hw4/c05.pdf

type MVAsh{F1 <: Function, F2 <: Function, R1 <: Range, R2 <: Range} <: AbstractAsh
    kernelx::F1
    kernely::F2
    rngx::R1
    rngy::R2
    mx::Int64
    my::Int64
    v::Matrix{Int}
    z::Matrix{Float64}
    nobs::Int64
end

function MVAsh(
        kernelx::Function, kernely::Function, rngx::Range, rngy::Range,
        mx::Int64, my::Int64, v::Matrix{Int}, z::Matrix{Float64}, nobs::Int64
        )
    @assert kernelx(0.1) > 0.0 && kernelx(-0.1) > 0.0 "kernelx must always be positive"
    @assert kernely(0.1) > 0.0 && kernely(-0.1) > 0.0 "kernely must always be positive"
    @assert length(rngx) > 1 "Need at least two bins"
    @assert length(rngy) > 1 "Need at least two bins"
    @assert length(rngx) == size(v, 2) == size(z, 2)
    @assert length(rngy) == size(v, 1) == size(z, 1)
    @assert mx > 0 "Smoothing parameter mx must be positive"
    @assert my > 0 "Smoothing parameter my must be positive"
    MVAsh{typeof(kernelx), typeof(kernely), typeof(rngx), typeof(rngy)}(
        kernelx, kernely, rngx, rngy, mx, my, v, z, nobs
    )
end
function Base.show(io::IO, o::MVAsh)
    println(io, typeof(o))
    println(io, "  > X")
    println(io, "    - m      : ", o.mx)
    println(io, "    - kernel : ", o.kernelx)
    println(io, "    - edges  : ", o.rngx)
    println(io, "  > Y")
    println(io, "    - m      : ", o.my)
    println(io, "    - kernel : ", o.kernely)
    println(io, "    - edges  : ", o.rngy)
end

function ash(
        x::Vector, y::Vector,
        rngx::Range = extendrange(x, .4, 100), rngy::Range = extendrange(y, .4, 100);
        mx::Int = 5, my::Int = 5,
        kernelx::Function = Kernels.biweight, kernely::Function = Kernels.biweight,
        warnout::Bool = true
    )
    ncol, nrow = length(rngx), length(rngy)
    o = MVAsh(
        kernelx, kernely, rngx, rngy, mx, my,
        zeros(Int64, nrow, ncol), zeros(nrow, ncol), 0
    )
    update!(o, x, y)
    ash!(o; warnout = warnout)
 end

function update!(o::MVAsh, x::AbstractVector, y::AbstractVector)
    n = length(x)
    @assert n == length(y)
    δx, δy = step(o.rngx), step(o.rngy)
    ax, ay = o.rngx[1], o.rngy[1]
    nbinx, nbiny = length(o.rngx), length(o.rngy)

    for i in 1:n
        kx = floor(Int, (x[i] - ax) / δx + 1.5)
        ky = floor(Int, (y[i] - ay) / δy + 1.5)
        if 1 <= kx <= nbinx && 1 <= ky <= nbiny
            o.v[ky, kx] += 1
        end
    end
    o.nobs += n
    o
end

function ash!(o::MVAsh;
        mx::Int64 = o.mx, my::Int64 = o.my,
        kernelx::Function = o.kernelx, kernely::Function = o.kernely,
        warnout::Bool = true
    )
    o.mx, o.my, o.kernelx, o.kernely = mx, my, kernelx, kernely
    δx, δy = step(o.rngx), step(o.rngy)
    nbinx, nbiny = length(o.rngx), length(o.rngy)

    w = zeros(2 * my - 1, 2 * mx - 1)
    for i = (1 - my):(my - 1), j = (1 - mx):(mx - 1)
        w[i + my, j + mx] = kernely(i / my) * kernelx(i / mx)
    end

    v, z = o.v, o.z
    for k = 1:nbiny
       for l = 1:nbinx
           if v[k, l] == 0
                l += 1
            else
                for i = max(1, k - my + 1):min(nbiny, k + my - 1)
                    for j = max(1, l - mx + 1):min(nbinx, l + mx - 1)
                        z[i, j] += v[k, l] * w[my + i - k, mx + j - l]
                    end
                end
            end
        end
    end

    # make density integrate to 1
    denom = 1.0 / (sum(o.z) * δx * δy)
    for i in eachindex(z)
        @inbounds z[i] *= denom
    end

    if warnout
        any(z[end, :] .!= 0)    && warn("nonzero density outside bounds at minimum(rngy)")
        any(z[1, :] .!= 0)      && warn("nonzero density outside bounds at maximum(rngy)")
        any(z[:, 1] .!= 0)      && warn("nonzero density outside bounds at minimum(rngx)")
        any(z[:, end] .!= 0)    && warn("nonzero density outside bounds at maximum(rngx)")
    end
    o
end

xyz(o::MVAsh) = (collect(o.rngx), collect(o.rngy), copy(o.z))
function Base.mean(o::MVAsh)
    meanx = mean(o.rngx, StatsBase.WeightVec(vec(sum(o.z, 1))))
    meany = mean(o.rngy, StatsBase.WeightVec(vec(sum(o.z, 2))))
    [meanx; meany]
end
function Base.var(o::MVAsh)
    varx = var(o.rngx, StatsBase.WeightVec(vec(sum(o.z, 1))))
    vary = var(o.rngy, StatsBase.WeightVec(vec(sum(o.z, 2))))
    [varx; vary]
end
Base.std(o::MVAsh) = sqrt(var(o))
function StatsBase.fit!(o::MVAsh, x::AbstractVector, y::AbstractVector; kw...)
    update!(o, x, y)
    ash!(o; kw...)
end

@recipe function f(o::MVAsh)
    seriestype --> :heatmap
    o.rngx, o.rngy, o.z
end


text:
module TrackFiles

include("trk.jl")

end

ext: jl


text:
using HDF5

function move_onehot(move::Point)
    board = BitArray(N, N)
    fill!(board, false)
    board[linearindex(move)] = true
    board
end

"Generate training data of form [(features, move), ...] for a single sgf file"
function generate_training_data(filename::AbstractString;
                                features::Vector{Function}=DEFAULT_FEATURES)
    sgf = load_sgf(filename)
    examples = Vector{Tuple{BitArray, BitArray}}()
    if sgf == nothing
        return examples
    end
    board = Board()
    while true
        move = get_next_move(sgf)
        if move == nothing
            break
        end
        if move != PASS_MOVE
            computed_features = get_features(board, features=features)
            answer = move_onehot(move)
            push!(examples, (computed_features, answer))
        end
        play_move(board, move)
    end
    return examples
end

"Generate training data for many sgf files"
function generate_training_data(filenames::Vector{AbstractString};
                                progress_update=100,
                                features::Vector{Function} = DEFAULT_FEATURES)
    examples = Vector{Tuple{BitArray, BitArray}}()
    tm = time()
    for (i, filename) in enumerate(filenames)
        if i % progress_update == 0
            println(STDERR, "$(i)/$(length(filenames)): $(time() - tm)")
        end
	game = generate_training_data(filename, features=features)
        if length(game) > 0
            push!(examples, game...)
        end
    end
    examples
end

"Extract training data directly to an hdf5 file"
function extract_to_hdf5(hf5_filename::AbstractString, filenames::Vector{AbstractString}, features::Vector{Function}; progress_update=1000, nlabels = 81)
    example_data = get_features(Board(), features=features)
    data_dim = size(example_data)

    hf5 = h5open(hf5_filename, "w")
    data = d_create(hf5, "data", datatype(UInt8),
                 ((data_dim..., 1), (data_dim..., -1)),
                 "chunk", (data_dim..., 1))
    label = d_create(hf5, "label", datatype(UInt8),
                     ((nlabels, 1), (nlabels, -1)),
                     "chunk", (nlabels, 1))
    tm = time()
    nexamples = 0
    for (i, filename) in enumerate(filenames)
        if i % progress_update == 0
            println(STDERR, "$(i)/$(length(filenames)): $(time() - tm)")
        end
        game = generate_training_data(filename, features=features)
        if length(game) > 0
            new_N = nexamples + length(game)
            set_dims!(data, (data_dim..., new_N))
            set_dims!(label, (nlabels, new_N))
            data_batch = batch_training_examples([x[1] for x in game])
            label_batch = batch_training_examples([x[2][:] for x in game])
            # Assumes data is 3 dims per example and labels are 1
            data[:, :, :, nexamples+1:new_N] = data_batch
            label[:, nexamples+1:new_N] = label_batch
            nexamples = new_N
        end
    end
    close(hf5)
end

"""
Create a vector by appending the given arrays in memory.
Reshape such that they are concatenated along the last+1 dimension

Wrote this because
X = cat(4,[x[1] for x in examples]...);
is slow for some reason
"""
function batch_training_examples(arrs)
    sizeper = length(arrs[1])
    out = similar(arrs[1], sizeper * length(arrs))
    for (i, arr) in enumerate(arrs)
        start = sizeper * (i - 1) + 1
        finish = sizeper * (i)
        out[start:finish] = arrs[i][:]
    end
    reshape(out, size(arrs[1])..., length(arrs))
end


text:
using BayesPlot
using Test

@testset "BayesPlot.jl" begin
    # Write your tests here.
end

max_stars_repo_name: kskyten/BayesPlot

max_issues_repo_name: kskyten/BayesPlot

max_forks_repo_name: kskyten/BayesPlot


text:
import BEAST; BE = BEAST;
using CompScienceMeshes
using Base.Test

T = Float64
P = SVector{3,T}
tol = eps(T)*10^3

# s1 = int d/R^3
# s2 = int 1/R
# s3 = int 1
# s4 = int R
# s5 = int R^2
# s6 = int R^3
# v1 = int (rho'-rho)/R^3
# v2 = int (rho'-rho)/R
# v3 = int (rho'-rho)
# v4 = int (rho'-rho)*R
# v5 = int (rho'-rho)*R
# g1 = int grad' (1/R)
# g2 = int grad' R
# g3 = int grad' R^2
# g4 = int grad' R^3
function withquadrules(triangle, r, n)

	u, w = BE.trgauss(n)

    n = (triangle[1] - triangle[3]) × (triangle[2] - triangle[3])
    n /= norm(n)
    d = (r - triangle[1]) ⋅ n
    ρ = r - d*n

    scalar = zeros(T,6)
    vector = zeros(P,6)
    gradgr = zeros(P,4)

	m = simplex(triangle[1], triangle[2], triangle[3])

    for i in 1 : length(w)
        q = neighborhood(m, u[:,i])
        s = cartesian(q)
        dq = w[i] * jacobian(q)

        R = norm(s - r)
        scalar[1] += d / R^3 * dq
        scalar[2] += 1 / R * dq
        scalar[3] += 1 * dq
        scalar[4] += R * dq
        scalar[5] += R^2 * dq
        scalar[6] += R^3 * dq
        vector[1] += (s - ρ) / R^3 * dq
        vector[2] += (s - ρ) / R * dq
        vector[3] += (s - ρ) * dq
        vector[4] += (s - ρ) * R * dq
        vector[5] += (s - ρ) * R^2 * dq
    end

    gradgr[1] = -1( vector[1] - scalar[1] * n )
    gradgr[2] = +1( vector[2] - d*scalar[2]*n )
    gradgr[3] = +2( vector[3] - d*scalar[3]*n )
    gradgr[4] = +3( vector[4] - d*scalar[4]*n )

    return scalar, vector, gradgr
end

tr = [
	point(0.0, 0.0, 0.0),
	point(2.0, 0.0, 0.0),
	point(0.0, 3.0, 0.0)]

pts = [
    (tr[1]+2*tr[2]-1*tr[3]),
    (tr[1]+tr[2])/2 + point(0.0, 0.0, 0.3),
    point(3.0, 0.0, 0.0),
    (tr[1]+tr[2]+tr[3])/3 + point(0.0, 0.0, 0.3),
    tr[2] + point(0.0, 0.0, 0.3),
    (tr[1]+tr[2]+tr[3])/3 + point(0.0, 0.0, 100.0)
]

for p in pts
    s1, v1, g1 = BE.wiltonints(tr[1],tr[2],tr[3],p)
    s2, v2, g2 = withquadrules(tr,p,13)

    for i in length(s1)
        @test norm(s1[i]-s2[i]) < 1.0e-5
    end

    for i in length(v1)
        @test norm(v1[i]-v2[i]) < 1.0e-5
    end

    for i in 1:length(g1)
        for j in 1:3
            @test norm(g1[i][j]-g2[i][j]) < 1.0e-5
        end
    end
end

Γ = readmesh(joinpath(dirname(@__FILE__),"assets","sphere2.in"))
nc = numcells(Γ)
t = chart(Γ, first(cells(Γ)))
s = chart(Γ, last(cells(Γ)))

X = BEAST.raviartthomas(Γ)
x = BEAST.refspace(X)

κ = 1.0
op = BEAST.MWSingleLayer3D(κ)

n = BE.numfunctions(x)
z1 = zeros(Complex128, n, n)
z2 = zeros(z1)

BE = BEAST

tqd = BE.quadpoints(x, [t], (12,))
bqd = BE.quadpoints(x, [s], (13,))

DQ_strategy = BE.DoubleQuadStrategy(tqd[1,1], bqd[1,1])
BEAST.momintegrals!(op, x, x, t, s, z1, DQ_strategy)

SE_strategy = BE.WiltonSEStrategy(
  tqd[1,1],
  BE.DoubleQuadStrategy(
	tqd[1,1],
	bqd[1,1],
  ),
)
BEAST.momintegrals!(op, x, x, t, s, z2, SE_strategy)

@test norm(z1-z2) < 1.0e-7


text:
using AlgebraicTraits
using Base.Test

@test ismonoid(String, *)
@test !ismonoid(String, +)

@test ismonoid(Int, +)
@test ismonoid(Int, *)


text:
@service S3

bucket_name = "aws-jl-test-issues---" * _now_formatted()
S3.create_bucket(bucket_name)

@testset "issue 223" begin
    # https://github.com/JuliaCloud/AWS.jl/issues/223
    body = "Hello World!"
    file_name = "contains spaces"

    try
        S3.put_object(bucket_name, file_name, Dict("body" => body))
        resp = S3.get_object(bucket_name, file_name)

        @test String(resp) == body
    finally
        S3.delete_bucket(bucket_name)
    end
end


text:
# Fast sodium current

# Rate of change of gating variables
function d_mna(m_na, vm)
    mα = 0.32 * (vm + 47.13) / (-expm1(-0.1 * (vm + 47.13)))
    mβ = 0.08 * exp(-vm / 11.0)
    d_m_na = mα - m_na * (mα + mβ)
end

function d_hna(h_na, vm)
    if vm >= -40
        hα = 0.0
        hβ = inv(0.13 * (1 + exp((vm + 10.66) / -11.1)))
    else
        hα = 0.135 * exp(-(vm + 80.0) / 6.8)
        hβ = 3.56 * exp(0.079 * vm) + 3.1e5 * exp(0.35 * vm)
    end
    d_h_na = hα - h_na * (hα + hβ)
end

function d_jna(j_na, vm)
    if vm >= -40
        jα = 0.0
        jβ = 0.3 * exp(-2.535e-7 * vm) / (1 + exp(-0.1 * (vm + 32)))
    else
        jα = ((-127140 * exp(0.2444 * vm) - 3.471e-5 * exp(-0.04391 * vm))
              * (vm + 37.78) / (1 + exp(0.311 * (vm + 79.23))))
        jβ = 0.1212 * exp(-0.01052 * vm) / (1 + exp(-0.1378 * (vm + 40.14)))
    end
    d_j_na = jα - j_na * (jα + jβ)
end

# Fast sodium current (iNa)
ina(vm, eNa, m_na, h_na, j_na, G_NA) = G_NA * (vm - eNa) * m_na^3 * h_na * j_na


text:
@testset "Simple shader execution" begin
    instance = Instance([], [])

    physical_device = find_physical_device(
        instance,
        check = props ->
            contains(props.device_name, get(ENV, "VKCOMPUTEUTILS_GPU", "SwiftShader")),
    )

    @test physical_device isa PhysicalDevice

    @test_throws ErrorException find_queue_family_idx(physical_device, check = _ -> false)

    qfam_idx = find_queue_family_idx(
        physical_device,
        check = qfp -> hasbits(qfp.queue_flags, QUEUE_TRANSFER_BIT | QUEUE_COMPUTE_BIT),
        score = qfp -> qfp.queue_count,
    )

    device = Device(physical_device, [DeviceQueueCreateInfo(qfam_idx, [1.0])], [], [])

    @test device isa Device

    memorytype_local_visible = find_memory_type_idx(
        physical_device,
        check = (mt, heap) -> hasbits(
            mt.property_flags,
            MEMORY_PROPERTY_HOST_VISIBLE_BIT | MEMORY_PROPERTY_DEVICE_LOCAL_BIT,
        ),
        score = (mt, heap) -> heap.size,
    )

    items = 100
    mem_size = sizeof(Float32) * items

    buffer = unwrap(
        create_buffer(
            device,
            mem_size,
            BUFFER_USAGE_STORAGE_BUFFER_BIT,
            SHARING_MODE_EXCLUSIVE,
            [qfam_idx],
        ),
    )
    @test buffer isa Buffer

    mem = unwrap(allocate_memory(device, mem_size, memorytype_local_visible))
    @test mem isa DeviceMemory

    memptr = map_memory(device, mem, 0, mem_size)
    data = unsafe_wrap(
        Vector{Float32},
        convert(Ptr{Float32}, unwrap(memptr)),
        items,
        own = false,
    )
    data .= 0
    flush_mapped_memory_ranges(device, [MappedMemoryRange(mem, 0, mem_size)])

    bind_buffer_memory(device, buffer, mem, 0)

    struct shader_push_consts
        val::Float32
        n::UInt32
    end

    shader_code = """
    #version 430

    layout (local_size_x_id = 0) in;
    layout (local_size_x = 1, local_size_y = 1, local_size_z = 1) in;

    layout(constant_id = 0) const uint blocksize = 1; // copy of local_size_x

    layout (push_constant) uniform Params
    {
        float val;
        uint n;
    } params;

    layout(std430, binding=0) buffer databuf
    {
        float data[];
    };

    void
    main()
    {
        uint i = gl_GlobalInvocationID.x;
        if(i < params.n) data[i] = params.val * i;
    }
    """

    const_local_size_x = UInt32(32)
    cpl = ComputeShaderPipeline(
        device,
        shader_code,
        1,
        (const_local_size_x,),
        shader_push_consts,
    )
    @test cpl isa ComputeShaderPipeline{shader_push_consts,1}

    cmdpool = unwrap(create_command_pool(device, qfam_idx))
    cbufs = unwrap(
        allocate_command_buffers(
            device,
            CommandBufferAllocateInfo(cmdpool, COMMAND_BUFFER_LEVEL_PRIMARY, 1),
        ),
    )
    cbuf = first(cbufs)

    compute_q = get_device_queue(device, qfam_idx, 0)

    write_descriptor_set_buffers(device, cpl, [buffer])

    begin_command_buffer(
        cbuf,
        CommandBufferBeginInfo(flags = COMMAND_BUFFER_USAGE_ONE_TIME_SUBMIT_BIT),
    )

    some_val = 1.2345

    cmd_bind_dispatch(
        cbuf,
        cpl,
        shader_push_consts(some_val, items),
        n_blocks(items, const_local_size_x),
        1,
        1,
    )
    end_command_buffer(cbuf)

    queue_submit(compute_q, [SubmitInfo([], [], [cbuf], [])])
    GC.@preserve cpl begin
        unwrap(queue_wait_idle(compute_q))
    end

    invalidate_mapped_memory_ranges(device, [MappedMemoryRange(mem, 0, mem_size)])

    @test isapprox(data, some_val .* (0:items-1))
end


text:
@doc raw"""
`SDE`: Stratonovich Stochastic Differential Equation

Defines a stochastic differential initial value problem
```math
\begin{aligned}
dq (t) &= v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} ,
\end{aligned}
```
with drift vector field ``v``, diffusion matrix ``B``,
initial conditions ``q_{0}``, the dynamical variable ``q``
taking values in ``\mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `vType <: Function`: type of `v`
* `BType <: Function`: type of `B`
* `pType <: Union{NamedTuple,Nothing}`: parameters type

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector field ``v``
* `m`:  dimension of the Wiener process
* `ns`: number of sample paths
* `v`:  function computing the deterministic vector field
* `B`:  function computing the d x m diffusion matrix
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)
* `parameters`: either a `NamedTuple` containing the equations parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

The functions `v` and `B`, providing the drift vector field and diffusion matrix,
`v(t, q, v)` and `B(t, q, B, col=0)`, where `t` is the current time, `q` is the
current solution vector, and `v` and `B` are the variables which hold the result
of evaluating the vector field ``v`` and the matrix ``B`` on `t` and `q` (if col==0),
or the column col of the matrix B (if col>0).

### Constructors

```julia
SDE(m, ns, v, B, t₀, q₀; parameters=NullParameters(), periodicity=zero(q₀[begin]))
SDE(m, ns, v, B, q₀::StateVector; kwargs...) = SDE(m, ns, v, B, 0.0, q₀; kwargs...)
SDE(m, ns, v, B, t₀, q₀::State; kwargs...) = SDE(m, ns, v, B, t₀, [q₀]; kwargs...)
SDE(m, ns, v, B, q₀::State; kwargs...) = SDE(m, ns, v, B, 0.0, q₀; kwargs...)
```

### Example

```julia
    function v(t, q, v, p)
        λ = p[:λ]
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(t, q, B, p, col=0)
        μ = p[:μ]
        if col==0 #whole matrix
            B[1,1] = μ*q[1]
            B[2,1] = μ*q[2]
        elseif col==1
            #just first column
        end
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.
    p = (λ=λ, μ=μ)

    sde = SDE(v, B, t₀, q₀; parameters=p)
```
"""
struct SDE{dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
           vType <: Function, BType <: Function,
           pType <: OptionalParameters} <: AbstractEquationSDE{Nothing,Nothing,Nothing}

    d::Int
    m::Int
    ns::Int
    v::vType
    B::BType
    t₀::tType
    q₀::Vector{arrayType}
    parameters::pType
    periodicity::arrayType

    function SDE(m, ns, v::vType, B::BType, t₀::tType, q₀::Vector{arrayType};
                 parameters=NullParameters(), periodicity=zero(q₀[begin])) where {
                        dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
                        vType <: Function, BType <: Function}

        d  = length(q₀[begin])
        ni = length(q₀)

        @assert all(length(q) == d for q in q₀)

        @assert ns ≥ 1
        @assert ni == 1 || ns == 1
        # either multiple deterministic initial conditions and one sample path
        # or one deterministic initial condition and multiple sample paths

        new{dType,tType,arrayType,vType,BType,typeof(parameters)}(d, m, ns, v, B, t₀, q₀, parameters, periodicity)
    end
end

SDE(m, ns, v, B, q₀::StateVector; kwargs...) = SDE(m, ns, v, B, 0.0, q₀; kwargs...)
SDE(m, ns, v, B, t₀, q₀::State; kwargs...) = SDE(m, ns, v, B, t₀, [q₀]; kwargs...)
SDE(m, ns, v, B, q₀::State; kwargs...) = SDE(m, ns, v, B, 0.0, q₀; kwargs...)

# const SDEHT{HT,DT,TT,AT,VT,BT,PT} = SDE{DT,TT,AT,VT,BT,HT,PT} # type alias for dispatch on Hamiltonian type parameter
# const SDEPT{PT,DT,TT,AT,VT,BT,HT} = SDE{DT,TT,AT,VT,BT,HT,PT} # type alias for dispatch on parameters type parameter
const SDEPT{PT,DT,TT,AT,VT,BT} = SDE{DT,TT,AT,VT,BT,PT} # type alias for dispatch on parameters type parameter

Base.hash(sde::SDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.ns,
                               hash(sde.v, hash(sde.B, hash(sde.t₀, hash(sde.q₀,
                               hash(sde.parameters, hash(sde.periodicity, h)))))))))

Base.:(==)(sde1::SDE, sde2::SDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.ns == sde2.ns
                             && sde1.v == sde2.v
                             && sde1.B == sde2.B
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.parameters == sde2.parameters
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(equ::SDE, t₀::Real, q₀::StateVector, ns::Int=equ.ns)
    @assert all([length(q) == ndims(equ) for q in q₀])
    SDE(equ.m, ns, equ.v, equ.B, t₀, q₀; parameters=equ.parameters, periodicity=equ.periodicity)
end

Base.similar(equ::SDE, t₀::Real, q₀::State, ns::Int=equ.ns; kwargs...) = similar(equ, t₀, [q₀], ns; kwargs...)
Base.similar(equ::SDE, q₀::Union{State,StateVector}, ns::Int=equ.ns; kwargs...) = similar(equ, equ.t₀, q₀, ns; kwargs...)

Base.ndims(sde::SDE) = sde.d
Base.axes(equ::SDE) = axes(equ.q₀[begin])
GeometricBase.nsamples(equ::SDE) = length(equ.q₀)
GeometricBase.periodicity(equ::SDE) = equ.periodicity

initial_conditions(equ::SDE) = (equ.t₀, equ.q₀)

# hashamiltonian(::SDEHT{<:Nothing}) = false
# hashamiltonian(::SDEHT{<:Function}) = true

hasparameters(::SDEPT{<:NullParameters}) = false
hasparameters(::SDEPT{<:NamedTuple}) = true

_get_v(equ::SDE) = hasparameters(equ) ? (t,q,v) -> equ.v(t, q, v, equ.parameters) : equ.v
_get_B(equ::SDE) = hasparameters(equ) ? (t,q,B,col=0) -> equ.B(t, q, B, equ.parameters, col) : equ.B
# _get_h(equ::SDE) = hasparameters(equ) ? (t,q) -> equ.h(t, q, equ.parameters) : equ.h

function functions(equ::SDE)
    names = (:v,:B)
    equs  = (_get_v(equ), _get_B(equ))

    # if hashamiltonian(equ)
    #     names = (names..., :h)
    #     equs  = (equs..., _get_h(equ))
    # end

    NamedTuple{names}(equs)
end


text:
export BidirectionalChannel, take!, put!, inChannel, outChannel, isready, isfull, isopen, eltype

struct BidirectionalChannel{T}
    in::Channel{T}
    out::Channel{T}
end

function BidirectionalChannel{T}(sz::Integer = 0) where {T}
    in = Channel{T}(sz)
    out = Channel{T}(sz)
    return BidirectionalChannel(in, out)
end

function BidirectionalChannel{T}(biChannel::BidirectionalChannel{T}) where {T}
    in = biChannel.out
    out = biChannel.in
    return BidirectionalChannel(in, out)
end

inChannel(biChannel::BidirectionalChannel) = biChannel.in
outChannel(biChannel::BidirectionalChannel) = biChannel.out

take!(biChannel::BidirectionalChannel{T}) where {T} = take!(inChannel(biChannel))
put!(biChannel::BidirectionalChannel{T}, value::T) where {T} = put!(outChannel(biChannel), value)

isready(biChannel::BidirectionalChannel) = isready(inChannel(biChannel))
isfull(biChannel::BidirectionalChannel) = length(outChannel(biChannel).data) >= outChannel(biChannel).sz_max
isopen(biChannel::BidirectionalChannel) = isopen(inChannel(biChannel)) && isopen(outChannel(biChannel))
function close(biChannel::BidirectionalChannel)
    close(biChannel.in)
    close(biChannel.out)
end

eltype(::Type{BidirectionalChannel{T}}) where {T} = T

wait(biChannel::BidirectionalChannel{T}) where {T} = wait(inChannel(biChannel))


text:
const DEFAULT_RESOLUTION = Ref((1920, 1080))

if Sys.iswindows()
    function primary_resolution()
        dc = ccall((:GetDC, :user32), Ptr{Cvoid}, (Ptr{Cvoid},), C_NULL)
        ntuple(2) do i
            Int(ccall((:GetDeviceCaps, :gdi32), Cint, (Ptr{Cvoid}, Cint), dc, (2 - i) + 117))
        end
    end
elseif Sys.isapple()
    const _CoreGraphics = "CoreGraphics.framework/CoreGraphics"
    function primary_resolution()
        dispid = ccall((:CGMainDisplayID, _CoreGraphics), UInt32,())
        height = ccall((:CGDisplayPixelsHigh,_CoreGraphics), Int, (UInt32,), dispid)
        width = ccall((:CGDisplayPixelsWide,_CoreGraphics), Int, (UInt32,), dispid)
        return (width, height)
    end
else
    # TODO implement linux
    primary_resolution() = DEFAULT_RESOLUTION[]
end

"""
Returns the resolution of the primary monitor.
If the primary monitor can't be accessed, returns (1920, 1080) (full hd)
"""
function primary_resolution end

"""
Returns a reasonable resolution for the main monitor.
(right now just half the resolution of the main monitor)
"""
reasonable_resolution() = primary_resolution() .÷ 2

#=
Conservative 7-color palette from Points of view: Color blindness, Bang Wong - Nature Methods
https://www.nature.com/articles/nmeth.1618?WT.ec_id=NMETH-201106
=#

const wong_colors = [
    RGB(230/255, 159/255, 0/255),
    RGB(86/255, 180/255, 233/255),
    RGB(0/255, 158/255, 115/255),
    RGB(240/255, 228/255, 66/255),
    RGB(0/255, 114/255, 178/255),
    RGB(213/255, 94/255, 0/255),
    RGB(204/255, 121/255, 167/255),
]

const default_palettes = Attributes(
    color = wong_colors,
    marker = [:circle, :xcross, :utriangle, :diamond, :dtriangle, :star8, :pentagon, :rect],
    linestyle = [nothing, :dash, :dot, :dashdot, :dashdotdot],
    side = [:left, :right]
)

const minimal_default = Attributes(
    palette = default_palettes,
    font = "Dejavu Sans",
    backgroundcolor = :white,
    color = :black,
    colormap = :viridis,
    marker = Circle,
    markersize = 0.1,
    linestyle = nothing,
    resolution = reasonable_resolution(),
    visible = true,
    clear = true,
    show_axis = true,
    show_legend = false,
    scale_plot = true,
    center = true,
    update_limits = true,
    axis = Attributes(),
    axis2d = Attributes(),
    axis3d = Attributes(),
    legend = Attributes(),
    axis_type = automatic,
    camera = automatic,
    limits = automatic,
    padding = Vec3f0(0.05),
    raw = false,
    SSAO = Attributes(
        # enable = false,
        bias = 0.025f0,       # z threshhold for occlusion
        radius = 0.5f0,       # range of sample positions (in world space)
        blur = Int32(2),      # A (2blur+1) by (2blur+1) range is used for blurring
        # N_samples = 64,       # number of samples (requires shader reload)
    ),
)

const _current_default_theme = deepcopy(minimal_default)

function current_default_theme(; kw_args...)
    return merge!(Attributes(kw_args), _current_default_theme)
end

function set_theme!(new_theme::Attributes)
    empty!(_current_default_theme)
    new_theme = merge!(new_theme, deepcopy(minimal_default))
    merge!(_current_default_theme, new_theme)
    return
end

function set_theme!(;kw_args...)
    set_theme!(Attributes(; kw_args...))
end


text:
using Underscores
using Utils

function load()
    open(joinpath(@__DIR__, "03_tree-map.txt"), "r") do f
        permutedims(reduce(hcat, collect.(eachline(f))))
    end
end

function count_trees(map, right, down)
    length, wrap = size(map)
    col, row = 1, 1
    trees = 0
    while row <= length
        if map[row, mod1(col, wrap)] == '#'
            trees += 1
        end
        col += right
        row += down
    end
    trees
end

function tree_product(map, slopes)
    @_ prod(count_trees(map, _...), slopes)
end

trees = load()
slopes = [(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)]

println("Trees for 3 and 1: $(count_trees(trees, 3, 1))")
println("Product of tree counts for $slopes: $(tree_product(trees, slopes))")

using Test

@testset "Toboggan Trajectory" begin
    @test size(trees, 1) == 323
    @test size(trees, 2) == 31

    @testset "example input" begin
        example_input = """
            ..##.......
            #...#...#..
            .#....#..#.
            ..#.#...#.#
            .#...##..#.
            ..#.##.....
            .#.#.#....#
            .#........#
            #.##...#...
            #...##....#
            .#..#...#.#
            """
        example = permutedims(reduce(hcat, collect.(eachline(IOBuffer(example_input)))))

        @test size(example) == (11, 11)

        @testset "part one" begin
            @test count_trees(example, 3, 1) == 7
        end

        @testset "part two" begin
            @test tree_product(example, slopes) == 336
        end
    end

end

max_stars_repo_name: CmdQ/AoC


text:
# Useful utilities to support internal implementation

macro check_argdims(cond)
    quote
        ($(esc(cond))) || throw(DimensionMismatch("Inconsistent argument dimensions."))
    end
end

_rcopy!(r::StridedVecOrMat, x::StridedVecOrMat) = (r === x || copyto!(r, x); r)

function _addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real)
    if c == one(c)
        for i = 1:length(a)
            @inbounds r[i] = a[i] + b[i]
        end
    else
        for i = 1:length(a)
            @inbounds r[i] = a[i] + b[i] * c
        end
    end
    return r
end

function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::Real)
    n = size(a, 1)
    for i = 1:n
        @inbounds a[i,i] += v
    end
    return a
end

function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::Vector, c::Real)
    n = size(a, 1)
    @check_argdims length(v) == n
    if c == one(c)
        for i = 1:n
            @inbounds a[i,i] += v[i]
        end
    else
        for i = 1:n
            @inbounds a[i,i] += v[i] * c
        end
    end
    return a
end

_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::Real) = _adddiag!(copy(a), v)
_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::Vector, c::Real) = _adddiag!(copy(a), v, c)
_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::Vector{T}) where {T<:Real} = _adddiag!(copy(a), v, one(T))

function wsumsq(w::AbstractVector, a::AbstractVector)
    @check_argdims(length(a) == length(w))
    s = zero(promote_type(eltype(w), eltype(a)))
    for i = 1:length(a)
        @inbounds s += abs2(a[i]) * w[i]
    end
    return s
end

function colwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    n = length(r)
    @check_argdims n == size(a, 2) == size(b, 2) && size(a, 1) == size(b, 1)
    for j = 1:n
        v = zero(promote_type(eltype(a), eltype(b)))
        @simd for i = 1:size(a, 1)
            @inbounds v += a[i, j]*b[i, j]
        end
        r[j] = v
    end
    return r
end

function colwise_sumsq!(r::AbstractArray, a::AbstractMatrix, c::Real)
    n = length(r)
    @check_argdims n == size(a, 2)
    for j = 1:n
        v = zero(eltype(a))
        @simd for i = 1:size(a, 1)
            @inbounds v += abs2(a[i, j])
        end
        r[j] = v*c
    end
    return r
end


text:
@testset "Edit a descriptor" begin

    f1 = Field("width")
    f1.typed = "integer"
    f1.constraints.required = true
    f2 = Field("name")
    f2.typed = "string"
    f2.constraints.required = false

    @testset "Create and set field properties" begin
        s = Schema()
        TableSchema.add_field(s, f1)
        TableSchema.add_field(s, f2)
        @test length(s.fields) == 2
        @test s.fields[1].constraints.required
        @test !s.fields[2].constraints.required
    end

    @testset "Modify primary and foreign keys" begin
        s = Schema("../data/schema_valid_full.json")
        @test length(s.primary_key) == 4
        @test_throws SchemaError TableSchema.set_primary_key(s, "invalid")
        TableSchema.set_primary_key(s, "home_location")
        @test length(s.primary_key) == 4
        TableSchema.set_primary_key(s, "gender")
        @test length(s.primary_key) == 5
        TableSchema.set_foreign_key(s, ["home_location"], "country", ["name"])
        @test length(s.foreign_keys) == 2
        @test_throws SchemaError TableSchema.set_foreign_key(s, ["invalid"], "", ["name"])
    end

    @testset "Validate changes to the descriptor" begin
        s = Schema("../data/schema_valid_full.json")
        validate(s, true)
        @test TableSchema.is_valid(s)
        TableSchema.remove_field(s, "position_title")
        TableSchema.remove_field(s, "first_name")
        @test_throws SchemaError validate(s, true)
    end

end


text:
# Color Definitions
PALETTE_FILE = joinpath(
    dirname(dirname(pathof(PowerGraphics))),
    "report_templates",
    "color-palette.yaml",
)

struct PaletteColor
    category::AbstractString
    RGB::AbstractString
    color::RGBA{Float64}
    order::Int64
end

function PaletteColor(category::String, RGB::String, order::Int64)
    rgba =
        parse.(Int64, strip.(split(strip(RGB, ['r', 'g', 'b', 'a', '(', ')', ' ']), ",")))
    color = Colors.RGBA(rgba[1] / 288, rgba[2] / 288, rgba[3] / 288, rgba[4])
    return PaletteColor(category, RGB, color, order)
end

function get_palette(file = nothing)
    file = isnothing(file) ? PALETTE_FILE : file
    palette_config = YAML.load_file(file)
    palette_colors = []
    for (k, v) in palette_config
        push!(palette_colors, PaletteColor(k, v["RGB"], v["order"]))
    end
    sort!(palette_colors, by = x -> x.order)
    return palette_colors
end
function get_default_palette()
    default_palette = []
    palette = get_palette()
    default_order = [6, 52, 14, 1, 32, 7, 18, 20, 27, 53, 17] # the default order from the color palette #
    for i in default_order
        for p in palette
            if p.order == i
                push!(default_palette, p)
            end
        end
    end
    return default_palette
end

# Recursively find all subtypes: useful for categorizing variables
function all_subtypes(t::Type)
    st = [t]
    for t in st
        union!(st, InteractiveUtils.subtypes(t))
    end
    return [split(string(s), ".")[end] for s in st]
end

GR_DEFAULT = permutedims(getfield.(get_default_palette(), :color))
FUEL_DEFAULT = getfield.(get_palette(), :color)
PLOTLY_DEFAULT = getfield.(get_default_palette(), :RGB)
PLOTLY_FUEL_DEFAULT = getfield.(get_palette(), :RGB)
CATEGORY_DEFAULT = getfield.(get_palette(), :category)

SUPPORTED_EXTRA_PLOT_KWARGS = [:linestyle, :linewidth]
SUPPORTED_PLOTLY_SAVE_KWARGS =
    [:autoplay, :post_script, :full_html, :animation_opts, :default_width, :default_height]

NEGATIVE_PARAMETERS = [PSY.PowerLoad]
SUPPORTED_CURTAILMENT_PARAMETERS = [PSI.ActivePowerTimeSeriesParameter]

SUPPORTED_CURTAILMENT_VARIABLES = [PSI.ActivePowerVariable]
SUPPORTED_LOAD_VARIABLES = [PSI.ActivePowerVariable]
SUPPORTED_STORAGE_VARIABLES = [PSI.ActivePowerInVariable, PSI.ActivePowerOutVariable]
SUPPORTED_SERVICE_VARIABLES = [PSI.ActivePowerReserveVariable]

SUPPORTED_OVERGENERATION_VARIABLE = PSI.SystemBalanceSlackDown
SUPPORTED_UNSERVEDENERGY_VARIABLES = PSI.SystemBalanceSlackUp
BALANCE_SLACKVARS = Dict(
    SUPPORTED_OVERGENERATION_VARIABLE => "Over Generation",
    SUPPORTED_UNSERVEDENERGY_VARIABLES => "Unserved Energy",
)

LOAD_RENAMING = Dict(
    :ActivePowerTimeSeriesParameter__PowerLoad => :Load,
    :ActivePowerVariable__PowerLoad => :Dispatchable_Load,
)

GENERATOR_MAPPING_FILE = joinpath(
    dirname(dirname(pathof(PowerGraphics))),
    "report_templates",
    "generator_mapping.yaml",
)

function match_fuel_colors(data::DataFrames.DataFrame, backend)
    if backend == Plots.PlotlyJSBackend()
        color_range = PLOTLY_FUEL_DEFAULT
    else
        color_range = FUEL_DEFAULT
    end
    color_fuel = DataFrames.DataFrame(fuels = CATEGORY_DEFAULT, colors = color_range)
    names = DataFrames.names(data)
    default =
        [(color_fuel[findall(in(["$(names[1])"]), color_fuel.fuels), :][:, :colors])[1]]
    for i in 2:length(names)
        @debug names[i] (color_fuel[findall(in(["$(names[i])"]), color_fuel.fuels), :][
            :,
            :colors,
        ])
        specific_color =
            (color_fuel[findall(in(["$(names[i])"]), color_fuel.fuels), :][:, :colors])[1]
        default = hcat(default, specific_color)
    end
    return default
end

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
# Julia wrapper for header: /usr/include/scip/cons_integral.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrIntegral(scip)
    ccall((:SCIPincludeConshdlrIntegral, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

module ClassifyController

using Genie.Renderer
using Genie.Router
using ClassifyModel

function anno()
    query = ClassifyModel.getquery(ClassifyModel.DataSet)
    if query === nothing
        ClassifyModel.dumpdata(ClassifyModel.DataSet)
        redirect_to(:get)
    else
        html!(:classify, :anno, query = "Q: " * query, labelset = ClassifyModel.Labels)
    end
end

function submit()
    @params(:label) |> println
    # ClassifyModel.update!(ClassifyModel.DataSet, ClassifyModel.Labels)
    ClassifyModel.update!(ClassifyModel.DataSet, @params(:label))
    redirect_to(:get_classify)
end

function save()
    ClassifyModel.dumpdata(ClassifyModel.DataSet)
    redirect_to(:get)
end

end


text: import Base.@deprecate_binding


text:
using Interpolations

"""Return the reduced BCS gap deltar = delta(T)/delta(T=0).

Parameters
----------
temp : float
    The reduced temperature temp = T/Tc, where Tc is the superconducting critical temperature.

Returns
-------
dr : float
    The reduced superconducting energy gap delta(T)/delta(T=0).

Note
----
Tabulated data from Muhlschlegel (1959).
Low-temperature analytic formula from Gao (2008).
"""
function deltar_bcs(temp)

    #The BCS gap factor
    const bcs = 1.764

    #These data points run from t=0.18 to t=1
    #in steps of 0.02 from Muhlschlegel (1959)
    delta_calc = [1.0, 0.9999, 0.9997, 0.9994, 0.9989,
           0.9982, 0.9971, 0.9957, 0.9938, 0.9915,
           0.9885, 0.985,  0.9809, 0.976,  0.9704,
           0.9641, 0.9569, 0.9488, 0.9399, 0.9299,
           0.919,  0.907,  0.8939, 0.8796, 0.864,
           0.8471, 0.8288, 0.8089, 0.7874, 0.764,
           0.7386, 0.711,  0.681,  0.648,  0.6117,
           0.5715, 0.5263, 0.4749, 0.4148, 0.3416,
           0.2436, 0.0]

    if 1.0 > temp >= 0.3
        #interpolate data from table
        #Get temperature "index"
        tindx = (temp-0.18)/0.02+1

        #Build interpolating object and get value
        interp_dr = interpolate(delta_calc, BSpline(Cubic(Flat())), OnGrid())
        dr = interp_dr[tindx]

    elseif temp >= 1.0 || temp < 0.0
        dr = 0.0
    elseif temp < 0.005
        dr = 1.0
    else
        #This expression does a nice job of smoothely connecting the table to zero temp
        #Taken from Gao 2008
        dr = exp(-sqrt(2*bcs*temp)*exp(-bcs/temp))
    end

    return dr
end

"""Return the approximated reduced BCS gap deltar = delta(T)/delta(T=0).

Parameters
----------
temp : float
    The reduced temperature temp = T/Tc, where Tc is the superconducting critical temperature.

Returns
-------
dr : float
    The reduced superconducting energy gap delta(T)/delta(T=0).

Note
----
Cosine approximation for gap taken from Pöpel (1989)."""
function deltar_cos(temp)
    if temp < 1
        return sqrt(cos(0.5*π*temp^2))
    else
        return 0.0
    end
end


text:
struct Spectrum_29
    cyclefrequency
    spectralamplitude
end
function Spectrum_29(spect::Spectrogram; unitfreq = 1, unitdens = 1)
    binwidth = spect.freq.multiplier
    amplitudes = sqrt.( 2 * spect.power * binwidth) * unitdens
    # amplitudes can be a matrix if the sample size was larger than the window length
    avgamplitudes = mean(amplitudes, dims = 2)
    Spectrum_29(spect.freq * unitfreq, avgamplitudes)
end

function sampledspectrum_29(v, sampledistance; minimumcyclefrequency = NaN, maximumwavelength = NaN)
    welldefined = xor(isnan(minimumcyclefrequency), isnan(maximumwavelength))
    @assert welldefined "either minimumcyclefrequency or maximumwavelength must be given"
    minimumcyclefrequency = isnan(minimumcyclefrequency) ? 1 / maximumwavelength : minimumcyclefrequency
    @assert dimension(sampledistance * minimumcyclefrequency) == NoDims
    unitλ =  oneunit(sampledistance)
    unitv = oneunit(eltype(v))
    maxcyclesamples =  round(Int, upreferred( 1 / (minimumcyclefrequency * sampledistance)))
    windowlength =  min(nextfastfft(maxcyclesamples * 8), length(v))
    # Calculate parameter for a Tukey window with 10 % of window
    # length outside of flat region
    x_flat = 0.5 * 0.9
    α = 1 - 2x_flat
    wind = tukey(windowlength, α)
    onsid = !(eltype(v) <: ComplexQuantity)
    spectul = spectrogram(v / unitv, windowlength,
        onesided = onsid,
        fs = 1 / (sampledistance / unitλ),
        window = wind
        )
    Spectrum_29(spectul, unitfreq = 1 / unitλ, unitdens = unitv)
end

function spectrumbars_29(spect::Spectrum_29)
    iter = zip(spect.cyclefrequency[2:end], spect.spectralamplitude[2:end])
    pts = Vector{Point}()
    maxa = maximum(spect.spectralamplitude[2:end])
    for (f, a) in iter
        λ = 1 / f
        y = a * 2EM * scale_pt_to_unit(m) / maxa
        push!(pts,  Point(λ, 0y))
        push!(pts,  Point(λ, y))
        push!(pts,  Point(λ, 0y))
    end
    pts
end

function draw_spectrumbars_29(origo, spect::Spectrum_29)
    spectpoints = spectrumbars_29(spect)
    prettypoly(origo .+ spectpoints, :stroke, vertexlabels = (n,l ) -> begin
        mod(n, 3) == 2 && circle(O, 0.006m, :stroke)
    end)
end

function draw_spectrum_29(origo, spect::Spectrum_29)
    draw_spectrumbars_29(origo, spect)
    λ_max = 1 / spect.cyclefrequency[2]
    @layer begin
        sethue(PALETTE[3])
        arrow(origo, origo + (0.0,  -2EM))
        settext("Amplitude", origo + (0.0,  -2EM ), markup=true)
        arrow(origo, origo + (min(1.05λ_max, 0.95 * scale_pt_to_unit(m) * WI), 0.0m))
        settext("Wave length, λ [m]", origo + ( 0.8 * scale_pt_to_unit(m) * WI, 0.0m), markup=true)
    end
end

islocalmax(v, i) = i > 1 && i < length(v) && v[i - 1] < v[i] > v[i + 1]

"""
   two_largest_maxima_amplitude_wavelength(spect::Spectrum_29)

   Returns two tuples, (λ_1, a_1), (λ_2, a_2)

where
   λ_1 is the wavelength of the peak amplitude in the spectrum
   a_1 is the peak amplitude in the spectrum
   (λ_2, a_2) is the second largest peak
"""
function two_largest_maxima_amplitude_wavelength(spect::Spectrum_29)
    amps = spect.spectralamplitude
    ua = eltype(amps)
    freqs = spect.cyclefrequency
    uλ = typeof(1 / freqs[1])

    local_maxima_indices = [i for i= 1:length(amps) if islocalmax(amps, i)]
    local_maxima_values = [amps[i] for i in local_maxima_indices]
    sorted_maxima_indices = sortperm(local_maxima_values, rev = true )
    ind1, ind2 = local_maxima_indices[sorted_maxima_indices[1]], local_maxima_indices[sorted_maxima_indices[2]]

    round3(typ, x) = typ <: Quantity ? round(typ, x, digits = 3) : round(x, digits = 3)

    λ_1 = round3(uλ, 1 / freqs[ind1])
    a_1 = round3(ua, amps[ind1])
    λ_2 = round3(uλ, 1 / freqs[ind2])
    a_2 = round3(ua, amps[ind2])

    return (λ_1, a_1), (λ_2, a_2)
end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

@time begin
    import CSV
    import MAT
    using DataFrames
    import SparseArrays
    using Serialization

    # -------------------------------------------------------------------
    # run add "https://github.com/josePereiro/Chemostat" in the
    # julia Pkg REPL for installing the package
    import Chemostat
    const ChU = Chemostat.Utils
    const ChSS = Chemostat.SteadyState
    const ChLP = Chemostat.LP

    # -------------------------------------------------------------------
    # Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
    # package, then you must activate the package enviroment (see README)
    import Chemostat_Heerden2013
    const ChH = Chemostat_Heerden2013
    const Hd  = ChH.HeerdenData;
    const Bd  = ChH.BegData
    const iJR = ChH.iJR904
end

## -------------------------------------------------------------------
# Tools
function partial_test(model, title  = "PARTIAL TEST")
    withcost = iJR.COST_IDER in model.rxns
    fbaout = withcost ?
        ChLP.fba(model, iJR.BIOMASS_IDER, iJR.COST_IDER; btol = 1e-2) :
        ChLP.fba(model, iJR.BIOMASS_IDER)
    ChU.tagprintln_inmw(title,
        "\nsize:             ", size(model),
        "\nobj_ider:         " , iJR.BIOMASS_IDER,
        "\nfba obj_val:      ", ChU.av(model, fbaout, iJR.BIOMASS_IDER),
        "\nmax exp obj_val:  ", maximum(Hd.val("D")),
        "\ncost_ider:        ", withcost ? iJR.COST_IDER : "not found",
        "\nfba cost_val:     ", withcost ? ChU.av(model, fbaout, iJR.COST_IDER) : "not found",
    )
end

## -------------------------------------------------------------------
# BASE MODEL
# Load Mat file
println("Original .mat model")
src_file = iJR.rawdir("iJR904.mat")
mat_model = MAT.matread(src_file)["model"]
model = ChU.MetNet(mat_model; reshape = true);

## -------------------------------------------------------------------
# the maximal experimental growth rate in Heerden2013 is ~0.2 1/h
# The raw model present a growth rate bigger than that, so it is ok
# to use it directly as base model
partial_test(model)

## -------------------------------------------------------------------
# Gene intactivations
# Heerden2013 describes some gene modifications
# https://doi.org/10.1186/1475-2859-12-80. Table 2
ChU.tagprintln_inmw("GENE MODIFICATIONS")
for (name, iders) in iJR.load_Hd_to_inactivate_map()

    name in ["Aspartate aminotransferase"] && continue ## Closing this one are lethal

    ChU.bounds!.([model], iders, 0.0, 0.0)
    @info("Closing gene modificated", iders)
    partial_test(model)
    println()
end

## -------------------------------------------------------------------
# Set bounds
# The abs maximum bounds will be set to iJR.MAX_ABS_BOUND
ChU.tagprintln_inmw("CLAMP BOUNDS",
    "\nabs max bound: ", iJR.MAX_ABS_BOUND
)
foreach(model.rxns) do ider
        ChU.isfixxed(model, ider) && return # fixxed reaction are untouched

        old_ub = ChU.ub(model, ider)
        new_ub = old_ub == 0.0 ? 0.0 : iJR.MAX_ABS_BOUND
        ChU.ub!(model, ider, new_ub)

        old_lb = ChU.lb(model, ider)
        new_lb = old_lb == 0.0 ? 0.0 : -iJR.MAX_ABS_BOUND
        ChU.lb!(model, ider, new_lb)
end
partial_test(model)

## -------------------------------------------------------------------
exchs = ChU.exchanges(model)
ChU.tagprintln_inmw("CLOSE EXCANGES",
    "\nChU.exchanges: ", exchs |> length
)
# Close, for now, all ChU.exchanges for avoiding it to be in revs
# The reversible reactions will be splited for modeling cost
# Exchanges have not associated cost, so, we do not split them
foreach(exchs) do idx
    ChU.ub!(model, idx, 0.0) # Closing all outtakes
    ChU.lb!(model, idx, 0.0) # Closing all intakes
end

## -------------------------------------------------------------------
# Enzymatic cost info
# The cost will be introduced as a reaction, we follow the same cost models as
# Beg et al. (2007): https://doi.org/10.1073/pnas.0609845104.
# A new balance equations is then added:
#        Σ(rᵢ*costᵢ) + tot_cost = 0
#    Because the cost coefficients (costᵢ) < 0 (it resamble a reactant), the system must allocate
#    the fluxes (rᵢ) so that Σ(rᵢ*costᵢ) = tot_cost, and tot_cost
#    are usually bounded [0.0, 1.0]
cost_info = Dict()
fwd_ider(rxn) = string(rxn, ChU.FWD_SUFFIX);
bkwd_ider(rxn) = string(rxn, ChU.BKWD_SUFFIX);
for rxn in model.rxns
    # The ChU.exchanges, the atpm and the biomass are synthetic reactions, so,
    # they have should not have an associated enzimatic cost
    any(startswith.(rxn, ["EX_", "DM_"])) && continue
    rxn == iJR.BIOMASS_IDER && continue
    rxn == iJR.ATPM_IDER && continue

    # Only the internal, non reversible reactions have an associated cost
    # We will split the rev reactions, so we save the cost for both versions (fwd, bkwd)
    if ChU.isrev(model, rxn)
        cost_info[fwd_ider(rxn)] = -iJR.beg_enz_cost(rxn)
        cost_info[bkwd_ider(rxn)] = -iJR.beg_enz_cost(rxn)
    else
        cost_info[rxn] = -iJR.beg_enz_cost(rxn)
    end
end

## -------------------------------------------------------------------
# Spliting revs
ChU.tagprintln_inmw("SPLITING REVS",
    "\nfwd_suffix:      ", ChU.FWD_SUFFIX,
    "\nbkwd_suffix:     ", ChU.BKWD_SUFFIX,
)

model = ChU.split_revs(model;
    get_fwd_ider = fwd_ider,
    get_bkwd_ider = bkwd_ider,
);

## -------------------------------------------------------------------
# Adding cost raction
cost_met_id = "cost"
cost_exch_id = iJR.COST_IDER
ChU.tagprintln_inmw("ADDING COST",
    "\ncosts to add: ", cost_info |> length,
    "\nmin abs coe:  ", cost_info |> values .|> abs |> minimum,
    "\nmax abs coe:  ", cost_info |> values .|> abs |> maximum,
    "\ncost met id:  ", cost_met_id,
    "\ncost exch id: ", cost_exch_id
)

M, N = size(model)
cost_met = ChU.Met(cost_met_id, S = collect(values(cost_info)),
    rxns = collect(keys(cost_info)), b = 0.0)
model = ChU.expanded_model(model, M + 1, N + 1)
ChU.set_met!(model, ChU.findempty(model, :mets), cost_met)
cost_exch = ChU.Rxn(cost_exch_id, S = [1.0], mets = [cost_met_id],
    lb = -iJR.MAX_ABS_BOUND, ub = 0.0, c = 0.0)
ChU.set_rxn!(model, ChU.findempty(model, :rxns), cost_exch);

## -------------------------------------------------------------------
# Set base ChU.exchanges
ChU.tagprintln_inmw("SETTING EXCHANGES")
# To control the intakes just the metabolites defined in the
# base_intake_info (The minimum medium) will be opened.
# The base model will be constraint as in a cultivation with xi = 1.0
# see Cossios paper (see README)
let
    ξ = 1.0
    base_intake_info = iJR.load_base_intake_info()

    # println("Minimum medium: ", iJR.base_intake_info)
    foreach(exchs) do idx
        ChU.ub!(model, idx, iJR.MAX_ABS_BOUND) # Opening all outakes
        ChU.lb!(model, idx, 0.0) # Closing all intakes
    end

    # see Cossios paper (see README) for details in the Chemostat bound constraint
    ChSS.apply_bound!(model, ξ, base_intake_info);
    # tot_cost is the exchange that controls the bounds of the
    # enzimatic cost contraint, we bound it to [0, 1.0]
    ChU.lb!(model, cost_exch_id, 0.0);
    ChU.ub!(model, cost_exch_id, 1.0);
end

## -------------------------------------------------------------------
function scale_model(model, scale_factor)
    base_nzabs_range = ChU.nzabs_range(model.S)
    base_size = size(model)

    # Scale model (reduce S ill-condition)
    model = ChU.well_scaled_model(model, scale_factor)

    scl_size = size(model)
    scl_nzabs_range = ChU.nzabs_range(model.S)

    @info("Model", exp,
        base_size, base_nzabs_range,
        scl_size, scl_nzabs_range
    ); println()
    return model
end

## -------------------------------------------------------------------
# FVA PREPROCESSING
compressed(model) = model |> ChU.struct_to_dict |> ChU.compressed_copy
MODEL_FILE = iJR.procdir("base_models.bson")
BASE_MODELS = isfile(MODEL_FILE) ?
    ChU.load_data(MODEL_FILE) :
    Dict("load_model" => compressed(model))
cGLCs = Hd.val("cGLC")
for (exp, cGLC) in enumerate(cGLCs)

    D = get!(BASE_MODELS, "fva_models", Dict())
    ChU.tagprintln_inmw("DOING FVA",
        "\nexp:             ", exp,
        "\ncGLC:            ", cGLC,
        "\ncProgress:       ", length(D),
        "\n"
    )
    haskey(D, exp) && continue # cached

    ## -------------------------------------------------------------------
    # prepare model
    # Scale model (reduce S ill-condition)
    scale_factor = 1000.0
    model0 = scale_model(deepcopy(model), scale_factor)

    intake_info = iJR.load_base_intake_info()
    # The only difference between experiments is the feed medium
    # concentration.
    ξ = Hd.val("xi", exp)
    intake_info[iJR.EX_GLC_IDER]["c"] = cGLC
    ChSS.apply_bound!(model0, ξ, intake_info;
        emptyfirst = true, ignore_miss = true
    )

    ## -------------------------------------------------------------------
    # fva
    partial_test(model0)
    fva_model = ChLP.fva_preprocess(model0,
        batchlen = 50,
        check_obj = iJR.BIOMASS_IDER,
        verbose = true
    )
    partial_test(fva_model)

    # storing
    D[exp] = compressed(fva_model)
end

## -------------------------------------------------------------------
# MAX MODEL
let
    # This model is bounded by the maximum rates found for EColi.
    # Data From:
    # Varma, (1993): 2465–73. https://doi.org/10.1128/AEM.59.8.2465-2473.1993.
    # Extract max exchages from FIG 3 to form the maximum polytope

    ChU.tagprintln_inmw("DOING MAX MODEL",
        "\n"
    )

    scale_factor = 1000.0
    max_model = scale_model(deepcopy(model), scale_factor)
    # max_model = deepcopy(model)
    # Biomass
    # 2.2 1/ h
    ChU.bounds!(max_model, iJR.BIOMASS_IDER, 0.0, 2.2)

    Hd_rxns_map = iJR.load_rxns_map()
    # 40 mmol / gDW h
    ChU.bounds!(max_model, Hd_rxns_map["GLC"], -40.0, 0.0)
    # 45 mmol/ gDW
    ChU.bounds!(max_model, Hd_rxns_map["AC"], 0.0, 40.0)
    # 55 mmol/ gDW h
    ChU.bounds!(max_model, Hd_rxns_map["FORM"], 0.0, 55.0)
    # 20 mmol/ gDW h
    ChU.bounds!(max_model, Hd_rxns_map["O2"], -20.0, 0.0)

    # fva
    max_model = ChLP.fva_preprocess(max_model,
        check_obj = iJR.BIOMASS_IDER,
        verbose = true
    );

    ## -------------------------------------------------------------------
    test_model = deepcopy(max_model)
    for (exp, D) in Hd.val(:D) |> enumerate
        cgD_X = Hd.cval(:GLC, exp) * Hd.val(:D, exp) / Hd.val(:DCW, exp)
        ChU.lb!(test_model, iJR.EX_GLC_IDER, -cgD_X)
        fbaout = ChLP.fba(test_model, iJR.BIOMASS_IDER, iJR.COST_IDER)
        biom = ChU.av(test_model, fbaout, iJR.BIOMASS_IDER)
        cost = ChU.av(test_model, fbaout, iJR.COST_IDER)
        @info("Test", exp, cgD_X, D, biom, cost); println()
    end

    ## -------------------------------------------------------------------
    # saving
    BASE_MODELS["max_model"] = compressed(max_model)
end;

## -------------------------------------------------------------------
# SAVING
ChU.save_data(MODEL_FILE, BASE_MODELS);


text:
using OpenAIGym
using PyCall
using Test

"""
`function time_steps(env::GymEnv{T}, num_eps::Int) where T`

run through num_eps eps, recording the time taken for each step and
how many steps were made. Doesn't time the `reset!` or the first step of each
episode (since higher chance that it's slower/faster than the rest, and we want
to compare the average time taken for each step as fairly as possible)
"""
function time_steps(env::GymEnv, num_eps::Int)
    t = 0.0
    steps = 0
    for i in 1:num_eps
        reset!(env)
        # step!(env, rand(env.actions)) # ignore the first step - it might be slow?
        t += (@elapsed steps += epstep(env))
    end
    steps, t
end

"""
Steps through an episode until it's `done`
assumes env has been `reset!`
"""
function epstep(env::GymEnv)
    steps = 0
    while true
        steps += 1
        r, s′ = step!(env, rand(env.actions))
        finished(env, s′) && break
    end
    steps
end

@testset "Gym Basics" begin

    pong = GymEnv(:Pong, :v4)
    pongnf = GymEnv(:PongNoFrameskip, :v4)
    pacman = GymEnv(:MsPacman, :v4)
    pacmannf = GymEnv(:MsPacmanNoFrameskip, :v4)
    cartpole = GymEnv(:CartPole)
    bj = GymEnv(:Blackjack)

    allenvs = [pong, pongnf, pacman, pacmannf, cartpole, bj]
    eps2trial = Dict(pong=>2, pongnf=>1, pacman=>2, pacmannf=>1, cartpole=>400, bj=>30000)
    atarienvs = [pong, pongnf, pacman, pacmannf]
    envs = allenvs

    @testset "string constructor" begin
        for name ∈ ("Pong-v4", "PongNoFrameskip-v4", "MsPacman-v4",
                    "MsPacmanNoFrameskip-v4", "CartPole-v0", "Blackjack-v0")
            env = GymEnv(name)
            @test !PyCall.ispynull(env.pyenv)
        end
    end

    @testset "envs load" begin
        # check they all work - no errors == no worries
        println("------------------------------ Check envs load ------------------------------")
        for (i, env) in enumerate(envs)
            a = rand(env.actions) |> OpenAIGym.pyaction
            action_type = a |> PyObject |> pytypeof
            println("env.pyenv: $(env.pyenv) action_type: $action_type (e.g. $a)")
            time_steps(env, 1)
            @test !ispynull(env.pyenv)
            println("------------------------------")
        end
    end

    @testset "julia speed test" begin
        println("------------------------------ Begin Julia Speed Check ------------------------------")
        for env in envs
            num_eps = eps2trial[env]
            steps, t = time_steps(env, num_eps)
            println("env.pyenv: $(env.pyenv) num_eps: $num_eps t: $t steps: $steps")
            println("microsecs/step (lower is better): ", t*1e6/steps)
            close(env)
            println("------------------------------")
        end
        println("------------------------------ End Julia Speed Check ------------------------------\n")
    end

    @testset "python speed test" begin
        println("------------------------------ Begin Python Speed Check ------------------------------")
        py"""
        import gym
        import numpy as np

        pong = gym.make("Pong-v4")
        pongnf = gym.make("PongNoFrameskip-v4")
        pacman = gym.make("MsPacman-v4");
        pacmannf = gym.make("MsPacmanNoFrameskip-v4");
        cartpole = gym.make("CartPole-v0")
        bj = gym.make("Blackjack-v0")

        allenvs = [pong, pongnf, pacman, pacmannf, cartpole, bj]
        eps2trial = {pong: 2, pongnf: 1, pacman: 2, pacmannf: 1, cartpole: 400, bj: 30000}
        atarienvs = [pong, pongnf, pacman, pacmannf];

        envs = allenvs

        import time
        class Timer(object):
            elapsed = 0.0
            def __init__(self, name=None):
                self.name = name

            def __enter__(self):
                self.tstart = time.time()

            def __exit__(self, type, value, traceback):
                Timer.elapsed = time.time() - self.tstart

        def time_steps(env, num_eps):
            t = 0.0
            steps = 0
            for i in range(num_eps):
                env.reset()
                with Timer():
                    steps += epstep(env)
                t += Timer.elapsed
            return steps, t

        def epstep(env):
            steps = 0
            while True:
                steps += 1
                action = env.action_space.sample()
                state, reward, done, info = env.step(action)
                if done == True:
                    break
            return steps

        for env in envs:
            num_eps = eps2trial[env]
            with Timer():
                steps, s = time_steps(env, num_eps)
            t = Timer.elapsed
            print("{env} num_eps: {num_eps} t: {t} steps: {steps} \n"
                  "microsecs/step (lower is better): {time}".format(
                    env=env, num_eps=num_eps, t=t, steps=steps,
                    time=t*1e6/steps))
            print("------------------------------")
        """
        println("------------------------------ End Python Speed Check ------------------------------")
    end  # @testset "python speed test"

    @testset "Base.show" begin
      let
        io = IOBuffer()
        env = GymEnv(:MsPacman, :v4)
        show(io, env)
        @test String(take!(io)) == "GymEnv MsPacman-v4\n" *
                                   "  TimeLimit\n" *
                                   "  r  = 0.0\n" *
                                   "  ∑r = 0.0"
      end

      let
        io = IOBuffer()
        env = GymEnv(:Blackjack)
        show(io, env)
        @test String(take!(io)) == "GymEnv Blackjack-v0\n" *
                                   "  r  = 0.0\n" *
                                   "  ∑r = 0.0"
      end
    end  # @testset "Base.show"
end


text:
using ApproxFun, Base.Test, Compat
    import ApproxFun:Multiplication,InterlaceOperator

import Compat.view

# test fast copy is consistent with getindex

C=ToeplitzOperator([1.,2.,3.],[4.,5.,6.])

@test_approx_eq full(C[1:5,1:5])    [4.0 5.0 6.0 0.0 0.0
                                     1.0 4.0 5.0 6.0 0.0
                                     2.0 1.0 4.0 5.0 6.0
                                     3.0 2.0 1.0 4.0 5.0
                                     0.0 3.0 2.0 1.0 4.0]

C=Conversion(Ultraspherical(1),Ultraspherical(2))

@test_approx_eq full(C[1:5,1:5])     [1.0 0.0 -0.3333333333333333 0.0  0.0
                                      0.0 0.5  0.0               -0.25 0.0
                                      0.0 0.0  0.3333333333333333 0.0 -0.2
                                      0.0 0.0  0.0                0.25 0.0
                                      0.0 0.0  0.0                0.0  0.2]

for M in (ToeplitzOperator([1.,2.,3.],[4.,5.,6.]),
                HankelOperator([1.,2.,3.,4.,5.,6.,7.]),
                Conversion(Ultraspherical(1),Ultraspherical(2)),
                Multiplication(Fun([1.,2.,3.],Chebyshev()),Chebyshev()))
    B=M[1:5,1:5]

    for k=1:5,j=1:5
        @test_approx_eq B[k,j] M[k,j]
    end

    @test_approx_eq M[1:5,1:5][2:5,1:5] M[2:5,1:5]
    @test_approx_eq M[1:5,2:5] M[1:5,1:5][:,2:end]
end

d=Interval(-10.,5.);
S=Chebyshev(d)

@test norm(Fun(Fun(Fun(exp,S),Ultraspherical(2,d)),S)-Fun(exp,S)) < 100eps()

@test_approx_eq copy(@compat view(Derivative(Ultraspherical(1)),1:2,1:2))[1,2] Derivative(Ultraspherical(1))[1,2]
@test_approx_eq exp(0.1) (Derivative()*Fun(exp,Ultraspherical(1)))(0.1)

T=Float64
f=Fun(exp)
d=domain(f)
D=Derivative(d)

Q=integrate(d)

@test norm((Q+I)*f-(integrate(f)+f)) < 100eps()
@test norm((Q)*f-(integrate(f))) < 100eps()

x=Fun(identity)
X=Multiplication(x,space(x))
d=Interval()

A=Conversion(Chebyshev(d),Ultraspherical(2,d))

@test norm(A\Fun(x.*f,rangespace(A))-(x.*f)) < 100eps()

@test norm((Conversion(Chebyshev(d),Ultraspherical(2,d))\(D^2*f))-f'') < 100eps()

@test norm(X*f-(x.*f)) < 100eps()

A=Conversion(Chebyshev(d),Ultraspherical(2,d))*X
@test norm((A*f.coefficients).coefficients-coefficients(x.*f,rangespace(A))) < 100eps()

## Special functions

x=Fun(identity)
@test norm(cos(x)-Fun(cos))<10eps()
@test norm(sin(x)-Fun(sin))<10eps()
@test norm(exp(x)-Fun(exp))<10eps()
@test norm(sin(x)./x-Fun(x->sinc(x/π)))<100eps()

P=ApproxFun.PermutationOperator([2,1])
@test_approx_eq P[1:4,1:4] [0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]

## Periodic

d=PeriodicInterval(0.,2π)
a=Fun(t-> 1+sin(cos(10t)),d)
D=Derivative(d)
L=D+a
f=Fun(t->exp(sin(t)),d)
u=L\f

@test norm(L*u-f) < 100eps()

d=PeriodicInterval(0.,2π)
a1=Fun(t->sin(cos(t/2)^2),d)
a0=Fun(t->cos(12sin(t)),d)
D=Derivative(d)
L=D^2+a1*D+a0

f=Fun([1,2,3,4,5],space(a1))

@test_approx_eq (Multiplication(a0,Fourier([0.,2π]))*f)(0.1)  (a0(0.1)*f(0.1))
@test_approx_eq ((Multiplication(a1,Fourier([0.,2π]))*D)*f)(0.1)  (a1(0.1)*f'(0.1))
@test_approx_eq (L.ops[1]*f)(0.1) f''(0.1)
@test_approx_eq (L.ops[2]*f)(0.1) a1(0.1)*f'(0.1)
@test_approx_eq (L.ops[3]*f)(0.1) a0(0.1)*f(0.1)
@test_approx_eq (L*f)(0.1) f''(0.1)+a1(0.1)*f'(0.1)+a0(0.1)*f(0.1)

f=Fun(t->exp(cos(2t)),d)
u=L\f

@test norm(L*u-f) < 1000eps()

## Check mixed

d=Interval()
D=Derivative(d)
x=Fun(identity,d)
A=D*(x*D)
B=D+x*D^2
C=x*D^2+D

f=Fun(exp)
@test_approx_eq (A.ops[end]*f)(0.1) f'(0.1)
@test_approx_eq ((x*D)*f)(0.1) 0.1*f'(0.1)
@test_approx_eq (A*f)(0.1) f'(0.1)+0.1*f''(0.1)
@test_approx_eq (B*f)(0.1) f'(0.1)+0.1*f''(0.1)
@test_approx_eq (C*f)(0.1) f'(0.1)+0.1*f''(0.1)

@test norm((A-B)[1:10,1:10]|>full)<eps()
@test norm((B-A)[1:10,1:10]|>full)<eps()
@test norm((A-C)[1:10,1:10]|>full)<eps()
@test norm((C-A)[1:10,1:10]|>full)<eps()
@test norm((C-B)[1:10,1:10]|>full)<eps()
@test norm((B-C)[1:10,1:10]|>full)<eps()

## Cached operator
@test cache(Derivative(Chebyshev(),2))[1,1]==0

S=Chebyshev()
io=ApproxFun.InterlaceOperator([InterlaceOperator(dirichlet(S));Derivative(Chebyshev());lneumann(S)])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]

S=Chebyshev()
io=ApproxFun.InterlaceOperator([InterlaceOperator(dirichlet(S));Derivative(Chebyshev())+Fun(cos);lneumann(S)])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]

S=Chebyshev()
io=ApproxFun.InterlaceOperator([InterlaceOperator(dirichlet(S));Derivative(Chebyshev())])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]

S=Chebyshev()
io=ApproxFun.InterlaceOperator([InterlaceOperator(dirichlet(S));Derivative(Chebyshev())+Fun(cos)])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]

S=Chebyshev()
io=ApproxFun.InterlaceOperator([Derivative(Chebyshev());InterlaceOperator(dirichlet(S))])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]

S=Chebyshev()
io=ApproxFun.InterlaceOperator([Derivative(Chebyshev())+Fun(cos);InterlaceOperator(dirichlet(S))])
co=cache(io)
@test co[1:100,1:100] == io[1:100,1:100]
@test co[1:100,1:100] == io[1:100,1:100]
@test co[200:300,200:300] == io[1:300,1:300][200:300,200:300]


text:
struct LsqFitResult{P, R, J, W <: AbstractArray}
    param::P
    resid::R
    jacobian::J
    converged::Bool
    wt::W
end

StatsBase.coef(lfr::LsqFitResult) = lfr.param
StatsBase.dof(lfr::LsqFitResult) = nobs(lfr) - length(coef(lfr))
StatsBase.nobs(lfr::LsqFitResult) = length(lfr.resid)
StatsBase.rss(lfr::LsqFitResult) = sum(abs2, lfr.resid)
StatsBase.weights(lfr::LsqFitResult) = lfr.wt
StatsBase.residuals(lfr::LsqFitResult) = lfr.resid
mse(lfr::LsqFitResult) = rss(lfr)/dof(lfr)

function check_data_health(xdata, ydata)
    if any(ismissing, xdata) || any(ismissing, ydata)
        error("Data contains `missing` values and a fit cannot be performed")
    end
    if any(isinf, xdata) || any(isinf, ydata) || any(isnan, xdata) || any(isnan, ydata)
        error("Data contains `Inf` or `NaN` values and a fit cannot be performed")
    end
end

# provide a method for those who have their own Jacobian function
function lmfit(f, g, p0::AbstractArray, wt::AbstractArray; kwargs...)
    r = f(p0)
    R = OnceDifferentiable(f, g, p0, copy(r); inplace=false)
    lmfit(R, p0, wt; kwargs...)
end

#for inplace f and inplace g
function lmfit(f!, g!, p0::AbstractArray, wt::AbstractArray, r::AbstractArray; kwargs...)
    R = OnceDifferentiable(f!, g!, p0, copy(r); inplace = true)
    lmfit(R, p0, wt; kwargs...)
end

#for inplace f only
function lmfit(f, p0::AbstractArray, wt::AbstractArray, r::AbstractArray; autodiff = :finite, kwargs...)
    R = OnceDifferentiable(f, p0, copy(r); inplace = true, autodiff = autodiff)
    lmfit(R, p0, wt; kwargs...)
end

function lmfit(f, p0::AbstractArray, wt::AbstractArray; autodiff = :finite, kwargs...)
    # this is a convenience function for the curve_fit() methods
    # which assume f(p) is the cost functionj i.e. the residual of a
    # model where
    #   model(xpts, params...) = ydata + error (noise)

    # this minimizes f(p) using a least squares sum of squared error:
    #   rss = sum(f(p)^2)
    #
    # returns p, f(p), g(p) where
    #   p    : best fit parameters
    #   f(p) : function evaluated at best fit p, (weighted) residuals
    #   g(p) : estimated Jacobian at p (Jacobian with respect to p)

    # construct Jacobian function, which uses finite difference method
    r = f(p0)
    autodiff = autodiff == :forwarddiff ? :forward : autodiff
    R = OnceDifferentiable(f, p0, copy(r); inplace = false, autodiff = autodiff)
    lmfit(R, p0, wt; kwargs...)
end

function lmfit(R::OnceDifferentiable, p0::AbstractArray, wt::AbstractArray; autodiff = :finite, kwargs...)
    results = levenberg_marquardt(R, p0; kwargs...)
    p = minimizer(results)
    return LsqFitResult(p, value!(R, p), jacobian!(R, p), converged(results), wt)
end

"""
    curve_fit(model, xdata, ydata, p0) -> fit
    curve_fit(model, xdata, ydata, wt, p0) -> fit

Fit data to a non-linear `model`. `p0` is an initial model parameter guess (see Example),
and `wt` is an optional array of weights.
The return object is a composite type (`LsqFitResult`), with some interesting values:

* `fit.resid` : residuals = vector of residuals
* `fit.jacobian` : estimated Jacobian at solution

additionally, it is possible to quiry the degrees of freedom with

* `dof(fit)`
* `coef(fit)`

## Example
```julia
# a two-parameter exponential model
# x: array of independent variables
# p: array of model parameters
model(x, p) = p[1]*exp.(-x.*p[2])

# some example data
# xdata: independent variables
# ydata: dependent variable
xdata = range(0, stop=10, length=20)
ydata = model(xdata, [1.0 2.0]) + 0.01*randn(length(xdata))
p0 = [0.5, 0.5]

fit = curve_fit(model, xdata, ydata, p0)
```
"""
function curve_fit end

function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, p0::AbstractArray; inplace = false, kwargs...)
    check_data_health(xdata, ydata)
    # construct the cost function
    T = eltype(ydata)

    if inplace
        f! = (F,p)  -> (model(F,xdata,p); @. F = F - ydata)
        lmfit(f!, p0, T[], ydata; kwargs...)
    else
        f = (p) -> model(xdata, p) - ydata
        lmfit(f, p0, T[]; kwargs...)
    end
end

function curve_fit(model, jacobian_model,
            xdata::AbstractArray, ydata::AbstractArray, p0::AbstractArray; inplace = false, kwargs...)
    check_data_health(xdata, ydata)

    T = eltype(ydata)

    if inplace
        f! = (F,p) -> (model(F,xdata,p); @. F = F - ydata)
        g! = (G,p)  -> jacobian_model(G, xdata, p)
        lmfit(f!, g!, p0, T[], copy(ydata); kwargs...)
    else
        f = (p) -> model(xdata, p) - ydata
        g = (p) -> jacobian_model(xdata, p)
        lmfit(f, g, p0, T[]; kwargs...)
    end
end

function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, wt::AbstractArray{T}, p0::AbstractArray; inplace = false, kwargs...) where T
    check_data_health(xdata, ydata)
    # construct a weighted cost function, with a vector weight for each ydata
    # for example, this might be wt = 1/sigma where sigma is some error term
    u = sqrt.(wt) # to be consistant with the matrix form

    if inplace
        f! = (F,p) -> (model(F,xdata,p); @. F = u*(F - ydata))
        lmfit(f!, p0, wt, ydata; kwargs...)
    else
        f = (p)  -> u .* ( model(xdata, p) - ydata )
        lmfit(f,p0,wt; kwargs...)
    end
end

function curve_fit(model, jacobian_model,
            xdata::AbstractArray, ydata::AbstractArray, wt::AbstractArray{T}, p0; inplace = false, kwargs...) where T
    check_data_health(xdata, ydata)
    u = sqrt.(wt) # to be consistant with the matrix form

    if inplace
        f! = (F,p) -> (model(F,xdata,p); @. F = u*(F - ydata))
        g! = (G,p) -> (jacobian_model(G, xdata, p); @. G = u*G )
        lmfit(f!, g!, p0, wt, ydata; kwargs...)
    else
        f = (p) -> u .* ( model(xdata, p) - ydata )
        g = (p) -> u .* ( jacobian_model(xdata, p) )
        lmfit(f, g, p0, wt; kwargs...)
    end
end

function curve_fit(model, xdata::AbstractArray, ydata::AbstractArray, wt::AbstractArray{T,2}, p0; kwargs...) where T
    check_data_health(xdata, ydata)

    # as before, construct a weighted cost function with where this
    # method uses a matrix weight.
    # for example: an inverse_covariance matrix

    # Cholesky is effectively a sqrt of a matrix, which is what we want
    # to minimize in the least-squares of levenberg_marquardt()
    # This requires the matrix to be positive definite
    u = cholesky(wt).U

    f(p) = u * ( model(xdata, p) - ydata )
    lmfit(f,p0,wt; kwargs...)
end

function curve_fit(model, jacobian_model,
            xdata::AbstractArray, ydata::AbstractArray, wt::AbstractArray{T,2}, p0; kwargs...) where T
    check_data_health(xdata, ydata)

    u = cholesky(wt).U

    f(p) = u * ( model(xdata, p) - ydata )
    g(p) = u * ( jacobian_model(xdata, p) )
    lmfit(f, g, p0, wt; kwargs...)
end

function estimate_covar(fit::LsqFitResult)
    # computes covariance matrix of fit parameters
    J = fit.jacobian

    if isempty(fit.wt)
        r = fit.resid

        # compute the covariance matrix from the QR decomposition
        Q, R = qr(J)
        Rinv = inv(R)
        covar = Rinv*Rinv'*mse(fit)
    else
        covar = inv(J'*J)
    end

    return covar
end

function StatsBase.stderror(fit::LsqFitResult; rtol::Real=NaN, atol::Real=0)
    # computes standard error of estimates from
    #   fit   : a LsqFitResult from a curve_fit()
    #   atol  : absolute tolerance for approximate comparisson to 0.0 in negativity check
    #   rtol  : relative tolerance for approximate comparisson to 0.0 in negativity check
    covar = estimate_covar(fit)
    # then the standard errors are given by the sqrt of the diagonal
    vars = diag(covar)
    vratio = minimum(vars)/maximum(vars)
    if !isapprox(vratio, 0.0, atol=atol, rtol=isnan(rtol) ? Base.rtoldefault(vratio, 0.0, 0) : rtol) && vratio < 0.0
        error("Covariance matrix is negative for atol=$atol and rtol=$rtol")
    end
    return sqrt.(abs.(vars))
end

function margin_error(fit::LsqFitResult, alpha=0.05; rtol::Real=NaN, atol::Real=0)
    # computes margin of error at alpha significance level from
    #   fit   : a LsqFitResult from a curve_fit()
    #   alpha : significance level, e.g. alpha=0.05 for 95% confidence
    #   atol  : absolute tolerance for approximate comparisson to 0.0 in negativity check
    #   rtol  : relative tolerance for approximate comparisson to 0.0 in negativity check
    std_errors = stderror(fit; rtol=rtol, atol=atol)
    dist = TDist(dof(fit))
    critical_values = quantile(dist, 1 - alpha/2)
    # scale standard errors by quantile of the student-t distribution (critical values)
    return std_errors * critical_values
end

function confidence_interval(fit::LsqFitResult, alpha=0.05; rtol::Real=NaN, atol::Real=0)
    # computes confidence intervals at alpha significance level from
    #   fit   : a LsqFitResult from a curve_fit()
    #   alpha : significance level, e.g. alpha=0.05 for 95% confidence
    #   atol  : absolute tolerance for approximate comparisson to 0.0 in negativity check
    #   rtol  : relative tolerance for approximate comparisson to 0.0 in negativity check
    std_errors = stderror(fit; rtol=rtol, atol=atol)
    margin_of_errors = margin_error(fit, alpha; rtol=rtol, atol=atol)
    confidence_intervals = collect(zip(coef(fit) - margin_of_errors, coef(fit) + margin_of_errors))
end

@deprecate standard_errors(args...; kwargs...) stderror(args...; kwargs...)
@deprecate estimate_errors(fit::LsqFitResult, confidence=0.95; rtol::Real=NaN, atol::Real=0) margin_error(fit, 1-confidence; rtol=rtol, atol=atol)


text:
module Utils

using Interpolations
using HDF5
using SpecialFunctions
using LinearAlgebra # pinv

export build_interpolation,
	   _readh5_file,
	   multiple_reflections,
	   find_closest,
	   blackbody_radiation,
	   wavelength2rgb,
	   simple_moving_average,
	   exponential_moving_average,
	   savitzky_golay,
	   flatten_arrays,
	   array_arrays,
	   average_polarisation,
	   gaussian,
	   lorentzian,
	   voigtian,
	   _oscillators_input,
	   unfold_bnd,
	   Info

function Info()
    tmp1 = "\n " *
	"\n " *
    "\n Available functions from Utils module:" *
	"\n " *
	"\n     build_interpolation(X)" *
	"\n     multiple_reflections(n)" *
	"\n     find_closest(x, x0)" *
	"\n     blackbody_radiation(λ, T)" *
	"\n     wavelength2rgb(λ)" *
	"\n     simple_moving_average(x)" *
	"\n     exponential_moving_average(x)" *
	"\n     flatten_arrays(x)" *
	"\n     array_arrays(x, y)" *
	"\n     average_polarisation(pol, Xp, Xs)" *
	"\n     gaussian(x, p)" *
	"\n     lorentzian(x, p)" *
	"\n     voigtian(x, p)" *
	"\n "
	"\n     To use any of these functions type: ?Utils.function" *
	"\n "
	return println(tmp1)
end

"""

	Find the index of the closest element in an array of reals x to a given real x0.

		idx = Utils.find_closest(x, x0)

			x: Data array
			x0: Element to search in x

			idx: Index

"""
function find_closest(x::AbstractArray{T1,N1}, x0::T1) where {T1<:Real, N1}
	idmin = findmin(abs.(x .- x0))[2][1]
	return idmin
end

"""Read the HDF5 file with the database."""
function _readh5_file(x::String, y::Symbol)
    datapath = joinpath(@__DIR__, "..", "data/")
	if isequal(y,:RI)
		fid = "RefractiveIndicesDB.h5"
	elseif isequal(y,:SP)
		fid = "ExampleSpectraDB.h5"
	end
    readf = h5open(joinpath(datapath,fid),"r") do file
        read(file,x)
    end
    return readf
end

"""

	Build interpolation objects from Matrix. Uses a Gridded(Linear()) grid.

		itp = Utils.build_interpolation(X)

			X: Matrix array with first column as independent variable and
               second column as dependent one

			itp: interpolation object

"""
function build_interpolation(X::Array{T1,2}) where {T1<:Float64}
	# Sort by ascending wavelength
	X = X[sortperm(X[:,1]),:]
    knots = (sort(vec(X[:,1])),)
	itp = interpolate(knots, vec(X[:,2]), Gridded(Linear()))
    return itp
end

"""

	Computes the reflection of a sequence of indices of refraction without
    interference and absorption effects.

		y = Utils.multiple_reflections(n)

			n: Vector with the indices of refraction of the involved layers

			y: Reflection

	author: Cuchu

"""
function multiple_reflections(n::AbstractArray{T1}) where {T1<:Real}
	r1::Float64 = 0.0
	r::Float64 = 0.0
	for i = 1:length(n)-1
		r2 = ((n[i] - n[i+1])/(n[i] + n[i+1]))^2
		r = (r1 + r2 - 2.0*r1*r2)/(1.0 - r1*r2)
		r1 = r
	end
	return r
end

"""

	Returns the spectral radiance (W·sr^-1·m^-3), given the wavelength (m)
    and the absolute temperature (K).

		B = Utils.blackbody_radiation(λ, T)

			λ: range of wavelengths [m]
			T: range of temperatures [K]

			B: spectral radiance [W·sr^-1·m^-3]

	Source: https://en.wikipedia.org/wiki/Planck%27s_law

"""
function blackbody_radiation(λ::AbstractArray{T1}, T::AbstractArray{T1}) where {T1<:Real}
	h::Float64 = 6.626070040e-34 # Planck's constant, J·s
	c::Float64 =  299792458.0 # speed of light, m·s
	kB::Float64 = 1.38064852e-23 # Boltzmann constant, J·K^-1
	R = Array{Float64,2}(undef,length(λ),length(T))
	num::Array = 2.0*h*c*c./λ.^5
	den::Array = h*c./(kB.*λ)
	for i in eachindex(T)
		R[:,i] .= num./expm1.(den./T[i])
	end
	return R
end

"""

	Returns an RGB matrix-color based on the wavelength input in nm.

		RBG = Utils.wavelength2rgb(λ)

			λ: range of wavelengths [nm]

	Wavelength range: 380nm and 780nm, black otherwise.

	Source: https://stackoverflow.com/questions/3407942/rgb-values-of-visible-spectrum

"""
function wavelength2rgb(λ::AbstractArray{<:Real})
    Γ::Float64 = 0.8
    # warm up
    R = Array{Float64, 1}(undef, length(λ))
    G = similar(R)
    B = similar(R)
    intensityCorrection = similar(R)
    R[380.0 .<= λ .<= 440.0] .= -(λ[380.0 .<= λ .<= 440.0] - 440.)./(440. - 380.)
    G[380.0 .<= λ .<= 440.0] .= 0.0
    B[380.0 .<= λ .<= 440.0] .= 1.0 .+ (0. * λ[380.0 .<= λ .<= 440.0])
    R[440.0 .< λ .<= 490.0] .= 0.0
    G[440.0 .< λ .<= 490.0] .= ((λ[440.0 .< λ .<= 490.0] - 440.)./(490. - 440.))
    B[440.0 .< λ .<= 490.0] .= 1.0
    R[490.0 .< λ .<= 510.0] .= 0.0
    G[490.0 .< λ .<= 510.0] .= 1.0
    B[490.0 .< λ .<= 510.0] .= (-(λ[490.0 .< λ .<= 510.0] - 510.)./(510. - 490.))
    R[510.0 .< λ .<= 580.0] .= ((λ[510.0 .< λ .<= 580.0] - 510.)./(580. - 510.))
    G[510.0 .< λ .<= 580.0] .= 1.0
    B[510.0 .< λ .<= 580.0] .= 0.0
    R[580.0 .< λ .<= 645.0] .= 1.0
    G[580.0 .< λ .<= 645.0] .= (-(λ[580.0 .< λ .<= 645.0] - 645.)./(645. - 580.))
    B[580.0 .< λ .<= 645.0] .= 0.0
    R[645.0 .< λ .<= 780.0] .= 1.0
    G[645.0 .< λ .<= 780.0] .= 0.0
    B[645.0 .< λ .<= 780.0] .= 0.0
    R[λ .> 780] .= 0.0
    R[λ .< 380] .= 0.0
    G[λ .> 780] .= 0.0
    G[λ .< 380] .= 0.0
    B[λ .> 780] .= 0.0
    B[λ .< 380] .= 0.0
    # LET THE INTENSITY SSS FALL OFF NEAR THE VISION LIMITS
    intensityCorrection[λ .> 700.0] .= 0.3 .+ 0.7.*(780.0 .- λ[λ .> 700.0])./(780. - 700.)
    intensityCorrection[λ .< 420.0] .= 0.3 .+ 0.7.*(λ[λ .< 420.0] - 380.0)./(420. - 380.)
    intensityCorrection[420.0 .< λ .< 700.0] .= 1.0.*λ[420.0 .< λ .< 700.0]
    # GAMMA ADJUST
    R .= (intensityCorrection.*R).^Γ
    G .= (intensityCorrection.*G).^Γ
    B .= (intensityCorrection.*B).^Γ
    # # Multiply by 255
    # R = R * 255
    # G = G * 255
    # B = B * 255
    return [R G B]
end

"""

	Smooth vector with simple moving average model.

		y = Utils.simple_moving_average(x; w=1)

			x: Vector to be smoothed out
			w: Integer with the window size

			y: Smoothed output
	Ref: https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average_(boxcar_filter)

"""
function simple_moving_average(y::AbstractVector{T1}; w::T2=1) where {T1<:Real, T2<:Int64}
    _y = [y[1]*ones(w); y; y[end]*ones(w)]
    m = float(2*w + 1)
    w = Int(w + 1)
    y_ = similar(y)
    for i in w:length(_y)-(w-1)
        y_[i-w+1] = sum(_y[i-w+1:i+w-1])
    end
    return y_ ./ m
end

"""

	Smooth vector with exponential moving average model.

		y = Utils.exponential_moving_average(x; α=1.0)

			x: Vector to be smoothed out
			α: Weight

			y: Smoothed output
	Ref: https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average

"""
function exponential_moving_average(y::AbstractVector{T1}; α::T1=1.0) where {T1<:Real}
    y_ = similar(y)
    y_[1] = y[1]
    for i in 2:length(y)
        y_[i] = α*y[i] + (1.0 - α)*y_[i-1]
    end
    return y_
end

"""

    Convert Array{Array{Real,1},1} to an Array{Float64,1}.

"""
function flatten_arrays(xin::Array{Array{T1,1},1}) where {T1<:Real}
    xout = []
    # Determine the total number of elements of xin
    for i in eachindex(1:length(xin)), j in eachindex(1:length(xin[i]))
        push!(xout, xin[i][j])
    end
    return float.(xout)
end

"""

    Convert Array{Real,1} to Array{Array{Float64,1},1}.

"""
function array_arrays(x::Array{T1,1}, _x::Array{Array{T1,1},1}) where{T1<:Real, T2<:Real}
    xfinal = deepcopy(_x)
    ki::Int64 = 1
    kf::Int64 = 0
    for i in eachindex(collect(1:length(xfinal)))
        kf += length(xfinal[i])
        xfinal[i][1] = x[ki]
        xfinal[i][2:end] = x[ki+1:kf]
        ki += length(xfinal[i])
    end
    return float.(xfinal)
end

"""

    Return the polarisation averaged spectrum.

        X = Utils.average_polarisation(pol, Xp, Xs)

            pol: indicates the polarisation
                 (pol=1.0 => p, pol=0.0 => s, 0.0<=pol<=1.0 => average)
            Xp: p/TM polarisation quantity (e.g. Rp)
            Xs: s/TE polarisation quantity (e.g. Rs)

            Xavg: unpolarized quantity = pol*Xp + (1.0 - pol)*Xs

"""
function average_polarisation(pol::T1, Xp::Array{T1}, Xs::Array{T1}) where {T1<:Float64}
	X = @. pol*Xp + (1.0 - pol)*Xs
    return X
end

"""

	Multipeak Gaussian PDF curve.

		modelg = Utils.gaussian(x, p)

			x = vector with data of the x-axis
			p = [
				  [p0], # offset
				  [A1, μ1, σ1], # first peak
				  [A2, μ2, σ2], # second peak
				  ..., # so on
				]
				p0: offset of the total curve
				Ai: amplitude of peak i
				μi: position of peak i (mean or expectation of the distribution)
				σi: standard deviation of peak i
				i = 1,...,n

	Source: https://en.wikipedia.org/wiki/Normal_distribution

"""
function gaussian(x::AbstractArray{T1}, p::Array{Array{T1,1},1}) where {T1<:Real}
    isequal(length(flatten_arrays(p)[2:end]), 3*(length(p)-1)) || throw("The inputs are not correct. For the Gaussian curve, each peak has three parameters plus the offset.")
    y = ones(length(x)).*p[1][1] # offset
    for i = 2:length(p)
        @. y += p[i][1]*exp(-((x - p[i][2])/p[i][3])^2)
    end
    return y
end

"""

	Multipeak Cauchy-Lorentz PDF curve.

		modell = Utils.lorentzian(x, p)

			x = vector with data of the x-axis
			p = [
				  [p0], # offset
				  [A1, μ1, Γ1], # first peak
				  [A2, μ2, Γ2], # second peak
				  ..., # so on
				]
				p0: offset of the total curve
				Ai: amplitude of peak i
				μi: position of peak i (location parameter)
				Γi: half-width at half-maximum (HWHM) of peak i
				i = 1,...,n

	Source: https://en.wikipedia.org/wiki/Cauchy_distribution

"""
function lorentzian(x::AbstractArray{T1}, p::Array{Array{T1,1},1}) where {T1<:Real}
    isequal(length(flatten_arrays(p)[2:end]), 3*(length(p)-1)) || throw("The inputs are not correct. For the Lorentzian curve, each peak has three parameters plus the offset.")
    y = ones(length(x)).*p[1][1] # offset
    for i = 2:length(p)
        @. y += p[i][1]/(1.0 + ((x - p[i][2])/p[i][3])^2)
    end
    return y
end

"""

	Single peak Voigt PDF curve. The function is computed using the Faddeeva's
    function through the scaled complementary error function of x.

		modelv = Utils.voigtian(x, p)

			x = vector with data of the x-axis
			p = [
				  [p0], # offset
				  [A1, μ1, σ1, Γ1], # first peak
				  [A2, μ2, σ2, Γ2], # second peak
				  ..., # so on
				]
				p0: offset of the total curve
				Ai: amplitude of the curve
				μi: position of peak i (location parameter)
				σi: Gaussian standard deviation of peak i
				Γi: Lorentzian half-width at half-maximum (HWHM) of peak i
				i = 1,...,n

		The μ parameter has been included in a custom fashion to allow the
        distribution to be uncentered.

	Source: https://en.wikipedia.org/wiki/Voigt_profile

"""
function voigtian(x::AbstractArray{T1}, p::Array{Array{T1,1},1}) where {T1<:Real}
    isequal(length(flatten_arrays(p)[2:end]), 4*(length(p)-1)) || throw("The inputs are not correct. For the Voigtian curve, each peak has four parameters plus the offset.")
    y = ones(length(x)).*p[1][1] # offset
    for i = 2:length(p)
        @. y += p[i][1]*real(erfcx(-im*((x - p[i][2] + im*p[i][3])/sqrt(2)/p[i][4])))
    end
    return y
end

"""

	Build the input for oscillators models. Each index of refraction associated with
	Lorentz can be modelled with arbitrary numbers of oscillators. This function provides a
	way to construct a valid input for them.

	The input for the following models can be constructed are:
		RI.drudelorentz, RI.tauclorentz, RI.forouhibloomer,
		RI.lorentzplasmon, RI.codylorentz

		x = _oscillators_input(argin, n)

			argin: Array with two elements.
				   The first element is an array with the parameters not associated to
				   oscillators, while the second element contains an array with all the
				   parameters of the oscillators.
			n: number of parameters in each oscillator
			x: Array to input in RI modelled with Lorentz.

	For instance, the Tauc-Lorentz model accepts inputs with two parameters
	(length(argin[1]) == 2) not associated with the oscillator and three (n=3) forming
	each oscillator. Say you want to pass a three oscillator model input to the
	Tauc-Lorentz:
		argin = [
			[4.0, 3.5],  # ϵinf, Eg
			[15.0, 6.0, 1.05, # osc 1, A1, E01, Γ1
			 8.41, 4.5, 0.2, # osc 2, A2, E02, Γ2
			 1.5, 4.6, 0.1, # osc 3, A3, E03, Γ3
			],
		]
		x = _oscillators_input(argin, 3)
		x = [
			  [4.0, 3.5],  # ϵinf, Eg
			  [15.0, 6.0, 1.05], # osc 1, A1, E01, Γ1
			  [8.41, 4.5, 0.2], # osc 2, A2, E02, Γ2
			  [1.5, 4.6, 0.1], # osc 3, A3, E03, Γ3
		]
		RI.tauclorentz(x, λ)

	Notice: Careful with the number of parameteres in each model.

"""
## Build the input for oscillators models.
# Each element of x has the parameters for the corresponding oscillator.
# function _oscillators_input(argin, n::Int64)
#     x = argin[1]
#     numel::Int64 = Int(length(argin[2:end])/n + 1)
#     k::Int64 = 0
#     for i = 2:numel
#         push!(x, argin[i+k:n+i-1+k])
#         k += n-1
#     end
#     return x
# end
function _oscillators_input(argin, n::Int64)
    x = [argin[1]]
    numel::Int64 = Int(length(argin[2])/n)
    k::Int64 = 0
    for i = 1:numel
        push!(x, argin[2][i+k:n+i-1+k])
        k += n-1
    end
    return x
end

"""

	unfold_bnd(lb, ub)

Returns two arrays with the lower and upper bounds in accepted format for BBO package.

"""
function unfold_bnd(lb, ub)
    return (collect(zip(flatten_arrays(lb), flatten_arrays(ub))))
end

end # Utils


text:
immutable Hyperboloid{T} <: Model
    a::T
    α::T
    μ::T
end

Base.show(io::IO, ::Type{Hyperboloid}) = print(io, "Classical Hyperboloid")

function Base.show(io::IO, m::Hyperboloid)
    @printf(
        io,
        "Classical Hyperboloid (a = %.3f, α = %.3f, μ = %.3f)",
        m.a,
        m.α,
        m.μ,
    )
end

function initial_parameters(::Type{Hyperboloid})
    return Float64[log(1.0), log(1.0), log(1.0)]
end

canonical_parameters(m::Hyperboloid) = (m.a, m.α, m.μ)

function Hyperboloid(parameters::Vector)
    a = exp(parameters[1])
    α = exp(parameters[2])
    μ = exp(parameters[3])
    return Hyperboloid(a, α, μ)
end

@inline function predict(
    m::Hyperboloid,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, α, μ = m.a, m.α, m.μ
    u1 = x1 * (1 + α * t1)^(-μ)
    u2 = x2 * (1 + α * t2)^(-μ)
    z = a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::Hyperboloid,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, α, μ = m.a, m.α, m.μ
    u1 = x1 * (1 + α * t1)^(-μ)
    u2 = x2 * (1 + α * t2)^(-μ)
    gr = (
        (u2 - u1) * a,
        a * (x2 * -μ * (1 + α * t2)^(-μ - 1) * t2 - x1 * -μ * (1 + α * t1)^(-μ - 1) * t1) * α,
        a * (-x2 * (1 + α * t2)^(-μ) * log(1 + α * t2) + x1 * (1 + α * t1)^(-μ) * log(1 + α * t1)) * μ,
    )
    return gr
end


text:
"""
    sliceview(A, code)
    slicecopy(A, code)

Slice array `A` according to `code`, a tuple of length `ndims(A)`,
in which `:` indicates a dimension of the slices, and `*` a dimension separating them.
For example if `code = (:,*,:)` then slices are either `view(A, :,i,:)`
or `A[:,i,:]` with `i=1:size(A,2)`.
"""
function sliceview(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
    [ view(A, i...) for i in iter ]
end

@doc @doc(sliceview)
function slicecopy(A::AbstractArray{T,N}, code::Tuple) where {T,N}
    N == length(code) || throw(ArgumentError("wrong code length"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(A,d) : Ref(:), Val(N))...)
    [ A[i...] for i in iter ]
end

sliceview(A::AbstractArray{T,N}) where {T,N} = sliceview(A, ntuple(i-> i==N ? (*) : (:), N))
slicecopy(A::AbstractArray{T,N}) where {T,N} = slicecopy(A, ntuple(i-> i==N ? (*) : (:), N))

"""
    copy_glue(A, code)

This is the inverse of `sliceview`, kept around only for ZygoteRules...
"""
@inline function copy_glue(A::AbstractArray{IT,N}, code::Tuple) where {IT,N}
    gluecodecheck(A, code)
    B = Array{eltype(first(A))}(undef, gluedsize(A, code))
    glue!(B, A, code)
end

@doc @doc(glue)
function glue!(B::AbstractArray{T,N}, A::AbstractArray{IT,ON}, code::Tuple) where {T,N,IT,ON}
    gluecodecheck(A, code)
    N == _ndims(A) + _ndims(first(A))  || throw(DimensionMismatch("wrong size target"))
    iter = Iterators.product(ntuple(d -> (code[d]==*) ? axes(B,d) : Ref(:), Val(N))...)
    for ind in iter
        decolon = filter(j -> !(j isa Colon), ind)
        copyto!(view(B, ind...), A[decolon...])
    end
    B
end

function gluecodecheck(A::AbstractArray, code::Tuple)
    colons = countcolons(code)
    inner = _ndims(first(A))
    outer = _ndims(A)
    colons == inner || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(first(A)) == $inner not cannot be glued with code = $(pretty(code))"))
    length(code) - colons == outer || throw(DimensionMismatch("wrong number of dimensions: " *
        "ndims(A) == $outer cannot be glued with code = $(pretty(code))"))
end

countcolons(code::Tuple{}) = 0
countcolons(code::Tuple) = Int(first(code) isa Colon) + countcolons(Base.tail(code))

_ndims(A) = ndims(A)
_ndims(A::Tuple) = 1

@generated function gluedsize(A::AbstractArray{IT, N}, code::Tuple) where {IT, N}
    list = Any[]
    dout = 1
    din = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(first(A),$din) ))
            din += 1
        else
            push!(list, :( size(A,$dout) ))
            dout += 1
        end
    end
    :( ($(list...),) )
end

using ChainRulesCore # TODO add tests?

# Rules moved from SliceView.jl

ChainRulesCore.rrule(::typeof(sliceview), A::AbstractArray, code::Tuple) =
    sliceview(A, code), Δ -> (NoTangent(), copy_glue(Δ, code), NoTangent())

ChainRulesCore.rrule(::typeof(copy_glue), A::AbstractArray, code::Tuple) =
    copy_glue(A, code), Δ -> (NoTangent(), sliceview(Δ, code), NoTangent())