text:
module Maths
import FiniteDifferences
import LinearAlgebra: Tridiagonal, mul!, ldiv!
import SpecialFunctions: erf, erfc, gamma
import StaticArrays: SVector
import Random: AbstractRNG, randn, GLOBAL_RNG
import FFTW
import Luna
import Luna.Utils: saveFFTwisdom, loadFFTwisdom
import Roots: fzero
import Dierckx
import Peaks

#= Pre-created finite difference methods for speed.
   Use (order+6)th order central finite differences with 2 adaptive steps up to
   11th order. Above that creates overflow errors, so we cap it. =#
FDMs = [FiniteDifferences.central_fdm(min(order+6, 11), order, adapt=1) for order=1:7]

"""
    derivative(f, x, order::Integer)

Calculate derivative of function f(x) at value x using finite differences
"""
function derivative(f, x, order::Integer)
    if order == 0
        return f(x)
    else
        scale = abs(x) > 0 ? x : 1.0
        FDMs[order](y->f(y*scale), x/scale)/scale^order
    end
end

"""
    gauss(x, σ; x0=0, power=2)

Gaussian or hypergaussian function over `x` with std dev `σ`, power in the exponent of `power`,
centred at `x0`.
"""
gauss(x, σ; x0=0, power=2) = exp(-1/2 * ((x-x0)/σ)^power)

"""
    gauss(x; fwhm, x0=0, power=2)

Gaussian or hypergaussian function over `x` with FWHM `fwhm`, power in the exponent of `power`,
centred at `x0`.
"""
gauss(x; x0=0, power=2, fwhm) = gauss(x, fwhm_to_σ(fwhm; power=power), x0 = x0, power=power)

fwhm_to_σ(fwhm; power=2) = fwhm / (2 * (2 * log(2))^(1 / power))

"""
    gaussnorm(σ; power=2)
    gaussnorm(;fwhm, power=2)

Norm of a Gaussian with standard deviation `σ` or FWHM `fwhm` and power in the exponent
`power`.
"""
gaussnorm(σ; power=2) = 2σ*2^(1/power)*gamma((power+1)/power)
gaussnorm(;fwhm, power=2) = gaussnorm(fwhm_to_σ(fwhm; power=power); power=power)

"""
    randgauss([rng=GLOBAL_RNG], μ, σ, args...)

Generate random numbers with normal distribution centred on `μ` with standard deviation `σ`.

Any additional `args` are passed onto `randn` and can be used to specify the output dimensions.
"""
randgauss(μ, σ, args...) = randgauss(GLOBAL_RNG, μ, σ, args...)
randgauss(rng::AbstractRNG, μ, σ, args...) = σ*randn(rng, args...) .+ μ

"""
    moment(x::Vector, y::Vector, n = 1)

Calculate the nth moment of the vector y.
"""
function moment(x::Vector, y::Vector, n=1)
    if length(x) ≠ length(y)
        throw(DomainError(x, "x and y must have same length"))
    end
    return sum(x.^n .* y) / sum(y)
end

"""
    moment(x::Vector, y, n=1; dim=1)

Calculate the nth moment of multi-dimensional array y along dimension dim.
"""
function moment(x::Vector, y, n=1; dim=1)
    if size(y, dim) ≠ length(x)
        throw(DomainError(y, "y must be of same length as x along dim"))
    end
    xshape = ones(Integer, ndims(y))
    xshape[dim] = length(x)
    return sum(reshape(x, Tuple(xshape)).^n .* y, dims=dim) ./ sum(y, dims=dim)
end

"""
    rms_width(x::Vector, y; dim=1)

RMS width of distribution y on axis x
"""
function rms_width(x::Vector, y::Vector; dim=1)
    return sqrt(moment(x, y, 2) - moment(x, y, 1)^2)
end

function rms_width(x::Vector, y; dim=1)
    return sqrt.(moment(x, y, 2, dim=dim) - moment(x, y, 1, dim=dim).^2)
end

"""
    fwhm(x, y; method=:linear, baseline=false, minmax=:min, level=0.5)

Calculate the full width at half maximum (FWHM) of `y` on the axis `x`.

`method` can be `:spline` or `:nearest`. `:spline` uses a [`CSpline`](@ref), whereas
`:nearest` finds the closest values either side of the crossing point and interpolates linearly.

If `baseline` is true, the width is not taken at
`level * maximum(y)`, but at half of the span of `y`, `level * (maximum(y) - minimum(y))`.

`minmax` determines whether the FWHM is taken at the narrowest (`:min`) or the widest (`:max`)
point of y.

The default `level=0.5` requests the full width at half maximum. Setting `level` to something
different computes the corresponding width. E.g. `level=0.1` for the 10% width.
"""
function fwhm(x, y::AbstractVector; kwargs...)
    left, right = level_xings(x, y; kwargs...)
    return abs(right-left)
end

function fwhm(x, y::AbstractArray; dim=1, kwargs...)
    idxlo = CartesianIndices(size(y)[1:dim-1])
    idxhi = CartesianIndices(size(y)[dim+1:end])
    outshape = collect(size(y))
    outshape[dim] = 1
    out = zeros(eltype(y), Tuple(outshape))
    for hi in idxhi
        for lo in idxlo
            out[lo, 1, hi] = fwhm(x, y[lo, :, hi]; kwargs...)
        end
    end
    out
end

"""
    level_xings(x, y; method=:linear, baseline=false, minmax=:min)

Find crossings of the curve `y` on the axis `x` with the value `level * maximum(y)`.

`method` can be `:spline`, `:linear` or `:nearest`. `:spline` uses a [`CSpline`](@ref).
`:linear` finds the closest values either side of the crossing point and interpolates linearly.
`:nearest` just uses the closest values either side of the crossing point.

If `baseline` is true, the width is not taken at
`level * maximum(y)`, but at half of the span of `y`, `level * (maximum(y) - minimum(y))`.

`minmax` determines whether the crossings are taken at the narrowest (`:min`) or the widest (`:max`)
point of y.
"""
function level_xings(x::AbstractVector, y::AbstractVector;
                     method=:linear, baseline=false, minmax=:min, level=0.5)
    minmax in (:min, :max) || error("minmax has to be :min or :max")
    if baseline
        val = minimum(y) + level*(maximum(y) - minimum(y))
    else
        val = level*maximum(y)
    end
    if !(method in (:spline, :linear, :nearest))
        error("Unknown FWHM method $method")
    end
    maxidx = argmax(y)
    xmax = x[maxidx]
    if method in (:nearest, :linear)
        try
            if minmax == :min
                lefti = findlast((x .< xmax) .& (y .< val))
                righti = findfirst((x .> xmax) .& (y .< val))
                (method == :nearest) && return (x[lefti], x[righti])

                left = linterpx(x[lefti], x[lefti+1], y[lefti], y[lefti+1], val)
                right = linterpx(x[righti-1], x[righti], y[righti-1], y[righti], val)
            else
                lefti = findfirst((x .< xmax) .& (y .> val))
                righti = findlast((x .> xmax) .& (y .> val))
                (method == :nearest) && return (x[lefti], x[righti])

                left = linterpx(x[lefti-1], x[lefti], y[lefti-1], y[lefti], val)
                right = linterpx(x[righti], x[righti+1], y[righti], y[righti+1], val)
            end
            return left, right
        catch
            return NaN, NaN
        end
    elseif method == :spline
        #spline method
        try
            spl = BSpline(x, y .- val)
            r = roots(spl; maxn=10000)
            rleft = r[r .< xmax]
            rright = r[r .> xmax]
            if minmax == :min
                return maximum(rleft), minimum(rright)
            else
                return minimum(rleft), maximum(rright)
            end
        catch e
            return NaN, NaN
        end
    else
        error("Unknown level_xings method $method")
    end
end

"""
    linterpx(x1, x2, y1, y2, val)

Given two points on a straight line, `(x1, y1)` and `(x2, y2)`, calculate the value of `x`
at which this straight line intercepts with `val`.

# Examples
```jldoctest
julia> x1 = 0; x2 = 1; y1 = 0; y2 = 2; # y = 2x
julia> linterpx(x1, x2, y1, y2, 0.5)
0.25
```
"""
function linterpx(x1, x2, y1, y2, val)
    slope = (y2-y1)/(x2-x1)
    icpt = y1 - slope*x1
    (val-icpt)/slope
end

"""
    hwhm(f, x0=0; direction=:fwd)

Find the value `x` where the function `f(x)` drops to half of its maximum, which is located
at `x0`. For `direction==:fwd`, search in the region `x > x0`, for :bwd, search in `x < x0`.
"""
function hwhm(f, x0=0; direction=:fwd)
    m = f(x0)
    fhalf(x) = f(x) - 0.5*m
    if direction == :fwd
        xhw = fzero(fhalf, x0, Inf)
    elseif direction == :bwd
        xhw = fzero(fhalf, -Inf, x0)
    end
    return abs(xhw - x0)
end

"""
    cumtrapz!([out, ] y, x; dim=1)

Trapezoidal integration for multi-dimensional arrays or vectors, in-place or with output array.

If `out` is omitted, `y` is integrated in place. Otherwise the result is placed into `out`.

`x` can be an array (the x axis) or a number (the x axis spacing).
"""
function cumtrapz! end

function cumtrapz!(y, x; dim=1)
    idxlo = CartesianIndices(size(y)[1:dim-1])
    idxhi = CartesianIndices(size(y)[dim+1:end])
    _cumtrapz!(y, x, idxlo, idxhi)
end

"""
    _cumtrapz!([out, ] y, x, idxlo, idxhi)

Inner function for multi-dimensional `cumtrapz!` - uses 1-D routine internally
"""
function _cumtrapz!(y, x, idxlo, idxhi)
    for lo in idxlo
        for hi in idxhi
            cumtrapz!(view(y, lo, :, hi), x)
        end
    end
end

function cumtrapz!(y::T, x) where T <: Union{SubArray, Vector}
    tmp = y[1]
    y[1] = 0
    for i in 2:length(y)
        tmp2 = y[i]
        y[i] = y[i-1] + 1//2 * (tmp + tmp2) * _dx(x, i)
        tmp = tmp2
    end
end

function cumtrapz!(out, y, x; dim=1)
    idxlo = CartesianIndices(size(y)[1:dim-1])
    idxhi = CartesianIndices(size(y)[dim+1:end])
    _cumtrapz!(out, y, x, idxlo, idxhi)
end

function _cumtrapz!(out, y, x, idxlo, idxhi)
    for lo in idxlo
        for hi in idxhi
            cumtrapz!(view(out, lo, :, hi), view(y, lo, :, hi), x)
        end
    end
end

function cumtrapz!(out, y::Union{SubArray, Vector}, x)
    out[1] = 0
    for i in 2:length(y)
        out[i] = out[i-1]+ 1//2*(y[i-1] + y[i])*_dx(x, i)
    end
end

"""
    _dx(x, i)

Calculate the axis spacing at index `i` given an axis `x`. If `x` is a number, interpret this
as `δx` directly
"""
_dx(x, i) = x[i] - x[i-1]
_dx(δx::Number, i) = δx

"""
    cumtrapz(y, x; dim=1)

Calculate the cumulative trapezoidal integral of `y`.

`x` can be an array (the x axis) or a number (the x axis spacing).
"""
function cumtrapz(y, x; dim=1)
    out = similar(y)
    cumtrapz!(out, y, x; dim=dim)
    return out
end

"""
    normbymax(x, dims)

Return an array normalised by its maximum value
"""
function normbymax(x, dims)
    return x ./ maximum(x; dims=dims)
end

function normbymax(x)
    return x ./ maximum(x)
end

"""
    log10_norm(x)

Normalised log10 i.e. maximum of output is 0
"""
function log10_norm(x)
    return log10.(normbymax(x))
end

function log10_norm(x, dims)
    return log10.(normbymax(x, dims))
end

"""
    errfun_window(x, xmin, xmax, width)
    errfun_window(x, xmin, xmax, width_left, width_right)

Filtering window based on the error function. If `width_left` and `width_right` are given
separately, the smoothing width is different on each side.
"""
function errfun_window(x, xmin, xmax, width)
    return @. 0.5 * (erf((x - xmin) / width) + erfc((x - xmax) / width) - 1)
end

function errfun_window(x, xmin, xmax, width_left, width_right)
    return @. 0.5 * (erf((x - xmin) / width_left) + erfc((x - xmax) / width_right) - 1)
end

"""
    planck_taper(x, xmin, xmax, ε)

Planck taper window as defined in the paper (https://arxiv.org/pdf/1003.2939.pdf eq(7)).

#Arguments
-`xmin` : lower limit (window is 0 here)
-`xmax` : upper limit (window is 0 here)
-`ε` : fraction of window width over which to increase from 0 to 1
"""
function planck_taper(x, xmin, xmax, ε)
    x0 = (xmax + xmin) / 2
    xc = x .- x0
    X = (xmax - xmin)
    x1  = -X / 2
    x2 = -X / 2 * (1 - 2ε)
    x3 = X / 2 * (1 - 2ε)
    x4 = X / 2
    return _taper.(xc, x1, x2, x3, x4)
end

"""
    planck_taper(x, left0, left1, right1, right0)

Planck taper window, but finding the taper width by defining 4 points:
The window increases from 0 to 1 between `left0` and `left1`, and then drops again
to 0 between `right1` and `right0`.
"""
function planck_taper(x, left0, left1, right1, right0)
    x0 = (right0 + left0) / 2
    xc = x .- x0
    X = right0 - left0
    εleft = abs(left1 - left0) / X
    εright = abs(right0 - right1) / X
    x1  = -X / 2
    x2 = -X / 2 * (1 - 2εleft)
    x3 = X / 2 * (1 - 2εright)
    x4 = X / 2
    return _taper.(xc, x1, x2, x3, x4)
