text:
module TimeSchemes

"""
  function Crank-nicolson(M, K, F, v0, t; outputNodes = [])

Solve a linear first order in time system of ODES using the crank-nicolson discretization scheme.

The form of the system is as follows:
    M v'(t) + K v(t) = F(t)
    v(0) = v_0
"""
function crank_nicolson(M, K, F, v0, t; outputNodes = [])
    N = length(v0) # number of equations in the system
    T = length(t)  # number of discretization times

    # ASSUMPTION : t is a equidistributed vector
    dt = t[2] - t[1]
    V = []
    if length(outputNodes) != 0
        V = zeros(length(outputNodes), T)
    else
        V = zeros(N, T) # solution to return
        outputNodes = 1:N
    end
    cV = v0
    V[outputNodes, 1] = cV[outputNodes]
    MI = M + dt/2 * K
    MO = M - dt/2 * K
    for n = 2:T
        cV = MI \ convert(Array{Float64, 1}, dt*(F[:, n-1]+F[:,n])/2 + MO* cV)
        V[outputNodes, n] = cV[outputNodes]
    end
    V, cV
end

"""
  function implicit_euler(M, K, F, v0, t; outputNodes = [])

Solve a linear first order in time system of ODES using an implicit Euler discretization scheme.

The form of the system is as follows:
    M v'(t) + K v(t) = F(t)
    v(0) = v_0
"""
function implicit_euler(M, K, F, v0, t; outputNodes = [])
    N = length(v0) # number of equations in the system
    T = length(t)  # number of discretization times

    # ASSUMPTION : t is a equidistributed vector
    dt = t[2] - t[1]
    V = []
    if length(outputNodes) != 0
        V = zeros(length(outputNodes), T)
    else
        V = zeros(N, T) # solution to return
        outputNodes = 1:N
    end
    cV = v0
    V[outputNodes, 1] = cV[outputNodes]
    MI = M + dt * K
    for n = 2:T
        cV = MI \ convert(Array{Float64, 1}, dt*F[:, n] + M * cV)
        V[outputNodes, n] = cV[outputNodes]
    end
    V, cV
end

"""
  function explicit_euler(M, K, F, v0, t; outputNodes = [])

Solve a linear first order in time system of ODES using an implicit Euler discretization scheme.

The form of the system is as follows:
    M v'(t) + K v(t) = F(t)
    v(0) = v_0
"""
function explicit_euler(M, K, F, v0, t; outputNodes = [])
    N = length(v0) # number of equations in the system
    T = length(t)  # number of discretization times

    # ASSUMPTION : t is a equidistributed vector
    dt = t[2] - t[1]
    V = []
    if length(outputNodes) != 0
        V = zeros(length(outputNodes), T)
    else
        V = zeros(N, T) # solution to return
        outputNodes = 1:N
    end
    cV = v0
    V[outputNodes, 1] = cV[outputNodes]
    MI = M - dt * K
    for n = 2:T
        cV = M \ convert(Array{Float64, 1}, dt*F[:, n-1] + MI * cV)
        V[outputNodes, n] = cV[outputNodes]
    end
    V, cV
end

"""
  function newmark(M, B, K, F, v0, v1, t; β = 0.25, γ = 0.5, outputNodes = [])

Solve a second order in time system of ODES using a Newmark discretization scheme.

The form of the system is as follows:
    v(0) = v_0,   v'(0) = v_1
"""
function newmark(M, B, K, F, v0, v1, t; β = 0.25, γ = 0.5, outputNodes = [])
  # β and γ are the parameters of the scheme

  N = length(v0) # number of equations in the system
  T = length(t)  # number of discretization times

  # ASSUMPTION : t is a equidistributed vector
  dt = t[2] - t[1]
  V = []
  if length(outputNodes) != 0
    V = zeros(length(outputNodes), T)
  else
    V = zeros(2 * N, T) # solution to return
    outputNodes = 1:(2*N)
  end

  p = 1:N               # first N components of the solution v
  q = (N + 1):(2 * N)   # last  N components of the solution v'

  cV = zeros(2 * N, 1)
  # initialization with initial contion
  cV[p] = v0                       # initial position
  cV[q] = v1                       # initial velocity
  a = M \ (F[:,1] - B * v1 - K * v0) # initial acceleration
  V[:, 1] = cV[outputNodes]

  MI = M + γ * dt * B + β * dt^2 * K

  # for details about this scheme see [Hughes, p. 490]
  for n = 2:T
    d = cV[p] + dt * cV[q] + dt^2 / 2 * (1 - 2 * β) * a
    v = cV[q] + (1 - γ) * dt * a
    a = MI \ convert(Array{Float64, 1}, F[:, n] - B * v - K * d)
    cV[p] = d + β * dt^2 * a
    cV[q] = v + γ * dt * a
    V[:, n] = cV[outputNodes]
  end
  V, cV
end

end # module


text:
"""
Chemical reactions.

* [`hydrolysis_ocs`](@ref)
"""
module Reactions

include("../docstring_style.jl")

using Canopy.Water: water_dissoc

export hydrolysis_ocs

"""
Calcualte the rate constant of carbonyl sulfide (OCS) abiotic hydrolysis [s^-1]
in natural water from temperature [K] and pH. Applicable range: temperature
5--30 °C and pH 4--10.

# Arguments

* `temp`: Temperature [K].
* `pH` (_optional_): pH value. Default is 7.
* `seawater` (_kwarg_): If `true`, calculate the OCS hydrolysis rate in the
  seawater (Radford-Knoery & Cutter, 1993). If `false` (default), calculate the
  OCS hydrolysis rate in the fresh water (Elliott et al., 1989).

# References

* Elliott, S., Lu, E., & Rowland, F. S. (1989). Rates and mechanisms for the
  hydrolysis of carbonyl sulfide in natural waters. _Environmental Science &
  Technology_, _23_(4), 458--461. <https://doi.org/10.1021/es00181a011>
* Radford-Knoery, J., & Cutter, G. A. (1993). Determination of carbonyl sulfide
  and hydrogen sulfide species in natural waters using specialized collection
  procedures and gas chromatography with flame photometric detection.
  _Analytical Chemistry_, _65_(8), 976--982.
  <https://doi.org/10.1021/ac00056a005>

# Examples

```jldoctest
julia> hydrolysis_ocs(288.15)
6.6018369330935074e-6

julia> hydrolysis_ocs(288.15, 9.)
3.648990496059149e-5

julia> hydrolysis_ocs(288.15, 8.2, seawater=true)
1.4069474819575507e-5
```
"""
function hydrolysis_ocs(temp, pH = 7.0; seawater::Bool = false)
    temp_ref = 298.15
    c_OH = 10.0^(pH - water_dissoc(temp))
    t_diff = 1.0 / temp - 1.0 / temp_ref

    # note: the parameter values used here differ from those in the original
    # publications due to a small correction for water pK_w
    if seawater
        # RC93 equation
        return 1.63838819502e-5 * exp(-6444.02904777 * t_diff) +
               11.7115101829 * exp(-2427.27401921 * t_diff) * c_OH
    else
        # E89 equation
        return 2.11834513803e-5 * exp(-10418.3722377 * t_diff) +
               14.16881179 * exp(-6469.11889197 * t_diff) * c_OH
    end
end

end  # module

ext: jl


text:
instructions = readlines("2a.txt")

struct Pos
    x::Any
    depth::Any
    aim::Any
end

forward(p::Pos, amount) = Pos(p.x + amount, p.depth + p.aim * amount, p.aim)
up(p::Pos, amount) = Pos(p.x, p.depth, p.aim - amount)
down(p::Pos, amount) = Pos(p.x, p.depth, p.aim + amount)

function fn_named(name)
    getfield(Main, Symbol(name))
end

pos = Pos(0, 0, 0)
for instruction in instructions
    direction, amount = split(instruction)
    amount = parse(Int, amount)
    pos = fn_named(direction)(pos, amount)
end

pos.x * pos.depth


text:
function combinepoints(Points)
    allpoints = Points[1]
    numsets  = size(Points,1)
    for i=2:numsets
        allpoints = [allpoints Points[i]]
    end
    return allpoints
end

"""
Map high dimensional points to lower dimension using classicalMDS
    get_low_dim_points(Points, dim)
    :param Points: [d,n,numsets] array of sets of points
    :param dim: the dimension to be mapped to
"""
function get_low_dim_points(Points, dim)
    allpoints = combinepoints(Points)
    dXX = calc_vec_dist(allpoints)
    lowdimpoints = transform(fit(MDS,dXX,maxoutdim=dim, distances=true))
    allpoints = combinepoints(Points)
    dXX = calc_vec_dist(allpoints)
    lowdimpoints = transform(fit(MDS,dXX,maxoutdim=dim, distances=true))
    numsets = size(Points,1)
    numpoints = size(Points[1],2)
    dimpoints = size(lowdimpoints,1)
    newpoints = zeros(dimpoints, numpoints , numsets)
    newpoints = assignpoints!(lowdimpoints, newpoints)

    return newpoints
end

function assignpoints!(lowdimpoints,newpoints)
    numsets = size(newpoints,3)
    numpoints = size(newpoints,2)
    for i=1:numsets
        newpoints[:,:,i] = lowdimpoints[:,1:numpoints]
        lowdimpoints = lowdimpoints[:,numpoints+1:end]
    end
    return newpoints
end

#bootstrap sampling from points
function bootstrapsamples(points)
    n = size(points)
    r = rand(1:n[2],n[2])
    samplepoints = zeros(n[1],n[2])
    for i=1:n[2]
        samplepoints[:,i] = points[:,r[i]]
    end
    return samplepoints
end

"""
Perform two-sample hypothesis testing
    dfr_hypothesistest(Points, dim)
    :param Points: (d,n) array of points
    :param dim: the dimension to be mapped to
"""
function fisherraotest(pdf1, pdf2, n1, n2, dFR)
    numsamples=200
    dfr_samples = ones(numsamples)
    r = rand(numsamples) #use this to select which distribution to sample from
    for i=1:numsamples
        if (r[i] > 0.5)
            sample_points1, ind = sample(pdf1,n1) #sample points from the pdf
            sample_points2, ind = sample(pdf1,n1) #sample points from the pdf
            sample_pdf1 = kde!(sample_points1) #estimate a pdf from the sample points
            sample_pdf2 = kde!(sample_points2) #estimate a pdf from the sample points
            sample_pdf1 = kde!(sample_points1)
            sample_pdf2 = kde!(sample_points2)
            #samplepoints1 = bootstrapsamples(newpoints1)
            #samplepoints2 = bootstrapsamples(newpoints1)
            dfr_samples[i] = fisherraodistance(sample_pdf1, sample_pdf2, sample_points1, sample_points2)
        else
            #samplepoints1 = bootstrapsamples(newpoints2)
            #samplepoints2 = bootstrapsamples(newpoints2)
            sample_points1, ind = sample(pdf2,n2) #sample points from the pdf
            sample_points2, ind = sample(pdf2,n2) #sample points from the pdf
            sample_pdf1 = kde!(sample_points1) #estimate a pdf from the sample points
            sample_pdf2 = kde!(sample_points2) #estimate a pdf from the sample points
            sample_pdf1 = kde!(sample_points1)
            sample_pdf2 = kde!(sample_points2)
            sample_pdf1 = kde!(sample_points1)
            sample_pdf2 = kde!(sample_points2)
            dfr_samples[i] = fisherraodistance(sample_pdf1, sample_pdf2, sample_points1, sample_points2)
        end
    end
    pvalue = length(dfr_samples[dfr_samples.>dFR]) / numsamples
    return [pvalue, dfr_samples]
end

"""
Calculate Euclidean pairwise distances between a set of points
     calc_vec_dist(allpoints)
    :param allpoints: (d,n) array of points
"""
function calc_vec_dist(allpoints)
    N = size(allpoints, 2)
    dXX = zeros(N,N)
    for i = 1:N
        v1 = allpoints[:,i]
        for j=1:N
            v2 = allpoints[:,j]
            dXX[i,j] = sqrt((sum(v1-v2).^2))
        end
    end
    return dXX
end

"""
Estimates the Fisher-Rao Distance between two densities
     fisherraodistance(pdf1, pdf2, points1, points2)
    :params pdf1, pdf2: Ball Density estimates returned from kde!
    :params points1, points2: two [d,n] sets of points being compared
"""
function fisherraodistance(pdf1, pdf2, points1, points2)
    n1 = size(points1,2)
    likelihood_points1_pdf1 = evaluateDualTree(pdf1, points1)
    likelihood_points1_pdf2 = evaluateDualTree(pdf2, points1)
    IP = sum( sqrt.(likelihood_points1_pdf2 ./ likelihood_points1_pdf1) )/n1
    IP = IP > 1 ? 1 : IP
    dfr1 = acos( IP )

    n2 = size(points2,2)
    likelihood_points2_pdf2 = evaluateDualTree(pdf2, points2)
    likelihood_points2_pdf1 = evaluateDualTree(pdf1, points2)
    IP = sum( sqrt.(likelihood_points2_pdf1 ./ likelihood_points2_pdf2) )/n2
    IP = IP > 1 ? 1 : IP
    dfr2 = acos( IP)

    return fisherRaoDistance = (dfr1 + dfr2)/2
end


text:
export last, LastNotFoundException

import Base: last
import Base: show

struct LastNotFoundException <: Exception end

"""
    last(; default = nothing)

Creates a last operator, which returns an Observable that emits only
the last item emitted by the source Observable.
Sends `LastNotFoundException` error message if a given source completes without emitting a single value.

# Arguments
- `default`: an optional default value to provide if no values were emitted

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> last(), logger())
;

# output

[LogActor] Data: 3
[LogActor] Completed

```

```jldoctest
using Rocket

source = from(Int[])
subscribe!(source |> last() |> catch_error((err, obs) -> of(1)), logger())
;

# output

[LogActor] Data: 1
[LogActor] Completed
```

```jldoctest
using Rocket

source = from(Int[])
subscribe!(source |> last(default = 1), logger())
;

# output

[LogActor] Data: 1
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
last(; default::D = nothing) where D = LastOperator{D}(default)

struct LastOperator{D} <: InferableOperator
    default :: D
end

function on_call!(::Type{L}, ::Type{R}, operator::LastOperator, source) where { L, R }
    return proxy(R, source, LastProxy{R}(operator.default))
end

operator_right(operator::LastOperator{Nothing}, ::Type{L}) where L        = L
operator_right(operator::LastOperator{D},       ::Type{L}) where { L, D } = Union{L, D}

struct LastProxy{L} <: ActorProxy
    default :: Union{L, Nothing}
end

actor_proxy!(::Type, proxy::LastProxy{L}, actor::A) where { L, A } = LastActor{L, A}(proxy.default, actor)

mutable struct LastActor{L, A} <: Actor{L}
    last   :: Union{L, Nothing}
    actor  :: A
end

function on_next!(actor::LastActor{L}, data::L) where L
    actor.last = data
end

function on_error!(actor::LastActor, err)
    error!(actor.actor, err)
end

function on_complete!(actor::LastActor)
    if actor.last !== nothing
        next!(actor.actor, actor.last)
        complete!(actor.actor)
    else
        error!(actor.actor, LastNotFoundException())
    end
end

Base.show(io::IO, ::LastOperator)         = print(io, "LastOperator()")
Base.show(io::IO, ::LastProxy{L}) where L = print(io, "LastProxy($L)")
Base.show(io::IO, ::LastActor{L}) where L = print(io, "LastActor($L)")


text:
using REPL: REPL, LineEdit

using REPL
using Unitful: s
using FileIO: load, save, loadstreaming, savestreaming
import LibSndFile
using PortAudio
using SampledSignals
using EllipsisNotation

const START_TIME = Ref{Int}(1)

function get_song()
    return load(joinpath(pkgdir(MinionsREPL), "minions.ogg"))
end

const MINIONS_SONG = get_song()

function press(step = 40_000)
    stream = PortAudioStream(0, 2)
    start = START_TIME[]
    try
        if start+step > 1357248
            START_TIME[] = stop = mod1(start+step, 1357248)
            write(stream, MINIONS_SONG[start:end, :])
            write(stream, MINIONS_SONG[1:stop, :])
        else
            START_TIME[] += step
            write(stream, MINIONS_SONG[start:start+step, :])
        end
    finally
        close(stream)
    end
end

function REPL.LineEdit.match_input(k::Dict{Char}, s::Union{Nothing,LineEdit.MIState}, term::Union{LineEdit.AbstractTerminal,IOBuffer}=LineEdit.terminal(s), cs::Vector{Char}=Char[], keymap::Dict{Char} = k)
    # if we run out of characters to match before resolving an action,
    # return an empty keymap function
    eof(term) && return (s, p) -> :abort
    c = read(term, Char)
    press()
    # Ignore any `wildcard` as this is used as a
    # placeholder for the wildcard (see normalize_key("*"))
    c == LineEdit.wildcard && return (s, p) -> :ok
    push!(cs, c)
    key = haskey(k, c) ? c : LineEdit.wildcard

    # if we don't match on the key, look for a default action then fallback on 'nothing' to ignore
    return LineEdit.match_input(get(k, key, nothing), s, term, cs, keymap)
end


text:
mutable struct SVM{T<:Real} <: Model
    C::T
    state

    SVM(; C::Real = 1, T = Float32) = new{T}(T(C))
end

get_kernel(::Type{Linear}) = LIBSVM.Kernel.Linear
get_kernel(::Type{Gaussian}) = LIBSVM.Kernel.RadialBasis
parameters(m::SVM) = (C = m.C, )
extract_scores(model::SVM, K) = LIBSVM.svmpredict(model.state, K)[2][2, :]


text:
"""
`Plot(graphics; nargs...)` → `Plot`

Description
============

Decoration for objects that are `GraphicsArea` (or `Canvas`).
It is used to surround the inner `GraphicsArea` object with
additional information such as a title, border, and axis labels.

Usage
======

    Plot(graphics; title = "", xlabel = "", ylabel = "", border = :solid, margin = 3, padding = 1, labels = true)

    Plot(x, y, canvastype; title = "", width = 40, height = 15, border = :solid, xlim = [0, 0], ylim = [0, 0], margin = 3, padding = 1, labels = true, grid = true)

Arguments
==========

- **`graphics`** : The `GraphicsArea` (e.g. a subtype of `Canvas`) that the plot should decorate

- **`x`** : The horizontal dimension for each point.

- **`y`** : The vertical dimension for each point.

- **`canvastype`** : The type of canvas that should be used for drawing.

- **`title`** : Text to display on the top of the plot.

- **`xlabel`** : Text to display on the x axis of the plot

- **`ylabel`** : Text to display on the y axis of the plot

- **`width`** : Number of characters per row that should be used for plotting.

- **`height`** : Number of rows that should be used for plotting. Not applicable to `barplot`.

- **`border`** : The style of the bounding box of the plot. Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`xlim`** : Plotting range for the x coordinate. `[0, 0]` stands for automatic.

- **`ylim`** : Plotting range for the y coordinate. `[0, 0]` stands for automatic.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

- **`grid`** : Can be used to hide the gridlines at the origin

Methods
========

- `title!(plot::Plot, title::String)`

- `xlabel!(plot::Plot, xlabel::String)`

- `ylabel!(plot::Plot, xlabel::String)`

- `annotate!(plot::Plot, where::Symbol, value::String)`

- `annotate!(plot::Plot, where::Symbol, row::Int, value::String)`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

see also
=========

`scatterplot`, `lineplot`, `BarplotGraphics`, `BrailleCanvas`, `BlockCanvas`, `AsciiCanvas`
"""
type Plot{T<:GraphicsArea}
    graphics::T
    title::String
    xlabel::String
    ylabel::String
    margin::Int
    padding::Int
    border::Symbol
    labels_left::Dict{Int,String}
    colors_left::Dict{Int,Symbol}
    labels_right::Dict{Int,String}
    colors_right::Dict{Int,Symbol}
    decorations::Dict{Symbol,String}
    colors_deco::Dict{Symbol,Symbol}
    show_labels::Bool
    autocolor::Int
end

function Plot{T<:GraphicsArea}(
        graphics::T;
        title::AbstractString = "",
        xlabel::AbstractString = "",
        ylabel::AbstractString = "",
        border::Symbol = :solid,
        margin::Int = 3,
        padding::Int = 1,
        labels = true)
    rows = nrows(graphics)
    cols = ncols(graphics)
    labels_left = Dict{Int,String}()
    colors_left = Dict{Int,Symbol}()
    labels_right = Dict{Int,String}()
    colors_right = Dict{Int,Symbol}()
    decorations = Dict{Symbol,String}()
    colors_deco = Dict{Symbol,Symbol}()
    Plot{T}(graphics, title, xlabel, ylabel,
            margin, padding, border,
            labels_left, colors_left, labels_right, colors_right,
            decorations, colors_deco, labels, 0)
end

function Plot{C<:Canvas, F<:AbstractFloat}(
        X::Vector{F}, Y::Vector{F}, ::Type{C} = BrailleCanvas;
        title::AbstractString = "",
        width::Int = 40,
        height::Int = 15,
        border::Symbol = :solid,
        xlim::Vector = [0.,0.],
        ylim::Vector = [0.,0.],
        margin::Int = 3,
        padding::Int = 1,
        labels::Bool = true,
        grid::Bool = true)
    length(xlim) == length(ylim) == 2 || throw(ArgumentError("xlim and ylim must only be vectors of length 2"))
    margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
    length(X) == length(Y) || throw(DimensionMismatch("X and Y must be the same length"))
    width = max(width, 5)
    height = max(height, 2)

    min_x, max_x = extend_limits(X, xlim)
    min_y, max_y = extend_limits(Y, ylim)
    origin_x = min_x
    origin_y = min_y
    p_width = max_x - origin_x
    p_height = max_y - origin_y

    canvas = C(width, height,
               origin_x = origin_x, origin_y = origin_y,
               width = p_width, height = p_height)
    new_plot = Plot(canvas, title = title, margin = margin,
                    padding = padding, border = border, labels = labels)

    min_x_str = string(roundable(min_x) ? round(Int, min_x, RoundNearestTiesUp) : min_x)
    max_x_str = string(roundable(max_x) ? round(Int, max_x, RoundNearestTiesUp) : max_x)
    min_y_str = string(roundable(min_y) ? round(Int, min_y, RoundNearestTiesUp) : min_y)
    max_y_str = string(roundable(max_y) ? round(Int, max_y, RoundNearestTiesUp) : max_y)
    annotate!(new_plot, :l, 1, max_y_str)
    annotate!(new_plot, :l, height, min_y_str)
    annotate!(new_plot, :bl, min_x_str)
    annotate!(new_plot, :br, max_x_str)
    if grid
        if min_y < 0 < max_y
            for i in linspace(min_x, max_x, width * x_pixel_per_char(typeof(canvas)))
                points!(new_plot, i, 0., :white)
            end
        end
        if min_x < 0 < max_x
            for i in linspace(min_y, max_y, height * y_pixel_per_char(typeof(canvas)))
                points!(new_plot, 0., i, :white)
            end
        end
    end
    new_plot
end

function next_color!{T<:GraphicsArea}(plot::Plot{T})
    cur_color = color_cycle[plot.autocolor + 1]
    plot.autocolor = ((plot.autocolor + 1) % length(color_cycle))
    cur_color
end

"""
`title(plot) →  String`

Returns the current title of the given plot.
Alternatively, the title can be changed with `title!`.
"""
function title{T<:GraphicsArea}(plot::Plot{T})
    plot.title
end

"""
`title!(plot, newtitle) →  Plot`

Sets a new title for the given plot.
Alternatively, the current title can be
queried using `title`.
"""
function title!{T<:GraphicsArea}(plot::Plot{T}, title::AbstractString)
    plot.title = title
    plot
end

"""
`xlabel(plot) →  String`

Returns the current label for the x-axis.
Alternatively, the x-label can be changed with `xlabel!`
"""
function xlabel{T<:GraphicsArea}(plot::Plot{T})
    plot.xlabel
end

"""
`xlabel!(plot, newlabel) →  Plot`

Sets a new x-label for the given plot.
Alternatively, the current label can be
queried using `xlabel`
"""
function xlabel!{T<:GraphicsArea}(plot::Plot{T}, xlabel::AbstractString)
    plot.xlabel = xlabel
    plot
end

"""
`ylabel(plot) →  String`

Returns the current label for the y-axis.
Alternatively, the y-label can be changed with `ylabel!`
"""
function ylabel{T<:GraphicsArea}(plot::Plot{T})
    plot.ylabel
end

"""
`ylabel!(plot, newlabel) →  Plot`

Sets a new y-label for the given plot.
Alternatively, the current label can be
queried using `ylabel`
"""
function ylabel!{T<:GraphicsArea}(plot::Plot{T}, ylabel::AbstractString)
    plot.ylabel = ylabel
    plot
end

"""
`annotate!(plot, where, value[, color])`