end

function _taper(xc, x1, x2, x3, x4)
    if xc < x1
        return zero(xc)
    elseif x1 < xc < x2
        z12 = (x2 - x1) / (xc - x1) + (x2 - x1) / (xc - x2)
        return 1 / (1 + exp(z12))
    elseif x2 <= xc <= x3
        return one(xc)
    elseif x3 < xc < x4
        z34 = (x3 - x4) / (xc - x3) + (x3 - x4) / (xc - x4)
        return 1 / (1 + exp(z34))
    else
        return zero(xc)
    end
end

"""
    hypergauss_window(x, xmin, xmax, power = 10)

Hypergaussian window function.
"""
function hypergauss_window(x, xmin, xmax, power = 10)
    fw = xmax - xmin
    x0 = (xmax + xmin) / 2
    return gauss(x, x0 = x0, fwhm = fw, power = power)
end

"""
    gabor(t, A, ts, fw)

Compute the Gabor transform (aka spectrogram or time-gated Fourier transform) of the vector
`A`, sampled on axis `t`, with windows centred at `ts` and a window FWHM of `fw`.
"""
function gabor(t::Vector, A::Vector, ts, fw)
    tmp = Array{eltype(A), 2}(undef, (length(t), length(ts)))
    for (ii, ti) in enumerate(ts)
        tmp[:, ii] = A .* Maths.gauss.(t; x0=ti, fwhm=fw)
    end
    _gaborFT(tmp)
end

_gaborFT(x::Array{T, 2}) where T <: Real = FFTW.rfft(x, 1)
_gaborFT(x::Array{T, 2}) where T <: Complex = FFTW.fft(x, 1)

"""
    hilbert(x; dim=1)

Compute the Hilbert transform, i.e. find the analytic signal from a real signal.
"""
function hilbert(x::Array{T,N}; dim = 1) where T <: Real where N
    xf = FFTW.fft(x, dim)
    n1 = size(xf, dim)÷2
    n2 = size(xf, dim)
    idxlo = CartesianIndices(size(xf)[1:dim - 1])
    idxhi = CartesianIndices(size(xf)[dim + 1:end])
    xf[idxlo, 2:n1, idxhi] .*= 2
    xf[idxlo, (n1+1):n2, idxhi] .= 0
    return FFTW.ifft(xf, dim)
end

"""
    plan_hilbert!(x; dim=1)

Pre-plan a Hilbert transform.

Returns a closure `hilbert!(out, x)` which places the Hilbert transform of `x` in `out`.
"""
function plan_hilbert!(x; dim=1)
    loadFFTwisdom()
    FT = FFTW.plan_fft(copy(x), dim, flags=Luna.settings["fftw_flag"])
    saveFFTwisdom()
    xf = Array{ComplexF64}(undef, size(FT))
    idxlo = CartesianIndices(size(xf)[1:dim - 1])
    idxhi = CartesianIndices(size(xf)[dim + 1:end])
    n1 = size(xf, dim)÷2
    n2 = size(xf, dim)
    xc = complex(x)
    function hilbert!(out, x)
        copyto!(xc, x)
        mul!(xf, FT, xc)
        xf[idxlo, 2:n1, idxhi] .*= 2
        xf[idxlo, (n1+1):n2, idxhi] .= 0
        ldiv!(out, FT, xf)
    end
    return hilbert!
end

"""
    plan_hilbert(x; dim=1)

Pre-plan a Hilbert transform.

Returns a closure `hilbert(x)` which returns the Hilbert transform of `x` without allocation.

!!! warning
    The closure returned always returns a reference to the same array buffer, which could lead
    to unexpected results if it is called from more than one location. To avoid this the array
    should either: (i) only be used in the same code segment; (ii) only be used transiently
    as part of a larger computation; (iii) be copied.
"""
function plan_hilbert(x; dim=1)
    out = complex(x)
    hilbert! = plan_hilbert!(x, dim=dim)
    function hilbert(x)
        hilbert!(out, x)
    end
    return hilbert
end

"""
    oversample(t, x; factor::Int=4, dim=1)

Oversample (smooth) an array by 0-padding in the frequency domain. Works for both real-valued
and complex-valued arrays.
"""
function oversample(t, x::Array{<:Real, N}; factor::Int=4, dim=1) where N
    if factor == 1
        return t, x
    end
    xf = FFTW.rfft(x, dim)

    len = size(xf, dim)
    newlen_t = factor * length(t)
    if iseven(newlen_t)
        newlen_ω = Int(newlen_t / 2 + 1)
    else
        newlen_ω = Int((newlen_t + 1) / 2)
    end
    δt = t[2] - t[1]
    δto = δt / factor
    Nto = collect(0:newlen_t - 1)
    to = t[1] .+ Nto .* δto

    shape = collect(size(xf))
    shape[dim] = newlen_ω
    xfo = Array{eltype(xf), N}(undef, Tuple(shape))
    fill!(xfo, 0.0)
    idxlo = CartesianIndices(size(xfo)[1:dim - 1])
    idxhi = CartesianIndices(size(xfo)[dim + 1:end])
    xfo[idxlo, 1:len, idxhi] .= factor .* xf
    return to, FFTW.irfft(xfo, newlen_t, dim)
end