`annotate!(plot, where, row, value[, color])`

This method is responsible for the setting
all the textual decorations of a plot.

Note that `where` can be any of: `:tl` (top-left),
`:t` (top-center), `:tr` (top-right),
`:bl` (bottom-left), `:b` (bottom-center),
`:br` (bottom-right), `:l` (left), `:r` (right).

If `where` is either `:l`, or `:r`, then `row`
can be between 1 and the number of character rows
of the plots canvas.
"""
function annotate!{T<:GraphicsArea}(plot::Plot{T}, where::Symbol, value::AbstractString, color::Symbol)
    where == :t || where == :b || where == :l || where == :r || where == :tl || where == :tr || where == :bl || where == :br || throw(ArgumentError("Unknown location: try one of these :tl :t :tr :bl :b :br"))
    if where == :l || where == :r
        for row = 1:nrows(plot.graphics)
            if where == :l
                if(!haskey(plot.labels_left, row) || plot.labels_left[row] == "")
                    plot.labels_left[row] = value
                    plot.colors_left[row] = color
                    return plot
                end
            elseif where == :r
                if(!haskey(plot.labels_right, row) || plot.labels_right[row] == "")
                    plot.labels_right[row] = value
                    plot.colors_right[row] = color
                    return plot
                end
            end
        end
    else
        plot.decorations[where] = value
        plot.colors_deco[where] = color
        return plot
    end
end

function annotate!{T<:GraphicsArea}(plot::Plot{T}, where::Symbol, value::AbstractString; color::Symbol=:white)
    annotate!(plot, where, value, color)
end

function annotate!{T<:GraphicsArea}(plot::Plot{T}, where::Symbol, row::Int, value::AbstractString, color::Symbol)
    if where == :l
        plot.labels_left[row] = value
        plot.colors_left[row] = color
    elseif where == :r
        plot.labels_right[row] = value
        plot.colors_right[row] = color
    else
        throw(ArgumentError("Unknown location: try one of these :l :r"))
    end
    plot
end

function annotate!{T<:GraphicsArea}(plot::Plot{T}, where::Symbol, row::Int, value::AbstractString; color::Symbol=:white)
    annotate!(plot, where, row, value, color)
end

function lines!{T<:Canvas}(plot::Plot{T}, args...; vars...)
    lines!(plot.graphics, args...; vars...)
    plot
end

function pixel!{T<:Canvas}(plot::Plot{T}, args...; vars...)
    pixel!(plot.graphics, args...; vars...)
    plot
end

function points!{T<:Canvas}(plot::Plot{T}, args...; vars...)
    points!(plot.graphics, args...; vars...)
    plot
end

function print_title(io::IO, padding::AbstractString, title::AbstractString; p_width::Int = 0)
    if title != ""
        offset = round(Int, p_width / 2 - length(title) / 2, RoundNearestTiesUp)
        offset = offset > 0 ? offset: 0
        tpad = repeat(" ", offset)
        print_with_color(:white, io, padding, tpad, title, "\n")
    end
end

function print_border_top(io::IO, padding::AbstractString, length::Int, border::Symbol = :solid)
    b = bordermap[border]
    border == :none || print_with_color(:white, io, padding, b[:tl], repeat(b[:t], length), b[:tr])
end

function print_border_bottom(io::IO, padding::AbstractString, length::Int, border::Symbol = :solid)
    b = bordermap[border]
    border == :none || print_with_color(:white, io, padding, b[:bl], repeat(b[:b], length), b[:br])
end

function Base.show(io::IO, p::Plot)
    b = UnicodePlots.bordermap[p.border]
    c = p.graphics
    border_length = ncols(c)

    # get length of largest strings to the left and right
    max_len_l = p.show_labels && !isempty(p.labels_left)  ? maximum([length(string(l)) for l in values(p.labels_left)]) : 0
    max_len_r = p.show_labels && !isempty(p.labels_right) ? maximum([length(string(l)) for l in values(p.labels_right)]) : 0
    if p.show_labels && p.ylabel != ""
        max_len_l += length(p.ylabel) + 1
    end

    # offset where the plot (incl border) begins
    plot_offset = max_len_l + p.margin + p.padding

    # padding-string from left to border
    plot_padding = repeat(" ", p.padding)

    # padding-string between labels and border
    border_padding = repeat(" ", plot_offset)

    # plot the title and the top border
    print_title(io, border_padding, p.title, p_width = border_length)
    if p.show_labels
        topleft_str  = get(p.decorations, :tl, "")
        topleft_col  = get(p.colors_deco, :tl, :white)
        topmid_str   = get(p.decorations, :t, "")
        topmid_col   = get(p.colors_deco, :t, :white)
        topright_str = get(p.decorations, :tr, "")
        topright_col = get(p.colors_deco, :tr, :white)
        if topleft_str != "" || topright_str != "" || topmid_str != ""
            topleft_len  = length(topleft_str)
            topmid_len   = length(topmid_str)
            topright_len = length(topright_str)
            print_with_color(topleft_col, io, border_padding, topleft_str)
            cnt = round(Int, border_length / 2 - topmid_len / 2 - topleft_len, RoundNearestTiesUp)
            pad = cnt > 0 ? repeat(" ", cnt) : ""
            print_with_color(topmid_col, io, pad, topmid_str)
            cnt = border_length - topright_len - topleft_len - topmid_len + 2 - cnt
            pad = cnt > 0 ? repeat(" ", cnt) : ""
            print_with_color(topright_col, io, pad, topright_str, "\n")
        end
    end
    print_border_top(io, border_padding, border_length, p.border)
    print(io, repeat(" ", max_len_r), plot_padding, "\n")

    # compute position of ylabel
    ylabRow = round(nrows(c) / 2, RoundNearestTiesUp)

    # plot all rows
    for row in 1:nrows(c)
        # Current labels to left and right of the row and their length
        left_str  = get(p.labels_left,  row, "")
        left_col  = get(p.colors_left,  row, :white)
        right_str = get(p.labels_right, row, "")
        right_col = get(p.colors_right, row, :white)
        left_len  = length(left_str)
        right_len = length(right_str)
        # print left annotations
        print(io, repeat(" ", p.margin))
        if p.show_labels
            if row == ylabRow
                # print ylabel
                print_with_color(:white, io, p.ylabel)
                print(io, repeat(" ", max_len_l - length(p.ylabel) - left_len))
            else
                # print padding to fill ylabel length
                print(io, repeat(" ", max_len_l - left_len))
            end
            # print the left annotation
            print_with_color(left_col, io, left_str)
        end
        # print left border
        print_with_color(:white, io, plot_padding, b[:l])
        # print canvas row
        printrow(io, c, row)
        #print right label and padding
        print_with_color(:white, io, b[:r])
        if p.show_labels
            print(io, plot_padding)
            print_with_color(right_col, io, right_str)
            print(io, repeat(" ", max_len_r - right_len))
        end
        print(io, "\n")
    end

    # draw bottom border and bottom labels
    print_border_bottom(io, border_padding, border_length, p.border)
    print(io, repeat(" ", max_len_r), plot_padding, "\n")
    if p.show_labels
        botleft_str  = get(p.decorations, :bl, "")
        botleft_col  = get(p.colors_deco, :bl, :white)
        botmid_str   = get(p.decorations, :b, "")
        botmid_col   = get(p.colors_deco, :b, :white)
        botright_str = get(p.decorations, :br, "")
        botright_col = get(p.colors_deco, :br, :white)
        if botleft_str != "" || botright_str != "" || botmid_str != ""
            botleft_len  = length(botleft_str)
            botmid_len   = length(botmid_str)
            botright_len = length(botright_str)
            print_with_color(botleft_col, io, border_padding, botleft_str)
            cnt = round(Int, border_length / 2 - botmid_len / 2 - botleft_len, RoundNearestTiesUp)
            pad = cnt > 0 ? repeat(" ", cnt) : ""
            print_with_color(botmid_col, io, pad, botmid_str)
            cnt = border_length - botright_len - botleft_len - botmid_len + 2 - cnt
            pad = cnt > 0 ? repeat(" ", cnt) : ""
            print_with_color(botright_col, io, pad, botright_str, "\n")
        end
        # abuse the print_title function to print the xlabel. maybe refactor this
        print_title(io, border_padding, p.xlabel, p_width = border_length)
    end
end


text:
using PTChain
using MHChain

function init_table(col_sums)
    n = length(col_sums)
    out = zeros(n, n)
    for i in 1:n
        out[1:col_sums[i], i] = ones(col_sums[i])
    end
    out[:]
end

function score_table(table, target_vector)
    out = rearrange_vector(table)
    return sum(abs.(sum(out, 2) .- target_vector))
end

function rearrange_vector(table)
    n = Int(sqrt(length(table)))
    out = reshape(table, n, n)
    out
end

function swap_columns(table)
    out = rearrange_vector(table)
    n = size(out, 1)
    col_n = rand(1:n)
    one_indx = rand(find(!iszero, out[:, col_n]))
    zero_indx = rand(find(iszero, out[:, col_n]))
    out[one_indx, col_n], out[zero_indx, col_n] = out[zero_indx, col_n], out[one_indx, col_n]
    return out[:]
end

function testWorkers(iterations, matrix_size)
	target_vector = fill(2, matrix_size)
	start_x = init_table(target_vector)[:]
	temps = [0.01, 1, 2.5, 5]
	out = zeros(iterations)
	f(x) = score_table(x, target_vector)
	for i = 1:iterations
		out[i], x = sampleTempering(start_x, f, temps, swap_columns, 15)
	end
	out
end

function simAnneal(n, T0, reps)
    scores = [1, 2, 2, 2, 3, 4]
    temps = T0 * (0.90).^(0:5:75)
    f(x) = score_table(x, scores)
    min_val, x, stats = simulatedAnnealing(init_table(scores), f, temps, swap_columns, fill(reps, length(temps)))
end

function repSimAnneal(n, T0, reps)
    a = 100
    c = 0
    while a > 0
        val, x = simAnneal(n, T0, reps)
        if val < a
            a = val
        end
    c = c + 1
    end
    c
end

function gradDescend()
    scores = [1, 2, 2, 2, 3, 4]
    f(x) = score_table(x, scores)
    min_val, x, stats = runChain(init_table(scores), f, 0.0, swap_columns, 1000)
end

function testMany(largest_size)
	iterations = 5
	min_size = 2
	c = zeros(largest_size, iterations)
	for i = min_size:largest_size
		c[i, :] = testWorkers(iterations, i)
		print(i)
	end
	c
end

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@userplot HScatter

@recipe function f(hs::HScatter; lags=nothing, tol=1e-1, distance=Euclidean())
  # get inputs
  spatialdata = hs.args[1]
  var₁ = hs.args[2]
  var₂ = length(hs.args) == 3 ? hs.args[3] : var₁

  # lookup valid data
  X₁, z₁ = valid(spatialdata, var₁)
  X₂, z₂ = valid(spatialdata, var₂)

  # compute pairwise distance
  m, n = length(z₁), length(z₂)
  pairs = [(i,j) for j in 1:n for i in j:m]
  ds = [evaluate(distance, view(X₁,:,i), view(X₂,:,j)) for (i,j) in pairs]

  # use quartiles by default
  lags == nothing && (lags = quantile(ds, [.00,.25,.50,.75]))

  xguide --> var₁
  yguide --> var₂
  legend --> false
  aspect_ratio --> :equal
  layout --> (1, length(lags))

  for (i, lag) in enumerate(lags)
    # find indices with given lag
    match = findall(abs.(ds .- lag) .< tol)

    if isempty(match)
      @warn "no points were found with lag = $lag, skipping..."
      continue
    end

    # scatter plot coordinates
    mpairs = view(pairs, match)
    x = z₁[first.(mpairs)]
    y = z₂[last.(mpairs)]

    # plot identity line
    @series begin
      subplot := i
      seriestype := :path
      primary := false
      linestyle := :dash
      color := :black

      xmin, xmax = extrema(x)
      ymin, ymax = extrema(y)
      vmin = min(xmin, ymin)
      vmax = max(xmax, ymax)

      [vmin, vmax], [vmin, vmax]
    end

    # plot h-scatter
    @series begin
      subplot := i
      seriestype := :scatter
      title --> @sprintf "lag = %.1f, corr = %.2f" lag cor(x, y)

      x, y
    end
  end
end


text:
using ARMA, Polynomials
using LinearAlgebra
using Test
using ToeplitzMatrices

@testset "ARMA" begin

# 1) Test padded_length: rounds up to convenient size for FFT
@testset "padded_length" begin
    @test ARMA.padded_length(1000) == 1024
    @test ARMA.padded_length(16) == 16
    @test ARMA.padded_length(12) == 16
    @test ARMA.padded_length(11) == 12
    @test ARMA.padded_length(10) == 12
    @test ARMA.padded_length(9) == 10
    @test ARMA.padded_length(8) == 8
    for i=1:10
        j = rand(1:1025)
        @test ARMA.padded_length(j) >= j
    end
end

# 2) Test estimate_covariance
@testset "estimate_covariance" begin
    @test estimate_covariance([0,2,0,-2]) == [2.0, 0.0, -2.0, 0.0]
    u = randn(2+2^13)
    r = u[3:end] + u[1:end-2] + 2*u[2:end-1]
    cv = estimate_covariance(r, 20)
    @test abs.(cv[1] - 6) < 1
    @test abs.(cv[2] - 4) < 1
    @test abs.(cv[3] - 1) < 1
    for lag = 4:19
        @test abs.(cv[lag]) < 1
    end

    @test estimate_covariance(u, 10) == estimate_covariance(u, 10, div(length(u), div(length(u)-1+150,150)))
    @test length(estimate_covariance(u)) == length(u)
end

# 3) Basic tests of ARMAModel constructors
function similar_list(a::Vector, b::Vector, eps)
    @assert length(a) == length(b)
    for a1 in a
        if all(abs2.(b .- a1) .> eps^2)
            return false
        end
    end
    for b1 in b
        if all(abs2.(a .- b1) .> eps^2)
            return false
        end
    end
    true
end

@testset "ARMA_constructors" begin
    p,q = 3,3
    rs = 1 .+ (randn(q) .^ 2)
    ps = 1 .+ (randn(p) .^ 2)
    variance = 1.0
    m = ARMAModel(rs, ps, variance)
    @test m.p == p
    @test m.q == q
    n = ARMAModel(m.thetacoef, m.phicoef)
    @test m.thetacoef == n.thetacoef
    @test m.phicoef == n.phicoef
    @test similar_list(m.roots_, n.roots_, 1e-7)
    @test similar_list(m.poles, n.poles, 1e-7)
end

# 4) Now complete tests of several models that have been worked out carefully
# on paper, as well as several that are randomly created.

@testset "ARMA_representations" begin
    # Generate 6 models of fixed parameters and order (2,0), (0,2), (1,1), (1,2), (2,1), (2,2)
    thetas=Dict('A'=>[2], 'B'=>[2,2.6,.8], 'C'=>[2,1.6], 'D'=>[2,2.6,.8], 'E'=>[2,1.6], 'F'=>[2,2.6,.8])
    phis = Dict('A'=>[1,-.3,-.4], 'B'=>[1], 'C'=>[1,-.8], 'D'=>[1,-.8], 'E'=>[1,-.3,-.4], 'F'=>[1,-.3,-.4])
    EPSILON = 2e-4

    # And generate 6 models of random order and random parameters
    for model in "GHIJKL"
        # Order will be 0<=p,q <=6.
        # Use rand^(-.3) for roots/poles. Negative power ensures abs.(r)>1, and
        # the 0.3 power concentrates the values near the unit circle.
        p = rand(0:6)
        q = rand(0:6)
        if p+q==0; p=q=5; end  # Don't test ARMA(0,0) model!
        roots_ = rand(q) .^ (-.3)
        poles = rand(p) .^ (-.3)

        # Want one negative pole, if p>=3
        if p>2
            poles[end] *= -1
        end

        # Half the time, on larger-order models, make one pair roots and/or poles complex.
        if p>2 && rand(0:1) == 1
            poles = complex(poles)
            poles[1] = complex(real(poles[1]),real(poles[2]))
            poles[2] = conj(poles[1])
        end

        if q>2 && rand(0:1) == 1
            roots_ = complex(roots_)
            roots_[1] = complex(real(roots_[1]),real(roots_[2]))
            roots_[2] = conj(roots_[1])
        end

        # Scale theta by 0.7 to avoid lucky cancellations in the tests.
        thetas[model] = ARMA.polynomial_from_roots(roots_) * 0.7
        phis[model] = ARMA.polynomial_from_roots(poles)
        phis[model] *= 1.0/phis[model][1]
    end

    # Loop over all the models specified by their rational function representation
    # in thetas[] and phis[]. For each model, construct it all 3 ways (theta,phi;
    # roots, poles, and variance; or sum-of-exponentials). Verify that the resulting
    # model has the same covariance and other key properties.

    for model in "ABCDEFGHIJKL"
        thcoef = thetas[model]
        phcoef = phis[model]
        if phcoef[1] != 1.0
            thcoef /= phcoef[1]
            phcoef /= phcoef[1]
        end
        @assert phcoef[1] == 1.0
        p = length(phcoef)-1
        q = length(thcoef)-1
        # println("Testing model $model of order ARMA($p,$q).")

        m1 = ARMAModel(thcoef, phcoef)

        roots_ = roots(Polynomial(thcoef))
        poles = roots(Polynomial(phcoef))
        expbases = 1.0 ./ poles

        # We'll be working with q+1 equations to find initial values
        # of psi: the Taylor expansion coefficients of theta(z)/phi(z).
        # See BD (3.3.3) for the q+1 initial equations and  (3.3.4) for the
        # homogeneous equations beyond the first q+1. Careful with the sign conventions,
        # b/c BD uses phi(z) = 1 - (phi1*z + phi2*z^2 + ...), while I prefer a + sign.
        phpad = fill(0.0, q+1)
        if q>p
            phpad[1:p+1] = phcoef
        else
            phpad[1:q+1] = phcoef[1:q+1]
        end

        psi = copy(thcoef)
        for j=1:q
            for k=1:j
                psi[j+1] -= phpad[k+1]*psi[1+j-k]
            end
        end

        # We have to solve for the first N=max(p,q)+1 values of covariance at once.
        # For these, see BD equation (3.3.8) for the first q+1 equations and (3.3.9)
        # for the remaining (p-q), if any.
        N = 1+max(p,q)
        phN = copy(phcoef)
        if q>p
            append!(phN, fill(0.0, q-p))
        end
        A = fill(0.0, N, N)
        for i=1:N
            for j=1:N
                col = 1+abs(j-i)
                A[i,col] += phN[j]
            end
        end
        rhs3_3_8 = fill(0.0, N)
        for k=1:q+1 # here j,k are both 1 larger than in BD 3.3.8.
            for j=k:q+1
                rhs3_3_8[k] += thcoef[j]*psi[1+j-k]
            end
        end
        gamma = A \ rhs3_3_8

        m2 = ARMAModel(roots_, poles, gamma[1])

        if q<p
            covarIV = Float64[]
        else
            covarIV = gamma[1:1+q-p]
        end
        B = Array{ComplexF64}(undef, p, p)
        for r=1:p
            for c=1:p
                B[r,c] = expbases[c]^(N-p+r-1)
            end
        end
        expampls = B \ gamma[N-p+1:N]

        # If p-q>1, then you can't really work from sum-of-exponentials to ARMA,
        # because you have to assume q=p-1, leading to infinte roots, etc etc.
        # As a hack, when this degenerate case is reached, skip the sum-exp representation.
        if p-q > 1
            m3 = m2
        else
            m3 = ARMAModel(expbases, expampls, covarIV)
        end

        # A) Check that model orders are equivalent
        # Take care with model m3, b/c it never sets q<p-1 when constructing.
        @test p == m1.p
        @test p == m2.p
        @test p == m3.p
        @test q == m1.q
        @test q == m2.q
        @test q == m3.q

        # B) Check that model covariance matches
        c1 = model_covariance(m1, 100)
        c2 = model_covariance(m2, 100)
        c3 = model_covariance(m3, 100)
        c0 = c1[1]
        @test all(abs.(c1 .- c2) .< EPSILON*c0)
        @test all(abs.(c1 .- c3) .< EPSILON*c0)

        # C) Check that the initial covariances match
        # While this should be redundant with above test, let's just be sure
        NIV = max(0,q-p+1)
        if NIV>0
            @test all(abs.(m1.covarIV[1:NIV].-m2.covarIV[1:NIV]) .< EPSILON*c0)
            @test all(abs.(m1.covarIV[1:NIV].-m3.covarIV[1:NIV]) .< EPSILON*c0)
        end

        # D) Check that the model rational function representation matches.
        if m1.q > 0
            maxcoef = maximum(abs.(m1.thetacoef))
            @test all(abs.(m1.thetacoef .- m2.thetacoef) .< EPSILON*maxcoef)
            # At this point, the m3 theta polynomial is not at all guaranteed to match
            # the others, so omit that test for now. If the model_covariance matches,
            # this test is not critical, but we'll think over how it can be improved.
        end

        maxcoef = maximum(abs.(m1.phicoef))
        @test all(abs.(m1.phicoef.-m2.phicoef) .< EPSILON*maxcoef)
        @test all(abs.(m1.phicoef.-m3.phicoef) .< EPSILON*maxcoef)

        # E) Test model_psd. This isn't easy to see how to test, other than re-implement
        # the model_psd code itself!
        N = 50
        freq = collect(range(0, stop=0.5, length=N))
        z = exp.(-2im*pi *freq)
        numer = m1.thetacoef[1] .+ fill(0.0im, N)
        for i=1:m1.q
            numer .+= m1.thetacoef[i+1] * (z.^i)
        end
        denom = m1.phicoef[1] .+ fill(0.0im, N)
        for i=1:m1.p
            denom .+= m1.phicoef[i+1] * (z.^i)
        end
        psd = abs2.(numer ./ denom)
        threshold = 1e-3 * maximum(abs.(psd[1]))
        mpsd = model_psd(m1, N)
        @test size(psd) == size(mpsd)
        @test all(abs.(psd .- mpsd) .< threshold)
        # @test all(abs.(psd - model_psd(m2, N)) .< threshold)
        # @test all(abs.(psd - model_psd(m3, N)) .< threshold)

        mpsd = model_psd(m1, freq)
        @test size(psd) == size(mpsd)
        @test all(abs.(psd .- mpsd) .< threshold)
        # @test all(abs.(psd - model_psd(m2, freq)) .< threshold)
        # @test all(abs.(psd - model_psd(m3, freq)) .< threshold)

    end
end

# 5) Test fitting data to a sum-of-exponentials representation
# and an ARMA model of order (p, q=p)
function test_sum_exp(ampls::Vector, bases::Vector, N::Integer)

    # First, make sure that ARMA.exponential_model does the right thing.
    t = 0:(N-1)
    signal = ARMA.exponential_model(t, ampls, bases)
    signal2 = zero(signal)
    for (b,a) in zip(bases, ampls)
        signal2 .+= real(a*(b.^(0:N-1)))
    end
    @test all(abs.(signal2-signal) .< 1e-6*minimum(abs.(ampls)))

    # Now add a tiny bit of noise, fit exponentials and see what happens.
    # Rather than testing the fit, test the model that it generates.
    noise_level = 1e-4
    signal .+= randn(N)*noise_level
    NB = length(bases)
    afit, bfit = fit_exponentials(signal, pmin=NB, pmax=NB)
    cmodel = ARMA.exponential_model(t, afit, bfit)
    @test all(abs.(cmodel .- signal) .< .1*maximum(signal))

    # Now test the full fitARMA function, with 0 and then 1 exceptional value.
    p = length(bases)
    model = fitARMA(signal, p, p-1, deltar=noise_level, pmin=p-2)
    cmodel = model_covariance(model, N)
    @test all(abs.(cmodel .- signal) .< .1*maximum(signal))

    signal[1] *= 2
    model = fitARMA(signal, p, p, deltar=noise_level, pmin=p-2)
    cmodel = model_covariance(model, N)
    @test all(abs.(cmodel .- signal) .< .1*maximum(signal))
end

@testset "exponential_fits" begin
    ampls=[5.0,4,3-1im,3+1im]
    bases=[.999,.98,.7+.1im,.7-.1im]
    test_sum_exp(ampls, bases, 400)

    ampls=[7.0,5,3-1im,3+1im]
    bases=[.99,.9,.1+.8im,.1-.8im]
    test_sum_exp(ampls, bases, 400)

    ampls=[1,2,3,4,5]
    bases=[.999, .95, .9, .8, .5]
    test_sum_exp(ampls, bases, 400)
end

# 6) Test toeplitz_whiten and toeplitz_whiten! with an ARMA(2,2) and 5 random vectors
@testset "toeplitz_whiten" begin
    r=[3,-3]
    poles = [1.25,-2]
    model = ARMAModel(r, poles, 10.0)
    N = 50
    Phi = fill(0.0, N, N)
    for i=1:model.p+1
        for col=1:N+1-i
            Phi[col+i-1,col] = model.phicoef[i]
        end
    end
    Theta = fill(0.0, N, N)
    for i=1:model.q+1
        for col=1:N+1-i
            Theta[col+i-1,col] = model.thetacoef[i]
        end
    end
    for i=1:5
        v = randn(N)
        correct_tw = Theta \ (Phi * v)
        tw = toeplitz_whiten(model, v)
        @test all(abs.(tw - correct_tw) .< 1e-6)
    end
end

arrays_similar(v::AbstractArray, w::AbstractArray, eps=1e-10) = all(abs.(v-w) .< eps)

# 7) Test internals used by whiten, unwhiten, solve_covariance, mult_covariance
@testset "whiten" begin
    for i=1:5
        N = 50
        v = randn(N)
        vx = copy(v)
        vx[2:end] += 0.8*v[1:end-1]
        vy = copy(v)
        vy[2:end] -= 0.3*v[1:end-1]
        vy[3:end] -= 0.4*v[1:end-2]

        @test arrays_similar( ARMA.convolve_same(v, [1, 0.8]), vx)
        @test arrays_similar( ARMA.deconvolve_same(vx, [1, 0.8]), v)
        @test arrays_similar( ARMA.convolve_same(v, [1, -.3, -.4]), vy)
        @test arrays_similar( ARMA.deconvolve_same(vy, [1, -.3, -.4]), v)
    end

    # Test whiten, unwhiten, solve_covariance, mult_covariance
    model23 = ARMAModel([1.2,1.1,1.02], [1.25, -2], 10)
    model32 = ARMAModel([1.25,-2], [1.2,1.1,1.02], 10)
    model52 = ARMAModel([1.25,-2], [6,2.5,1.2,1.1,1.02], 10)
    model25 = ARMAModel([6,2.5,1.2,1.1,1.02], [1.25, -2], 10)
    for model in (model23, model32, model25, model52)
        N = 16
        gamma = model_covariance(model, N)
        R = SymmetricToeplitz(gamma)
        L = cholesky(R).L
        x = fill(0.0, N)
        y = fill(0.0, N)
        x[1:model.p+1] = model.phicoef
        y[1] = x[1]
        Phi = Toeplitz(x, y)
        # must force symmetry, or numerical precision will make cholesky() fail.
        RR = Hermitian(Phi*R*Phi')
        LL = cholesky(RR).L

        solver = ARMASolver(model, N)
        for N in [2,4,6,8,10,14]
            v = randn(N)
            @test arrays_similar(LL[1:N,1:N]\v, solver.LL[1:N,1:N]\v, 1e-7)
            @test arrays_similar(LL[1:N,1:N]*v, solver.LL[1:N,1:N]*v, 1e-7)
            @test arrays_similar(L[1:N,1:N]\v, whiten(solver, v), 1e-6)
            @test arrays_similar(L[1:N,1:N]*v, unwhiten(solver, v), 1e-7)
            @test arrays_similar(R[1:N,1:N]*v, mult_covariance(solver, v), 1e-7)
            @test arrays_similar(R[1:N,1:N]\v, solve_covariance(solver, v), 1e-3)
            Rinv = inverse_covariance(solver, N)
            @test arrays_similar(R[1:N,1:N]*Rinv, Matrix{Float64}(I, N, N), 1e-7)
            @test arrays_similar(Rinv*R[1:N,1:N], Matrix{Float64}(I, N, N), 1e-7)

            # Test that they can be applied to matrices as well as vectors
            M = randn(N, 4)
            @test arrays_similar(L[1:N,1:N]\M, whiten(solver, M), 1e-5)
            @test arrays_similar(L[1:N,1:N]*M, unwhiten(solver, M), 1e-7)
            @test arrays_similar(R[1:N,1:N]*M, mult_covariance(solver, M), 1e-7)
            @test arrays_similar(R[1:N,1:N]\M, solve_covariance(solver, M), 1e-3)
        end

    end
end

include("hdf5test.jl")
end


text:
"""
    rasterio!(dataset::AbstractDataset, buffer::Array{<:Any, 3},
        bands; <keyword arguments>)
    rasterio!(dataset::AbstractDataset, buffer::Array{<:Any, 3}, bands, rows,
        cols; <keyword arguments>)
    rasterio!(rasterband::AbstractRasterBand, buffer::Matrix{<:Any};
        <keyword arguments>)
    rasterio!(rasterband::AbstractRasterBand, buffer::Matrix{<:Any}, rows,
        cols; <keyword arguments>)

Read/write a region of image data from multiple bands.

This method allows reading a region of one or more `RasterBand`s from this
dataset into a buffer, or writing data from a buffer into a region of the
`RasterBand`s. It automatically takes care of data type translation if the
element type (`<:Any`) of the buffer is different than that of the
`RasterBand`. The method also takes care of image decimation / replication
if the buffer size (`xsz × ysz`) is different than the size of the
region being accessed (`xsize × ysize`).

The `pxspace`, `linespace` and `bandspace` parameters allow reading into
or writing from various organization of buffers.

For highest performance full resolution data access, read and write on \"block
boundaries\" as returned by `blocksize()`, or use the `readblock!()` and
`writeblock!()` methods.

### Parameters
* `rows`        A continuous range of rows expressed as a
                `UnitRange{<:Integer}`, such as 2:9.
* `cols`        A continuous range of columns expressed as a
                `UnitRange{<:Integer}`, such as 2:9.
* `access`      Either `GF_Read` to read a region of data, or
                `GF_Write` to write a region of data.
* `xoffset`     The pixel offset to the top left corner of the region to be
                accessed. It will be `0` (default) to start from the left.
* `yoffset`     The line offset to the top left corner of the region to be
                accessed. It will be `0` (default) to start from the top.
* `xsize`       The width of the region of the band to be accessed in pixels.
* `ysize`       The height of the region of the band to be accessed in lines.
* `buffer`      The buffer into which the data should be read, or from which it
                should be written. It must contain `≥ xsz * ysz * <# of bands>`
                words of type `eltype(buffer)`. It is organized in left to
                right, top to bottom pixel order. Spacing is controlled by the
                `pxspace`, and `linespace` parameters