function oversample(t, x::Array{<:Complex,N}; factor::Int=4, dim=1) where N
    if factor == 1
        return t, x
    end
    xf = FFTW.fftshift(FFTW.fft(x, dim), dim)

    len = size(xf, dim)
    newlen = factor * length(t)
    δt = t[2] - t[1]
    δto = δt / factor
    Nto = collect(0:newlen - 1)
    to = t[1] .+ Nto .* δto

    sidx  = (newlen - len)//2 + 1
    iseven(newlen) || (sidx -= 1//2)
    iseven(len) || (sidx += 1//2)
    startidx = Int(sidx)
    endidx = startidx+len-1

    shape = collect(size(xf))
    shape[dim] = newlen
    xfo = Array{eltype(xf), N}(undef, Tuple(shape))
    fill!(xfo, 0.0)
    idxlo = CartesianIndices(size(xfo)[1:dim - 1])
    idxhi = CartesianIndices(size(xfo)[dim + 1:end])
    xfo[idxlo, startidx:endidx, idxhi] .= factor .* xf
    return to, FFTW.ifft(FFTW.ifftshift(xfo, dim), dim)
end

"""
    aitken_accelerate(f, x0; n0 = 0, rtol = 1e-6, maxiter = 10000)

Find limit of a series by Aitken acceleration.
"""
function aitken_accelerate(f, x0; n0 = 0, rtol = 1e-6, maxiter = 10000)
    n = n0
    x0 = f(x0, n)
    x1 = f(x0, n + 1)
    x2 = f(x1, n + 2)
    Ax = aitken(x0, x1, x2)
    success = false
    while ~success && n < maxiter
        n += 1
        Axprev = Ax
        x0 = x1
        x1 = x2
        x2 = f(x2, n + 2)
        Ax = aitken(x0, x1, x2)

        if 2 * abs(Ax - Axprev) / abs(Ax + Axprev) < rtol
            success = true
        end
    end
    return Ax, success, n
end

function aitken(x0, x1, x2)
    den = (x0 - x1) - (x1 - x2)
    return x0 - (x1 - x0)^2 / den
end

"""
    converge_series(f, x0; n0 = 0, rtol = 1e-6, maxiter = 10000)

Find limit of series by brute force. The iteration is stopped when the relative change in
accumulated value is smaller than `rtol`.
"""
function converge_series(f, x0; n0 = 0, rtol = 1e-6, maxiter = 10000)
    n = n0
    x1 = x0
    success = false
    while ~success && n < maxiter
        x1 = f(x0, n)

        if 2 * abs(x1 - x0) / abs(x1 + x0) < rtol
            success = true
        end

        n += 1
        x0 = x1
    end
    return x1, success, n
end

"""
    check_spline_args(x, y)

Ensure that the x array contains unique and sorted values (while preserving the
relaionship between the x and y values).
"""
function check_spline_args(x, y)
    if any(diff(x) .== 0)
        error("entries in x must be unique")
    end
    if !issorted(x)
        idcs = sortperm(x)
        x = x[idcs]
        y = y[idcs]
    end
    x, y
end

"""
     make_spline_ifun(x, ifun)

If `ifun != nothing` then `ifun(x0)` should return the index of the first element in x
which is bigger than x0. Otherwise, it defaults two one of two options:
1. If `x` is uniformly spaced, the index is calculated based on the spacing of `x`
2. If `x` is not uniformly spaced, a `FastFinder` is used instead.
"""
function make_spline_ifun(x, ifun)
    if ifun === nothing
        δx = x[2] - x[1]
        if all(diff(x) .≈ δx)
            # x is uniformly spaced - use fast lookup
            xmax = maximum(x)
            xmin = minimum(x)
            N = length(x)
            ffast(x0) = x0 <= xmin ? 2 :
                        x0 >= xmax ? N :
                        ceil(Int, (x0-xmin)/(xmax-xmin)*(N-1))+1
            ifun = ffast
        else
            # x is not uniformly spaced - use brute-force lookup
            ifun = FastFinder(x)
        end
    end
    ifun
end

"""
    CSpline

Simple cubic spline, see e.g.:
http://mathworld.wolfram.com/CubicSpline.html Boundary
conditions extrapolate with initially constant gradient
"""
struct CSpline{Tx,Ty,Vx<:AbstractVector{Tx},Vy<:AbstractVector{Ty}, fT}
    x::Vx
    y::Vy
    D::Vy
    ifun::fT
    bounds_error::Bool
end

# make broadcast like a scalar
Broadcast.broadcastable(c::CSpline) = Ref(c)

"""
    CSpline(x, y [, ifun]; bounds_error=false)

Construct a `CSpline` to interpolate the values `y` on axis `x`.

If given, `ifun(x0)` should return the index of the first element in x which is bigger
than x0. Otherwise, it defaults two one of two options:
1. If `x` is uniformly spaced, the index is calculated based on the spacing of `x`
2. If `x` is not uniformly spaced, a `FastFinder` is used instead.
"""
function CSpline(x, y, ifun=nothing; bounds_error=false)
    x, y = check_spline_args(x, y)
    R = similar(y)
    R[1] = y[2] - y[1]
    for i in 2:(length(y)-1)
        R[i] = y[i+1] - y[i-1]
    end
    R[end] = y[end] - y[end - 1]
    @. R *= 3
    d = fill(4.0, size(y))
    d[1] = 2.0
    d[end] = 2.0
    dl = fill(1.0, length(y) - 1)
    M = Tridiagonal(dl, d, dl)
    D = M \ R
    CSpline(x, y, D, make_spline_ifun(x, ifun), bounds_error)
end

"""
    (c::CSpline)(x0)

Evaluate the `CSpline` at coordinate `x0`
"""
function (c::CSpline)(x0)
    if c.bounds_error
        x0 < c.x[1] && throw(DomainError("CSpline evaulated out of bounds, $x0 < $(c.x[1])"))
        x0 > c.x[end] && throw(DomainError("CSpline evaulated out of bounds, $x0 > $(c.x[end])"))
    end
    i = c.ifun(x0)
    x0 == c.x[i] && return c.y[i]
    x0 == c.x[i-1] && return c.y[i-1]
    t = (x0 - c.x[i - 1])/(c.x[i] - c.x[i - 1])
    (c.y[i - 1]
        + c.D[i - 1]*t
        + (3*(c.y[i] - c.y[i - 1]) - 2*c.D[i - 1] - c.D[i])*t^2
        + (2*(c.y[i - 1] - c.y[i]) + c.D[i - 1] + c.D[i])*t^3)
end

"""
    linterp(x, x1, y1, x2, y2)

Linear interpolation of `y` at position `xp`, between points `(x1, y1)` and `(x2, y2)`.
For `xp` outside interval `[x1, x2]` this corresponds to linear extrapolation.
"""
function linterp(xp, x1, y1, x2, y2)
    t = (xp - x1) / (x2 - x1)
    (1 - t) * y1 + t * y2
end

"""
    LinTerp

Linear interpolation.
"""
struct LinTerp{Tx,Ty,Vx<:AbstractVector{Tx},Vy<:AbstractVector{Ty}, fT}
    x::Vx
    y::Vy
    ifun::fT
end

# make  broadcast like a scalar
Broadcast.broadcastable(l::LinTerp) = Ref(l)

"""
    LinTerp(xp, xs, ys)

Construct a linear interpolator over `xs`, `ys`.
"""
function LinTerp(xs, ys)
    LinTerp(xs, ys, make_spline_ifun(xs, nothing))
end

"""
    (l::LinTerp)(x)

Evaluate a linear interpolator `LinTerp` at point `x`.
"""
function (l::LinTerp)(x)
    i = l.ifun(x)
    linterp(x, l.x[i - 1], l.y[i - 1], l.x[i], l.y[i])
end

"""
    fftfreq(x)

Calculate frequency vector k from samples x for complex-valued FFT
"""
function fftfreq(x)
    Dx = abs(x[2] - x[1])
    all(diff(x) .≈ Dx) || error("x must be spaced uniformly")
    N = length(x)
    n = collect(range(0, length=N))
    (n .- N/2) .* 2π/(N*Dx)
end

"""
    rfftfreq(x)

Calculate frequency vector k from samples x for rFFT
"""
function rfftfreq(x)
    Dx = abs(x[2] - x[1])
    all(diff(x) .≈ Dx) || error("x must be spaced uniformly")
    Nx = length(x)
    Nf = length(x)÷2 + 1
    n = collect(range(0, length=Nf))
    n .* 2π/(Nx*Dx)
end

rfftnorm(δt) = 2*fftnorm(δt)
fftnorm(δt) = δt/sqrt(2π)

"""
    FastFinder

Callable type which accelerates index finding for the case where inputs are usually in order.
"""
mutable struct FastFinder{xT, xeT}
    x::xT
    mi::xeT
    ma::xeT
    N::Int
    ilast::Int
    xlast::xeT
end

"""
    FastFinder(x)

Construct a `FastFinder` to find indices in the array `x`.

!!! warning
    `x` must be sorted in ascending order for `FastFinder` to work.
"""
function FastFinder(x::AbstractArray)
    issorted(x) || error("Input array for FastFinder must be sorted in ascending order.")
    if any(diff(x) .== 0)
        error("Entries in input array for FastFinder must be unique.")
    end
    FastFinder(x, x[1], x[end], length(x), 0, typemin(x[1]))
end

"""
    (f::FastFinder)(x0)

Find the first index in `f.x` which is larger than `x0`.

This is similar to [`findfirst`](@ref), but it starts at the index which was last used.
If the new value `x0` is close to the previous `x0`, this is much faster than `findfirst`.
"""
function (f::FastFinder)(x0::Number)
    # Default cases if we're out of bounds
    if x0 <= f.x[1]
        f.xlast = x0
        f.ilast = 1
        return 2
    elseif x0 >= f.x[end]
        f.xlast = x0
        f.ilast = f.N
        return f.N
    end
    if f.ilast == 0 # first call -- f.xlast is not set properly so comparisons won't work
        # return using brute-force method instead
        f.ilast = findfirst(x -> x>x0, f.x)
        return f.ilast
    end
    if x0 == f.xlast # same value as before - no work to be done
        return f.ilast
    elseif x0 < f.xlast # smaller than previous value - go through array backwards
        f.xlast = x0
        for i = f.ilast:-1:1
            if f.x[i] < x0
                f.ilast = i+1 # found last idx where x < x0 -> at i+1, x > x0
                return i+1
            end
        end
        # we only get to this point if we haven't found x0 - return the lower bound (2)
        f.ilast = 1
        return 2
    else # larger than previous value - just pick up where we left off
        f.xlast = x0
        for i = f.ilast:f.N
            if f.x[i] > x0
                f.ilast = i
                return i
            end
        end
        # we only get to this point if we haven't found x0 - return the upper bound (N)
        f.ilast = f.N
        return f.N
    end
end

struct RealBSpline{sT,hT,fT}
    rspl::sT
    h::hT
    hh::hT
    ifun::fT
end

struct CmplxBSpline{sT}
    rspl::sT
    ispl::sT
end

Broadcast.broadcastable(rs::RealBSpline) = Ref(rs)
Broadcast.broadcastable(cs::CmplxBSpline) = Ref(cs)

"""
    BSpline(x, y)

Construct a `RealBSpline` or `CmplxSpline` of given `order` (1 to 5, default=3)
to interpolate the values `y` on axis `x`.

If given, `ifun(x0)` should return the index of the first element in x which is bigger
than x0. Otherwise, it defaults two one of two options:
1. If `x` is uniformly spaced, the index is calculated based on the spacing of `x`
2. If `x` is not uniformly spaced, a `FastFinder` is used instead.

# Note
For accurate derivatives use 5th order splines
"""
function BSpline(x::AbstractVector{Tx}, y::AbstractVector{T}; ifun = nothing, order=3) where {Tx <: Real, T <: Real}
    check_spline_args(x, y)
    rspl = Dierckx.Spline1D(x, y, bc="extrapolate", k=order, s=0.0)
    h = zeros(eltype(y), order + 1)
    hh = similar(h)
    RealBSpline(rspl, h, hh, make_spline_ifun(Dierckx.get_knots(rspl), ifun))
end

function BSpline(x::AbstractVector{Tx}, y::AbstractVector{T}; kwargs...) where {Tx <: Real, T <: Complex}
    CmplxBSpline(BSpline(x, real(y); kwargs...),
                 BSpline(x, imag(y); kwargs...))
end

"""
    (cs::RealBSpline)(x)

Evaluate the `RealBSpline` at coordinate(s) `x`
"""
function (rs::RealBSpline)(x)
    splev!(rs.h, rs.hh, rs.rspl.t, rs.rspl.c, rs.rspl.k, x, rs.ifun)
end

"""
    (cs::CmplxSpline)(x)

Evaluate the `CmplxBSpline` at coordinate(s) `x`
"""
function (cs::CmplxBSpline)(x)
    complex(cs.rspl(x), cs.ispl(x))
end

"""
    derivative(rs::RealBSpline, x, order::Integer)

Calculate derivative of the spline `rs`. For `1 <= order < k` (where `k` is the spline order)
this uses an optimised routine.
For `order >= k` this falls back to the generic finite difference based method.
"""
function derivative(rs::RealBSpline, x, order::Integer)
    if order == 0
        return rs(x)
    elseif order < rs.rspl.k
        return Dierckx.derivative(rs.rspl, x, order)
    else
        invoke(derivative, Tuple{Any,Any,Integer}, rs, x, order)
    end
end

"""
    derivative(cs::CmplxBSpline, x, order::Integer)

Calculate derivative of the spline `cs`. For `1 <= order < k` (where `k` is the spline order)
this uses an optimised routine.
For `order >= k` this falls back to the generic finite difference based method.
"""
function derivative(cs::CmplxBSpline, x, order::Integer)
    complex(derivative(cs.rspl, x, order), derivative(cs.ispl, x, order))
end

"""
    differentiate_spline(rs::RealBSpline; order::Integer=1)

Return a new spline which is the derivative of `rs` and has the same support.
Higher orders are obtained by repeated spline differentiation.
"""
function differentiate_spline(rs::RealBSpline, order::Integer)
    x = Dierckx.get_knots(rs.rspl)
    dspl = rs
    for i = 1:order
        dy = derivative.(dspl, x, 1)
        dspl = BSpline(x, dy, order=rs.rspl.k)
    end
    dspl
end

"""
    roots(rs::RealBSpline)

Find the roots of the spline `rs`.
"""
function roots(rs::RealBSpline; maxn::Int=128)
    Dierckx.roots(rs.rspl; maxn=maxn)
end

"""
    splev!(h, hh, t, c, k, x, ifun)

Evaluate a spline s(x) of degree k, given in its b-spline representation.

# Arguments
- `h::ArrayPT<:Real,1}`: work space
- `hh::ArrayPT<:Real,1}`: work space
- `t::Array{T<:Real,1}`: the positions of the knots
- `c::Array{T<:Real,1}`: the b-spline coefficients
- `k::Integer`: the degree of s(x)
- `x::Real`: the point to evaluate at
- `ifun::Function`: a function to find the index `i` s.t. `t[i] > x`

Adapted from splev.f from Dierckx
http://www.netlib.org/dierckx/index.html
"""
function splev!(h, hh, t, c, k, x, ifun)
    k1 = k + 1
    k2 = k1 + 1
    nk1 = length(t) - k1
    tb = t[k1]
    te = t[nk1 + 1]
    l = k1
    l1 = l + 1
    # search for knot interval t(l) <= arg < t(l+1)
    l = ifun(x) - 1 + k
    # evaluate the non-zero b-splines at x.
    fpbspl!(h, hh, t, k, x, l)
    # find the value of s(x) at x.
    sp = 0.0
    ll = l - k1
    for j = 1:k1
        ll = ll + 1
        sp += c[ll]*h[j]
    end
    sp
end

"""
    fpbspl!(h, hh, t, k, x, l)

Evaluate the (k+1) non-zero b-splines of
degree k at t[l] <= x < t[l+1] using the stable recurrence
relation of de boor and cox.

# Arguments
- `h::ArrayPT<:Real,1}`: work space
- `hh::ArrayPT<:Real,1}`: work space
- `t::Array{T<:Real,1}`: the positions of the knots
- `k::Integer`: the degree of s(x)
- `x::Real`: the point to evaluate at
- `l::Integer`: the active knot location: `t[l] <= x < t[l+1]`

Adapted from fpbspl.f from Dierckx
http://www.netlib.org/dierckx/index.html
"""
function fpbspl!(h, hh, t, k, x, l)
    h[1] = one(x)
    for j = 1:k
        for i = 1:j
            hh[i] = h[i]
        end
        h[1] = 0.0
        for i = 1:j
            li = l+i
            lj = li-j
            if t[li] != t[lj]
                f = hh[i]/(t[li] - t[lj])
                h[i] = h[i] + f*(t[li] - x)
                h[i+1] = f*(x - t[lj])
            else
                h[i+1] = 0.0
            end
        end
    end
    return h
end

"""
    pkfw(x, y, pki; level=0.5, skipnonmono=true, closest=5)

Find the full width of a peak in `y` over `x` centred at index `pki`.

The default `level=0.5` requests the full width at half maximum. Setting `level` to something
different computes the corresponding width. E.g. `level=0.1` for the 10% width.

`skipnonmono=true` skips peaks which are not monotonically increaing/decreasing before/after the peak.

`closest=5` sets the minimum number of indices for the full width.
"""
function pkfw(x, y, pki; level=0.5, skipnonmono=true, closest=5)
    val = level*y[pki]
    iup = findnext(x -> x < val, y, pki)
    if iup == nothing
        iup = length(x)
    end
    idn = findprev(x -> x < val, y, pki)
    if idn == nothing
        idn = 1
    end
    if skipnonmono
        if any(diff(y[pki:iup]) .> 0)
            return missing
        end
        if any(diff(y[idn:pki]) .< 0)
            return missing
        end
    end
    if (iup - idn) < closest
        return missing
    end
    up = Maths.linterpx(x[iup - 1], x[iup], y[iup - 1], y[iup], val)
    dn = Maths.linterpx(x[idn], x[idn + 1], y[idn], y[idn + 1], val)
    return up - dn
end

"""
    findpeaks(x, y; threshold=0.0, filterfw=true)

Find isolated peaks in a signal `y` over `x` and return their value, FWHM and index.
`threshold=0.0` allows filtering peaks above a threshold value.
If `filterfw=true` then only peaks with a clean FWHM are returned.
"""
function findpeaks(x, y; threshold=0.0, filterfw=true)
    pkis = Peaks.argmaxima(y, 10)
    proms = Peaks.peakproms(pkis, y)
    pks = [(peak=y[pki], fw=pkfw(x, y, pki), position=x[pki], index=pki) for pki in pkis]
    # filter out peaks with missing fws
    if filterfw
        pks = filter(x -> !(x.fw === missing), pks)
    end
    # filter out peaks below threshold
    pks = filter(x -> x.peak > threshold, pks)
end

end


text:
# # Univariate Gaussian Process

# This package implements a univariate Bayesian Optimization with the Gaussian Process used as the surrogate model
# and visualisation of the optimization process using Plots.

# This package could be a great tool to visualize Gaussian processes to get a better understanding of them
# and/or to experiment with Bayesian optimization with different acquisition functions, covariance functions
# and other hyperparameters.

# ## Bayesian Optimization - brief introduction

# Bayesian optimization is an optimization technique which cleverly selects evaluation points
# from the objective function domain to optimize it while minimizing the number of its evaluations.

# The main advantages of Bayesian optimization is that it can optimize black-box
# and it is effective for optimizing functions which are computationally expensive to evaluate.

# Bayesian optimization works as follows:
# 1. Select the next evaluation point by maximizing the acquisition function.
# 2. Evaluate the objective function at the selected point.
# 3. Augment the data with the new evaluation.
# 4. Update the surrogate model based on the data.
# 5. Repeat.

# Terms:
# - **objective function** - The to-be-optimized function. Can be blackbox, usually is expensive to evaluate.
# - **acquisition function** - A function which to each point of the domain assigns a fitness
#                           for being the point of the next evaluation.
# - **data** - Contains the information from all previous objective function evaluations.
# - **surrogate model** - A probabilistic model estimating the objective function.

# Despite its name, Bayesian optimization can also be used to approximate functions
# instead of optimizing them by using a different acquisition functions.

# ## Gaussian Process - brief introduction

# Gaussian process can be used as the surrogate model for the Bayesian optimization.

# Gaussian process represents the estimated value of the objective function at each point of the domain
# as a one-dimensional normal distribution. This means that the approximation of the whole objective function
# (assuming its continuous) is an infinite-dimensional normal distribution over all possible forms
# the objective function could have.

# ## Example usage

# Let's assume that we wish to maximize a black-box objective function. That means that we can query
# the objective function for values at any point of its domain but we do not know its analytical form.
# Furthermore, let's assume that we are limited in the number of evaluations we can perform because this
# unknow function takes a very long time to evaluate. Let's say we can only perform 10 evaluations.

# Let's define this function and the bounds within which we want maximize it.
# Furthermore, let's assume that we are limited to 5 evaluations.

obj_func(x::Real) = sin(x)

bounds = (-π, π);
evals = 10;

# Here, the objective function is simply sinus, but pretend that we do not know that.

# Now let's maximize our objective function using the UnivariateGaussianProcess package.

# First load the package.

using Revise
using UnivariateGaussianProcess

# Now we have to choose the parameters of the Bayesian optimization
# as well as the Gaussian process surrogate model.

cov_func = matern3;
acq_func = (x, data) -> 3 * ei(x, data) + 1/2;

# We will use the `matern1` covariance function for the Gaussian process
# and the `expected improvement` acquisition function for the Bayesian optimization.

# (Note that the coeficients in the acquisition function definition have no effect on the optimization process.
# They are there just to make the generated plots more readable.)

# Now we can run the algorithm.

data = gaussian_process(obj_func, bounds, cov_func, acq_func; evals, plot_only_final=false);

# You can see the optimization process on the plots.
# (You can set the keyword argument `plot_only_final` to true to reduce the output length.)

# To get the best function value found so far by the algorithm use the following function.

GP_res = get_max(data)[2]

# ## Comparison to BlackBoxOptim

# Now let's try to maximize our objective function using the BlackBoxOptim package
# which can optimize black-box function but is not specialized to minimize the number of evaluations.

# We will perform 10 runs as this optimization algorithm is stochastic.

using BlackBoxOptim
using Random
using Statistics
Random.seed!(555)

obj_func(x::AbstractArray{<:Real}) = - obj_func(first(x))

fitness = Float64[]
for i in 1:10
    res = bboptimize(obj_func; NumDimensions=1, SearchRange=bounds, MaxFuncEvals=evals, TraceMode=:silent);
    append!(fitness, - best_fitness(res))
end

# Now let's see the average result of the BlackBoxOptim package.

BBO_res = mean(fitness)

# We can see that the average result is significantly worse than the result of the Bayseian optimization as expected.

# This of course does not prove that the Bayesian optimization performs better, but it is a simple example
# showing that Bayesian optimization can work well with only few objective function evaluations.

ext: jl


text:
using CSV
using DataFrames

"""
Perform primitive check for smoothness of derivative.
If derivative û is smooth, elements of diff(û) should be small.
Use parameter-sets from MATLAB examples. =#
"""
function _test_demo_large(data, scale, precond, diff_kernel)
    @test sum(
        abs.(
            diff(
                tvdiff(
                    data,
                    40,
                    1e-1;
                    ε=1e-8,
                    scale=scale,
                    precond=precond,
                    diff_kernel=diff_kernel,
                ),
            ),
        ),
    ) < 1e3
end

function _testset_demo_large(scale, preconds, diff_kernels)
    # Load data
    file = CSV.File("./data/demo_large.csv")
    data = DataFrame(file).largescaledata

    @testset "Demo large" begin
        for precond in preconds
            @testset "$precond" begin
                for diff_kernel in diff_kernels
                    @testset "$diff_kernel" begin
                        _test_demo_large(data, scale, precond, diff_kernel)
                    end
                end
            end
        end
    end
end


text:
# diffusion map calculation

function get_adjacency(data::Matrix, kernel::T) where T <: AbstractKernel
    n = size(data)[1]
    A = zeros(n, n)
    Threads.@threads for j in 1:n
        @views for i in j+1:n
            A[i,j] = A[j,i] = similarity(kernel, data[i, :], data[j, :])
        end
    end
    return A
end

function thresholding!(A::Matrix, threshold::Int)
    if threshold < size(A,1)
        n = size(A,1)
        cutOff = rand(size(A,1))
        Threads.@threads for i in 1:n
            cutOff[i] = sort(A[:,i], rev=true)[threshold]
        end
        Threads.@threads for j in 1:n
            @views for i in j:n
                    A[i,j] = A[j,i] = (A[i,j] < cutOff[i]) && (A[i,j] < cutOff[j]) ? 0.0 : A[i,j]
            end
        end
    end
end

function get_laplacian(A::Matrix, method::RowNormalizedLaplacian)
    D = Diagonal(sum(A, dims=2)[:])
    L = D - A
    return inv(D) * L
end

function get_laplacian(A::Matrix, method::NormalizedAdjacency)
    D = Diagonal(sum(A, dims=2)[:])
    return inv(D) * A
end

function get_laplacian(A::Matrix, method::Adjacency)
    return A
end

function get_laplacian(A::Matrix, method::SymmetricLaplacian)
    D = Diagonal(sum(A, dims=2)[:])
    L = D - A
    return D^(-1/2) * L * D^(-1/2)
end

function get_laplacian(A::Matrix, method::RegularLaplacian)
    D = Diagonal(sum(A, dims=2)[:])
    return D - A
end

function eigen_solve(L::Matrix, eigenSolver::FullEigen, sort_for::Symbol)
    λ, ϕ = eigen(L)
    iter = Dict(:SR => 1:size(L,1), :LR => size(L,1):-1:1)
    ϕ_vec = [ϕ[:,i] for i in iter[sort_for]]
    return λ, ϕ_vec
end

function eigen_solve(L::Matrix, eigenSolver::ArpackEigen, sort_for::Symbol)
    λ, ϕ = eigs(sparse(L), nev=eigenSolver.n_first, which=sort_for)
    ϕ_vec = [ϕ[:,i] for i in 1:length(λ)]
    return real.(λ), real.(ϕ_vec)
end

function eigen_solve(L::Matrix, eigenSolver::KrylovEigen, sort_for::Symbol)
    λ, ϕ = eigsolve(sparse(L), eigenSolver.n_first, sort_for)
    return real.(λ), real.(ϕ)
end

# solve the whole Diffusionmap
function solve(dm::Diffusionmap; eigenSolver::T=FullEigen()) where T <: AbstractEigenSolver
    A = get_adjacency(dm.data, dm.kernel)
    thresholding!(A, dm.threshold)
    L = get_laplacian(A, dm.laplace_type)
    sort_for = typeof(dm.laplace_type)==Adjacency || typeof(dm.laplace_type)==NormalizedAdjacency ? :LR : :SR
    λ, ϕ = eigen_solve(L, eigenSolver, sort_for)
    return λ, ϕ
end


text:
"""
    CutDelta

A representation of the change in a plan's cut edges resulting from a flip step.
The `cut_edges_before` set represents cut edges that were cut edges before the flip step,
but are not necessarily cut edges after the flip step; the `cut_edges_after` set represents
edges that were not necessarily cut edges before the flip step but are cut edges after the
flip step. Note that the union of these sets is not guaranteed to be {∅}, as it is possible
to have a cut edge that _stays_ a cut edge after a flip step, albeit with different
assignments. The `Δ` field stores the net change in the number of cut edges as a result
of the flip and is easily computable from the `cut_edges_before` and `cut_edges_after`
fields; it is included merely for convenience.

The `neighbors` field stores the neighbors of the node being flipped that shared a district
assignment with the node before the flip step, but not after the flip step. This is used
for contiguity checks.
"""
struct CutDelta{T<:Int}
    cut_edges_before::BitSet
    cut_edges_after::BitSet
    neighbors::BitSet
    Δ::T
end

"""
    Flip

Represents a single state change in a chain run. Flips are intended to be passed to
constraints to determine whether the proposal is valid. If the proposal is accepted,
the plan is updated and the flip step is stored for the purpose of computing statistics.

In general, chains have fairly tight population constraints in addition to the implicit
contiguity constraint and other constraints (e.g. maximum cut edges). It is computationally
cheaper to first compute _only_ the population delta for a proposed flip and abort the flip
if population constraints are not met, rather than also computing a cut edges delta
(necessary for contiguity checks) that is discarded if the population constraint is not met.
For this reason, the `cut_delta` field can be `Nothing`. The `add_cut_delta` function
should be called to generate an augmented `Flip` once population constraints have been
checked.
"""
abstract type AbstractFlip end

struct Flip <: AbstractFlip
    nodes::Array{Int}
    populations::Array{Int}
    left_district::Int
    right_district::Int
    left_pop::Int
    right_pop::Int
    old_assignments::Array{Int}
    new_assignments::Array{Int}
    cut_delta::Union{CutDelta, Missing}
    balanced_cuts::Int
    # TODO: weight/latency
end

struct DummyFlip <: AbstractFlip
    reason::AbstractString
end

function DummyFlip()::DummyFlip
    return DummyFlip("")
end

"""
    CutDelta(flip, plan, graph)

Return the `CutDelta` determined by a plan and a proposed flip step.
"""
function CutDelta(flip::Flip, plan::Plan, graph::IndexedGraph)::CutDelta
    altered_nodes = Dict{Int, Int}(node => index
                                   for (index, node) in enumerate(flip.nodes))
    cut_edges_before = BitSet(Int[])
    cut_edges_after = BitSet(Int[])
    neighbors = BitSet(Int[])
    for node in union(plan.district_nodes[flip.left_district],
                      plan.district_nodes[flip.right_district])
        @inbounds for index in 1:graph.neighbors_per_node[node]
            neighbor = graph.node_neighbors[index, node]
            edge_index = graph.src_dst_to_edge[node, neighbor]
            if ((plan.assignment[node] == flip.left_district &&
                 plan.assignment[neighbor] == flip.right_district) ||
                (plan.assignment[node] == flip.right_district &&
                 plan.assignment[neighbor] == flip.left_district))
                push!(cut_edges_before, edge_index)
            end

            if node in keys(altered_nodes)
                node_index = altered_nodes[node]
                if neighbor in keys(altered_nodes)
                    neighbor_index = altered_nodes[neighbor]
                    if (flip.new_assignments[node_index] !=
                        flip.new_assignments[neighbor_index])
                        push!(cut_edges_after, edge_index)
                    end
                else
                    if flip.new_assignments[node_index] != plan.assignment[neighbor]
                        push!(cut_edges_after, edge_index)
                    end
                end
            end
        end
    end
    Δ = length(cut_edges_after) - length(cut_edges_before)
    return CutDelta(cut_edges_before, cut_edges_after, neighbors, Δ)
end

"""
    add_cut_delta(flip, graph, plan)

Return an augmented version of a `Flip` with the `cut_delta` field populated.
"""
function add_cut_delta(flip::Flip, graph::IndexedGraph, plan::Plan)::Flip
    cut_delta = CutDelta(flip, plan, graph)
    return Flip(flip.nodes, flip.populations, flip.left_district,
                flip.right_district, flip.left_pop, flip.right_pop,
                flip.old_assignments, flip.new_assignments, cut_delta,
                flip.balanced_cuts)
end


text:
using FactCheck
using ValueDispatch

facts("Value dispatch fns") do

    context("@on with `=`") do

        @dispatch fizzbuzz(n::Int) = (n % 3 == 0, n % 5 == 0)

        @fact_throws fizzbuzz(15)

        @on (true,true)  fizzbuzz(n) = "fizzbuzz"
        @on (true,false) fizzbuzz(n) = "fizz"
        @on (false,true) fizzbuzz(n) = "buzz"
        @on :default     fizzbuzz(n) = n

        @fact fizzbuzz(15) --> "fizzbuzz"
        @fact fizzbuzz(3)  --> "fizz"
        @fact fizzbuzz(5)  --> "buzz"
        @fact fizzbuzz(4)  --> 4

    end

    context("@on with `function`") do

        @dispatch function bangclang(n::Int)
            (n % 3 == 0, n % 5 == 0)
        end

        @fact_throws bangclang(15)

        @on (true,true) function bangclang(n)
            "bangclang"
        end
        @on (true,false) function bangclang(n)
            "bang"
        end
        @on (false,true) function bangclang(n)
            "clang"
        end
        @on :default function bangclang(n)
            n
        end

        @fact bangclang(15) --> "bangclang"
        @fact bangclang(3)  --> "bang"
        @fact bangclang(5)  --> "clang"
        @fact bangclang(4)  --> 4

    end

    context("using `register`") do

        @dispatch bizzbazz(n::Int) = (n % 3 == 0, n % 5 == 0)

        register(bizzbazz, (true,true)) do n
            "bizzbazz"
        end
        register(bizzbazz, (true,false)) do n
            "bizz"
        end
        register(n -> "bazz", bizzbazz, (false,true))
        register(n -> n, bizzbazz, :default)

        @fact bizzbazz(15) --> "bizzbazz"
        @fact bizzbazz(3)  --> "bizz"
        @fact bizzbazz(5)  --> "bazz"
        @fact bizzbazz(4)  --> 4

    end

end


text:
@testset "Plots recipe" begin
    mtg = read_mtg("files/simple_plant.mtg", Dict);

    recipe = RecipesBase.apply_recipe(Dict{Symbol,Any}(), mtg)

    @test length(recipe) == length(mtg)
    # ?NB: both should be equal because 6 lines (all except first point have an edge),
    # and 1 scatter.

    @test recipe[1].plotattributes == Dict{Symbol,Any}(:label => "", :seriescolor => :black, :seriestype => :line)
    @test recipe[1].args == ([0.0, 0.0], [0.0, 0.1]) # coordinates of the two vertex for the edge
    @test recipe[6].args == ([0.0, 0.7071067811865474], [1.3, 2.0071067811865477]) # same, last one

    df_coordinates = MultiScaleTreeGraph.mtg_coordinates_df(mtg, force = true)

    @test recipe[7].args[1] == df_coordinates.XX
    @test recipe[7].plotattributes == Dict{Symbol,Any}(
            :label => "",
            :seriescolor => :viridis,
            :hover => [
                "name: `node_1`, link: `/`, symbol: `\$`, index: `0`",
                "name: `node_2`, link: `/`, symbol: `Individual`, index: `0`",
                "name: `node_3`, link: `/`, symbol: `Axis`, index: `0`",
                "name: `node_4`, link: `/`, symbol: `Internode`, index: `0`",
                "name: `node_5`, link: `+`, symbol: `Leaf`, index: `0`",
                "name: `node_6`, link: `<`, symbol: `Internode`, index: `1`",
                "name: `node_7`, link: `+`, symbol: `Leaf`, index: `0`"
            ],
            :seriestype => :scatter,
            :colorbar_entry => false,
            :marker_z => Any[1, 1, 1, 1, 2, 1, 2]
        )
end


text:
module MarkTrapz
    export mark_trapz!
    function mark_trapz!(robot, num_steps::Integer)

        line = interface_line(robot.move!)
        trajectories = interface_trajectories(robot)

        putmarkers!(side, num_steps) = line.movements!(robot.putmarker!, side, num_steps)

        trajectories.snake!(Ost, Nord) do fold_direct
            (fold_direct==Ost) && robot.putmarker!()
            putmarkers!(fold_direct, num_steps)
            (fold_direct==West) && (num_steps -= 2)
            return true
        end
    end

    function mark_trapz!(robot)
        robot = interface_protected_robot(robot)
        line = interface_line(robot.move!)

        num_steps = line.get_num_movements!(Ost)
        line.movements!(West)
        mark_trapz!(robot, num_steps)

    end

end

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


text:
#MDDatasets DataHR defninitions
#-------------------------------------------------------------------------------

#==Main types
===============================================================================#

#Hyper-rectangle -representation of data:
#-------------------------------------------------------------------------------
mutable struct DataHR{T} <: DataMD
	sweeps::Vector{PSweep}
	elem::Array{T}

	function DataHR{T}(sweeps::Vector{PSweep}, elem::Array{T}) where {T}
		if !elemallowed(DataMD, eltype(elem))
			msg = "Can only create DataHR{T} for T ∈ {DataF1, DataFloat, DataInt, DataComplex}"
			throw(ArgumentError(msg))
		elseif size(DataHR, sweeps) != size(elem)
			throw(ArgumentError("Number of sweeps must match dimensionality of elem"))
		end
		return new(sweeps, elem)
	end
end

#Shorthand (because default (non-parameterized) constructor was overwritten):
DataHR(sweeps::Vector{PSweep}, a::Array{T}) where T = DataHR{T}(sweeps, a)

#Construct DataHR from Vector{PSweep}:
(::Type{DataHR{T}})(sweeps::Vector{PSweep}) where T = DataHR{T}(sweeps, Array{T}(undef, size(DataHR, sweeps)...))

#Construct DataHR{DataF1} from DataHR{Number}
#Collapse inner-most sweep (last dimension), by default:
#TODO: use convert(...) instead?
function (::Type{DataHR{DataF1}})(d::DataHR{T}) where T<:Number
	sweeps = d.sweeps[1:end-1]
	x = d.sweeps[end].v

	#Collapsed last sweep:
	if length(sweeps) < 1
		y = d.elem
		return DataF1(x, reshape(y, length(y)))
	end

	#Reduced DataHR structure:
	result = DataHR{DataF1}(sweeps) #Construct empty results
	_sub = collect(subscripts(result))
	for inds in _sub
		y = d.elem[inds...,:]
		result.elem[inds...] = DataF1(x, reshape(y, length(y)))
	end
	return result
end

#Relay function, so people can blindly convert to DataHR{DataF1} using any DataHR:
(::Type{DataHR{DataF1}})(d::DataHR{DataF1}) = d

#==Type promotions
===============================================================================#
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:DataHR, T2<:Number} = DataHR

#==Accessor functions
===============================================================================#
#Compute the size of a DataHR array from a Vector{PSweep}:
function Base.size(::Type{DataHR}, sweeps::Vector{PSweep})
	dims = Int[]
	for s in sweeps
		push!(dims, length(s.v))
	end
	if 0 == length(dims) #Without sweeps, you can still have a single subset
		push!(dims, 1)
	end
	return tuple(dims...)
end

#Returns the dimension corresponding to the given parameter name:
function dimension(::Type{DataHR}, sweeps::Vector{PSweep}, id::String)
	dim = findfirst((s)->(id==s.id), sweeps)
	ensure(dim!=nothing, ArgumentError("Sweep not found: $id."))
	return dim
end
dimension(d::DataHR, id::String) = dimension(DataHR, d.sweeps, id)

#Returns an element subscripts iterator for a DataHR corresponding to Vector{PSweep}.
subscripts(::Type{DataHR}, sweeps::Vector{PSweep}) =
	subscripts(size(DataHR, sweeps))
subscripts(d::DataHR) = subscripts(d.elem)

#Dimensionality of DataHR array:
Base.ndims(::Type{DataHR}, sweeps::Vector{PSweep}) = max(1, length(sweeps))
Base.ndims(d::DataHR) = ndims(DataHR, d.sweeps)

#Obtain sweep info
#-------------------------------------------------------------------------------
sweeps(d::DataHR) = d.sweeps
sweep(d::DataHR, dim::Int) = d.sweeps[dim].v
sweep(d::DataHR, dim::String) = d.sweeps[dim].v

#Returns parameter sweep coordinates corresponding to given subscript:
function coordinates(d::DataHR, subscr::Tuple=0)
	result = []
	if length(d.sweeps) > 0
		for i in 1:length(subscr)
			push!(result, sweep(d, i)[subscr[i]])
		end
	end
	return result
end

"""
    paramlist(d::Vector{PSweep})

Return a list of parameter in a list of PSweep.
"""
paramlist(sl::Vector{PSweep}) = String[s.id for s in sl]

"""
    paramlist(d::DataHR)

Return a list of parameter values being swept.
"""
paramlist(d::DataHR) = paramlist(d.sweeps)

#==Help with construction
===============================================================================#

#Implement "fill(DataHR, ...) do sweepval" syntax:
function Base.fill!(fn::Function, d::DataHR)
	for inds in subscripts(d)
		d.elem[inds...] = fn(coordinates(d, inds)...)
	end
	return d
end
Base.fill(fn::Function, ::Type{DataHR{T}}, swlist::Vector{PSweep}) where T =
	fill!(fn, DataHR{T}(swlist))
Base.fill(fn::Function, ::Type{DataHR{T}}, sweep::PSweep) where T = fill(fn, DataHR{T}, PSweep[sweep])
#Default: DataHR -> creates DataHR{DataF1}??? Is this a bad idea?
Base.fill(fn::Function, ::Type{DataHR}, swlist::Vector{PSweep}) = fill(fn, DataHR{DataF1}, swlist)

#Fill with a specified value:
Base.fill(elem_0::Number, ::Type{DataHR}, swlist::Vector{PSweep}) =
	DataHR(swlist, fill(elem_0, size(DataHR, swlist)))
Base.fill(elem_0::DataF1, ::Type{DataHR}, swlist::Vector{PSweep}) =
	DataHR(swlist, fill!(Array{DataF1}(undef, size(DataHR, swlist)), elem_0)) #Must be array of abstract DatF1[]
Base.fill(elem_0::DF1_Num, ::Type{DataHR}, sweep::PSweep) = fill(elem_0, DataHR, PSweep[sweep])

#zeros & ones:
Base.zeros(::Type{DataHR}, swlist::Vector{PSweep}) = fill(0.0, DataHR, swlist)
Base.zeros(::Type{DataHR}, sweep::PSweep) = zeros(DataHR, PSweep[sweep])
Base.ones(::Type{DataHR}, swlist::Vector{PSweep}) = fill(1.0, DataHR, swlist)
Base.ones(::Type{DataHR}, sweep::PSweep) = ones(DataHR, PSweep[sweep])

#==Data generation
===============================================================================#
#Generate a DataHR object containing the value of a given swept parameter:
function parameter(::Type{DataHR}, sweeps::Vector{PSweep}, sweepno::Int)
	sw = sweeps[sweepno].v #Sweep of interest
	T = eltype(sw)
	result = DataHR{T}(sweeps)
	for inds in subscripts(result)
		result.elem[inds...] = sw[inds[sweepno]]
	end
	return result
end
parameter(::Type{DataHR}, sweeps::Vector{PSweep}, id::String) =
	parameter(DataHR, sweeps, dimension(DataHR, sweeps, id))
parameter(d::DataHR, id::String) = parameter(DataHR, d.sweeps, id)

#==Dataset reductions
===============================================================================#

#Like getindex(A, inds...), but with DataHR:
#TODO: Support array "view"s instead of copying with "getindex" for efficiency??
function getsubarrayind(d::DataHR{T}, inds...) where T
	sweeps = PSweep[]
	idx = 1
	for rng in inds
		sw = d.sweeps[idx]

		#Only provide a sweep if user selects a range of more than one element:
		addsweep = Colon == typeof(rng) || length(rng)>1
		if addsweep
			push!(sweeps, PSweep(sw.id, sw.v[rng]))
		end
		idx +=1
	end
	return DataHR{T}(sweeps, reshape(getindex(d.elem, inds...), size(DataHR, sweeps)))
end

#getindex(DataHR, inds...), using key/value pairs:
function getsubarraykw(d::DataHR{T}; kwargs...) where T
	sweeps = PSweep[]
	indlist = Vector{Int}[]
	for sweep in d.sweeps
		keepsweep = true
		arg = getkwarg(kwargs, Symbol(sweep.id))
		if arg != nothing
			inds = indices(sweep, arg)
			push!(indlist, inds)
			if length(inds) > 1
				keepsweep = false
				push!(sweeps, PSweep(sweep.id, sweep.v[inds...]))
			end
		else #Keep sweep untouched:
			push!(indlist, 1:length(sweep.v))
			push!(sweeps, sweep)
		end
	end
	return DataHR{T}(sweeps, reshape(getindex(d.elem, indlist...), size(DataHR, sweeps)))
end

#TODO: Support array "view"s instead of copying with "getindex" for efficiency??
#NOTE: Used to be called "sub", which was replaced by "view".
function getsubarray(d::DataHR{T}, args...; kwargs...) where T
	if length(kwargs) > 0
		return getsubarraykw(d, args...; kwargs...)
	else
		return getsubarrayind(d, args...)
	end
end

#==User-friendly show functions
===============================================================================#
#Also changes string():
Base.print(io::IO, ::Type{DataHR{DataF1}}) = print(io, "DataHR{DataF1}")
Base.print(io::IO, ::Type{DataHR{DataFloat}}) = print(io, "DataHR{DataFloat}")
Base.print(io::IO, ::Type{DataHR{DataInt}}) = print(io, "DataHR{DataInt}")
Base.print(io::IO, ::Type{DataHR{DataComplex}}) = print(io, "DataHR{DataComplex}")

function Base.show(io::IO, ds::DataHR)
	szstr = string(size(ds.elem))
	typestr = string(typeof(ds))
	print(io, "$typestr$szstr[\n")
	for inds in subscripts(ds)
		if isassigned(ds.elem, inds...)
			subset = ds.elem[inds...]
			print(io, " $inds: "); show(io, subset); println(io)
		else
			println(io, " $inds: UNDEFINED")
		end
	end
	print(io, "]\n")
end

function Base.show(io::IO, ds::DataHR{T}) where T<:Number
	szstr = string(size(ds.elem))
	typestr = string(typeof(ds))
	print(io, "$typestr$szstr:\n")
	print(io, ds.elem)
end

#Last line


text:
### A Pluto.jl notebook ###
# v0.18.2

using Markdown
using InteractiveUtils

# ╔═╡ 88d07650-7a69-11eb-0a3c-4168940f5e2d
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
			Pkg.PackageSpec(name="PowerModelsDistributionStateEstimation",   version="0.6.0"),
			Pkg.PackageSpec(name="PowerModelsDistribution", version="0.11.4"),
		Pkg.PackageSpec(name="Ipopt", version="0.6.5"),
			Pkg.PackageSpec(name="Plots"), Pkg.PackageSpec(name="CSV")])
	 Pkg.build("Plots")
	using PowerModelsDistributionStateEstimation
	using PowerModelsDistribution
	using Ipopt
	using Plots
	using CSV
	nothing