* `xsz`         The width of the buffer into which the desired region is
                to be read, or from which it is to be written.
* `ysz`         The height of the buffer into which the desired region is
                to be read, or from which it is to be written.
* `bands`       The list of bands (`1`-based) to be read/written.
* `pxspace`     The byte offset from the start of a pixel value in the `buffer`
                to the start of the next pixel value within a scanline. By
                default (i.e., `0`) the size of `eltype(buffer)` will be used.
* `linespace`   The byte offset from the start of one scanline in pBuffer to
                the start of the next. By default (i.e., `0`) the value of
                `sizeof(eltype(buffer)) * xsz` will be used.
* `bandspace`   The byte offset from the start of one bands data to the start
                of the next. By default (`0`), it will be `linespace * ysz`
                implying band sequential organization of the buffer.

### Returns
`CE_Failure` if the access fails, otherwise `CE_None`.
"""
function rasterio! end

function rasterio!(
    dataset::AbstractDataset,
    buffer::T,
    bands,
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
    bandspace::Integer = 0,
)::T where {T<:Array{<:Any,3}}
    rasterio!(
        dataset,
        buffer,
        bands,
        0,
        0,
        size(buffer, 1),
        size(buffer, 2),
        access,
        pxspace,
        linespace,
        bandspace,
    )
    return buffer
end

function rasterio!(
    dataset::AbstractDataset,
    buffer::Array{T,3},
    bands,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
    bandspace::Integer = 0,
)::Array{T,3} where {T<:Any}
    xsize = cols[end] - cols[1] + 1
    xsize < 0 && error("invalid window width")
    ysize = rows[end] - rows[1] + 1
    ysize < 0 && error("invalid window height")
    rasterio!(
        dataset,
        buffer,
        bands,
        cols[1] - 1,
        rows[1] - 1,
        xsize,
        ysize,
        access,
        pxspace,
        linespace,
        bandspace,
    )
    return buffer
end

function rasterio!(
    dataset::AbstractDataset,
    buffer::Array{T,3},
    bands,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
    bandspace::Integer = 0,
    extraargs = Ptr{GDAL.GDALRasterIOExtraArg}(C_NULL),
)::Array{T,3} where {T<:Any}
    # `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg
    # structure with additional arguments to specify resampling and
    # progress callback, or `NULL` for default behaviour. The
    # `GDAL_RASTERIO_RESAMPLING` configuration option can also be
    # defined to override the default resampling to one of `BILINEAR`,
    # `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.
    (dataset == C_NULL) && error("Can't read NULL dataset")
    xbsize, ybsize, zbsize = size(buffer)
    nband = length(bands)
    bands = isa(bands, Vector{Cint}) ? bands : Cint.(collect(bands))
    @assert nband == zbsize
    result = GDAL.gdaldatasetrasterioex(
        dataset.ptr,
        access,
        xoffset,
        yoffset,
        xsize,
        ysize,
        pointer(buffer),
        xbsize,
        ybsize,
        convert(GDALDataType, T),
        nband,
        pointer(bands),
        pxspace,
        linespace,
        bandspace,
        extraargs,
    )
    @cplerr result "Access in DatasetRasterIO failed."
    return buffer
end

function rasterio!(
    rasterband::AbstractRasterBand,
    buffer::T,
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
)::T where {T<:Matrix{<:Any}}
    rasterio!(
        rasterband,
        buffer,
        0,
        0,
        width(rasterband),
        height(rasterband),
        access,
        pxspace,
        linespace,
    )
    return buffer
end

function rasterio!(
    rasterband::AbstractRasterBand,
    buffer::T,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
)::T where {T<:Matrix{<:Any}}
    xsize = length(cols)
    xsize < 1 && error("invalid window width")
    ysize = length(rows)
    ysize < 1 && error("invalid window height")
    rasterio!(
        rasterband,
        buffer,
        cols[1] - 1,
        rows[1] - 1,
        xsize,
        ysize,
        access,
        pxspace,
        linespace,
    )
    return buffer
end

function rasterio!(
    rasterband::AbstractRasterBand,
    buffer::Matrix{T},
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
    access::GDALRWFlag = GF_Read,
    pxspace::Integer = 0,
    linespace::Integer = 0,
    extraargs = Ptr{GDAL.GDALRasterIOExtraArg}(C_NULL),
)::Matrix{T} where {T<:Any}
    # `psExtraArg`  (new in GDAL 2.0) pointer to a GDALRasterIOExtraArg
    # structure with additional arguments to specify resampling and
    # progress callback, or `NULL` for default behaviour. The
    # `GDAL_RASTERIO_RESAMPLING` configuration option can also be
    # defined to override the default resampling to one of `BILINEAR`,
    # `CUBIC`, `CUBICSPLINE`, `LANCZOS`, `AVERAGE` or `MODE`.
    (rasterband == C_NULL) && error("Can't read NULL rasterband")
    xbsize, ybsize = size(buffer)
    result = GDAL.gdalrasterioex(
        rasterband.ptr,
        access,
        xoffset,
        yoffset,
        xsize,
        ysize,
        pointer(buffer),
        xbsize,
        ybsize,
        convert(GDALDataType, T),
        pxspace,
        linespace,
        extraargs,
    )
    @cplerr result "Access in RasterIO failed."
    return buffer
end

function read!(rb::AbstractRasterBand, buffer::T)::T where {T<:Matrix{<:Any}}
    rasterio!(rb, buffer, GF_Read)
    return buffer
end

function read!(
    rb::AbstractRasterBand,
    buffer::T,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:Matrix{<:Any}}
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize)
    return buffer
end

function read!(
    rb::AbstractRasterBand,
    buffer::T,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::T where {T<:Matrix{<:Any}}
    rasterio!(rb, buffer, rows, cols)
    return buffer
end

function read(rb::AbstractRasterBand)
    return rasterio!(rb, Array{pixeltype(rb)}(undef, width(rb), height(rb)))
end

function read(
    rb::AbstractRasterBand,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::Matrix{pixeltype(rb)}
    buffer = Matrix{pixeltype(rb)}(undef, xsize, ysize)
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize)
    return buffer
end

function read(
    rb::AbstractRasterBand,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::Matrix{pixeltype(rb)}
    buffer = Matrix{pixeltype(rb)}(undef, length(cols), length(rows))
    rasterio!(rb, buffer, rows, cols)
    return buffer
end

function write!(
    rb::AbstractRasterBand,
    buffer::T,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:Matrix{<:Any}}
    rasterio!(rb, buffer, xoffset, yoffset, xsize, ysize, GF_Write)
    return buffer
end

function write!(
    rb::AbstractRasterBand,
    buffer::T,
    rows::UnitRange{<:Integer} = 1:height(rb),
    cols::UnitRange{<:Integer} = 1:width(rb),
)::T where {T<:Matrix{<:Any}}
    rasterio!(rb, buffer, rows, cols, GF_Write)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    i::Integer,
)::T where {T<:Matrix{<:Any}}
    read!(getband(dataset, i), buffer)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    indices,
)::T where {T<:Array{<:Any,3}}
    rasterio!(dataset, buffer, indices, GF_Read)
    return buffer
end

function read!(dataset::AbstractDataset, buffer::T)::T where {T<:Array{<:Any,3}}
    nband = nraster(dataset)
    @assert size(buffer, 3) == nband
    rasterio!(dataset, buffer, collect(Cint, 1:nband), GF_Read)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    i::Integer,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:Matrix{<:Any}}
    read!(getband(dataset, i), buffer, xoffset, yoffset, xsize, ysize)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    indices,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:Array{<:Any,3}}
    rasterio!(dataset, buffer, indices, xoffset, yoffset, xsize, ysize)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    i::Integer,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::T where {T<:Matrix{<:Any}}
    read!(getband(dataset, i), buffer, rows, cols)
    return buffer
end

function read!(
    dataset::AbstractDataset,
    buffer::T,
    indices,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::T where {T<:Array{<:Any,3}}
    rasterio!(dataset, buffer, indices, rows, cols)
    return buffer
end

read(dataset::AbstractDataset, i::Integer) = read(getband(dataset, i))

function read(dataset::AbstractDataset, indices)::Array{pixeltype(dataset),3}
    buffer = Array{pixeltype(dataset)}(
        undef,
        width(dataset),
        height(dataset),
        length(indices),
    )
    rasterio!(dataset, buffer, indices)
    return buffer
end

function read(dataset::AbstractDataset)::Array{pixeltype(dataset),3}
    buffer = Array{pixeltype(dataset)}(
        undef,
        width(dataset),
        height(dataset),
        nraster(dataset),
    )
    read!(dataset, buffer)
    return buffer
end

function read(
    dataset::AbstractDataset,
    i::Integer,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::Matrix{pixeltype(dataset)}
    band = getband(dataset, i)
    buffer = Matrix{pixeltype(band)}(undef, xsize, ysize)
    read!(dataset, buffer, i, xoffset, yoffset, xsize, ysize)
    return buffer
end

function read(
    dataset::AbstractDataset,
    indices,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::Array{pixeltype(dataset),3}
    buffer = Array{pixeltype(dataset)}(undef, xsize, ysize, length(indices))
    rasterio!(dataset, buffer, indices, xsize, ysize, xoffset, yoffset)
    return buffer
end

function read(
    dataset::AbstractDataset,
    i::Integer,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::Matrix{pixeltype(dataset)}
    buffer = read(getband(dataset, i), rows, cols)
    return buffer
end

function read(
    dataset::AbstractDataset,
    indices,
    rows::UnitRange{<:Integer},
    cols::UnitRange{<:Integer},
)::Array{pixeltype(dataset),3}
    buffer = Array{pixeltype(dataset),3}(
        undef,
        length(cols),
        length(rows),
        length(indices),
    )
    rasterio!(dataset, buffer, indices, rows, cols)
    return buffer
end

function write!(
    dataset::T,
    buffer::Matrix{<:Any},
    i::Integer,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:AbstractDataset}
    write!(getband(dataset, i), buffer, xoffset, yoffset, xsize, ysize)
    return dataset
end

function write!(
    dataset::T,
    buffer::Array{<:Any,3},
    indices,
    xoffset::Integer,
    yoffset::Integer,
    xsize::Integer,
    ysize::Integer,
)::T where {T<:AbstractDataset}
    rasterio!(
        dataset,
        buffer,
        indices,
        xoffset,
        yoffset,
        xsize,
        ysize,
        GF_Write,
    )
    return dataset
end

function write!(
    dataset::T,
    buffer::Matrix{<:Any},
    i::Integer,
    rows::UnitRange{<:Integer} = 1:height(getband(dataset, i)),
    cols::UnitRange{<:Integer} = 1:width(getband(dataset, i)),
)::T where {T<:AbstractDataset}
    write!(getband(dataset, i), buffer, rows, cols)
    return dataset
end

function write!(
    dataset::T,
    buffer::Array{<:Any,3},
    indices,
    rows::UnitRange{<:Integer} = 1:height(dataset),
    cols::UnitRange{<:Integer} = 1:width(dataset),
)::T where {T<:AbstractDataset}
    rasterio!(dataset, buffer, indices, rows, cols, GF_Write)
    return dataset
end

"""
    readblock!(rb::AbstractRasterBand, xoffset::Integer, yoffset::Integer,
        buffer)