end

# ╔═╡ 39a2e490-25bf-11eb-2d3a-05fda8c6ecd4
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));

# ╔═╡ 84ed7c50-2270-11eb-035c-33bda27ed57e
md""" ## ⚡ Power System State Estimation ⚡    with _PowerModelsDistributionStateEstimation.jl_
"""

# ╔═╡ 44035512-2271-11eb-2ffb-2f3c2b8ea231
md"This tutorial is made with Pluto notebooks. You can find details on Pluto.jl (and hence on how to run this notebook) on its [github page](https://github.com/fonsp/Pluto.jl) or have a look at the [JuliaCon2020 presentation](https://www.youtube.com/watch?v=IAF8DjrQSSk&ab_channel=TheJuliaProgrammingLanguage). This notebook is available within the PowerModelsDistributionStateEstimation package, in the /example folder.

⚠⚠ Warning! ⚠⚠ This notebook refers to version 0.6.0 of the package. The package is in active development. The idea is that we improve functionalities in the future, without breaking the old ones. But you never know..
"

# ╔═╡ 94c4f342-2277-11eb-2932-89244e8ce3c9
md"________________________________________________________________________________"

# ╔═╡ 9ef21550-2277-11eb-070f-af3e43b62fef
md"### Let's get started: install and call packages"

# ╔═╡ fa6b2e80-226d-11eb-0500-b7b51a64a01b
md"The following packages are needed for your first state estimation script:
1) PowerModelsDistributionStateEstimation.jl to build a state estimation model/problem,
2) PowerModelsDistribution.jl, for power flow calculations and to parse network data,
3) A solver (in this case Ipopt), to solve the state estimation problem,
4) Plots.jl, to visualize some results,
5) CSV.jl, to read CSV files, we need it towards the end of the notebook.
Note, the packages are installed in a temporary environment (see cell below). This might take a while to compile if it is the first time you run this notebook.
"

# ╔═╡ 7626c130-7a91-11eb-1f0a-a14118b621b4
PMDSE = PowerModelsDistributionStateEstimation

# ╔═╡ 1cffc500-258e-11eb-217b-053f1ced96ae
md"_________________________________________________________________________________

### Let's get network and measurement data"

# ╔═╡ c8874710-258c-11eb-0447-49d7bc03eb39
md"Let's choose a simple network to design a state estimator for. There is one stored in PMDSE, in test/data/extra/network, we just need to find the path. PMDSE exports the 'BASE_DIR' constant, which allows you to figure out where your installation of the package is located. You can create a new cell and display the value of this variable. Otherwise, let's just define the network path:"