Read a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use RasterIO().

### Parameters
* `xoffset` the horizontal block offset, with zero indicating the left most
            block, 1 the next block and so forth.
* `yoffset` the vertical block offset, with zero indicating the top most block,
            1 the next block and so forth.
* `buffer`  the buffer into which the data will be read. The buffer must be
            large enough to hold GetBlockXSize()*GetBlockYSize() words of type
            GetRasterDataType().
"""
function readblock!(
    rb::AbstractRasterBand,
    xoffset::Integer,
    yoffset::Integer,
    buffer::T,
)::T where {T<:Any}
    result = GDAL.gdalreadblock(rb.ptr, xoffset, yoffset, buffer)
    @cplerr result "Failed to read block at ($xoffset,$yoffset)"
    return buffer
end

"""
    writeblock!(rb::AbstractRasterBand, xoffset::Integer, yoffset::Integer,
        buffer)

Write a block of image data efficiently.

This method accesses a "natural" block from the raster band without resampling,
or data type conversion. For a more generalized, but potentially less efficient
access use RasterIO().

### Parameters
* `xoffset` the horizontal block offset, with zero indicating the left most
            block, 1 the next block and so forth.
* `yoffset` the vertical block offset, with zero indicating the left most block,
            1 the next block and so forth.
* `buffer`  the buffer from which the data will be written. The buffer must be
            large enough to hold GetBlockXSize()*GetBlockYSize() words of type
            GetRasterDataType().