# ╔═╡ 1863f5be-258e-11eb-3b21-1d852ef2601e
ntw_path = joinpath(PMDSE.BASE_DIR, "test/data/extra/networks/case3_unbalanced.dss")

# ╔═╡ fcb32f32-258d-11eb-2d1b-d5a2629873cb
md"Similarly, you can access the measurement file relative to this network, in test/data/extra/measurements. The name of the file is `case3_meas.csv`. Call the path `msr_path`."

# ╔═╡ c0672120-258e-11eb-142c-b52038715ee9
msr_path = joinpath(PMDSE.BASE_DIR, "test/data/extra/measurements/case3_meas.csv")

# ╔═╡ cdaf7070-258f-11eb-0325-dfc2acfc139b
	md"Our network data is an OpenDSS file. Then, we can use the PowerModelsDistribution parser to obtain the data dictionary in the format required for our calculations. You should parse to a `MATHEMATICAL_MODEL`, because what we do next is adding the measurement data on top of the network data. This can't be done on the `ENGINEERING_MODEL`. Mathematical and Engineering models are standard input network data representations. If you are not familiar with them, you should probably read the relative documentation in PowerModelsDistribution, see [this](https://lanl-ansi.github.io/PowerModelsDistribution.jl/stable/manual/eng-data-model.html), and [this](https://lanl-ansi.github.io/PowerModelsDistribution.jl/stable/manual/math-model.html) and the sections thereafter."

# ╔═╡ fd38b2e2-258d-11eb-21b7-11f4d4d3d8bf
data = PowerModelsDistribution.parse_file(ntw_path; data_model=MATHEMATICAL)

# ╔═╡ 160082f0-2591-11eb-38f3-4b21df157408
md"Now we can add the measurement data to the data dictionary above. The information will be stored in `data[\"meas\"]`. The function we need is `add_measurements!()`. It takes as argument the network data dictionary `data`, the path to the measurement file `msr_path` and the parameter `actual_meas`. If the measurement data consists of actual meter measurements, set the latter as `true`. If they are not measurement data but, e.g, power flow results, set it to `false`. More details later :)."

# ╔═╡ fdb89140-258d-11eb-1299-f553c29a42b3
add_measurements!(data, msr_path, actual_meas = true)

# ╔═╡ 9378bcbe-25af-11eb-14c2-83e297d1ff10
md"Inspecting the data dictionary allows to see how the network looks like. Although a picture is handier:

![case3 unbalanced grid pic](https://raw.githubusercontent.com/Electa-Git/PowerModelsDistributionStateEstimation.jl/master/examples/assets/case3_pic.PNG)

We can see that the grid is unbalanced and has one three-phase generator, modelled like a slackbus.

"

# ╔═╡ 945f9c50-2599-11eb-385d-b16132a4fb00
md"_________________________________________________________________________________
### Let's finalize and solve our SE problem"

# ╔═╡ 1dff2810-2598-11eb-07d6-89579af705b0
md"We are almost ready to run a state estimation, we just need to state what the criteria of the state estimation will be. [SE criteria](https://electa-git.github.io/PowerModelsDistributionStateEstimation.jl/dev/se_criteria/) are `wls`, `wlav`, `rwlav`, etc. Finally, we set a rescaler. Depending on the solver and on the problem case, the low/high weights of some measurements might cause convergence issues. You can set a rescaler that multiplies the weight of all measurements."

# ╔═╡ 1b210e60-2598-11eb-1e18-978f114b214b
data["se_settings"] = Dict{String,Any}("criterion" => "rwlav", "rescaler" => 1)

# ╔═╡ b40a1b30-2598-11eb-1fef-3fc962eeb5e3
md"To solve the SE problem, we need a solver. We previously picked Ipopt, because it's free and can solve nonlinear problems. Let's call it and set some solver parameters:"

# ╔═╡ 0447f220-2599-11eb-043b-d33b03eec478
slv = PMDSE.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-8, "print_level"=>0)

# ╔═╡ 115d6710-2599-11eb-1d41-2326cd36aa9f
md"Finally, the last thing to decide is the power flow formulation that describes our problem. Let's go for a classic 'AC' formulation in polar coordinates. It can be chosen it in two different ways: either calling the generic `solve_mc_se(data, model_type, solver)` and passing the formulation, e.g. `_PMD.ACPUPowerModel` as argument `model_type` or by calling the function that directly involves that formulation, as follows:"

# ╔═╡ 490720c0-2599-11eb-0282-67d4a9d8605f
se_result = solve_acp_mc_se(data, slv)

# ╔═╡ 83950ea0-2599-11eb-302d-2759e4f60de4
md"_________________________________________________________________________________
### Is the estimation accurate?"

# ╔═╡ 71d07410-259a-11eb-16ed-01dc66e83598
md"Let's see what the load magnitudes are, according to the state estimator:"

# ╔═╡ bdf47d00-259a-11eb-2827-3d563dd59e37
se_result["solution"]["load"]

# ╔═╡ fd07f302-259a-11eb-1adf-237ce041edbd
md"Looking at the grid picture, this results seem perhaps even too accurate! Let's quickly have a look at the measurement data:"

# ╔═╡ 415489b0-259b-11eb-37b2-fd455d8aeb90
data["meas"]

# ╔═╡ 58065260-259b-11eb-3a8f-13216b7b7f1b
md"There's virtually no error on the loads! This makes me suspect that the those measurement data are actually not real measurement values, but just the results of a power flow. Let's run a power flow with the real loads and see what the voltage magnitude is, at each bus. We can run a three-phase unbalanced power flow using PowerModelsDistribution's function:"

# ╔═╡ 6ef98650-259a-11eb-3d7f-9dab68152f5d
begin
	pf_results = PowerModelsDistribution.solve_mc_pf(data, PowerModelsDistribution.ACPUPowerModel, slv)
	pf_results["solution"]["bus"]
end

# ╔═╡ d38df080-259d-11eb-14c8-13bd1dbb8d52
md"It does look like they are the same.. An easier way to compare power flow results (or anything that has the same dictionary structure!) to the SE result is by using the `calculate_voltage_magnitude_error` function:"

# ╔═╡ c5588ca2-259d-11eb-0fbc-718d01b1448b
delta, max_err, avg = calculate_voltage_magnitude_error(se_result, pf_results)

# ╔═╡ ba998ac2-28d2-11eb-384a-e748590fb9f7
md"We can even plot it:"

# ╔═╡ c51ff062-28d2-11eb-2732-dd00cba047e3
scatter(delta, xlabel="Index [-]", ylabel="ϵ [p.u.]", ylims = [0, max_err*1.05], title="Absolute voltage difference ϵ", legend = false)

# ╔═╡ a7e4eaf0-259e-11eb-0797-6d7b6db230c3
md"There is basically no error. This is too good to be true. That measurement file must have been created artificially from power flow results! How and why are explained in the next section.

⚠**NB**⚠: if it is the first time you run this and you get errors > e-6, you might want to manually re-run the cells above where ipopt is called and state estimation and power flow are run (you do this by clicking the triangle symbol at the bottom right of the cells). Sometimes something goes wrong in the solution process at the first go."

# ╔═╡ 8fa94df0-259e-11eb-1c60-8f8662091210
md"________________________________________________________________________________
### Creating artificial measurement data"

# ╔═╡ 65dcb0e0-259c-11eb-0df5-9781707c96a9
md"Probably that simple small network doesn't really exist. Or maybe it exists but there are no measurements, but the user wanted to explore state estimation designs regardless, for a future installation. If this is also your case, you can provide realistic power values to your load (or generators if you have any, except for the slackbus), and run a power flow. We already did this above, and got `pf_results`. Now let's use that to build artificial measurements."

# ╔═╡ 8ee68da0-25b3-11eb-02f7-5b3909402ce9
begin
	different_msr_path = pwd()*"//tmp_meas.csv"
	write_measurements!(PowerModelsDistribution.ACPUPowerModel, data, pf_results, different_msr_path)
end

# ╔═╡ dbdb3790-25b4-11eb-36e2-a3b2074002f4
md"Let's have a look at the csv file we just created"

# ╔═╡ c3334430-25b4-11eb-1c17-bb32bd1dfe59
meas_df = CSV.read(different_msr_path)

# ╔═╡ fbbe5d30-25b4-11eb-1319-ddf94502c717
md"The function `write_measurements!` takes the power flow results and produces a CSV file. The generated measurements are the voltage magnitude at each bus with a load or a generator, and active and reactive power of each load and generator. If fewer measurements are wished for, the user can delete the csv/dataframe rows. If different or more data is wished for, a customized function needs to be written, or the measurement data needs to be imported/created somewhere else. In column `par_1` of the csv, the values of a given variable from the power flow results are stored. `par_2` gives information on the uncertainty of the \"fake\" or actual meter which is placed there. Our meters here are fake and their error is assumed to follow a gaussian distribution. `par_2` is the standard deviation of the fake meter, and serves two scopes:

1) the inverse of its square is the weight of that measurement,
2) if the `actual_meas` argument in `add_measurements!` is set to `false`, a fake measurement is generated, sampling from the gaussian distribution which has as sigma `par_2` and mu `par_1`.
Let's see:
"

# ╔═╡ 97d180ce-25b5-11eb-1ae4-4f6681b31964
begin
	add_measurements!(data, msr_path, actual_meas = false)
	data["meas"]["5"]
end

# ╔═╡ 22f33860-25b7-11eb-201b-07199a206651
meas_df[5, :]

# ╔═╡ 406b0350-25b7-11eb-181b-b1af1beb6c5a
md"
____________________________________________________________________________________
### Exercise

It is left as excercise to the reader to see that the SE results with the new measurement data, which contain errors, is less accurate than the case in which there was no error on the measurements.
"

# ╔═╡ e2c20b70-25b8-11eb-1f8a-492df3e740b0
md"
__________________________________________________________________________________
### Is this it?
"

# ╔═╡ 260632d0-25b9-11eb-121c-b32d4595a5b1
md"No, but is the core of the package. All in all, there are a already quite some things mentioned above that you can play with, to design your estimator:
1) Choice of criterion
2) Choice of power flow equations
3) Choice of measurements, that can be \"built\" with the help of the shown functionalities
4) Rescaler

More information on the offered criteria and power flow equations can be found in the documentation.
Let's conclude with a bunch of additional handy functions:
`update_all_bounds!()` allows to add bounds on loads, generators, and buses. To have a subset, the following also exist: `update_load_bounds!()`, `update_generator_bounds!()`, `update_voltage_bounds!()`.

`convert_rectangular_to_polar!()` in postprocessing, to convert results/entries in rectangular coordinates, into polar coordinates.

`assign_start_to_variables!()` to give the solver an initial value.

"

# ╔═╡ 04839890-25ba-11eb-2ef5-13c22aa3e9a0
begin
   update_all_bounds!(data; v_min = 0.85, v_max = 1.15, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
   assign_start_to_variables!(data)
	data["bus"]
end

# ╔═╡ 1b5fe1c0-25bc-11eb-0a5a-d33b22f94dc4
md" Finally, PowerModelsDistributionStateEstimation is designed to be flexible and easy to expand. So, once you are familiar with the package, you can look into the optimization problem itself and add your own constraints, objectives, variables, functions.. 🚀🌕"

# ╔═╡ Cell order:
# ╟─39a2e490-25bf-11eb-2d3a-05fda8c6ecd4
# ╟─84ed7c50-2270-11eb-035c-33bda27ed57e
# ╟─44035512-2271-11eb-2ffb-2f3c2b8ea231
# ╟─94c4f342-2277-11eb-2932-89244e8ce3c9
# ╟─9ef21550-2277-11eb-070f-af3e43b62fef
# ╟─fa6b2e80-226d-11eb-0500-b7b51a64a01b
# ╠═88d07650-7a69-11eb-0a3c-4168940f5e2d
# ╟─7626c130-7a91-11eb-1f0a-a14118b621b4
# ╟─1cffc500-258e-11eb-217b-053f1ced96ae
# ╟─c8874710-258c-11eb-0447-49d7bc03eb39
# ╠═1863f5be-258e-11eb-3b21-1d852ef2601e
# ╟─fcb32f32-258d-11eb-2d1b-d5a2629873cb
# ╠═c0672120-258e-11eb-142c-b52038715ee9
# ╟─cdaf7070-258f-11eb-0325-dfc2acfc139b
# ╠═fd38b2e2-258d-11eb-21b7-11f4d4d3d8bf
# ╟─160082f0-2591-11eb-38f3-4b21df157408
# ╠═fdb89140-258d-11eb-1299-f553c29a42b3
# ╟─9378bcbe-25af-11eb-14c2-83e297d1ff10
# ╟─945f9c50-2599-11eb-385d-b16132a4fb00
# ╟─1dff2810-2598-11eb-07d6-89579af705b0
# ╠═1b210e60-2598-11eb-1e18-978f114b214b
# ╟─b40a1b30-2598-11eb-1fef-3fc962eeb5e3
# ╠═0447f220-2599-11eb-043b-d33b03eec478
# ╟─115d6710-2599-11eb-1d41-2326cd36aa9f
# ╠═490720c0-2599-11eb-0282-67d4a9d8605f
# ╟─83950ea0-2599-11eb-302d-2759e4f60de4
# ╟─71d07410-259a-11eb-16ed-01dc66e83598
# ╠═bdf47d00-259a-11eb-2827-3d563dd59e37
# ╟─fd07f302-259a-11eb-1adf-237ce041edbd
# ╠═415489b0-259b-11eb-37b2-fd455d8aeb90
# ╟─58065260-259b-11eb-3a8f-13216b7b7f1b
# ╠═6ef98650-259a-11eb-3d7f-9dab68152f5d
# ╟─d38df080-259d-11eb-14c8-13bd1dbb8d52
# ╠═c5588ca2-259d-11eb-0fbc-718d01b1448b
# ╟─ba998ac2-28d2-11eb-384a-e748590fb9f7
# ╟─c51ff062-28d2-11eb-2732-dd00cba047e3
# ╟─a7e4eaf0-259e-11eb-0797-6d7b6db230c3
# ╟─8fa94df0-259e-11eb-1c60-8f8662091210
# ╟─65dcb0e0-259c-11eb-0df5-9781707c96a9
# ╠═8ee68da0-25b3-11eb-02f7-5b3909402ce9
# ╟─dbdb3790-25b4-11eb-36e2-a3b2074002f4
# ╠═c3334430-25b4-11eb-1c17-bb32bd1dfe59
# ╟─fbbe5d30-25b4-11eb-1319-ddf94502c717
# ╠═97d180ce-25b5-11eb-1ae4-4f6681b31964
# ╠═22f33860-25b7-11eb-201b-07199a206651
# ╟─406b0350-25b7-11eb-181b-b1af1beb6c5a
# ╟─e2c20b70-25b8-11eb-1f8a-492df3e740b0
# ╟─260632d0-25b9-11eb-121c-b32d4595a5b1
# ╠═04839890-25ba-11eb-2ef5-13c22aa3e9a0
# ╟─1b5fe1c0-25bc-11eb-0a5a-d33b22f94dc4


text:
function gauss(x, w)
    1/sqrt(2.0*pi)./w.*exp(-x.^2/2.0./w.^2)
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
abstract type EvalMode end

struct Training <: EvalMode end

struct Rando <: EvalMode
	percenttest::AbstractFloat
end

struct Cross <: EvalMode
	numfolds::Integer
end

mutable struct Static <: EvalMode
	filename::AbstractString
	data::Matrix
	Static(filename) = new(filename)
end

Base.show(io::IO, ::Training) = print(io, "training")
Base.show(io::IO, ::Rando) = print(io, "random")
Base.show(io::IO, ::Cross) = print(io, "cross")
Base.show(io::IO, ::Static) = print(io, "static")

function EvalMode(args::Array)
	mode = popfirst!(args)
	if mode == "training"
		Training()
	elseif mode == "random"
		Rando(parse(Float64, popfirst!(args)))
	elseif mode == "cross"
		Cross(parse(Int, popfirst!(args)))
	elseif mode == "static"
		Static(popfirst!(args))
	end
end

struct ParseArgs
	normalize::Bool
	arff::AbstractString
	learner::AbstractString
	evaluation::EvalMode
	other::Dict{Symbol,AbstractString}
end

function parseargs(args::Array)
	normalize = false
	arff = ""
	learner = ""
	evaluation = nothing
	other = Dict()
	try
		while !isempty(args)
			option = popfirst!(args)
			if option == "-V"
				ENV["JULIA_DEBUG"] = "MLToolkit"
			elseif option == "-N"
				normalize = true
			elseif option == "-A"
				arff = popfirst!(args)
			elseif option == "-L"
				learner = popfirst!(args)
			elseif option == "-E"
				evaluation = EvalMode(args)
			elseif option == "--seed"
				Random.seed!(parse(popfirst!(args)))
			elseif startswith(option, "--")
				other[Symbol(option[3:end])] = popfirst!(args)
			else
				error("Invalid parameter: $option")
			end
		end
		if arff != "" && learner != "" && evaluation != nothing
			return ParseArgs(normalize, arff, learner, evaluation, other)
		end
	catch
	end
	@error join(["Usage:",
	"MLSystemManager -L learningAlgorithm -A ARFF_File -E evaluationMethod [extraParameters] [OPTIONS]",
	"OPTIONS:",
	"-V Print the confusion matrix and learner accuracy on individual class values",
	"-N Use normalized data",
	"",
	"Possible evaluation methods are:",
	"-E training",
	"-E static testARFF_File",
	"-E random training_percentage",
	"-E cross numOfFolds"], "\n")
	exit(1)
end

max_stars_repo_name: coreywoodfield/mltoolkit-julia

max_issues_repo_name: coreywoodfield/mltoolkit-julia

max_forks_repo_name: coreywoodfield/mltoolkit-julia


text:
@testset "rewrite" begin

    phs = Set([:x])
    @test matchex(:(_x ^ 2), :(a ^ 2)) == Dict(:_x => :a)
    @test matchex(:(x ^ 2), :(a ^ 2), phs=phs) == Dict(:x => :a)
    @test matchex(:(x ^ 2), :(a ^ 2)) == nothing

    phs = Set([:x, :y])
    @test rewrite(:(a + 2b), :(_x + 2_y), :(_x * _y)) == :(a * b)
    @test rewrite(:(a + 2b), :(x + 2y), :(x * y); phs=phs) == :(a * b)

    set_default_placeholders(Set([:x, :y]))
    @test rewrite(:(a + 2b), :(x + 2y), :(x * y)) == :(a * b)
    set_default_placeholders(Set(Symbol[]))

    @test tryrewrite(:(a + 2b), :(x + 2y), :(x * y); phs=[:x, :y]) == :(a * b)
    @test tryrewrite(:(a + 2b), :(x + 3y), :(x * y); phs=[:x, :y]) == nothing

    @test without(:(x * (m == n)), :(_i == _j)) == :x

    @test subs(:(x^n); n=2) == :(x^2)
    @test subs(:(x^n), Dict(:n => 2)) == :(x^2)
    @test subs(:(_mod._op); _mod=:Main, _op=:inc) == :(Main.inc)

    phs = Set([:x, :I])
    @test matchex(:(x[I...]), :(A[i,j,k]); phs=phs) == Dict(:x => :A, :I => [:i, :j, :k])
    @test matchex(:(foo(_args...)), :(foo(x, y))) == Dict(:_args => [:x, :y])

    ex = :(foo(bar(foo(A))))
    rules = [:(foo(x)) => :(quux(x)),
             :(bar(x)) => :(baz(x))]
    @test rewrite_all(ex, rules; phs=[:x]) == :(quux(baz(quux(A))))

    ex = :(foo(bar(foo(A))))
    pat = :(foo(x))
    rpat = :(quux(x))
    @test rewrite_all(ex, pat, rpat; phs=[:x]) == :(quux(bar(quux(A))))

    @test matchingex(:(_x + _y), :(a + a))
    @test !matchingex(:(_x + _y), :(a + a); exact=true)

    ex = :(foo(a, b, c))
    pat = :(foo(xs...))
    rpat = :(bar(xs...))
    @test rewrite(ex, pat, rpat; phs=[:xs]) == :(bar(a, b, c))
end


text:
import PlutoSliderServer
PlutoSliderServer.run_directory("."; static_export=true, config_toml_path="./PlutoDeployment.toml")

max_stars_repo_name: alvarmaciel/paginas-pluto

max_issues_repo_name: alvarmaciel/paginas-pluto

max_forks_repo_name: alvarmaciel/paginas-pluto


text:
#!/usr/bin/env julia

export pe9

pe9(n::Integer=1000)=[a*b*c for a in 1:n, b in 1:n, c in 1:n if a^2+b^2==c^2 && a+b+c==n && a<b<c]|>first

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe9()|>print
end

max_stars_repo_path: 9. Special Pythagorean triplet/9.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 9. Special Pythagorean triplet/9.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 9. Special Pythagorean triplet/9.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
###### update ####################

function update(ta::TimeArray{T,N,D}, tstamp::D, val::AbstractArray{T,N}) where {T,N,D}

    if length(ta) == 0
        throw(ArgumentError(
            "updating empty time arrays is not supported, " *
            "please use a scalable approach"))
    elseif tstamp < maximum(timestamp(ta))
        throw(ArgumentError("only appending operations supported"))
    else
        t    = vcat(timestamp(ta), tstamp)
        vals = vcat(values(ta), val)
        uta  = TimeArray(t, vals, colnames(ta), meta(ta))
    end
    uta
end

function update(ta::TimeArray{T,N,D}, tstamp::D, val::T) where {T,N,D}

    if length(ta) == 0
        throw(ArgumentError(
            "updating empty time arrays is not supported, " *
            "please use a scalable approach"))
    elseif tstamp < maximum(timestamp(ta))
        throw(ArgumentError("only appending operations supported"))
    else
        t    = vcat(timestamp(ta), tstamp)
        vals = vcat(values(ta), val)
        uta  = TimeArray(t, vals, colnames(ta), meta(ta))
    end
    uta
end

###### rename, rename! ####################

"""
    rename(ta::TimeArray, colnames::Vector{Symbol})
    rename(ta::TimeArray, colname::Symbol)
    rename(ta::TimeArray, orig => new, ...)
    rename(f::Base.Callable, ta, colnametyp)

Rename the columns of a `TimeArray`.

See also [`rename!`](@ref).

# Arguments
- `colnametyp` is the input type for the function `f`. The valid value is
  `Symbol` or `String`.
"""
function rename(ta::TimeArray, colnames::Vector{Symbol})
    length(colnames) == size(ta, 2) || throw(ArgumentError("Colnames length mismatch"))
    TimeArray(ta, colnames = colnames, unchecked = true)
end

"""
    rename!(ta::TimeArray, colnames::Vector{Symbol})
    rename!(ta::TimeArray, colname::Symbol)
    rename!(ta::TimeArray, orig => new, ...)
    rename!(f::Base.Callable, ta, colnametyp)

In-place rename the columns of a `TimeArray`.

See also [`rename`](@ref).

# Arguments
- `colnametyp` is the input type for the function `f`. The valid value is
  `Symbol` or `String`.
"""
function rename!(ta::TimeArray, colnames::Vector{Symbol})
    length(colnames) == size(ta, 2) ? (_colnames(ta)[:] = colnames) :
    length(colnames) > 0 && throw(ArgumentError("Colnames length mismatch"))
    ta
end

# shared interfaces
for f ∈ [:rename, :rename!]
    @eval begin
        $f(ta::TimeArray, colnames::Symbol) = $f(ta, [colnames])
        $f(ta::TimeArray) = throw(MethodError($f, (ta,)))
        $f(ta::TimeArray, pairs::Pair{Symbol,Symbol}...) =
            $f(ta, _mapcol(colnames(ta), pairs))
        $f(f::Base.Callable, ta::TimeArray, colnametyp::Type{Symbol} = Symbol) =
            $f(ta, map(f, colnames(ta)))
        $f(f::Base.Callable, ta::TimeArray, colnametyp::Type{String}) =
            $f(Symbol ∘ f ∘ string, ta)
    end
end

"""
utility function for `rename` and `rename!`
"""
function _mapcol(col, pairs)
    col′ = copy(col)
    for (key, val) ∈ pairs
        i = findfirst(isequal(key), col)
        if i ≡ nothing  # FIXME: `isnothing(i)` in Julia 1.1
            throw(ArgumentError("Unknown column `$key`"))
        end
        col′[i] = val
    end
    col′
end


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function mwg_main(par::ParSsm, h::Int64, nDraws::Array{Int64, 1}, burnin::Array{Int64, 1},
                  mwg_const::Array{Float64, 1}, par_ind::BoolParSsm, t=0::Int64, end_oos=0::Int64)

# ----------------------------------------------------------------------------------------------------------------------
# Metropolis-Within-Gibbs algorithm: mainframe
# ----------------------------------------------------------------------------------------------------------------------

     # -----------------------------------------------------------------------------------------------------------------
     # Initialise
     # -----------------------------------------------------------------------------------------------------------------

     # Dimensions
     k    = size(par.T)[1];
     n, m = size(par.y);

     # Target acceptance rate
     acc_target = 30;

     # No. parameters to be estimated
     par_size = SizeParSsm(sum(par_ind.d),
                                 sum(sum(par_ind.Z)),
                                 sum(sum(par_ind.R)),
                                 sum(par_ind.c),
                                 sum(sum(par_ind.T)),
                                 sum(sum(par_ind.Q)),
                                 sum(sum(par_ind.λ)),
                                 sum(sum(par_ind.ρ)),
                                 sum(par_ind.d) + sum(sum(par_ind.Z)) + sum(sum(par_ind.R)) +
                                    sum(par_ind.c) + sum(sum(par_ind.T)) + sum(sum(par_ind.Q)) +
                                    sum(sum(par_ind.λ)) + sum(sum(par_ind.ρ)));

     # -----------------------------------------------------------------------------------------------------------------
     # Lower and upper bounds
     # -----------------------------------------------------------------------------------------------------------------

     # Numerical constants
     xi        = 1e-3; # it must be small
     MIN_var   = 0;
     MIN_coeff = -Inf;
     MIN_λ     = xi;
     MIN_ρ     = xi;
     MAX_var   = Inf;
     MAX_coeff = Inf;
     MAX_λ     = pi;
     MAX_ρ     = 0.97;

     # MIN
     MIN = [MIN_var*ones(par_size.R);
                  MIN_coeff*ones(par_size.d + par_size.Z);
                  MIN_var*ones(par_size.Q);
                  MIN_coeff*ones(par_size.c + par_size.T);
                  MIN_λ*ones(par_size.λ);
                  MIN_ρ*ones(par_size.ρ)];

     # MAX
     MAX = [MAX_var*ones(par_size.R);
                  MAX_coeff*ones(par_size.d + par_size.Z);
                  MAX_var*ones(par_size.Q);
                  MAX_coeff*ones(par_size.c + par_size.T);
                  MAX_λ*ones(par_size.λ);
                  MAX_ρ*ones(par_size.ρ)];

     # -----------------------------------------------------------------------------------------------------------------
     # Set prior parameters
     # -----------------------------------------------------------------------------------------------------------------

     # Define prior distribution objects
     prior_opt = PriorOpt(Normal(0, 1/xi),
                                InverseGamma(3, 1), #TruncatedNormal(0, 1, 0, Inf),
                                par_size.λ*logpdf.(Uniform(MIN_λ, MAX_λ), MIN_λ),
                                par_size.ρ*logpdf.(Uniform(MIN_ρ, MAX_ρ), MIN_ρ));

     # Transformations: 1 natural logarithm, 2 no transformations, 3 generalized logit
     opt_transf = convert(Array{Int64, 1}, [1*ones(par_size.R);
                                                  2*ones(par_size.d + par_size.Z);
                                                  1*ones(par_size.Q);
                                                  2*ones(par_size.c + par_size.T);
                                                  3*ones(par_size.λ);
                                                  3*ones(par_size.ρ)]);

     # -----------------------------------------------------------------------------------------------------------------
     # Initialisation of the Metropolis-Within-Gibbs algorithm
     # -----------------------------------------------------------------------------------------------------------------

     # Initial guess for θ_ini_bound
     θ_ini_bound = [ones(par_size.R);
                          zeros(par_size.d);
                          ones(par_size.Z);
                          ones(par_size.Q);
                          zeros(par_size.c);
                          ones(par_size.T);
                          (2*pi/32)*ones(par_size.λ);
                          0.5*ones(par_size.ρ)];

     # Run init
     θ_ini_unb = get_par_unb(θ_ini_bound, MIN, MAX, opt_transf);

     algorithm_name = "Initialisation";
     acc_rate       = 0.0;

     chain_θ_unb_init = Array{Any}(undef, 1);
     Σ_init           = Array{Float64}(Matrix(I, length(θ_ini_unb), length(θ_ini_unb)));

     while acc_rate < 25.0 || acc_rate > 35.0

          chain_θ_unb_init, _, _, _, _, _ =
               mwg_run(θ_ini_unb, par, h, par_ind, par_size, prior_opt, MIN, MAX, opt_transf, nDraws[1], burnin[1],
                    mwg_const[1], algorithm_name, t, end_oos, Σ_init, 0);

          acc_rate = get_progress(chain_θ_unb_init, burnin[1]+1, nDraws[1], 0);

          if acc_rate < 25.0 || acc_rate > 35.0
               delta_acc    = acc_rate-acc_target;
               mwg_const[1] = mwg_const[1]*exp(delta_acc/100);
          end
     end

     θ_unb_init_clean = chain_θ_unb_init[:, burnin[1]+1:end]';
     θ_start_unb      = median(θ_unb_init_clean, dims=1)[:];
     Σ_mwg            = cov(θ_unb_init_clean, dims=1);

     # -----------------------------------------------------------------------------------------------------------------
     # Metropolis-Within-Gibbs algorithm
     # -----------------------------------------------------------------------------------------------------------------

     algorithm_name = "Metropolis-Within-Gibbs";
     acc_rate       = 0.0;

     chain_θ_unb   = Array{Any}(undef, 1);
     chain_θ_bound = Array{Any}(undef, 1);
     distr_α       = Array{Any}(undef, 1);
     distr_fcst    = Array{Any}(undef, 1);
     distr_par     = Array{Any}(undef, 1);

     while acc_rate < 25.0 || acc_rate > 35.0

          chain_θ_unb, chain_θ_bound, distr_α, distr_fcst, par, distr_par =
               mwg_run(θ_start_unb, par, h, par_ind, par_size, prior_opt, MIN, MAX, opt_transf, nDraws[2], burnin[2],
                    mwg_const[2], algorithm_name, t, end_oos, Σ_mwg, 1);

          acc_rate = get_progress(chain_θ_unb, burnin[2]+1, nDraws[2], 0);

          if acc_rate < 25.0 || acc_rate > 35.0
               delta_acc    = acc_rate-acc_target;
               mwg_const[2] = mwg_const[2]*exp(delta_acc/100);
          end
     end

     return distr_α, distr_fcst, chain_θ_unb, chain_θ_bound, mwg_const, acc_rate, par, par_size, distr_par;
end


text:
using SDELab2 # SDELab2.jl should be findable via LOAD_PATH
println("run_vdp2.jl")
type vdp_param
  alpha
  beta
  A
  B
end

function drift_f(u,p)
  return [u[2], p.alpha*u[1]+p.beta*u[2]-p.A*u[1]^3-p.B*u[1]^2*u[2]]
end
#
function diff_g_mat(u,dw,alpha)
  return [0;alpha*u[1]*dw]
end
#
function diff_g_vecs(alpha)
  m=1
  A=Array(Function,m)
  g1(u)=[0,u[1]*alpha]
  A[1]=g1
  return A
end
# dimensions
d=2
m=1
noise_type=3 # commutative
# params
p=vdp_param(-1,0.1,1.,1.)
sig=0.1
# initial data
u=[0.,0.001]
t0=0
tf=500
#
fcn=set_fcn(d,m,x->drift_f(x,p),(x,dw)->diff_g_mat(x,dw,sig),diff_g_vecs(sig), 3)
opt=set_opt(0.01,"SIBDF")
opt["do_plotting"]=true
opt["output_plot_type"]="phase_plot"
opt["output_sel"]=[1,2]
#using PyPlot
println("Start sde_strong_soln")
tic(); t,y,W=sde_strong_solution(fcn, linspace(t0,tf,10000), u, opt); run_time=toc()

#


text:
module BitsFields

export BitField, BitFields, ByRef, set!, refvalue

abstract type AbstractBitField  end
abstract type AbstractBitFields end

"""
    UBits

Unsigned types available to hold multiple bitfields.
"""
const UBits = Union{UInt128, UInt64, UInt32, UInt16, UInt8}

"""
    IBits

Signed types available to hold bitfield values.
"""
const IBits = Union{Int128, Int64, Int32, Int16, Int8}

"""
    bitsof

As sizeof is with bytes, bitsof is with bits.
"""
bitsof(::Type{T}) where {T} = sizeof(T) * 8

include("byref.jl")
include("bitfield.jl")
include("bitfields.jl")

Base.zero(bitfield::BitField{U}) where {U<:UBits} = zero(U)
Base.zero(bitfields::BitFields{N,U}) where {N, U<:UBits} = Ref(zero(eltype(bitfields[1])))
Base.zero(x::NamedTuple{Names, NTuple{N,BitField}}) where {Names, N} = Ref(zero(eltype(x[1])))

end # BitsFields


text:
## test subLocalGRaphCopy! function of cloudgraphs

# 0. Load required packages
using Caesar

# 1. Get the authentication and session information
include(joinpath(ENV["HOME"],"Documents","blandauthlocal.jl"))
cg, addrdict = standardcloudgraphsetup(addrdict=addrdict)
addrdict["sessionId"] = "Hackathon"
addrdict["robotId"] = "PixieBot"

# 2. Create a local empty factor graph object
fg = initfg(sessionname=addrdict["sessionId"], robotname=addrdict["robotId"], cloudgraph=cg)
# fg.stateless = true # so that we don't update local elements

# 3.a Fetch a portion of the graph
Caesar.subLocalGraphCopy!(fg, ["x1";"x3";], neighbors=1, reqbackendset=false, reqSolvable=false)

# see what is going on
Graphs.plot(fg.g)

# ===========ALTERNATIVE========================================================

fg = initfg(sessionname=addrdict["sessionId"], robotname=addrdict["robotId"], cloudgraph=cg)
# 3.b Or fetch the entire factor graph
Caesar.fullLocalGraphCopy!(fg, reqbackendset=false, reqSolvable=false)

#


text:
"""
	voronoiarea(x::Vector, y::Vector, rw) -> Vector

Compute the area of each Voronoi cell in the tesselation generated by `x` and `y`.

The vector `rw` specifies the boundary rectangle as `[xmin, xmax, ymin, ymax]`.
By default, `rw` is the unit rectangle.
"""
function voronoiarea(x::AbstractVector{Float64}, y::AbstractVector{Float64}, rw::Vector{Float64}=[0.0;1.0;0.0;1.0])
	pts, SCALEX, SCALEY = fit2boundingbox(x, y, rw)

	# Areas for scaled points
	C = voronoicells(pts)
	A = voronoiarea(C)

	# scale!(A, SCALEX*SCALEY)
	A .*= SCALEX*SCALEY

	return A
end

"""
	voronoiarea(C::Tessellation) -> Vector

Compute the area of each of the Voronoi cells in `C`.

Note that if the polygons of `C` are not ordered, they will be changed in-place.
"""
function voronoiarea(C::Tessellation)
	NC = length(C)
	A = Array{Float64}(undef, NC)

	for n in 1:NC
		A[n] = polyarea(C[n])
	end

	return A
end

"""
	polyarea(p::AbstractPoints2D)

Compute the area of the polygon with vertices `p` using the shoelace formula.
If the points in `p` are not sorted, they will be sorted **in-place**.
"""
function polyarea(pts::Vector{<:VoronoiDelaunay.AbstractPoint2D})
	issorted(pts) || sort!(pts)

	Np = length(pts)
	A = getx(pts[1])*( gety(pts[2]) - gety(pts[Np]) ) + getx(pts[Np])*( gety(pts[1]) - gety(pts[Np-1]) )

	for n in 2:Np-1
		A += getx(pts[n])*( gety(pts[n+1]) - gety(pts[n-1]) )
	end

	return 0.5*abs(A)
end

# Compute the average point of pts
function mean(pts::Vector{<:VoronoiDelaunay.AbstractPoint2D})
	# Average point
	ax = 0.0
	ay = 0.0

	for p in pts
		ax += getx(p)
		ay += gety(p)
	end

	Np = length(pts)
	VoronoiDelaunay.Point2D(ax/Np, ay/Np)
end

# Addition and subtraction for AbstractPoint2D
for op in [:+, :-]
	@eval begin
		Base.$op(p::VoronoiDelaunay.AbstractPoint2D, q::VoronoiDelaunay.AbstractPoint2D) = VoronoiDelaunay.Point2D( $op(getx(p), getx(q)), $op(gety(p), gety(q)) )
	end
end

Base.:*(a::Float64, p::VoronoiDelaunay.AbstractPoint2D) = VoronoiDelaunay.Point2D( a*getx(p), a*gety(p) )

# sorting for AbstractPoints2D
for name in [:sort!,:issorted]
	@eval begin
		function Base.$name(pts::Vector{T}) where T<:VoronoiDelaunay.AbstractPoint2D
			center = mean(pts)
			centralize = p -> p - center
			$name( pts, by=centralize )
		end
	end
end

# http://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order
function Base.isless(p::VoronoiDelaunay.AbstractPoint2D, q::VoronoiDelaunay.AbstractPoint2D)
	if getx(p) >= 0.0 && getx(q) < 0.0
		return true
	elseif getx(p) < 0.0 && getx(q) >= 0.0
		return false
	elseif getx(p) == getx(q) == 0.0
		if gety(p) >= 0.0 || gety(q) >= 0.0
			return gety(p) > gety(q)
		else
			return gety(p) < gety(q)
		end
	end

	det = getx(p)*gety(q) - getx(q)*gety(p)
	if det < 0.0
		return true
	elseif det > 0.0
		return false
	end

	# p and q are on the same line from the center; check which one is
	# closer to the origin
	origin = VoronoiDelaunay.Point2D(0.0, 0.0)
	dist_squared(p,origin) > dist_squared(q,origin)
end


text:
function build_ode(name,ics,pars,vars,tdomain)
  dsargs = ds[:args]()
  set_name(dsargs,name)
  set_ics(dsargs,ics)
  set_pars(dsargs,pars)
  set_vars(dsargs,vars)
  set_tdomain(dsargs,tdomain)
  dsargs
end

function build_ode(f::DiffEqBase.AbstractParameterizedFunction,u0,tspan,p)
  name = string(typeof(f))
  pars = Dict{String,Any}(); vars = Dict{String,Any}(); ics = Dict{String,Any}()
  for i in 1:length(f.params)
    pars[string(f.params[i])] = p[i]
  end
  for i in 1:length(f.syms)
    vars[string(f.syms[i])] = string(f.funcs[i])
    ics[string(f.syms[i])] = u0[i]
  end
  build_ode(name,ics,pars,vars,tspan)
end

function solve_ode(dsargs,alg=:Vode_ODEsystem,name="Default Name")
  DS = ds[:Generator][alg](dsargs)
  traj = DS[:compute](name)
  d = interpert_traj(traj)

  #=
  # Interpolations
  t = 5.4
  traj(t)[:coordarray]
  =#
end

function interpert_traj(traj)
  d = Dict{Symbol,Vector{Float64}}()
  d[Symbol(traj[:indepvarname])] = first(first(values(traj[:underlyingMesh]())))
  depvars = PyDict(traj[:sample]())
  for k in keys(depvars)
    d[Symbol(k)] = depvars[k]
  end
  d
end