"""
function writeblock!(
    rb::T,
    xoffset::Integer,
    yoffset::Integer,
    buffer,
)::T where {T<:AbstractRasterBand}
    result = GDAL.gdalwriteblock(rb.ptr, xoffset, yoffset, buffer)
    @cplerr result "Failed to write block at ($xoffset, $yoffset)"
    return rb
end


text:
using Documenter
using Literate

julia_files(dir) = sort(filter(endswith(".jl"), readdir(dir, join=true)))

function generate_markdowns(dir)
    for file in julia_files(dir)
        Literate.markdown(
            file,
            dir;
            documenter = true,
        )
    end
end

generate_markdowns(joinpath(@__DIR__, "src", "about"))
generate_markdowns(joinpath(@__DIR__, "src", "tutorial"))
generate_markdowns(joinpath(@__DIR__, "src", "howto"))

if get(ENV, "JULIA_DOCUMENTER_CI", "OFF") == "ON"
    using Pkg
    Pkg.develop(path=joinpath(dirname(@__DIR__), "generator"))
else
    push!(LOAD_PATH, joinpath(dirname(@__DIR__), "generator"))
end

using Vulkan
using VulkanGen

makedocs(;
    modules=[Vulkan, VulkanGen],
    format=Documenter.HTML(prettyurls = true),
    pages=[
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "About" => [
            "Motivations" => "about/motivations.md",
            "Interfacing with the C API" => "about/interfacing.md",
        ],
        "Tutorial" => [
            "Getting started" => "tutorial/getting_started.md",
            "Vulkan types" => "tutorial/types.md",
            "Vulkan functions" => "tutorial/functions.md",
            "Error handling" => "tutorial/error_handling.md",
            "Resource management" => "tutorial/resource_management.md",
            "Dispatch" => "tutorial/dispatch.md",
            "In-depth tutorial" => "tutorial/indepth.md",
        ],
        "How to" => [
            "howto/preferences.md",
            "howto/debugging.md",
        ],
        "API" => "api.md",
        "Utility" => "utility.md",
        "Troubleshooting" => "troubleshooting.md",
        "Developer documentation" => [
            "Overview" => "dev/overview.md",
            "Vulkan specification" => "dev/spec.md",
            "Generator" => "dev/gen.md",
        ],
    ],
    repo="https://github.com/JuliaGPU/Vulkan.jl/blob/{commit}{path}#L{line}",
    sitename="Vulkan.jl",
    authors="serenity4 <cedric.bel@hotmail.fr>",
    strict=false,
    doctest=false,
    checkdocs=:exports,
    linkcheck=:true,
)

deploydocs(
    repo = "github.com/JuliaGPU/Vulkan.jl.git",
)


text:
"""
`init_subspec!(m::PoolModel)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::PoolModel)
    if subspec(m) == "ss0"
        return
    elseif subspec(m) == "ss1"
        return ss1!(m)
    elseif subspec(m) == "ss2"
        return ss2!(m)
    elseif subspec(m) == "ss3"
        return ss3!(m)
    elseif subspec(m) == "ss4"
        return ss4!(m)
    elseif subspec(m) == "ss5"
        return ss5!(m)
    elseif subspec(m) == "ss6"
        return ss6!(m)
    elseif subspec(m) == "ss7"
        return ss7!(m)
    elseif subspec(m) == "ss8"
        return ss8!(m)
    elseif subspec(m) == "ss9"
        return ss9!(m)
    elseif subspec(m) == "ss10"
        return ss10!(m)
    else
        error("This subspec is not defined")
    end
end

"""
```
ss1!(m::PoolModel)
```

Initialize subspec 1 of `PoolModel` which uses the original Matlab predictive densities.
"""
function ss1!(m::PoolModel)
end

"""
```
ss2!(m::PoolModel)
```

Initialize subspec 2 of `PoolModel` which uses the corrected `Model805` predictive density
and the original `Model904` predictive density from Matlab.
"""
function ss2!(m::PoolModel)
end

"""
```
ss3!(m::PoolModel)
```

Initialize subspec 3 of `PoolModel` which uses the original Matlab predictive densities.
This uses Prior 2 in the archived Matlab code.

```
ρ ∼ Β(0.8,0.1),
μ ∼ N(0, Φ^{-1}(0.75))
σ^2 ∼ IG(1,4)
```
"""
function ss3!(m::PoolModel)
    m <= parameter(:ρ, 0.5, (1e-5, 0.999), (1e-5,0.999), SquareRoot(), Beta(0.8,0.1),
                   fixed = false,
                   description="ρ: persistence of AR processing underlying λ.",
                   tex_label="\\rho")
    m <= parameter(:μ, 0., (-Inf, Inf), (-Inf, Inf), Untransformed(), Normal(0., quantile(Normal(), 0.75)),
                   fixed = false,
                   description="μ: drift of AR processing underlying λ.",
                   tex_label="\\mu")
    m <= parameter(:σ, 1., (1e-5, Inf), (1e-5, Inf), ModelConstructors.Exponential(), RootInverseGamma(2., 2.),
                   fixed = false,
                   description="σ: volatility of AR processing underlying λ.",
                   tex_label="\\sigma")

end

"""
```
ss4!(m::PoolModel)
```

Initialize subspec 4 of `PoolModel` which uses the corrected `Model805` predictive density.
This uses Prior 2 in the archived Matlab code.

```
ρ ∼ Β(0.8,0.1),
μ ∼ N(0, Φ^{-1}(0.75))
σ^2 ∼ IG(1,4)
```
"""
function ss4!(m::PoolModel)
    ss3!(m)
end

"""
```
ss5!(m::PoolModel)
```

Initialize subspec 5 of `PoolModel` which uses the original Matlab predictive densities.
This uses Prior 3 in the archived Matlab code.

```
ρ ∼ Β(0.8,0.1),
σ^2 ∼ IG(1,4)
```
"""
function ss5!(m::PoolModel)
    m <= parameter(:ρ, 0.5, (1e-5, 0.999), (1e-5,0.999), SquareRoot(), Beta(0.8,0.1),
                   fixed = false,
                   description="ρ: persistence of AR processing underlying λ.",
                   tex_label="\\rho")
    m <= parameter(:σ, 1., (1e-5, Inf), (1e-5, Inf), ModelConstructors.Exponential(), RootInverseGamma(2., 2.),
                   fixed = false,
                   description="σ: volatility of AR processing underlying λ.",
                   tex_label="\\sigma")
end

"""
```
ss6!(m::PoolModel)
```

Initialize subspec 6 of `PoolModel` which uses the corrected `Model805` predictive density.
This uses Prior 3 in the archived Matlab code.

```
ρ ∼ Β(0.8,0.1),
σ^2 ∼ IG(1,4)
```
"""
function ss6!(m::PoolModel)
    ss5!(m)
end

"""
```
ss7!(m::PoolModel)
```

Initialize subspec 7 of `PoolModel` which uses the original Matlab predictive densities.
This uses Prior 2 in the paper:

```
ρ ∼ Β(0.8,0.1),
μ ∼ N(0, Φ^{-1}(0.75))
σ^2 ∼ IG(2,1)
```
"""
function ss7!(m::PoolModel)
    m <= parameter(:ρ, 0.5, (1e-5, 0.999), (1e-5,0.999), SquareRoot(), Beta(0.8,0.1),
                   fixed = false,
                   description="ρ: persistence of AR processing underlying λ.",
                   tex_label="\\rho")
    m <= parameter(:μ, 0., (-Inf, Inf), (-Inf, Inf), Untransformed(), Normal(0., quantile(Normal(), 0.75)),
                   fixed = false,
                   description="μ: drift of AR processing underlying λ.",
                   tex_label="\\mu")
    m <= parameter(:σ, 1., (1e-5, Inf), (1e-5, Inf), ModelConstructors.Exponential(), RootInverseGamma(4., 1 / sqrt(2.)),
                   fixed = false,
                   description="σ: volatility of AR processing underlying λ.",
                   tex_label="\\sigma")

end

"""
```
ss8!(m::PoolModel)
```

Initialize subspec 8 of `PoolModel` which uses the corrected `Model805` predictive density.
This uses Prior 2 in the paper:

```
ρ ∼ Β(0.8,0.1),
μ ∼ N(0, Φ^{-1}(0.75))
σ^2 ∼ IG(2,1)
```
"""
function ss8!(m::PoolModel)
    ss7!(m)
end

"""
```
ss9!(m::PoolModel)
```

Initialize subspec 9 of `PoolModel` which uses the original Matlab predictive densities.
This uses Prior 3 in the paper.

```
ρ ∼ Β(0.8,0.1),
σ^2 ∼ IG(2,1)
```
"""
function ss9!(m::PoolModel)
    m <= parameter(:ρ, 0.5, (1e-5, 0.999), (1e-5,0.999), SquareRoot(), Beta(0.8,0.1),
                   fixed = false,
                   description="ρ: persistence of AR processing underlying λ.",
                   tex_label="\\rho")
    m <= parameter(:σ, 1., (1e-5, Inf), (1e-5, Inf), ModelConstructors.Exponential(), RootInverseGamma(4., 1 / sqrt(2.)),
                   fixed = false,
                   description="σ: volatility of AR processing underlying λ.",
                   tex_label="\\sigma")
end

"""
```
ss10!(m::PoolModel)
```

Initialize subspec 10 of `PoolModel` which uses the corrected `Model805` predictive density.
This uses Prior 3 in the paper.

```
ρ ∼ Β(0.8,0.1),
σ^2 ∼ IG(2,1)
```
"""
function ss10!(m::PoolModel)
    ss9!(m)
end


text:
devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticTransformer),
    :TT => DeviceModel(TapTransformer, StaticTransformer),
)
services = Dict{Symbol, ServiceModel}()

@testset "Operation Model kwargs with CopperPlatePowerModel base" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    c_sys5_re = build_system("c_sys5_re")
    c_sys14 = build_system("c_sys14")

    @test_throws ArgumentError OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        bad_kwarg = 10,
    )
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    moi_tests(op_problem, true, 120, 0, 120, 120, 24, false)
    op_problem =
        OperationsProblem(TestOpProblem, template, c_sys14; optimizer = OSQP_optimizer)
    moi_tests(op_problem, false, 120, 0, 120, 120, 24, false)
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        use_forecast_data = false,
        optimizer = GLPK_optimizer,
    )
    moi_tests(op_problem, false, 5, 0, 5, 5, 1, false)
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        use_forecast_data = false,
        use_parameters = false,
        optimizer = GLPK_optimizer,
    )
    moi_tests(op_problem, false, 5, 0, 5, 5, 1, false)

    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        optimizer = GLPK_optimizer,
        slack_variables = true,
    )
    moi_tests(op_problem, false, 168, 0, 120, 120, 24, false)

end

@testset "Test optimization debugging functions" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    MOIU.attach_optimizer(op_problem.psi_container.JuMPmodel)
    con_index = collect(get_all_constraint_index(op_problem))
    length = size(con_index, 1)
    constraint =
        op_problem.psi_container.constraints[last(con_index)[1]].data[last(con_index)[2]]
    @test get_con_index(op_problem, length) == constraint
    @test get_con_index(op_problem, length + 1) == nothing
    var_index = collect(get_all_var_index(op_problem))
    length = size(var_index, 1)
    var = op_problem.psi_container.variables[last(var_index)[1]].data[last(var_index)[2]]
    @test get_var_index(op_problem, length) == var
    @test get_var_index(op_problem, length + 1) == nothing
end

@testset "Test passing custom JuMP model" begin
    my_model = JuMP.Model()
    my_model.ext[:PSI_Testing] = 1
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5 = build_system("c_sys5")
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5,
        my_model;
        optimizer = GLPK_optimizer,
        use_parameters = true,
    )
    @test haskey(op_problem.psi_container.JuMPmodel.ext, :PSI_Testing)
    @test (:params in keys(op_problem.psi_container.JuMPmodel.ext)) == true
end

@testset "Operation Model Constructors with Parameters" begin

    networks = [
        CopperPlatePowerModel,
        StandardPTDFModel,
        DCPPowerModel,
        NFAPowerModel,
        ACPPowerModel,
        ACRPowerModel,
        ACTPowerModel,
        DCPLLPowerModel,
        LPACCPowerModel,
        SOCWRPowerModel,
        QCRMPowerModel,
        QCLSPowerModel,
    ]

    thermal_gens = [
        ThermalStandardUnitCommitment,
        ThermalDispatch,
        ThermalRampLimited,
        ThermalDispatchNoMin,
    ]

    c_sys5 = build_system("c_sys5")
    c_sys5_re = build_system("c_sys5_re")
    c_sys5_bat = build_system("c_sys5_bat")
    systems = [c_sys5, c_sys5_re, c_sys5_bat]
    for net in networks, thermal in thermal_gens, system in systems, p in [true, false]
        @testset "Operation Model $(net) - $(thermal) - $(system)" begin
            devices = Dict{Symbol, DeviceModel}(
                :Generators => DeviceModel(ThermalStandard, thermal),
                :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
            )
            branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(Line, StaticLine))
            template = OperationsProblemTemplate(net, devices, branches, services)
            op_problem = OperationsProblem(
                TestOpProblem,
                template,
                system;
                use_parameters = p,
                PTDF = build_PTDF5(),
            )
            @test :nodal_balance_active in keys(op_problem.psi_container.expressions)
            @test (:params in keys(op_problem.psi_container.JuMPmodel.ext)) == p
        end
    end

    @testset "Operations template constructors" begin
        c_sys5 = build_system("c_sys5")
        op_problem_ed = PSI.EconomicDispatchProblem(c_sys5)
        op_problem_uc = PSI.UnitCommitmentProblem(c_sys5)
        moi_tests(op_problem_uc, false, 480, 0, 240, 120, 144, true)
        moi_tests(op_problem_ed, false, 120, 0, 168, 120, 24, false)
        ED = PSI.run_economic_dispatch(c_sys5; optimizer = fast_lp_optimizer)
        UC = PSI.run_unit_commitment(c_sys5; optimizer = fast_lp_optimizer)
        @test ED.optimizer_log[:primal_status] == MOI.FEASIBLE_POINT
        @test UC.optimizer_log[:primal_status] == MOI.FEASIBLE_POINT
    end
end


text:
using PyPlot
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../test/useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function displacement(alpha, x)
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function stress_field(lambda, mu, alpha, x)
    s11 =
        (lambda + 2mu) * alpha * pi * x[2] * cos(pi * x[1]) -
        lambda * alpha * pi * sin(pi * x[2])
    s22 =
        -(lambda + 2mu) * alpha * pi * sin(pi * x[2]) +
        lambda * alpha * pi * x[2] * cos(pi * x[1])
    s12 = alpha * mu * (3x[1]^2 + sin(pi * x[1]))
    return [s11, s22, s12]
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    penalty,
    alpha,
)

    L, W = CutCell.widths(cutmesh)
    lambda, mu = CutCell.lame_coefficients(stiffness, +1)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )
    displacementbc = CutCell.DisplacementCondition(
        x -> displacement(alpha, x),
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, L, W),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(
        sysmatrix,
        interfacecondition,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        displacementbc,
        cutmesh,
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement
end

function spatial_closest_points(refclosestpoints, refclosestcellids, mesh)
    dim, npts = size(refclosestpoints)
    spclosestpoints = zeros(dim, npts)
    for i = 1:npts
        cellmap = CutCell.cell_map(mesh, refclosestcellids[i])
        spclosestpoints[:, i] .= cellmap(refclosestpoints[:, i])
    end
    return spclosestpoints
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

function traction_force(stressvector, normal)
    return [
        stressvector[1] * normal[1] + stressvector[3] * normal[2],
        stressvector[3] * normal[1] + stressvector[2] * normal[2],
    ]
end

function traction_force_at_points(stresses, normals)
    npts = size(stresses)[2]
    @assert size(normals)[2] == npts

    return hcat(
        [traction_force(stresses[:, i], normals[:, i]) for i = 1:npts]...,
    )
end

K1, K2 = 247.0, 192.0    # Pa
mu1, mu2 = 126.0, 87.0   # Pa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 1.0
rho2 = 1.0
V1 = inv(rho1)
V2 = inv(rho2)

theta0 = 0.0
transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

width = 1.0
displacementscale = 0.01 * width
penaltyfactor = 1e1

polyorder = 3
numqp = required_quadrature_order(polyorder) + 2
nelmts = 9
center = [width / 2, width / 2]
inradius = width / 4

dx = width / nelmts
meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
penalty = penaltyfactor / dx * meanmoduli

basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
levelset = InterpolatingPolynomial(1, basis)
levelsetcoeffs = CutCell.levelset_coefficients(
    x -> -circle_distance_function(x, center, inradius),
    mesh,
)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
interfacequads =
    CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

nodaldisplacement = solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    penalty,
    displacementscale,
)

refseedpoints, spatialseedpoints, seedcellids =
    CutCell.seed_zero_levelset_with_interfacequads(interfacequads, cutmesh)
normals = CutCell.collect_interface_normals(interfacequads, cutmesh)

spatialpoints = spatialseedpoints
referencepoints = refseedpoints
referencecellids = seedcellids

relspatialpoints = spatialpoints .- center
angularposition = angular_position(relspatialpoints)
sortidx = sortperm(angularposition)

angularposition = angularposition[sortidx]
referencepoints = referencepoints[:, sortidx]
spatialpoints = spatialpoints[:, sortidx]
referencecellids = referencecellids[sortidx]
normals = normals[:, sortidx]

productdisplacement = CutCell.displacement_at_reference_points(
    referencepoints,
    referencecellids,
    +1,
    basis,
    nodaldisplacement,
    cutmesh,
)
parentdisplacement = CutCell.displacement_at_reference_points(
    referencepoints,
    referencecellids,
    -1,
    basis,
    nodaldisplacement,
    cutmesh,
)
exactdisplacement = hcat(
    [
        displacement(displacementscale, spatialpoints[:, i]) for
        i = 1:size(spatialpoints)[2]
    ]...,
)

# fig, ax = PyPlot.subplots(2, 1)
# ax[1].plot(angularposition, productdisplacement[1, :], label = "product u1")
# ax[1].plot(angularposition, parentdisplacement[1, :], label = "parent u1")
# ax[1].grid()
# ax[1].legend()
# ax[2].plot(angularposition, productdisplacement[2, :], label = "product u2")
# ax[2].plot(angularposition, parentdisplacement[2, :], label = "parent u2")
# ax[2].grid()
# ax[2].legend()
# fig.tight_layout()
# fig
#

productstress = CutCell.product_stress_at_reference_points(
    referencepoints,
    referencecellids,
    basis,
    stiffness,
    transfstress,
    theta0,
    nodaldisplacement,
    cutmesh,
)
parentstress = CutCell.parent_stress_at_reference_points(
    referencepoints,
    referencecellids,
    basis,
    stiffness,
    nodaldisplacement,
    cutmesh,
)

producttraction = traction_force_at_points(productstress,normals)
parenttraction = traction_force_at_points(parentstress,normals)

# fig, ax = PyPlot.subplots(2, 1)
# ax[1].plot(angularposition, producttraction[1, :], label = "product")
# ax[1].plot(angularposition, parenttraction[1, :], label = "parent")
# ax[1].set_title("t1")
# ax[1].legend()
# ax[1].grid()
# ax[2].plot(angularposition, producttraction[2, :], label = "product")
# ax[2].plot(angularposition, parenttraction[2, :], label = "parent")
# ax[2].set_title("t2")
# ax[2].legend()
# ax[2].grid()
# fig.tight_layout()
# fig

productpressure = CutCell.pressure_at_points(productstress)
parentpressure = CutCell.pressure_at_points(parentstress)

fig,ax = PyPlot.subplots()
ax.plot(angularposition,productpressure,label="product")
ax.plot(angularposition,parentpressure,label="parent")
ax.legend()
ax.grid()
ax.set_title("Pressure on interface")
fig

# productdevstress = CutCell.deviatoric_stress_at_points(productstress,productpressure)
# parentdevstress = CutCell.deviatoric_stress_at_points(parentstress,parentpressure)
#
# productdevnorm = CutCell.stress_inner_product_over_points(productdevstress)
# parentdevnorm = CutCell.stress_inner_product_over_points(parentdevstress)
#
# productnormaldevstress = CutCell.normal_stress_component_over_points(productdevstress,normals)
# parentnormaldevstress = CutCell.normal_stress_component_over_points(parentdevstress,normals)
#
# productspecificvolume = V1*(1.0 .- productpressure/K1)
# parentspecificvolume = V2*(1.0 .- parentpressure/K2)
#
#
# productp1 = productspecificvolume .* productpressure
# productp2 = V1/(2K1)*(productpressure.^2)
# productp3 = -(productspecificvolume .* productnormaldevstress)
# productp4 = V1/(4mu1)*productdevnorm
# productpotential = productp1+productp2+productp3+productp4
#
# parentp1 = parentspecificvolume .* parentpressure
# parentp2 = V2/(2K2)*(parentpressure.^2)
# parentp3 = -(parentspecificvolume .* parentnormaldevstress)
# parentp4 = V2/(4mu2)*parentdevnorm
# parentpotential = parentp1+parentp2+parentp3+parentp4
#
# potentialdifference = productpotential - parentpotential


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    RandomPath(domain)

A random path on a spatial `domain`.
"""
struct RandomPath{D<:AbstractDomain} <: AbstractPath{D}
  domain::D
  permut::Vector{Int}

  function RandomPath{D}(domain, permut) where {D<:AbstractDomain}
    @assert length(permut) == npoints(domain) "incorrect dimension"
    new(domain, permut)
  end
end

RandomPath(domain) = RandomPath{typeof(domain)}(domain, randperm(npoints(domain)))

Base.iterate(p::RandomPath, state=1) = state > npoints(p.domain) ? nothing : (p.permut[state], state + 1)


text:
# File holding the cross validation object

using Base.Threads
using ProgressMeter
using SparseGrids
include("RFR.jl")

check_random_state(seed::Int) = MersenneTwister(seed)
check_random_state(rng::AbstractRNG) = rng

# Cross Validation object
mutable struct cross_val

    # internal parameters
    mse_list::Vector{Float64}
    bias_list::Vector{Float64}
    variance_list::Vector{Float64}
    regressor_list::Vector{AbstractRegressor}
    dictionary_list::Vector{Dict}
    n_models::Union{Nothing, Int}

    # external parameters
    parameter_list::Dict
    random_state::Union{AbstractRNG, Int}

    # Initializing new object
    cross_val(parameter_list;
            random_state = 0) = new(
        [], [], [], [],
        make_single_dicts(parameter_list), nothing,
        parameter_list,
        check_random_state(random_state)
    )
end

# Function for making a single dictionary out of multiple parameters
# Needed when supplying multiple vectors for tuning parameters
# For example optimizing over α and max_depth
function make_single_dicts(d::Dict)

    names = collect(keys(d))
    # Get all combinations of values from the dictionary
    combs = combvec([float(v) for (k,v) in d])

    result_array = Array{Dict}(undef, length(combs))
    @inbounds for i in 1:length(combs)
        result_array[i] = Dict(zip(names, combs[i]))
    end

    return result_array
end

function initiate_RF!(cv::cross_val)
    cv.n_models = length(cv.dictionary_list)
    # println("Making ", cv.n_models, " different specifications")

    # Allocate space holding the results for each random forest
    cv.regressor_list = Vector{AbstractRegressor}(undef, cv.n_models)

    @inbounds for i in 1:cv.n_models
        rf = RFR(param_dict = cv.dictionary_list[i])
        cv.regressor_list[i] = rf
    end

    # Make room for validation results
    cv.mse_list = Float64[]
    cv.bias_list = Float64[]
    cv.variance_list = Float64[]

    return cv
end

# Function for getting the MSE
function get_mse(pred, y)
    bias = mean(pred .- y)
    variance = var(pred)
    mse = mean((pred .- mean(y)).^2)

    # Simlpe test to see if everything worked
    @assert mse == variance + bias^2 "mse != variance + bias^2"

    return bias, variance, mse
end

# Fitting the cross-validation object to data
function fit!(cv::cross_val, X::Matrix, Y::Matrix; nfolds::Int=3)

    # Initializing cross-validation object with random forests
    initiate_RF!(cv)

    incr = Int(floor(size(X, 1)/nfolds))
    println("Fitting ", nfolds, " folds")

    # Iterating over the list of dictionaries to fit the random forests
    for forest in cv.regressor_list

        result_matrix = Array{Float64}(undef, nfolds, 3)

        # Iterating through the folds
        ind = 1:1
        @inbounds for i = 1:nfolds
            # Increase interval for the next fold
            ind = (ind[end]:i*incr)

            # Splitting into hold-out sample
            train_cond = (1:size(X,1) .< ind[1]) + (1:size(X,1) .> ind[end]) .== 1
            test_cond = ind[1] .< 1:size(X, 1) .< ind[end]

            x_train = X[train_cond, :]
            x_test = X[test_cond, :]
            y_train = Y[train_cond, :]
            y_test = Y[test_cond, :]

            # Fit and predict given hold-out sample
            fit!(forest, x_train, y_train)
            pred = predict(forest, x_test)

            b, v, m = get_mse(pred, y_test)
            result_matrix[i, :] = [b, v, m]
        end

        # Assign values
        bias, variance, mse = mean(result_matrix, dims=2)
        validation_results!(cv, bias, variance, mse)
    end
end

# Add results to list which holds all results
function validation_results!(cv::cross_val, bias::Float64, variance::Float64, mse::Float64)
    push!(cv.bias_list, bias)
    push!(cv.variance_list, variance)
    push!(cv.mse_list, mse)
end

# Get best model out of all tried for given measure
function best_model(cv::cross_val, measure::String)

    ind = -1
    if measure == "mse"
        ind = cv.mse_list .== minimum(cv.mse_list)
    elseif measure == "bias"
        ind = cv.bias_list .== minimum(cv.bias_list)
    elseif measure == "variance"
        ind = cv.variance_list .== minimum(cv.variance_list)
    end

    return cv.regressor_list[ind][1]
end

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
@testset "pretty printing" begin
    P = sprint(show, MIME"text/plain"(), A)
    @test contains(P, "ClimArray")
    @test contains(P, "Lon")
    @test contains(P, "data")
    @test contains(P, "44")
end

@testset "General Statistics" begin
@testset "Dropping dimensions" begin
    dt = timemean(A)
    @test dt isa ClimArray
    @test !hasdim(dt, Time)
    @test hasdim(dt, Lon)
    @test hasdim(dt, Lat)
    dl = zonalmean(A)
    @test dl isa ClimArray
    @test !hasdim(dl, Lon)
    @test hasdim(dl, Time)
    @test hasdim(dl, Lat)
    da = latmean(A)
    @test da isa ClimArray
    @test !hasdim(da, Lat)
    @test hasdim(da, Lon)
    @test hasdim(da, Time)
    ds = spacemean(A)
    @test ds isa ClimArray
    @test !hasdim(ds, Lon)
    @test !hasdim(ds, Lat)
    @test hasdim(ds, Time)
    ds = dropagg(mean, A, Time)
    @test dt isa ClimArray
    @test !hasdim(dt, Time)
    @test hasdim(dt, Lon)
    @test hasdim(dt, Lat)
end

@testset "dropagg with weighting" begin
    # spatiotemporal weights:
    W = zero(A)
    W[Time(5)] .= 1.0
    res = dropagg(mean, A, Time, W)
    @test all(res .≈ A[Time(5)])
    @test dims(A, Lon).val == dims(res, Lon).val
    res = dropagg(std, A, Time, W)
    @test all(x -> isapprox(x, 0; atol = 1e-8), res)

    # just time weight
    w = zeros(length(t))
    w[5] = 1.0
    res = dropagg(mean, A, Time, w)
    @test all(res .≈ A[Time(5)])
    @test dims(A, Lon).val == dims(res, Lon).val
end

end # general statistics


text:
# Autogenerated wrapper script for Jansson_jll for x86_64-apple-darwin
export libjansson

JLLWrappers.@generate_wrapper_header("Jansson")
JLLWrappers.@declare_library_product(libjansson, "@rpath/libjansson.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libjansson,
        "lib/libjansson.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct CustomInterpreter <: AbstractInterpreter
    #= native =#

    native::NativeInterpreter
    optimize::Bool
    compress::Bool
    discard_trees::Bool

    #= debugging =#
    depth::Ref{Int}

    function CustomInterpreter(world         = get_world_counter();
                               inf_params    = gen_inf_params(),
                               opt_params    = gen_opt_params(),
                               optimize      = true,
                               compress      = false,
                               discard_trees = false,
                               # dummy kwargs so that kwargs for other functions can be passed on
                               __kwargs...)
        native = NativeInterpreter(world; inf_params, opt_params)
        return new(native,
                   optimize,
                   compress,
                   discard_trees,
                   Ref(0),
                   )
    end
end

function gen_inf_params(; # more constant prop, more correct reports ?
                          aggressive_constant_propagation::Bool = true,
                          # turn this off to get profiles on `throw` blocks, this might be good to default
                          # to `true` since `throw` calls themselves will be reported anyway
                          unoptimize_throw_blocks::Bool = true,
                          # dummy kwargs so that kwargs for other functions can be passed on
                          __kwargs...)
    return @static VERSION ≥ v"1.6.0-DEV.837" ?
           InferenceParams(; aggressive_constant_propagation,
                             unoptimize_throw_blocks,
                             ) :
           InferenceParams(; aggressive_constant_propagation,
                             )
end

function gen_opt_params(; # inlining should be disabled for `CustomInterpreter`, otherwise virtual stack frame
                          # traversing will fail for frames after optimizer runs on
                          inlining = false,
                          # dummy kwargs so that kwargs for other functions can be passed on
                          __kwargs...)
    return OptimizationParams(; inlining,
                                )
end

InferenceParams(interp::CustomInterpreter) = InferenceParams(interp.native)
OptimizationParams(interp::CustomInterpreter) = OptimizationParams(interp.native)
get_world_counter(interp::CustomInterpreter) = get_world_counter(interp.native)

lock_mi_inference(::CustomInterpreter, ::MethodInstance) = nothing
unlock_mi_inference(::CustomInterpreter, ::MethodInstance) = nothing

function add_remark!(interp::CustomInterpreter, sv, s)
    ret = add_remark!(interp.native, sv, s)
end

may_optimize(interp::CustomInterpreter) = interp.optimize
may_compress(interp::CustomInterpreter) = interp.compress
may_discard_trees(interp::CustomInterpreter) = interp.discard_trees

max_stars_repo_name: aviatesk/AbstractInterpreterPlayground

max_issues_repo_name: aviatesk/AbstractInterpreterPlayground

max_forks_repo_name: aviatesk/AbstractInterpreterPlayground