text:
"""
    ground_state_mincut(g::AGraph, h::Vector, J) -> σ, E

<<<<<<< HEAD:src/Ising/mincut.jl
Exactly compute the ground state of an Ising model with ferromagneting (i.e. nonnegative)
couplings using a minimum cut algorithm implemented in FatGraphs.jl package.
=======
Computes the ground state of a Random Field Ising Model using
a minimum cut algorithm implemented in Erdos.jl package.
>>>>>>> 3b8f1a4... FatGraphs -> Erdos:src/RFIM/mincut.jl

Coupling `J` can be in the form of a costant, a matrix, or a
vector of vectors (adjacency list style).

Returns a vector `σ` taking values ±1 and the ground state energy `E`.
"""
function ground_state_mincut(g::AGraph, h::Vector, J)
    N = nv(g)
    source = N + 1
    target = N + 2
    dg, c = net_capacity(g, h, J)
    f, F, labels = maximum_flow(dg, source, target, c)
    σ = 3 - labels[1:N]*2
    E = energy(g, σ, h, J)
    return σ, E
end

function net_capacity(g::AGraph, h::Vector{T}, J) where T
    N = nv(g)
    dg = digraph(g)
    add_vertices!(dg, 2)
    source = N+1
    target = N+2
    for i=1:N
        if h[i] > 0
            add_edge!(dg, source, i)
            add_edge!(dg, i, source)
        else
            add_edge!(dg, i, target)
            add_edge!(dg, target, i)
        end
    end

    C = promote_type(Jtype(J), T)
    c = Vector{Vector{C}}()
    for i=1:N
        neigs = neighbors(dg, i)
        push!(c, zeros(length(neigs)))
        for (k,j) in enumerate(neigs)
            if j <= N
                c[i][k] = getJ(J, i, j, k)
            elseif j == target
                c[i][k] = abs(h[i])
            end

        end
    end
    push!(c, zeros(degree(dg, source)))
    push!(c, zeros(degree(dg, target)))

    for (k,j) in enumerate(neighbors(dg, source))
        c[source][k] = abs(h[j])
    end

    return dg, c
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
function parse_allocation( filename :: String ) :: Dict
  n = 0
  key = "x"
  out=Dict{String,Tuple}()
  oldv = -1
  newv = -1
  for line in eachline( filename )
    n += 1
    r = (n % 3)
    if r==1
      key = line
      if key == "#"
        break;
      end
    elseif r == 2
      oldv = parse(Int64, line )
    else
      newv = parse(Int64, line )
      out[key]=(oldv,newv)
    end;
    print("$key\n");
  end # loop
  out
end # parse_allocation

ext: jl


text:
@testset "vectorizing parameters" begin
    A = Matrix(reshape(1:9, 3, :))
    @test vec_parameters(UpperTriangular(A)) == [1, 4, 5, 7, 8, 9]
    @test vec_parameters(LowerTriangular(A)) == [1, 2, 3, 5, 6, 9]
    @test vec_parameters(Symmetric(A)) == [1, 4, 5, 7, 8, 9] # may change
    @test vec_parameters(A) == Vector(1:9)
    @test vec_parameters(1.0) == [1.0]
    @test vec_parameters((1.0, [2, 3])) == [1.0, 2.0, 3.0]
end

@testset "fallback methods" begin
    msg(fname) = "You need to define `IndirectLikelihood.$(fname)` with this model type."

    output = @capture_err begin
        @test_throws MethodError MLE("something", "nonsensical")
    end
    @test contains(output, msg("MLE"))

    output = @capture_err begin
        @test_throws MethodError common_random(RNG, "undefined")
    end
    @test contains(output, msg("common_random"))

    output = @capture_err begin
        @test_throws MethodError loglikelihood("something", "nonsensical", 2)
    end
    @test contains(output, msg("loglikelihood"))

    output = @capture_err begin
        @test_throws MethodError simulate_data(RNG, "something", "nonsensical", 2)
    end
    @test contains(output, msg("simulate_data"))
end


text:
using DelayNetTools
using Test

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

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
#
# Calculates requests per second from client to server without response, but
# using only naked sockets.
# Results: rps: 20785.69948410947, time: 48.109999895095825
#
server = listen(2000)
sock   = accept(server)

function run()
  for i=1:1000000
    write(sock, Int(1))
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
module BinaryOps
using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedBinaryOperator, AbstractBinaryOp,
    valid_vec, juliaop, toGBType, GBType
import ..SuiteSparseGraphBLAS: juliaop
using ..libgb
export BinaryOp
function BinaryOp(name)
    if isGxB(name) || isGrB(name) #If it's a built-in drop the prefix
        simplifiedname = name[5:end]
    else
        simplifiedname = name
    end
    containername = Symbol(simplifiedname, "_T")
    exportedname = Symbol(simplifiedname)
    structquote = quote
        struct $containername <: AbstractBinaryOp
            typedops::Dict{Tuple{DataType, DataType}, TypedBinaryOperator}
            name::String
            $containername() = new(Dict{Tuple{DataType, DataType}, TypedBinaryOperator}(), $name)
        end
    end
    @eval($structquote) #eval container *type* into Types submodule
    constquote = quote
        const $exportedname = $containername()
        export $exportedname
    end
    @eval($constquote) #eval actual op into BinaryOps submodule
    return getproperty(BinaryOps, exportedname)
end
struct GenericBinaryOp <: AbstractBinaryOp
    typedops::Dict{Tuple{DataType, DataType}, TypedBinaryOperator}
    name::String
    GenericBinaryOp(name) = new(Dict{Tuple{DataType, DataType}, TypedBinaryOperator}(), name)
end

#This is adapted from the fork by cvdlab.
#Add a new GrB_BinaryOp to an AbstractBinaryOp
function _addbinaryop(
    op::AbstractBinaryOp,
    fn::Function,
    ztype::GBType{T},
    xtype::GBType{U},
    ytype::GBType{V}
) where {T,U,V}
    function binaryopfn(z, x, y)
        unsafe_store!(z, fn(x, y))
        return nothing
    end
    opref = Ref{libgb.GrB_BinaryOp}()
    binaryopfn_C = @cfunction($binaryopfn, Cvoid, (Ptr{T}, Ref{U}, Ref{V}))
    libgb.GB_BinaryOp_new(opref, binaryopfn_C, ztype, xtype, ytype, op.name)
    op.typedops[(U, V)] = TypedBinaryOperator{xtype, ytype, ztype}(opref[])
    return nothing
end

function BinaryOp(fn::Function; keep=true)
    @warn "Use built-in functions where possible, user defined functions are less performant.
        \nSee the documentation for a list of available built-in functions."
    name = string(fn)
    if keep
        op = BinaryOp(name)
        funcquote = quote
            BinaryOp(::typeof($fn)) = $op
            juliaop(::typeof($op)) = $fn
        end
        @eval($funcquote)
    else
        op = GenericBinaryOp(name)
    end
    return op
end

function BinaryOp(fn::Function, ztype, xtype, ytype; keep=true)
    op = BinaryOp(fn; keep)
    _addbinaryop(op, fn, toGBType(ztype), toGBType(xtype), toGBType(ytype))
    return op
end

#xtype == ytype == ztype
function BinaryOp(fn::Function, type::DataType; keep=true)
    return BinaryOp(fn, type, type, type; keep)
end

#Vectors of _type, add one function for each triple.
function BinaryOp(
    fn::Function,
    ztype::Vector{DataType},
    xtype::Vector{DataType},
    ytype::Vector{DataType};
    keep = true
)
    op = BinaryOp(fn; keep)
    length(ztype) == length(xtype) == length(ytype) ||
        throw(DimensionMismatch("Lengths of ztype, xtype, and ytype must match"))
    for i ∈ 1:length(ztype)
        _addbinaryop(op, fn, toGBType(ztype[i]), toGBType(xtype[i]), toGBType(ytype[i]))
    end
    return op
end

#Vector of type, xtype == ytype == ztype
function BinaryOp(fn::Function, type::Vector{DataType}; keep = true)
    return BinaryOp(fn, type, type, type; keep)
end
end
const BinaryUnion = Union{AbstractBinaryOp, TypedBinaryOperator}

function _createbinaryops()
    builtins = [
    "GrB_FIRST",
    "GrB_SECOND",
    "GxB_POW",
    "GrB_PLUS",
    "GrB_MINUS",
    "GrB_TIMES",
    "GrB_DIV",
    "GxB_RMINUS",
    "GxB_RDIV",
    "GxB_PAIR",
    "GxB_ANY",
    "GxB_ISEQ",
    "GxB_ISNE",
    "GxB_ISGT",
    "GxB_ISLT",
    "GxB_ISGE",
    "GxB_ISLE",
    "GrB_MIN",
    "GrB_MAX",
    "GxB_LOR",
    "GxB_LAND",
    "GxB_LXOR",
    "GxB_ATAN2",
    "GxB_HYPOT",
    "GxB_FMOD",
    "GxB_REMAINDER",
    "GxB_LDEXP",
    "GxB_COPYSIGN",
    "GrB_BOR",
    "GrB_BAND",
    "GrB_BXOR",
    "GrB_BXNOR",
    "GxB_BGET",
    "GxB_BSET",
    "GxB_BCLR",
    "GxB_BSHIFT",
    "GrB_EQ",
    "GrB_NE",
    "GrB_GT",
    "GrB_LT",
    "GrB_GE",
    "GrB_LE",
    "GxB_CMPLX",
    "GxB_FIRSTI",
    "GxB_FIRSTI1",
    "GxB_FIRSTJ",
    "GxB_FIRSTJ1",
    "GxB_SECONDI",
    "GxB_SECONDI1",
    "GxB_SECONDJ",
    "GxB_SECONDJ1",
    ]
    for name ∈ builtins
        BinaryOps.BinaryOp(name)
    end
end

function _load(binary::AbstractBinaryOp)
    booleans = [
        "GrB_FIRST",
        "GrB_SECOND",
        "GxB_POW",
        "GrB_PLUS",
        "GrB_MINUS",
        "GrB_TIMES",
        "GrB_DIV",
        "GxB_RMINUS",
        "GxB_RDIV",
        "GxB_PAIR",
        "GxB_ANY",
        "GxB_ISEQ",
        "GxB_ISNE",
        "GxB_ISGT",
        "GxB_ISLT",
        "GxB_ISGE",
        "GxB_ISLE",
        "GrB_MIN",
        "GrB_MAX",
        "GxB_LOR",
        "GxB_LAND",
        "GxB_LXOR",
        "GrB_EQ",
        "GrB_NE",
        "GrB_GT",
        "GrB_LT",
        "GrB_GE",
        "GrB_LE",
    ]
    integers = [
        "GrB_FIRST",
        "GrB_SECOND",
        "GxB_POW",
        "GrB_PLUS",
        "GrB_MINUS",
        "GrB_TIMES",
        "GrB_DIV",
        "GxB_RMINUS",
        "GxB_RDIV",
        "GxB_PAIR",
        "GxB_ANY",
        "GxB_ISEQ",
        "GxB_ISNE",
        "GxB_ISGT",
        "GxB_ISLT",
        "GxB_ISGE",
        "GxB_ISLE",
        "GrB_MIN",
        "GrB_MAX",
        "GxB_LOR",
        "GxB_LAND",
        "GxB_LXOR",
        "GrB_BOR",
        "GrB_BAND",
        "GrB_BXOR",
        "GrB_BXNOR",
        "GxB_BGET",
        "GxB_BSET",
        "GxB_BCLR",
        "GxB_BSHIFT",
        "GrB_EQ",
        "GrB_NE",
        "GrB_GT",
        "GrB_LT",
        "GrB_GE",
        "GrB_LE",
    ]
    unsignedintegers = [
        "GrB_FIRST",
        "GrB_SECOND",
        "GxB_POW",
        "GrB_PLUS",
        "GrB_MINUS",
        "GrB_TIMES",
        "GrB_DIV",
        "GxB_RMINUS",
        "GxB_RDIV",
        "GxB_PAIR",
        "GxB_ANY",
        "GxB_ISEQ",
        "GxB_ISNE",
        "GxB_ISGT",
        "GxB_ISLT",
        "GxB_ISGE",
        "GxB_ISLE",
        "GrB_MIN",
        "GrB_MAX",
        "GxB_LOR",
        "GxB_LAND",
        "GxB_LXOR",
        "GrB_BOR",
        "GrB_BAND",
        "GrB_BXOR",
        "GrB_BXNOR",
        "GxB_BGET",
        "GxB_BSET",
        "GxB_BCLR",
        "GxB_BSHIFT",
        "GrB_EQ",
        "GrB_NE",
        "GrB_GT",
        "GrB_LT",
        "GrB_GE",
        "GrB_LE",
    ]
    floats = [
        "GrB_FIRST",
        "GrB_SECOND",
        "GxB_POW",
        "GrB_PLUS",
        "GrB_MINUS",
        "GrB_TIMES",
        "GrB_DIV",
        "GxB_RMINUS",
        "GxB_RDIV",
        "GxB_PAIR",
        "GxB_ANY",
        "GxB_ISEQ",
        "GxB_ISNE",
        "GxB_ISGT",
        "GxB_ISLT",
        "GxB_ISGE",
        "GxB_ISLE",
        "GrB_MIN",
        "GrB_MAX",
        "GxB_LOR",
        "GxB_LAND",
        "GxB_LXOR",
        "GxB_ATAN2",
        "GxB_HYPOT",
        "GxB_FMOD",
        "GxB_REMAINDER",
        "GxB_LDEXP",
        "GxB_COPYSIGN",
        "GrB_EQ",
        "GrB_NE",
        "GrB_GT",
        "GrB_LT",
        "GrB_GE",
        "GrB_LE",
        "GxB_CMPLX",
    ]

    positionals = [
        "GxB_FIRSTI",
        "GxB_FIRSTI1",
        "GxB_FIRSTJ",
        "GxB_FIRSTJ1",
        "GxB_SECONDI",
        "GxB_SECONDI1",
        "GxB_SECONDJ",
        "GxB_SECONDJ1",
    ]

    complexes = [
        "GxB_FIRST",
        "GxB_SECOND",
        "GxB_POW",
        "GxB_PLUS",
        "GxB_MINUS",
        "GxB_TIMES",
        "GxB_DIV",
        "GxB_RMINUS",
        "GxB_RDIV",
        "GxB_PAIR",
        "GxB_ANY",
        "GxB_ISEQ",
        "GxB_ISNE",
        "GxB_EQ",
        "GxB_NE",
    ]
    name = binary.name
    if name ∈ booleans
        binary[Bool] = TypedBinaryOperator(load_global(name * "_BOOL", libgb.GrB_BinaryOp))
    end

    if name ∈ integers
        binary[Int8] = TypedBinaryOperator(load_global(name * "_INT8", libgb.GrB_BinaryOp))
        binary[Int16] = TypedBinaryOperator(load_global(name * "_INT16", libgb.GrB_BinaryOp))
        binary[Int32] = TypedBinaryOperator(load_global(name * "_INT32", libgb.GrB_BinaryOp))
        binary[Int64] = TypedBinaryOperator(load_global(name * "_INT64", libgb.GrB_BinaryOp))
    end

    if name ∈ unsignedintegers
        binary[UInt8] = TypedBinaryOperator(load_global(name * "_UINT8", libgb.GrB_BinaryOp))
        binary[UInt16] = TypedBinaryOperator(load_global(name * "_UINT16", libgb.GrB_BinaryOp))
        binary[UInt32] = TypedBinaryOperator(load_global(name * "_UINT32", libgb.GrB_BinaryOp))
        binary[UInt64] = TypedBinaryOperator(load_global(name * "_UINT64", libgb.GrB_BinaryOp))
    end

    if name ∈ floats
        binary[Float32] = TypedBinaryOperator(load_global(name * "_FP32", libgb.GrB_BinaryOp))
        binary[Float64] = TypedBinaryOperator(load_global(name * "_FP64", libgb.GrB_BinaryOp))
    end
    if name ∈ positionals
        binary[Any] = TypedBinaryOperator(load_global(name * "_INT64", libgb.GrB_BinaryOp))
    end
    name = "GxB_" * name[5:end]
    if name ∈ complexes
        binary[ComplexF32] = TypedBinaryOperator(load_global(name * "_FC32", libgb.GrB_BinaryOp))
        binary[ComplexF64] = TypedBinaryOperator(load_global(name * "_FC64", libgb.GrB_BinaryOp))
    end
end

ztype(::TypedBinaryOperator{X, Y, Z}) where {X, Y, Z} = Z
xtype(::TypedBinaryOperator{X, Y, Z}) where {X, Y, Z} = X
ytype(::TypedBinaryOperator{X, Y, Z}) where {X, Y, Z} = Y


text:
const ρ = 0.72
const ρϵ = 0.71

immutable ℕ end

Base.in(n::Integer,::Type{ℕ}) = n > 0
Base.in(n::Real,::Type{ℕ}) = (ν = round(Int,n); n == ν && ν ∈ ℕ)
Base.in(n::Complex,::Type{ℕ}) = imag(n) == 0 && real(n) ∈ ℕ
Base.in(n::Dual,::Type{ℕ}) = dualpart(n) == 0 && realpart(n) ∈ ℕ

immutable ℕ₀ end

Base.in(n::Integer,::Type{ℕ₀}) = n ≥ 0
Base.in(n::Real,::Type{ℕ₀}) = (ν = round(Int,n); n == ν && ν ∈ ℕ₀)
Base.in(n::Complex,::Type{ℕ₀}) = imag(n) == 0 && real(n) ∈ ℕ₀
Base.in(n::Dual,::Type{ℕ₀}) = dualpart(n) == 0 && realpart(n) ∈ ℕ₀

immutable ℤ end

Base.in(n::Integer,::Type{ℤ}) = true
Base.in(n::Real,::Type{ℤ}) = n == round(Int,n)
Base.in(n::Complex,::Type{ℤ}) = imag(n) == 0 && real(n) ∈ ℤ
Base.in(n::Dual,::Type{ℤ}) = dualpart(n) == 0 && realpart(n) ∈ ℤ

abeqcd(a,b,cd) = isequal(a,b) && isequal(b,cd)
abeqcd(a,b,c,d) = isequal(a,c) && isequal(b,d)

absarg(z) = abs(angle(z))

sqrtatanhsqrt(x) = x == 0 ? one(x) : (s = sqrt(-x); atan(s)/s)
sqrtasinsqrt(x) = x == 0 ? one(x) : (s = sqrt(x); asin(s)/s)
sinnasinsqrt(n,x) = x == 0 ? one(x) : (s = sqrt(x); sin(n*asin(s))/(n*s))
cosnasinsqrt(n,x) = cos(n*asin(sqrt(x)))
expnlog1pcoshatanhsqrt(n,x) = x == 0 ? one(x) : (s = sqrt(x); (exp(n*log1p(s))+exp(n*log1p(-s)))/2)
expnlog1psinhatanhsqrt(n,x) = x == 0 ? one(x) : (s = sqrt(x); (exp(n*log1p(s))-exp(n*log1p(-s)))/(2n*s))

sqrtatanhsqrt{T<:Real}(x::Union{T,Dual{T}}) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); atanh(s)/s) : (s = sqrt(-x); atan(s)/s)
sqrtasinsqrt{T<:Real}(x::Union{T,Dual{T}}) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); asin(s)/s) : (s = sqrt(-x); asinh(s)/s)
sinnasinsqrt{T<:Real}(n,x::Union{T,Dual{T}}) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); sin(n*asin(s))/(n*s)) : (s = sqrt(-x); sinh(n*asinh(s))/(n*s))
cosnasinsqrt{T<:Real}(n,x::Union{T,Dual{T}}) = x > 0 ? cos(n*asin(sqrt(x))) : cosh(n*asinh(sqrt(-x)))
expnlog1pcoshatanhsqrt{T<:Real}(n,x::Union{T,Dual{T}}) = x == 0 ? one(x) : x > 0 ? exp(n/2*log1p(-x))*cosh(n*atanh(sqrt(x))) : exp(n/2*log1p(-x))*cos(n*atan(sqrt(-x)))
expnlog1psinhatanhsqrt{T<:Real}(n,x::Union{T,Dual{T}}) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); exp(n/2*log1p(-x))*sinh(n*atanh(s))/(n*s)) : (s = sqrt(-x); exp(n/2*log1p(-x))*sin(n*atan(s))/(n*s))

expm1nlog1p(n,x) = x == 0 ? one(x) : expm1(n*log1p(x))/(n*x)

speciallog(x) = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); 3(atanh(s)-s)/s^3) : (s = sqrt(-x); 3(s-atan(s))/s^3)
function speciallog(x::Float64)
    if x > 0.2
        s = sqrt(x)
        3(atanh(s)-s)/s^3
    elseif x < -0.2
        s = sqrt(-x)
        3(s-atan(s))/s^3
    else
        speciallogseries(x)
    end
end
function speciallog(x::Complex128)
    if abs(x) > 0.2
        s = sqrt(-x)
        3(s-atan(s))/s^3
    else
        speciallogseries(x)
    end
end
# The Maclaurin series fails to be accurate to 1e-15 near x ≈ ±0.2. So we use a highly accurate Chebyshev expansion.
speciallogseries(x::Union{Float64,Dual128}) = @clenshaw(5.0x,1.0087391788544393911192,1.220474262857857637288e-01,8.7957928919918696061703e-03,6.9050958578444820505037e-04,5.7037120050065804396306e-05,4.8731405131379353370205e-06,4.2648797509486828820613e-07,3.800372208946157617901e-08,3.434168059359993493634e-09,3.1381484326392473547608e-10,2.8939845618385022798906e-11,2.6892186934806386106143e-12,2.5150879096374730760324e-13,2.3652490233687788117887e-14,2.2349973917002118259929e-15,2.120769988408948118084e-16)
speciallogseries(x::Union{Complex128,DualComplex256}) = @evalpoly(x,1.0000000000000000000000,5.9999999999999999999966e-01,4.2857142857142857142869e-01,3.3333333333333333333347e-01,2.7272727272727272727292e-01,2.3076923076923076923072e-01,1.9999999999999999999996e-01,1.7647058823529411764702e-01,1.5789473684210526315786e-01,1.4285714285714285714283e-01,1.3043478260869565217384e-01,1.2000000000000000000000e-01,1.1111111111111111111109e-01,1.0344827586206896551722e-01,9.6774193548387096774217e-02,9.0909090909090909090938e-02,8.5714285714285714285696e-02,8.1081081081081081081064e-02,7.6923076923076923076907e-02,7.3170731707317073170688e-02)

tanpi(z) = sinpi(z)/cospi(z)
@vectorize_1arg Number tanpi

const libm = Base.libm_name

unsafe_gamma(x::Float64) = ccall((:tgamma,libm),  Float64, (Float64,), x)
unsafe_gamma(x::Float32) = ccall((:tgammaf,libm),  Float32, (Float32,), x)
unsafe_gamma(x::Real) = unsafe_gamma(float(x))
function unsafe_gamma(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_gamma, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, Base.MPFR.ROUNDING_MODE[end])
    return z
end
unsafe_gamma(z::Complex) = gamma(z)
unsafe_gamma(z::Dual) = (r = realpart(z);w = unsafe_gamma(r); dual(w, w*digamma(r)*dualpart(z)))
@vectorize_1arg Number unsafe_gamma

"""
`@lanczosratio(z,ϵ,c₀,c...)`

Compute `∑_{i=1}^N cᵢ/(z-1+i)/(z-1+i+ϵ) / ( c₀ + ∑_{i=1}^N cᵢ/(z-1+i) )`
"""
macro lanczosratio(z, ϵ, c₀, c...)
    ex_num = :(zero(zm1))
    ex_den = esc(c₀)
    for i = 0:length(c)-1
        temp = :(inv(zm1+$i))
        ex_num = :(muladd($(esc(c[i+1])),$temp*inv(zm1pϵ+$i),$ex_num))
        ex_den = :(muladd($(esc(c[i+1])),$temp,$ex_den))
    end
    ex = :($ex_num/$ex_den)
    Expr(:block, :(zm1 = $(esc(z))), :(zm1pϵ = $(esc(z))+$(esc(ϵ))), ex)
end

Hsum(z::Union{Float64,Complex128,Dual128,DualComplex256},ϵ::Union{Float64,Complex128,Dual128,DualComplex256}) = @lanczosratio(z,ϵ,0.99999999999999709182,57.156235665862923517,-59.597960355475491248,14.136097974741747174,-0.49191381609762019978,0.33994649984811888699E-4,0.46523628927048575665E-4,-0.98374475304879564677E-4,0.15808870322491248884E-3,-0.21026444172410488319E-3,0.21743961811521264320E-3,-0.16431810653676389022E-3,0.84418223983852743293E-4,-0.26190838401581408670E-4,0.36899182659531622704E-5)

function H(z::Union{Float64,Complex128,Dual128,DualComplex256},ϵ::Union{Float64,Complex128,Dual128,DualComplex256})
    zm0p5 = z-0.5
    zpgm0p5 = zm0p5+4.7421875
    if real(z) ≥ 1/2
        if z == z+ϵ # ϵ is numerical 0
            zm0p5/zpgm0p5 + log(zpgm0p5) - 1. - Hsum(z,ϵ)
        else
            expm1( zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*Hsum(z,ϵ)) )/ϵ
        end
    else
        tpz = tanpi(z)
        if z == z+ϵ # ϵ is numerical 0
            H(1.-z,ϵ) - π/tpz
        else
            temp = (cospi(ϵ) + sinpi(ϵ)/tpz)*H(1.-z,-ϵ) + .5ϵ*(π*sinc(.5ϵ))^2 - π*sinc(ϵ)/tpz
            temp/(1.-ϵ*temp)
        end
    end
end

"""
Compute the function (1/Γ(z)-1/Γ(z+ϵ))/ϵ
"""
function G(z::Union{Float64,Complex128,Dual128,DualComplex256},ϵ::Union{Float64,Complex128,Dual128,DualComplex256})
    n,zpϵ = round(Int,real(z)),z+ϵ
    if abs(ϵ) > 0.1
        (inv(unsafe_gamma(z))-inv(unsafe_gamma(zpϵ)))/ϵ
    elseif z ≠ zpϵ
        m = round(Int,real(zpϵ))
        if z == n && n ≤ 0
            -inv(ϵ*unsafe_gamma(zpϵ))
        elseif zpϵ == m && m ≤ 0
            inv(ϵ*unsafe_gamma(z))
        elseif abs(z+abs(n)) < abs(zpϵ+abs(m))
            H(z,ϵ)/unsafe_gamma(zpϵ)
        else
            H(zpϵ,-ϵ)/unsafe_gamma(z)
        end
    else # ϵ is numerical 0
        if z == n && n ≤ 0
            (-1)^(n+1)*unsafe_gamma(1-n)
        else
            digamma(z)/unsafe_gamma(z)
        end
    end
end

G(z,ϵ) = ϵ == 0 ? digamma(z)/unsafe_gamma(z) : (inv(unsafe_gamma(z))-inv(unsafe_gamma(z+ϵ)))/ϵ

"""
Compute the function ((z+ϵ)ₘ-(z)ₘ)/ϵ
"""
function P(z::Number,ϵ::Number,m::Int)
    n₀ = -round(Int,real(z))
    if ϵ == 0
        if 0 ≤ n₀ < m
            ret1,ret2,n = one(z),zero(z),0
            while n < m
                n == n₀ && (n+=1; continue)
                ret1 *= z+n
                ret2 += inv(z+n)
                n+=1
            end
            ret1 + pochhammer(z,m)*ret2
        else
            ret = zero(z)
            for n=0:m-1
                ret += inv(z+n)
            end
            pochhammer(z,m)*ret
        end
    else
        if 0 ≤ n₀ < m
            zpϵ,ret1,ret2,n = z+ϵ,one(z),zero(z),0
            while n < m
                n == n₀ && (n+=1; continue)
                ret1 *= zpϵ+n
                ret2 += log1p(ϵ/(z+n))
                n+=1
            end
            ret1 + pochhammer(z,m)*expm1(ret2)/ϵ
        else
            ret = zero(z)
            for n=0:m-1
                ret += log1p(ϵ/(z+n))
            end
            pochhammer(z,m)*expm1(ret)/ϵ
        end
    end
end

E(z,ϵ) = ϵ == 0 ? z : expm1(ϵ*z)/ϵ

@vectorize_2arg Number E
@vectorize_2arg Number G
@vectorize_2arg Number H
G(z::AbstractVector{BigFloat},ϵ::BigFloat) = BigFloat[G(zi,ϵ) for zi in z]

# Transformation formula w = 1-z

reconeα₀(a,b,c,m::Int,ϵ) = ϵ == 0 ? (-1)^m*gamma(m)*gamma(c)/(gamma(a+m)*gamma(b+m)) : gamma(c)/(ϵ*gamma(1-m-ϵ)*gamma(a+m+ϵ)*gamma(b+m+ϵ))
reconeβ₀(a,b,c,w,m::Int,ϵ) = abs(ϵ) > 0.1 ? ( pochhammer(float(a),m)*pochhammer(b,m)/(gamma(1-ϵ)*gamma(a+m+ϵ)*gamma(b+m+ϵ)*gamma(m+1)) - w^ϵ/(gamma(a)*gamma(b)*gamma(m+1+ϵ)) )*gamma(c)*w^m/ϵ : ( (G(1.0,-ϵ)/gamma(m+1)+G(m+1.0,ϵ))/(gamma(a+m+ϵ)*gamma(b+m+ϵ)) - (G(float(a)+m,ϵ)/gamma(b+m+ϵ)+G(float(b)+m,ϵ)/gamma(a+m))/gamma(m+1+ϵ) - E(log(w),ϵ)/(gamma(a+m)*gamma(b+m)*gamma(m+1+ϵ)) )*gamma(c)*pochhammer(float(a),m)*pochhammer(b,m)*w^m
reconeγ₀(a,b,c,w,m::Int,ϵ) = gamma(c)*pochhammer(float(a),m)*pochhammer(b,m)*w^m/(gamma(a+m+ϵ)*gamma(b+m+ϵ)*gamma(m+1)*gamma(1-ϵ))

function Aone(a,b,c,w,m::Int,ϵ)
    αₙ = reconeα₀(a,b,c,m,ϵ)*one(w)
    ret = m ≤ 0 ? zero(w) : αₙ
    for n = 0:m-2
        αₙ *= (a+n)*(b+n)/((n+1)*(1-m-ϵ+n))*w
        ret += αₙ
    end
    ret
end

function Bone(a,b,c,w,m::Int,ϵ)
    βₙ,γₙ = reconeβ₀(a,b,c,w,m,ϵ)*one(w),reconeγ₀(a,b,c,w,m,ϵ)*w
    ret,err,n = βₙ,1.0,0
    while err > 10eps()
        βₙ = (a+m+n+ϵ)*(b+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(b+m+n)/(m+n+1) - (a+m+n) - (b+m+n) - ϵ + (a+m+n+ϵ)*(b+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
        ret += βₙ
        err = errcheck(βₙ,ret)
        γₙ *= (a+m+n)*(b+m+n)/((m+n+1)*(n+1-ϵ))*w
        n+=1
    end
    ret
end

function _₂F₁one(a,b,c,z)
    m = round(Int,real(c-(a+b)))
    ϵ = c-(a+b)-m
    w = 1-z
    (-1)^m/sinc(ϵ)*(Aone(a,b,c,w,m,ϵ) + Bone(a,b,c,w,m,ϵ))
end

# Transformation formula w = 1/z

recInfα₀(a,b,c,m::Int,ϵ) = ϵ == 0 ? (-1)^m*gamma(m)*gamma(c)/(gamma(a+m)*gamma(c-a)) : gamma(c)/(ϵ*gamma(1-m-ϵ)*gamma(a+m+ϵ)*gamma(c-a))
recInfβ₀(a,b,c,w,m::Int,ϵ) = abs(ϵ) > 0.1 ? ( pochhammer(float(a),m)*pochhammer(float(1-c+a),m)/(gamma(1-ϵ)*gamma(a+m+ϵ)*gamma(c-a)*gamma(m+1)) - (-w)^ϵ*pochhammer(float(1-c+a)+ϵ,m)/(gamma(a)*gamma(c-a-ϵ)*gamma(m+1+ϵ)) )*gamma(c)*w^m/ϵ : ( (pochhammer(float(1-c+a)+ϵ,m)*G(1.0,-ϵ)-P(1-c+a,ϵ,m)/gamma(1-ϵ))/(gamma(c-a)*gamma(a+m+ϵ)*gamma(m+1)) + pochhammer(float(1-c+a)+ϵ,m)*( (G(m+1.0,ϵ)/gamma(a+m+ϵ) - G(float(a)+m,ϵ)/gamma(m+1+ϵ))/gamma(c-a) - (G(float(c-a),-ϵ) - E(-log(-w),-ϵ)/gamma(c-a-ϵ))/(gamma(m+1+ϵ)*gamma(a+m)) ) )*gamma(c)*pochhammer(float(a),m)*w^m
recInfγ₀(a,b,c,w,m::Int,ϵ) = gamma(c)*pochhammer(float(a),m)*pochhammer(float(1-c+a),m)*w^m/(gamma(a+m+ϵ)*gamma(c-a)*gamma(m+1)*gamma(1-ϵ))

function AInf(a,b,c,w,m::Int,ϵ)
    αₙ = recInfα₀(a,b,c,m,ϵ)*one(w)
    ret = m ≤ 0 ? zero(w) : αₙ
    for n = 0:m-2
        αₙ *= (a+n)*(1-c+a+n)/((n+1)*(1-m-ϵ+n))*w
        ret += αₙ
    end
    ret
end

function BInf(a,b,c,w,m::Int,ϵ)
    βₙ,γₙ = recInfβ₀(a,b,c,w,m,ϵ)*one(w),recInfγ₀(a,b,c,w,m,ϵ)*w
    ret,err,n = βₙ,1.0,0
    while err > 10eps()
        βₙ = (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(1-c+a+m+n)/(m+n+1) - (a+m+n) - (1-c+a+m+n) - ϵ + (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
        ret += βₙ
        err = errcheck(βₙ,ret)
        γₙ *= (a+m+n)*(1-c+a+m+n)/((m+n+1)*(n+1-ϵ))*w
        n+=1
    end
    ret
end

function _₂F₁Inf(a,b,c,z)
    m = round(Int,real(b-a))
    ϵ = b-a-m
    w = inv(z)
    (-1)^m*(-w)^a/sinc(ϵ)*(AInf(a,b,c,w,m,ϵ) + BInf(a,b,c,w,m,ϵ))
end

function _₂F₁maclaurin(a::Number,b::Number,c::Number,z::Number)
    T = promote_type(typeof(a),typeof(b),typeof(c),typeof(z))
    S₀,S₁,err,j = one(T),one(T)+a*b*z/c,one(real(T)),1
    while err > 10eps(T)
        rⱼ = (a+j)/(j+1)*(b+j)/(c+j)
        S₀,S₁ = S₁,S₁+(S₁-S₀)*rⱼ*z
        err = errcheck(S₁-S₀,S₀)
        j+=1
    end
    return S₁
end

function _₂F₁maclaurinalt(a::Number,b::Number,c::Number,z::Number)
    T = promote_type(typeof(a),typeof(b),typeof(c),typeof(z))
    C,S,err,j = one(T),one(T),one(real(T)),0
    while err > 10eps(T)
        C *= (a+j)/(j+1)*(b+j)/(c+j)*z
        S += C
        err = errcheck(C,S)
        j+=1
    end
    return S
end

function _₂F₁continuation(s::Number,t::Number,c::Number,z₀::Number,z::Number)
    T = promote_type(typeof(s),typeof(t),typeof(c),typeof(z₀),typeof(z))
    izz₀,d0,d1 = inv(z-z₀),one(T),s/(2s-t+one(T))*((s+1)*(1-2z₀)+(t+1)*z₀-c)
    S₀,S₁,izz₀j,err,j = one(T),one(T)+d1*izz₀,izz₀,one(real(T)),2
    while err > 10eps(T)
        d0,d1,izz₀j = d1,(j+s-one(T))/j/(j+2s-t)*(((j+s)*(1-2z₀)+(t+1)*z₀-c)*d1 + z₀*(1-z₀)*(j+s-2)*d0),izz₀j*izz₀
        S₀,S₁ = S₁,S₁+d1*izz₀j
        err = errcheck(S₁-S₀,S₀)
        j+=1
    end
    return S₁
end

function _₂F₁continuationalt(a::Number,c::Number,z₀::Number,z::Number)
    T = promote_type(typeof(a),typeof(c),typeof(z₀),typeof(z))
    izz₀ = inv(z-z₀)
    e0,e1 = one(T),(a+one(T))*(one(T)-2z₀)+(2a+one(T))*z₀-c
    f0,f1 = zero(T),one(T)-2z₀
    cⱼ = log(z₀-z)+2digamma(one(T))-digamma(a)-digamma(c-a)
    S₀ = cⱼ
    cⱼ += 2/one(T)-one(T)/a
    C = a*izz₀
    S₁,err,j = S₀+(e1*cⱼ-f1)*C,one(real(T)),2
    while err > 10eps(T)
        f0,f1 = f1,(((j+a)*(1-2z₀)+(2a+1)*z₀-c)*f1+z₀*(1-z₀)*(j-1)*f0+(1-2z₀)*e1+2z₀*(1-z₀)*e0)/j
        e0,e1 = e1,(((j+a)*(1-2z₀)+(2a+1)*z₀-c)*e1+z₀*(1-z₀)*(j-1)*e0)/j
        C *= (a+j-1)*izz₀/j
        cⱼ += 2/T(j)-one(T)/(a+j-one(T))
        S₀,S₁ = S₁,S₁+(e1*cⱼ-f1)*C
        err = errcheck(S₁-S₀,S₀)
        j+=1
    end
    return S₁
end

function _₂F₁logsum(a::Number,b::Number,z::Number,w::Number,s::Int)
    T = promote_type(typeof(a),typeof(b),typeof(z),typeof(w))
    cⱼ = 2digamma(one(T))-digamma(a)-digamma(b)+s*log1p(-z)
    C,S,err,j = one(T),cⱼ,one(real(T)),0
    while err > 10eps(T)
        C *= (a+j)/(j+1)^2*(b+j)*w
        cⱼ += 2/(j+one(T))-one(T)/(a+j)-one(T)/(b+j)
        S += C*cⱼ
        err = errcheck(C,S)
        j+=1
    end
    return S
end

function _₂F₁logsumalt(a::Number,b::Number,z::Number,w::Number)
    T = promote_type(typeof(a),typeof(b),typeof(z),typeof(w))
    d,cⱼ = one(T)-b,2digamma(one(T))-digamma(a)-digamma(b)-log(-w)
    C,S,err,j = one(T),cⱼ,one(real(T)),0
    while err > 10eps(T)
        C *= (a+j)/(j+1)^2*(d+j)*w
        cⱼ += 2/(j+one(T))-one(T)/(a+j)+one(T)/(b-(j+one(T)))
        S += C*cⱼ
        err = errcheck(C,S)
        j+=1
    end
    return S
end

function _₂F₁taylor(a::Number,b::Number,c::Number,z::Number)
    T = promote_type(typeof(a),typeof(b),typeof(c),typeof(z))
    z₀ = abs(z) < 1 ? ρϵ*sign(z) : sign(z)/ρϵ
    q₀,q₁ = _₂F₁(a,b,c,z₀),a*b/c*_₂F₁(a+1,b+1,c+1,z₀)
    S₀,zz₀ = q₀,z-z₀
    S₁,err,zz₀j,j = S₀+q₁*zz₀,one(real(T)),zz₀,0
    while err > 10eps(T)
        q₀,q₁ = q₁,((j*(2z₀-one(T))-c+(a+b+one(T))*z₀)*q₁ + (a+j)*(b+j)/(j+one(T))*q₀)/(z₀*(one(T)-z₀)*(j+2))
        zz₀j *= zz₀
        S₀,S₁ = S₁,S₁+q₁*zz₀j
        err = errcheck(S₀-S₁,S₀)
        j+=1
    end
    return S₁
end

function _₃F₂maclaurin(a₁::Number,a₂::Number,a₃::Number,b₁::Number,b₂::Number,z::Number)
    T = promote_type(typeof(a₁),typeof(a₂),typeof(a₃),typeof(b₁),typeof(b₂),typeof(z))
    S₀,S₁,err,j = one(T),one(T)+(a₁*a₂*a₃*z)/(b₁*b₂),one(real(T)),1
    while err > 100eps(T)
        rⱼ = ((a₁+j)*(a₂+j)*(a₃+j))/((b₁+j)*(b₂+j)*(j+1))
        S₀,S₁ = S₁,S₁+(S₁-S₀)*rⱼ*z
        err = errcheck(S₁-S₀,S₀)
        j+=1
    end
    return S₁
end

function mFnmaclaurin{S<:Number,V<:Number}(a::AbstractVector{S},b::AbstractVector{V},z::Number)
    T = promote_type(S,V,typeof(z))
    S₀,S₁,err,j = one(T),one(T)+prod(a)*z/prod(b),one(real(T)),1
    while err > 100eps(T)
        rⱼ = inv(j+one(T))
        for i=1:length(a) rⱼ *= a[i]+j end
        for i=1:length(b) rⱼ /= b[i]+j end
        S₀,S₁ = S₁,S₁+(S₁-S₀)*rⱼ*z
        err = errcheck(S₁-S₀,S₀)
        j+=1
    end
    return S₁
end

hypot2(a,b) = hypot(a,b)
hypot2{T<:AbstractFloat}(a::Complex{T},b::Complex{T}) = hypot(hypot(a.re,b.re),hypot(a.im,b.im))
hypot2{T<:AbstractFloat}(a::T,b::Complex{T}) = hypot(hypot(a,b.re),b.im)
hypot2{T<:AbstractFloat}(a::Complex{T},b::T) = hypot(hypot(a.re,b),a.im)

errcheck(x,y) = abs(x/y)
errcheck(x::Dual,y::Dual) = hypot2(realpart(x),dualpart(x))/hypot2(realpart(y),dualpart(y))

# This is from Base because real(::Type{Dual{BigFloat}}) doesn't exist in the scop of rtoldefault, required by isapprox.

# isapprox: approximate equality of numbers
function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0)
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y)))
end

const ≈ = isapprox
≉(x,y) = !(x ≈ y)

# default tolerance arguments
rtoldefault{T<:AbstractFloat}(::Type{T}) = sqrt(eps(T))
rtoldefault{T<:Real}(::Type{T}) = 0
rtoldefault{T<:Real}(::Type{Dual{T}}) = rtoldefault(T)
rtoldefault{T<:Number,S<:Number}(x::Union{T,Type{T}}, y::Union{S,Type{S}}) = rtoldefault(promote_type(real(T),real(S)))


text:
using Distributions
using JLD
using ArgParse
using NCDatasets
using LinearAlgebra
# Import EnsembleKalmanProcesses modules
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributions

include(joinpath(@__DIR__, "helper_funcs.jl"))

"""
ek_update(iteration_::Int64)

Update CLIMAParameters ensemble using Ensemble Kalman Inversion,
return the ensemble and the parameter names.
"""
function ek_update(iteration_::Int64)
    # Recover versions from last iteration
    versions = readlines("versions_$(iteration_).txt")
    n_params = 0
    u_names = String[]
    open("$(versions[1]).output/$(versions[1])", "r") do io
        append!(
            u_names,
            [strip(string(split(line, "(")[1]), [' ']) for (index, line) in enumerate(eachline(io)) if index % 3 == 2],
        )
        n_params = length(u_names)
    end
    # Recover ensemble from last iteration, [N_ens, N_params]
    u = zeros(length(versions), n_params)
    for (ens_index, version_) in enumerate(versions)
        open("$(version_).output/$(version_)", "r") do io
            u[ens_index, :] = [parse(Float64, line) for (index, line) in enumerate(eachline(io)) if index % 3 == 0]
        end
    end
    u = Array(u')
    # Set averaging period for loss function
    t0_ = 0.0
    tf_ = 1800.0
    # Get observations (CliMA)
    yt = zeros(0)
    yt_var_list = []
    y_names = ["u", "v"]
    yt_, yt_var_ = get_clima_profile("truth_output", y_names, ti = t0_, tf = tf_, get_variance = true)
    # Add nugget to variance (regularization)
    yt_var_ = yt_var_ + Matrix(0.1I, size(yt_var_)[1], size(yt_var_)[2])
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)

    # Get outputs
    g_names = y_names
    g_ = get_clima_profile("$(versions[1]).output", g_names, ti = t0_, tf = tf_)
    g_ens = zeros(length(versions), length(g_))
    for (ens_index, version) in enumerate(versions)
        g_ens[ens_index, :] = get_clima_profile("$(version).output", g_names, ti = t0_, tf = tf_)
    end
    g_ens = Array(g_ens')
    # Construct EKP
    ekobj = EnsembleKalmanProcess(u, yt_, yt_var_, Inversion())
    # Advance EKP
    update_ensemble!(ekobj, g_ens)
    # Get new step
    u_new = get_u_final(ekobj)
    return u_new, u_names
end

# Read iteration number of ensemble to be recovered
s = ArgParseSettings()
@add_arg_table s begin
    "--iteration"
    help = "Calibration iteration number"
    arg_type = Int
    default = 1
end
parsed_args = parse_args(ARGS, s)
iteration_ = parsed_args["iteration"]

# Perform update
ens_new, param_names = ek_update(iteration_)
params_arr = [col[:] for col in eachcol(ens_new)]

# Generate new identifiers
versions = map(param -> generate_cm_params(param, param_names), params_arr)
open("versions_$(iteration_+1).txt", "w") do io
    for version in versions
        write(io, "clima_param_defs_$(version).jl\n")
    end
end


text:
"""
Container that mimics R vector behaviour.
Elements could be accessed either by indices as a normal vector,
or (optionally) by string keys as a dictionary.
"""
struct DictoVec{T}
    data::Vector{T}
    name2index::Dict{RString, Int}
    index2name::Vector{Union{RString, Nothing}}

    function DictoVec(data::AbstractVector{T}, names::AbstractVector{<:AbstractString} = Vector{RString}()) where T
        if !isempty(names) && length(data) != length(names)
            throw(DimensionMismatch("Lengths of data ($(length(data))) and element names ($(length(names))) differ"))
        end
        n2i = Dict{RString, Int}()
        i2n = fill!(similar(data, Union{RString, Nothing}), nothing)
        @inbounds for (i, k) in enumerate(names)
            if k != "" && k !== nothing
                n2i[k] = i
                i2n[i] = k
            else
                i2n[i] = nothing
            end
        end
        new{T}(data, n2i, i2n)
    end
end

Base.eltype(::Type{DictoVec{T}}) where T = T
Base.eltype(dict::DictoVec) = eltype(typeof(dict))
Base.length(dict::DictoVec) = length(dict.data)
Base.isempty(dict::DictoVec) = isempty(dict.data)

# key-based indexing
Base.haskey(dict::DictoVec, key) = haskey(dict.name2index, key)

function Base.setindex!(dict::DictoVec, value, key)
    ix = get(dict.name2index, key, 0)
    if ix > 0 # existing key
        setindex!(dict.data, value, ix)
    else # new key
        _key = convert(keytype(dict.name2index), key) # throws if key is not compatible
        push!(dict.data, value) # first try to add value (throws if it's not compatible)
        push!(dict.index2name, _key)
        dict.name2index[_key] = length(dict.data)
    end
end

function Base.getindex(dict::DictoVec, key)
    ix = get(dict.name2index, key, 0)
    if ix > 0
        return getindex(dict.data, ix)
    else
        throw(KeyError(key))
    end
end

# integer-based indexing (overrides key-based indexing)
Base.setindex!(dict::DictoVec, value, index::Integer) =
    setindex!(dict.data, value, index)
Base.setindex!(dict::DictoVec, values, indices::AbstractVector{<:Integer}) =
    setindex!(dict.data, values, indices)

Base.getindex(dict::DictoVec, index::Integer) = getindex(dict.data, index)
Base.getindex(dict::DictoVec, indices::AbstractVector{<:Integer}) = getindex(dict.data, indices)

Base.get(dict::DictoVec, index::Integer, default) = get(dict.data, index, default)
Base.get(f::Function, dict::DictoVec, index::Integer) = index ∈ eachindex(dict.data) ? dict.data[index] : f()

function Base.get(dict::DictoVec, key, default)
    ix = get(dict.name2index, key, 0)
    return ix > 0 ? dict.data[ix] : default
end

function Base.get(f::Function, dict::DictoVec, key)
    ix = get(dict.name2index, key, 0)
    return get(f, dict, ix)
end

function Base.delete!(dict::DictoVec, key)
    ix = pop!(dict.name2index, key, 0)
    if ix > 0
        deleteat!(dict.data, ix)
        deleteat!(dict.index2name, ix)
        if ix <= length(dict)
            # update indices
            for (k, oldix) in pairs(dict.name2index)
                if oldix > ix
                    dict.name2index[k] = oldix - 1
                end
            end
        end
    end
    return dict
end

Base.keys(dict::DictoVec) = keys(dict.name2index)

Base.values(dict::DictoVec) = dict.data

function Base.show(io::IO, dict::DictoVec)
    first = true
    print(io, typeof(dict), "(")
    n = 0
    for i in eachindex(dict.data)
        first || print(io, ',')
        first = false
        k = dict.index2name[i]
        if k !== nothing
            show(io, k)
            print(io, "=>")
        end
        show(io, dict.data[i])
        n += 1
        # limit && n >= 10 && (print(io, "…"); break)
    end
    print(io, ")")
end

function Base.convert(::Type{Dict{RString,Any}}, dv::DictoVec)
    res = Dict{RString,Any}()
    # elements that are only referenced by the index would not be in the result
    for (k, v) in dv.name2index
        res[k] = dv.data[v]
    end
    return res
end


text:
import StatsBase

"""
"""
function get_binary_thresholds(
        yscore::AbstractVector{<:Real};
        additionalthreshold::AbstractFloat = 0.5,
        )
    if !all(0 .<= yscore .<= 1)
        error("not all scores are in [0,1]")
    end
    result = sort(
        unique(
            vcat(
                0 - eps(),
                0,
                0 + eps(),
                1 - eps(),
                1,
                1 + eps(),
                additionalthreshold,
                unique(yscore),
                )
            );
        rev = false,
        )
    return result
end


text:
"""
    Structure holding the residual of the linearized residual r. This residual linearization is computed at
    a linearization point z0, θ0.
    The residual is decomposed into blocks to best exploit the structure of the problem.
         [rdyn ] -> Linearized dynamics constraints -> size = nx
    r  = |rrst | -> Rest of the linearized constraints -> size = ny
         [rbil ] -> Bilinear complementarity constraints -> size = ny
    The linear parts:
        [rdyn]
		|rrst| = [r0 + rz0*(z-z0) + rθ0*(θ-θ0)
		[....]
    The bilinear part:
        rbil = y1 .* y2 .- κ
"""
mutable struct RLin{T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}
    # Reference residual
    rdyn0::SVector{nx,T}
    rrst0::SVector{ny,T}
    rbil0::SVector{ny,T}

    # Residual
    rdyn::SVector{nx,T}
    rrst::SVector{ny,T}
    rbil::SVector{ny,T}
    rbil_cache::Vector{T}

    # Reference residual jacobian rz0
    Dx::SMatrix{nx,nx,T,nxx}
    Dy1::SMatrix{nx,ny,T,nxy}
    Rx::SMatrix{ny,nx,T,nxy}
    Ry1::SMatrix{ny,ny,T,nyy}
    Ry2::SVector{ny,T}

    # Reference residual jacobian rθ0
    rθdyn::SMatrix{nx,nθ,T,nxθ}
    rθrst::SMatrix{ny,nθ,T,nyθ}
    rθbil::SMatrix{ny,nθ,T,nyθ}

    # Reference z0 and θ0
    x0::SVector{nx,T}
    y10::SVector{ny,T}
    y20::SVector{ny,T}
    θ0::SVector{nθ,T}
    # Values of z and θ
    x::SVector{nx,T}
    y1::SVector{ny,T}
    y2::SVector{ny,T}
    θ::SVector{nθ,T}

    # Indices
    nz::Int
    nθ::Int
    ix::SVector{nx,Int}
    iy1::SVector{ny,Int}
    iy2::SVector{ny,Int}
    iθ::SVector{nθ,Int}
    idyn::SVector{nx,Int}
    irst::SVector{ny,Int}
    ibil::SVector{ny,Int}
    ialt::SVector{nc,Int}

    # Altitude
    alt::SVector{nc,T}
    alt_zeros::SVector{nn,T}
end

function RLin(s::Simulation, z0::AbstractVector{T}, θ0::AbstractVector{T},
        r0::AbstractVector{T}, rz0::AbstractMatrix{T}, rθ0::AbstractMatrix{T}) where {T}

    model = s.model
    env = s.env

    nq = model.nq
    nc = model.nc
    nb = nc * friction_dim(env)
    nz = num_var(model, env)
    nθ = num_data(model)
    nx = nq
    ny = 2nc + nb
    nn = nc + nb

    # Terms
    off = 0
    idyn, irst, ibil, ialt = linearization_term_index(model, env)

    # Vars
    off = 0
    ix, iy1, iy2 = linearization_var_index(model, env)
    iθ = Vector(1:nθ)

    # Residual
    rdyn0 = SVector{nx,T}(r0[idyn])
    rrst0 = SVector{ny,T}(r0[irst])
    rbil0 = SVector{ny,T}(r0[ibil])
    rdyn = zeros(SVector{nx,T})
    rrst = zeros(SVector{ny,T})
    rbil = zeros(SVector{ny,T})
    rbil_cache = zeros(ny)

    # Matrices
    Dx = SMatrix{nx,nx,T,nx^2}(rz0[idyn,ix])
    Dy1 = SMatrix{nx,ny,T,nx*ny}(rz0[idyn,iy1])
    Rx = SMatrix{ny,nx,T,nx*ny}(rz0[irst,ix])
    Ry1 = SMatrix{ny,ny,T,ny^2}(rz0[irst,iy1])
    Ry2 = SVector{ny,T}(diag(rz0[irst,iy2]))
    y1 = SVector{ny,T}(diag(rz0[ibil,iy2]))
    y2 = SVector{ny,T}(diag(rz0[ibil,iy1]))
    rθdyn = SMatrix{nx,nθ,T,nx*nθ}(rθ0[idyn,iθ])
    rθrst = SMatrix{ny,nθ,T,ny*nθ}(rθ0[irst,iθ])
    rθbil = SMatrix{ny,nθ,T,ny*nθ}(rθ0[ibil,iθ])

    # Vectors
    x0  = SVector{nx,T}(z0[ix])
    y10 = SVector{ny,T}(z0[iy1])
    y20 = SVector{ny,T}(z0[iy2])
    θ0  = SVector{nθ,T}(θ0)
    x  = zeros(SVector{nx,T})
    y1 = zeros(SVector{ny,T})
    y2 = zeros(SVector{ny,T})
    θ  = zeros(SVector{nθ,T})

    return RLin{T,nx,ny,nθ,nx^2,nx*ny,ny*ny,nx*nθ,ny*nθ,nc,nn}(
        rdyn0,
        rrst0,
        rbil0,
        rdyn,
        rrst,
        rbil,
        rbil_cache,

        Dx,
        Dy1,
        Rx,
        Ry1,
        Ry2,
        rθdyn,
        rθrst,
        rθbil,

        x0,
        y10,
        y20,
        θ0,
        x,
        y1,
        y2,
        θ,

        nz,
        nθ,
        SVector{nx,Int}(ix),
        SVector{ny,Int}(iy1),
        SVector{ny,Int}(iy2),
        SVector{nθ,Int}(iθ),
        SVector{nx,Int}(idyn),
        SVector{ny,Int}(irst),
        SVector{ny,Int}(ibil),
        SVector{nc,Int}(ialt),
        zeros(SVector{nc,T}),
        zeros(SVector{nc + nb,T}))
end

"""
    Structure holding the Jacobian of the linearized residual rz. This Jacobian linearization is computed at
    a linearization point z0, θ0.
    The Jacobian is decomposed into blocks to best exploit its sparsity and structure.
         [Dx  Dy1       0        ] -> Linearized dynamics constraints -> size = nx
    rz = |Rx  Ry1       Diag(Ry2)| -> Rest of the linearized constraints -> size = ny
         [0   Diag(y2)  Diag(y1) ] -> Bilinear complementarity constraints -> size = ny
          |      |         |
		  |		 |	       -> 1st set of variables associated with the bilinear constraints [γ1, b1, ψ] -> size = ny
          |      -> 2nd set of variables associated with the bilinear constraints [s1, η, s2] -> size = ny
		  -> Rest of the variables [q2] -> size = nx
    The constant parts are the ones that remains equals to the blocks in rz0:
        Dx, Dy1, Rx, Ry1, Ry2
    The variable parts:
        y1, y2.
	We define matrix A, it is used to initialize the interior point solver with
	a least square solution to the linearized constraints (i.e. excluding bilienar constraints).
	A = [Dx  Dy1       0        ] -> Linearized dynamics constraints -> size = nx
	    [Rx  Ry1       Diag(Ry2)] -> Rest of the linearized constraints -> size = ny
	δz = A' * inv(A * A') * δr
	We introduce the following notations:
	A' * inv(A * A') = [A1 A2]
					   |A3 A4|
	                   [A5 A6]
"""
mutable struct RZLin{T,nx,ny,nxx,nxy,nyy}
    # Reference residual jacobian rz0
    Dx::SMatrix{nx,nx,T,nxx}
    Dy1::SMatrix{nx,ny,T,nxy}
    Rx::SMatrix{ny,nx,T,nxy}
    Ry1::SMatrix{ny,ny,T,nyy}
    Ry2::SVector{ny,T}
    y1::SVector{ny,T}
    y2::SVector{ny,T}

    # Schur complement
    D::SMatrix{ny,ny,T,nyy}
    S::Schur{T,nx,ny,nxx,nxy,nyy}
    M::Matrix{T}
    M1::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}
    M2::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}
    M3::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}
    M4::SubArray{T,2,Matrix{T},Tuple{Vector{Int},Vector{Int}},false}

	# Least Squares
	A1::SMatrix{nx,nx,T,nxx}
	A2::SMatrix{nx,ny,T,nxy}
	A3::SMatrix{ny,nx,T,nxy}
	A4::SMatrix{ny,ny,T,nyy}
	A5::SMatrix{ny,nx,T,nxy}
	A6::SMatrix{ny,ny,T,nyy}

    # Indices
    ix::SVector{nx,Int}
    iy1::SVector{ny,Int}
    iy2::SVector{ny,Int}
    idyn::SVector{nx,Int}
    irst::SVector{ny,Int}
    ibil::SVector{ny,Int}
end

function RZLin(s::Simulation, rz0::AbstractMatrix{T}) where {T}
    model = s.model
    env = s.env

    nq = model.nq
    nc = model.nc
    nb = nc * friction_dim(env)

    nz = num_var(model, env)
    nx = nq
    ny = 2nc + nb

    # Terms
    off = 0
    idyn, irst, ibil, ialt = linearization_term_index(model, env)

    # Vars
    off = 0
    ix, iy1, iy2 = linearization_var_index(model, env)

    # Fill the matrix blocks rz0s
    Dx = SMatrix{nx,nx,T,nx^2}(rz0[idyn,ix])
    Dy1 = SMatrix{nx,ny,T,nx*ny}(rz0[idyn,iy1])
    Rx = SMatrix{ny,nx,T,nx*ny}(rz0[irst,ix])
    Ry1 = SMatrix{ny,ny,T,ny^2}(rz0[irst,iy1])
    Ry2 = SVector{ny,T}(diag(rz0[irst,iy2]))
    y1 = SVector{ny,T}(diag(rz0[ibil,iy2]))
    y2 = SVector{ny,T}(diag(rz0[ibil,iy1]))

    # Schur complement
    D = Ry1 - Diagonal(Ry2 .* y2 ./ y1)
    M = [Dx Dy1;
         Rx D  ]
    S = Schur(M; n=nx, m=ny)

    M = Array(M)
    M1 = view(M, collect(1:nx), collect(1:nx))
    M2 = view(M, collect(1:nx), collect(nx .+ (1:ny)))
    M3 = view(M, collect(nx .+ (1:ny)), collect(1:nx))
    M4 = view(M, collect(nx .+ (1:ny)), collect(nx .+ (1:ny)))

	# Least squares
	A = [Dx Dy1 zeros(nx, ny);
		 Rx Ry1 Diagonal(Ry2)]
	AA = A' * inv(A * A')
	A1 = SMatrix{nx,nx,T,nx*nx}(AA[1:nx,              1:nx])
	A2 = SMatrix{nx,ny,T,nx*ny}(AA[1:nx,              nx .+ (1:ny)])
	A3 = SMatrix{ny,nx,T,nx*ny}(AA[nx .+ (1:ny),      1:nx])
	A4 = SMatrix{ny,ny,T,ny*ny}(AA[nx .+ (1:ny),      nx .+ (1:ny)])
	A5 = SMatrix{ny,nx,T,nx*ny}(AA[nx + ny .+ (1:ny), 1:nx])
	A6 = SMatrix{ny,ny,T,ny*ny}(AA[nx + ny .+ (1:ny), nx .+ (1:ny)])

    return RZLin{T,nx,ny,nx^2,nx*ny,ny^2}(
        Dx,
        Dy1,
        Rx,
        Ry1,
        Ry2,
        y1,
        y2,
        D,
        S,
        Array(M),
        M1,
        M2,
        M3,
        M4,
		A1,
		A2,
		A3,
		A4,
		A5,
		A6,
        SVector{nx,Int}(ix),
        SVector{ny,Int}(iy1),
        SVector{ny,Int}(iy2),
        SVector{nx,Int}(idyn),
        SVector{ny,Int}(irst),
        SVector{ny,Int}(ibil),
        )
end

"""
    Structure holding the Jacobian of the linearized residual rθ. This Jacobian linearization is computed at
    a linearization point z0, θ0.
"""
mutable struct RθLin{T,nx,ny,nθ,nxθ,nyθ}
    # Reference residual jacobian rθ0
	rθ0::Matrix{T}
	rθdyn0::SMatrix{nx,nθ,T,nxθ}
	rθrst0::SMatrix{ny,nθ,T,nyθ}
    rθbil0::SMatrix{ny,nθ,T,nyθ}

    # Indices
    iθ::SVector{nθ,Int}
    idyn::SVector{nx,Int}
    irst::SVector{ny,Int}
    ibil::SVector{ny,Int}
end

function RθLin(s::Simulation, rθ0::AbstractMatrix{T}) where {T}
    model = s.model
    env = s.env

    nq = model.nq
    nc = model.nc
    nb = nc * friction_dim(env)

	nz = num_var(model, env)
    nθ = num_data(model)
    nx = nq
    ny = 2nc + nb

    # Terms
    off = 0

    idyn, irst, ibil, ialt = linearization_term_index(model, env)
    iθ = collect(1:nθ)

    # Fill the matrix blocks rθ0s
	rθdyn0 = SMatrix{nx,nθ,T,nx*nθ}(rθ0[idyn,:])
	rθrst0 = SMatrix{ny,nθ,T,ny*nθ}(rθ0[irst,:])
	rθbil0 = SMatrix{ny,nθ,T,ny*nθ}(rθ0[ibil,:])

    return RθLin{T,nx,ny,nθ,nx*nθ,ny*nθ}(
        Matrix(rθ0),
		rθdyn0,
		rθrst0,
		rθbil0,
        SVector{nθ,Int}(iθ),
        SVector{nx,Int}(idyn),
        SVector{ny,Int}(irst),
        SVector{ny,Int}(ibil),
        )
end

"""
	Update the residual r.
"""
function rlin!(r::RLin{T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}, z::Vector{T}, θ::Vector{T}, κ) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}
    r.x  = z[r.ix]
    r.y1 = z[r.iy1]
    r.y2 = z[r.iy2]
    r.θ  = θ[r.iθ]
    r.rdyn = r.rdyn0 + r.Dx*(r.x - r.x0) + r.Dy1*(r.y1 - r.y10)                         + r.rθdyn*(r.θ - r.θ0)
    r.rrst = r.rrst0 + r.Rx*(r.x - r.x0) + r.Ry1*(r.y1 - r.y10) + r.Ry2.*(r.y2 - r.y20) + r.rθrst*(r.θ - r.θ0) + SVector{ny}([r.alt; r.alt_zeros])
    r.rbil = r.y1 .* r.y2 .- κ[1]
    return nothing
end

# """
# 	Update the Jacobian rz, and update its Schur complement factorization.
# """
function rzlin!(rz::RZLin{T,nx,ny,nxx,nxy,nyy}, z::Vector{T}, θ::Vector{T}; reg = 0.0) where {T,nx,ny,nxx,nxy,nyy}
    # Unpack
    iy1 = rz.iy1
    iy2 = rz.iy2
    # id = rz.id
    Ry1 = rz.Ry1
    Ry2 = rz.Ry2

    # Update the matrix
    rz.y1 = z[iy1]
    rz.y2 = z[iy2]

	y1_reg = max.(rz.y1, reg)
	y2_reg = max.(rz.y2, reg)

    # update D in Schur complement
    rz.D = rz.Ry1 - Diagonal(rz.Ry2 .* y2_reg ./ y1_reg)

    # update Schur complement
    schur_factorize!(rz.S, rz.D)
    return nothing
end

function residual_violation(ip::InteriorPoint, r::RLin{T}; nquat::Int=0) where {T}
    e1 = norm(r.rdyn, Inf)
    e2 = norm(r.rrst, Inf)
    max(e1, e2)
end

function bilinear_violation(ip::InteriorPoint, r::RLin{T}; nquat::Int=0) where {T}
    norm(r.rbil, Inf)
end

function general_correction_term!(r::RLin{T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}, Δ::AbstractVector{T}, ortr::Vector{Int},
    socr::Vector{Int}, socri::Vector{Vector{Int}}, ortΔ::Vector{Vector{Int}}, socΔ::Vector{Vector{Vector{Int}}}) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}
    Δo1 = @views Δ[ortΔ[1]]
    Δo2 = @views Δ[ortΔ[2]]
    r.rbil_cache .= Δo1 .* Δo2
    r.rbil += SVector{ny}(r.rbil_cache)
    return nothing
end

"""
	Computes the search direction via a linear solve. The factorization has been done when updating rz,
	so this step should be very fast.
"""
function linear_solve!(Δ::Vector{T}, rz::RZLin{T,nx,ny,nxx,nxy,nyy},
        r::RLin{T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}; reg::T = 0.0) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}
    # unpack
    rdyn = r.rdyn
    rrst = r.rrst
    rbil = r.rbil
    Ry2 = rz.Ry2
    y1 = rz.y1
    y2 = rz.y2

	y1_reg = max.(reg, y1)
	y2_reg = max.(reg, y2)

    u = rdyn
    v = rrst - Ry2 .*rbil ./ y1_reg
    schur_solve!(rz.S, u, v)
    Δ[rz.ix]  = rz.S.x
    Δ[rz.iy1] = rz.S.y
    Δ[rz.iy2] = (rbil .- y2_reg .* Δ[rz.iy1]) ./ y1_reg
    return nothing
end

"""
	Computes the Jacobian of the solution with respect to the parameters θ
	via a linear solve. The factorization has been done when updating rz,
	so this step should be ~ fast.
"""
function linear_solve!(δz::Matrix{T}, rz::RZLin{T,nx,ny,nxx,nxy,nyy},
	rθ::RθLin{T,nx,ny,nθ,nxθ,nyθ}; reg::T = 0.0) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ}
    # unpack
	rθdyn = rθ.rθdyn0
	rθrst = rθ.rθrst0
	rθbil = rθ.rθbil0
	Ry2 = rz.Ry2
	y1  = rz.y1
	y2  = rz.y2
	ix  = rz.ix
	iy1 = rz.iy1
	iy2 = rz.iy2

	y1_reg = max.(reg, y1)
	y2_reg = max.(reg, y2)

	for i in eachindex(1:nθ)
		# We remark that rθbil is empty by construction of rθ.
		u = rθdyn[:,i]
		# v = rθrst[:,i] - Ry2 .* rθbil[:,i] ./ y1
		v = rθrst[:,i]
		schur_solve!(rz.S, u, v)
		@. δz[ix,i]  .= rz.S.x
		@. δz[iy1,i] .= rz.S.y
		δz[iy2,i] .= (rθbil[:,i] .- y2_reg .* δz[iy1,i]) ./ y1_reg
		# @. δz[iy2,i] .= .- y2 .* δz[iy1,i] ./ y1
	end
    return nothing
end

function linear_solve!(solver::EmptySolver, δz::Matrix{T}, rz::RZLin{T,nx,ny,nxx,nxy,nyy},
	rθ::RθLin{T,nx,ny,nθ,nxθ,nyθ}; reg::T = 0.0, fact::Bool = true) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ}
	# fact is useless here, since the factorization only happens once,
	# but necessary to be consistent with lu_solver
	linear_solve!(δz, rz, rθ, reg = reg)
	return nothing
end

function linear_solve!(solver::EmptySolver, Δ::Vector{T}, rz::RZLin{T,nx,ny,nxx,nxy,nyy},
        r::RLin{T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}; reg::T = 0.0, fact::Bool = true) where {T,nx,ny,nθ,nxx,nxy,nyy,nxθ,nyθ,nc,nn}
	# fact is useless here, since the factorization only happens once,
	# but necessary to be consistent with lu_solver
	linear_solve!(Δ, rz, r, reg = reg)
	return nothing
end

function update!(r::RLin{T}, z0::Vector{T}, θ0::Vector{T},
    r0::Vector{T}, rz0::Matrix{T}, rθ0::Matrix{T}) where {T}
    idyn = r.idyn
    irst = r.irst
    ibil = r.ibil
    iθ = r.iθ
    ix = r.ix
    iy1 = r.iy1
    iy2 = r.iy2

    # Reference residual
    r.rdyn0 = r0[idyn]
    r.rrst0 = r0[irst]
    r.rbil0 = r0[ibil]

    # Reference residual jacobian rz0
    r.Dx  = rz0[idyn, ix]
    r.Dy1 = rz0[idyn, iy1]
    r.Rx  = rz0[irst, ix]
    r.Ry1 = rz0[irst, iy1]
    r.Ry2 = diag(rz0[irst, iy2])

    # Reference residual jacobian rθ0
    r.rθdyn = rθ0[idyn, iθ]
    r.rθrst = rθ0[irst, iθ]
    r.rθbil = rθ0[ibil, iθ]

    # # Reference z0 and θ0
    r.x0  = z0[ix]
    r.y10 = z0[iy1]
    r.y20 = z0[iy2]
    r.θ0 = θ0[iθ]
    return nothing
end

function update!(rz::RZLin{T}, rz0::Matrix{T}) where T
	idyn = rz.idyn
	irst = rz.irst
	ibil = rz.ibil
	ix = rz.ix
	iy1 = rz.iy1
	iy2 = rz.iy2

    # Fill the matrix blocks rz0s
	rz.Dx  = rz0[idyn, ix]
	rz.Dy1 = rz0[idyn, iy1]
	rz.Rx  = rz0[irst, ix]
	rz.Ry1 = rz0[irst, iy1]
	rz.Ry2 = diag(rz0[irst, iy2])
	rz.y1  = diag(rz0[ibil, iy2])
	rz.y2  = diag(rz0[ibil, iy1])

	# Schur complement
	rz.D = rz.Ry1 - Diagonal(rz.Ry2 .* rz.y2 ./ rz.y1)
	# M = [rz.Dx rz.Dy1;
	#      rz.Rx rz.D  ]
 	# nx = length(ix)
	# ny = length(iy1)
	# rz.S = Schur(M; n=nx, m=ny)
    update!(rz.S, rz.Dx, rz.Dy1, rz.Rx, rz.D)
	return nothing
end

function update!(rθ::RθLin{T}, rθ0::Matrix{T}) where T
	rθ.rθdyn0 = rθ0[rθ.idyn, rθ.iθ]
	rθ.rθrst0 = rθ0[rθ.irst, rθ.iθ]
	rθ.rθbil0 = rθ0[rθ.ibil, rθ.iθ]
	return nothing
end

function rz!(ip::InteriorPoint, rz::RZLin{T}, z::AbstractVector{T},
    θ::AbstractVector{T}; reg::T = 0.0) where {T}
    rzlin!(rz, z, θ, reg = reg)
    return nothing
end

# function rz!(rz::RZLin{T}, z::AbstractVector{T},
# 		θ::AbstractVector{T}; reg::T = 0.0) where {T}
# 	rz!(rz, z, reg = reg)
# 	return nothing
# end

function rθ!(ip::InteriorPoint, rθ::RθLin{T}, z::AbstractVector{T},
    θ::AbstractVector{T}) where {T}
    return nothing
end

function rθ!(rθ::RθLin{T}, z::AbstractVector{T},
    θ::AbstractVector{T}) where {T}
    return nothing
end


text:
function link_epcdata(env,project,activeTiles)

    # Link EPC data
    DataSource = "EPC"
    selected_cols = [:OSG_UPRN,:Dwelling_Type,:Property_Type,:Historical,
        :Total_floor_area_m2_,:Prop_Summ_Total_Floor_Area,:Top_Storey,
        :Habitable_Room_Count,:Heated_Room_Count,
        :Wall_Summary,:Roof_Summary,:Floor_Summary,:Window_Summary,
        :Main_Heating_Summary,:Main_Heating_Controls_Summary,:Secondary_Heating_Summary,
        :Hot_Water_Summary,:Lighting_Summary,:Air_Tightness_Summary,
        :Part_1_Construction_Age_Band,:Multiple_Glazing_Type,
        ]
    project = ExtractData(project,DataSource,activeTiles,selected_cols;method="uprn")

    # Remove historical records (for dwellings with more than one entry)
    project.dat["EPC"] = project.dat["EPC"][project.dat["EPC"].Historical.==false,:]
    epcdata = project.dat["EPC"]

    # Find indices in EPC register that correspond to target UPRNs
    link_uprns(uprns,uprn) = findfirst(uprns.==uprn)
    function avoidmissings(vec)
        vec1 = fill!(Array{Int64}(undef,length(vec)),0)
        i_nonmissing = ismissing.(vec).==false
        vec1[i_nonmissing] .= vec[i_nonmissing]
        return vec1
    end
    epc_uprns = avoidmissings(project.dat["EPC"].OSG_UPRN)
    i_linked_uprns = link_uprns.((epc_uprns,),project.dat["master"].UPRN)

    # Find indices (in master data) which have an EPC
    i_epc = findall(isnothing.(i_linked_uprns).==false)

    # Add 'EPC_' suffix to all EPC column headings, except UPRN (i=1)
    colnames = names(project.dat["EPC"])
    suffix(str) = "EPC_$str"
    newcolnames = suffix.(colnames)

    # Create new EPC cols to append master stock data
    L = length(i_linked_uprns)
    function newcol(project,i_linked_uprns,i_epc,L,colname,newcolname)
        df = DataFrame(Symbol(newcolname) => fill!(Array{Union{Missing,Any}}(undef,L),missing))
        df[!,Symbol(newcolname)][i_epc] .= project.dat["EPC"][!,colname][i_linked_uprns[i_epc]]
        return df
    end
    newcols = newcol.((project,),(i_linked_uprns,),(i_epc,),(L,),colnames,newcolnames)
    newcols = hcat(newcols...)

    # Combine master and EPC data
    project.dat["master"] = hcat(project.dat["master"],newcols)

    # Area ratio, between EPC data and OS-evaluated. 1.0 assumed where EPC missing
    arearatio = project.dat["master"].EPC_Total_floor_area_m2_./project.dat["master"].OS_area_eval
    project.dat["master"][!,:arearatio] = ones(nrow(project.dat["master"]))
    project.dat["master"].arearatio[ismissing.(arearatio).==false] = arearatio[ismissing.(arearatio).==false]

    return project
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f1068fca-0331-11eb-20de-a98c65d5c2dc
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	using Optim, Plots
    using Distributions
	using GeneralQP
	using JuMP, Ipopt
end

# ╔═╡ b4764948-0330-11eb-3669-974d75ab1134
md"""
# Optimisation techniques
"""

# ╔═╡ 1b769f0c-0332-11eb-1efb-178c1985f3df
md"""
We will be using [`Optim`](https://julianlsolvers.github.io/Optim.jl/stable/) for several applications, both uni- and multivariate optimization.
"""

# ╔═╡ 165f35b0-0332-11eb-12e7-f7939d389e58
md"""
## Optimizing a function without gradient information
### Straightforward example
For a univariation function, you need to provide an upper and lower bound
```Julia
optimize(f, lower, upper)
```
Try to optimize $x^3 - 6x + x^2 +2$
* Compare the result between both methods (`Brent()` vs `GoldenSection()`).
* Intermediate results can be store by using the `store_trace=true` keyword argument. The type of the returned object is `Optim.UnivariateOptimizationResults`. The numerical value of each entry can be read by using `.value`.
* Illustrate the evolution of the estimated optimum.
"""

# ╔═╡ 77d89aa8-07d2-11eb-1bbd-a5c896e3ecfe
begin
	f(x) =  x .^3 - 6 * x + x.^2 + 2
	res_brent = optimize(f, 0, 10, Optim.Brent(), store_trace=true)
	res_golden = optimize(f, 0, 10, Optim.GoldenSection(), store_trace=true)
	res_brent.trace
end

# ╔═╡ f27797e5-b83f-4375-ab17-480c09fd7b7f
let
	# generating the illustration
	p1 = plot([v.iteration for v in res_brent.trace], [v.value for v in res_brent.trace], label="brent method", marker=:circle)
	plot!([v.iteration for v in res_golden.trace], [v.value for v in res_golden.trace], label="golden section method", marker=:circle, markeralpha=0.5)
	title!("function value")
	p2 = plot([v.iteration for v in res_brent.trace], [v.metadata["minimizer"] for v in res_brent.trace], label="brent method", marker=:circle)
	plot!([v.iteration for v in res_golden.trace], [v.metadata["minimizer"] for v in res_golden.trace], label="golden section method", marker=:circle, markeralpha=0.5)
	title!("minimizer")

	plot(p1,p2)
	xlabel!("Iterations")
end

# ╔═╡ 7bce2500-0332-11eb-2b63-87dc0d713825

# ╔═╡ e6294d6a-0334-11eb-3829-51ee2b8cadaf
md"""
### Data fitting
Suppose we have a random periodic signal with noise, i.e. $y_i = a \sin(x_i-b) + c + \epsilon_i$ and we wish to determine the values of $a,b$ and $c$ that minimize the difference between the “model” and the measured points.

* Define an error function
* Determine possible values for $a,b$ and $c$
* What is the effect of the initial values? Make an illustration of the error for different starting values.
"""

# ╔═╡ 66be5114-0335-11eb-01a9-c594b92937bf
# data generation (simulated sample)
begin
	# actual values
	a = 3; b = pi/3; c=10;
	# noise distribution
	e=Normal(0,0.1)
	# function
	F(x;a=a,b=b,c=c) = a*sin.(x .- b) .+ c
	# sample length
	n = 10;
	# domain
	xmin = 0; xmax = 20
	d=Uniform(xmin, xmax)
	# random data points (independent variable)
	x = sort(rand(d, n))
	# dependent variable with noise
	y = F.(x) .+ rand(e,n);
end

# ╔═╡ 15ef34dc-0336-11eb-0de5-b942e8871db8
# illustration
begin
	settings = Dict(:xlabel=>"x",:ylabel=>"y",:title=>"a sin(x - b) + c")
	X = range(xmin, xmax, length=50)
	scatter(x,y, label="sample")
	plot!(X,F.(X), label="ground truth"; settings...)
end

# ╔═╡ e949078a-07d3-11eb-0146-8115e335b2e9
begin
	"""
		errfun(v, x=x, y=y)

	où v = [a;b;c ] pour minimiser F(x;a=a,b=b,c=c)
	"""
	function errfun(v, x=x, y=y)
		a,b,c = v
		return sum((F.(x,a=a;b=b,c=c) .- y ) .^2)
	end

	res = optimize(errfun, Float64[1;1;1])
	a_opt, b_opt, c_opt = res.minimizer
	plot(X,F.(X), label="ground truth"; settings...)
	plot!(X, F.(X,a=a_opt, b=b_opt, c=c_opt), label="optimised"; settings...)
end

# ╔═╡ 901eed25-2aea-42a4-bd7e-68b919a8766c
res

# ╔═╡ 71d2bf30-0336-11eb-28ed-95518b9204a7

# ╔═╡ d0a304f2-0336-11eb-0d20-6de833d964b3

# ╔═╡ add5faba-03b8-11eb-0cc7-15f19eb1e0e2
md"""
## Optimisation with gradient information
Suppose we want to minimize a function ``\mathbb{R}^3 \mapsto \mathbb{R}``:

``\min g(\bar{x}) = x_1 ^2 + 2.5\sin(x_2) - x_1^2x_2^2x_3^2 ``

Compare the results (computation time) using
1. a zero order method (i.e. no gradients used)
2. the function and its gradient (both newton and bfgs method)
3. the function, its gradient and the hessian

You can evaluate the performance using the `@time` macro. For a more detailed and representative analysis, you can use the package [`BenchmarkTools`](https://github.com/JuliaCI/BenchmarkTools.jl) (we will go into detail in the session about performance)
"""

# ╔═╡ 9396ccae-03ba-11eb-27a9-1b88ee4fb45f
begin
	g(x) = x[1]^2 + 2.5*sin(x[2]) - x[1]^2*x[2]^2*x[3]^2
	initial_x = [-0.6;-1.2; 0.135];
end

# ╔═╡ 94abb7a2-03bb-11eb-1ceb-1dff8aa3cce7
begin
	# gradients
	function dg!(G, x)
		G[1] = 2*x[1] - 2*x[1]*x[2]^2*x[3]^2
		G[2] = 2.5*cos(x[2]) - 2*x[1]^2*x[2]*x[3]^2
		G[3] = -2*x[1]^2*x[2]^2*x[3]
	end

	function h!(H,x)
		H[1,1] = 2 - 2*x[2]^2*x[3]^2
		H[1,2] = -4*x[1]*x[2]*x[3]^2
		H[1,3] = -4*x[1]*x[2]^2*x[3]
		H[2,1] = -4*x[1]*x[2]*x[3]^2
		H[2,2] = -2.5*sin(x[2]) - 2*x[1]^2*x[3]^2
		H[2,3] = -4*x[1]^2*x[2]*x[3]
		H[3,1] = -4*x[1]*x[2]^2*x[3]
		H[3,2] = -4*x[1]^2*x[2]*x[3]
		H[3,3] = -2*x[1]^2*x[2]^2
	end
end

# ╔═╡ 6f552aa6-07d5-11eb-18a5-b32ed233c403
begin
	println("start method comparison")
	for _ in 1:2
		@time optimize(g, initial_x)
		@time optimize(g, dg!, initial_x, Newton())
		@time optimize(g, dg!, initial_x, BFGS())
		@time optimize(g, dg!, h!, initial_x)
	end
	println("finished method comparison")
end

# ╔═╡ 966b88dc-03bc-11eb-15a4-b5492ddf4ede
md"""
## Optimize the optimizer
You could study the influence of the optimization methods and try to optimize them as well (this is sometimes refered to as hyperparameter tuning).

Try to create a method that minimizes the amount of iterations by modifying the parameter $\eta$ from the `BFGS` method.

**Note:**
* Look at the documentation for possible values of $\eta$.
* This is merely as a proof of concept and will not come up with a significant improvement for this case.

"""

# ╔═╡ 68263aea-07d0-11eb-24f8-6383a3a1e09d
begin
	function optimme(η)
		res = optimize(g, dg!, initial_x,ConjugateGradient(eta=η))
		return res.iterations
	end

	optimize(optimme, 0, 20)
end

# ╔═╡ d5feb956-058b-11eb-1a7d-eb56768f00b1

# ╔═╡ ec264b44-03c2-11eb-1695-cbf638f8cea9
md"""
## Sequential Quadratic Programming
**Reminder**:  typical problem layout:

```math
\begin{align}
\min_{\vec{x}}\, & f\left(\vec{x}\right)\, \textrm{subject to}\, \begin{cases}
\vec{h}\left(\vec{x}\right)=\vec{0}\, \\
\vec{g}\left(\vec{x}\right)\leq\vec{0}
\end{cases}
\end{align}
```

where  ``f:\mathbb{R}^{n}\rightarrow\mathbb{R}``, ``\vec{h}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}``,
``m\leq n, \text{and } \vec{g}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}``.

This allows you to deal with non-linear constraints.

Try to solve the following problem:
```math
\begin{align}
\min_{\bar{x}} f(\bar{x}) = -x_1 - x_2 \\
\text{ST:} \begin{cases}-x_1^2 + x_2 \ge 0\\1- x_1^2 - x_2^2 \ge 0\end{cases}
\end{align}
```

1. Solve the problem with [NLopt](https://github.com/JuliaOpt/NLopt.jl) (placed in a seperate notebook)
2. Solve the problem with [JuMP](https://github.com/jump-dev/JuMP.jl) (combined with [Ipopt](https://github.com/jump-dev/Ipopt.jl))
"""

# ╔═╡ 02aedfec-d075-425b-8db6-e60d5d8a79c3
model

# ╔═╡ 053bae8a-087d-11eb-2e8c-73c41fb4e005
let
	model = Model(Ipopt.Optimizer)
	@variable(model, x[1:2])
	@objective(model, Min, -x[1] - x[2])
	@constraint(model, - x[1] ^2 + x[2] >= 0)
	@constraint(model, 1 - x[1] ^2 - x[2] ^2 >= 0)
	optimize!(model)
	println(termination_status(model))
	println("minimum: $(objective_value(model))")
	value.(x)
end

# ╔═╡ b5cf333b-5ffa-45d5-b9c0-00abc4b63196

# ╔═╡ Cell order:
# ╟─b4764948-0330-11eb-3669-974d75ab1134
# ╟─1b769f0c-0332-11eb-1efb-178c1985f3df
# ╠═f1068fca-0331-11eb-20de-a98c65d5c2dc
# ╟─165f35b0-0332-11eb-12e7-f7939d389e58
# ╠═77d89aa8-07d2-11eb-1bbd-a5c896e3ecfe
# ╠═f27797e5-b83f-4375-ab17-480c09fd7b7f
# ╟─7bce2500-0332-11eb-2b63-87dc0d713825
# ╟─e6294d6a-0334-11eb-3829-51ee2b8cadaf
# ╠═66be5114-0335-11eb-01a9-c594b92937bf
# ╠═15ef34dc-0336-11eb-0de5-b942e8871db8
# ╠═e949078a-07d3-11eb-0146-8115e335b2e9
# ╠═901eed25-2aea-42a4-bd7e-68b919a8766c
# ╟─71d2bf30-0336-11eb-28ed-95518b9204a7
# ╟─d0a304f2-0336-11eb-0d20-6de833d964b3
# ╟─add5faba-03b8-11eb-0cc7-15f19eb1e0e2
# ╠═9396ccae-03ba-11eb-27a9-1b88ee4fb45f
# ╠═94abb7a2-03bb-11eb-1ceb-1dff8aa3cce7
# ╠═6f552aa6-07d5-11eb-18a5-b32ed233c403
# ╟─966b88dc-03bc-11eb-15a4-b5492ddf4ede
# ╠═68263aea-07d0-11eb-24f8-6383a3a1e09d
# ╠═d5feb956-058b-11eb-1a7d-eb56768f00b1
# ╟─ec264b44-03c2-11eb-1695-cbf638f8cea9
# ╠═02aedfec-d075-425b-8db6-e60d5d8a79c3
# ╠═053bae8a-087d-11eb-2e8c-73c41fb4e005
# ╠═b5cf333b-5ffa-45d5-b9c0-00abc4b63196

max_stars_repo_path: Exercises/PS04 - Optimisation techniques I.jl


text:
using ..ChemistryFeaturization.AbstractType: AbstractCodec
import ..ChemistryFeaturization.encode
import ..ChemistryFeaturization.decode
using Flux: onecold

export encode, decode, output_shape

"""
    OneHotOneCold(categorical, bins)

AbstractCodec type which uses a dummy variable (as defined in statistical literature), i.e., which employs
one-hot encoding and a one-cold decoding scheme.
"""
struct OneHotOneCold <: AbstractCodec
    categorical::Bool
    bins::Vector
end

output_shape(ohoc::OneHotOneCold) =
    ohoc.categorical ? length(ohoc.bins) : length(ohoc.bins) - 1

"A flexible version of Flux.onehot that can handle both categorical and continuous-valued encoding."
function encode(ohoc::OneHotOneCold, val)
    local bin_index, onehot_vec
    if ohoc.categorical
        onehot_vec = [0.0 for i = 1:length(ohoc.bins)]
        bin_index = findfirst(isequal(val), ohoc.bins)
    else
        @assert eltype(ohoc.bins) <: Number "Your bins aren't numbers...are you sure you didn't mean for this feature to be categorical?"
        @assert ohoc.bins[1] <= val <= ohoc.bins[end] "The value $val is outside the range of bins $ohoc.bins"
        onehot_vec = [0.0 for i = 1:(length(ohoc.bins)-1)]
        bin_index = searchsorted(ohoc.bins, val).stop
        if bin_index == length(ohoc.bins) # got the max value
            bin_index = bin_index - 1
        elseif isapprox(val, ohoc.bins[1]) # sometimes will get 0 if this doesn't get checked
            bin_index = 1
        end
    end
    onehot_vec[bin_index] = 1.0
    return onehot_vec
end

encode(ohoc::OneHotOneCold, ::Missing) = Vector{Missing}(missing, output_shape(ohoc))

encode(ohoc::OneHotOneCold, vals::Vector) = encode.(Ref(ohoc), vals)

function encode(ohoc::OneHotOneCold, vals::Array)
    # output is one dimension larger
    output = Array{Union{Bool,Missing}}(missing, size(vals)..., output_shape(ohoc))
    for ind in CartesianIndices(vals)
        output[ind, :] = encode(ohoc, vals[ind])
    end
    return output
end

function decode(ohoc::OneHotOneCold, encoded::Vector)
    @assert length(encoded) == output_shape(ohoc)
    local decoded
    if ohoc.categorical # return value
        decoded = onecold(encoded, ohoc.bins)
    else # return bounds (TODO: should this be a tuple or a vector..? I like tuple for distinguishing from encoded vectors, but it doesn't play so nice with broadcasting...)
        decoded = (onecold(encoded, ohoc.bins[1:end-1]), onecold(encoded, ohoc.bins[2:end]))
    end
    return decoded
end

# this is separate from the Array case because at the moment there are different conventions of index ordering for ndims == 2 vs. ndims > 2...oops, we should probably change that, but it will require changes in AtomicGraphNets as well
function decode(ohoc::OneHotOneCold, encoded::Matrix)
    @assert size(encoded)[1] == output_shape(ohoc)
    local decoded
    decoded_length = size(encoded)[2]
    decoded_eltype = eltype(ohoc.bins)
    if ohoc.categorical
        decoded = Vector{Union{decoded_eltype,Missing}}(missing, decoded_length)
    else
        decoded =
            Vector{Union{Tuple{decoded_eltype,decoded_eltype},Missing}}(missing, decoded_length)
    end
    for i = 1:decoded_length
        vec = encoded[:, i] # this is the key difference from the Array dispatch below
        if !all(ismissing.(vec))
            decoded[i] = decode(ohoc, vec)
        end
    end
    return decoded
end

function decode(ohoc::OneHotOneCold, encoded::Array)
    @assert size(encoded)[end] == output_shape(ohoc)
    local decoded
    decoded_size = size(encoded)[1:end-1]
    decoded_eltype = eltype(ohoc.bins)
    if ohoc.categorical
        decoded = Array{Union{decoded_eltype,Missing}}(missing, decoded_size...)
    else
        decoded = Array{Union{Tuple{decoded_eltype,decoded_eltype},Missing}}(missing, decoded_size...)
    end
    for ind in CartesianIndices(decoded)
        vec = encoded[ind, :]
        if !all(ismissing.(vec))
            decoded[ind] = decode(ohoc, vec)
        end
    end
    return decoded
end


text:
function con_parse(x::AbstractString)
  word = "[[:alnum:]_\\.\\[\\]\"]+"
  pattern_lgl_left = "($word)( *)(?=[=<>≤≥\\!])" # logical operation
  pattern_lgl_right = "(?<=[=<>≤≥])( *)($word)" # logical operation
  pattern_op = r"([=<>≤≥\!]+)"

  k = str_extract_strip(x, pattern_lgl_left)
  v = str_extract_strip(x, pattern_lgl_right)
  op = str_extract_strip(x, pattern_op)
  k, v, op
end

function _con_dt_transform(dname::AbstractString, con::AbstractString)
  k, v, op = con_parse(con)
  k = "$dname.$k"
  op = ".$op"
  "($k $op $v)"
end

function _con_dt_transform(dname::AbstractString, con::Vector{<:AbstractString})
  map(x -> _con_dt_transform(dname, x), con)
end

function con_dt_transform(con::AbstractString; dname="d")
  cons_and = split(con, "&")
  cons_or = cons_and |> x -> split.(x, "|")

  con2 = map(con_or -> begin
      _con_dt_transform(dname, con_or) |> x -> join(x, " .| ")
    end, cons_or) |> x -> join(x, " .& ")
  "$dname[$con2, :]"
end

export con_split, con_parse, con_dt_transform,
  paste


text:
# For better tree printing and navigation
## Things we need to define
function AbstractTrees.children(node::TrapezoidalSearchNode)
    if node.left === nothing
        if node.right === nothing
            return ()
        else
            return (node.right, )
        end
    else
        if node.right === nothing
            return (node.left, )
        else
            return (node.left, node.right)
        end
    end
end

## Things that make printing prettier
AbstractTrees.printnode(io::IO, node::TrapezoidalSearchNode) = print(io, node)

function buildmega(table)
    s = 0
    for x in table.geometry
        s += length(x.points) - length(x.parts)
    end
    res = Vector{EData{typeof(table.geometry[1].points[1].x)}}(undef, s)
    offset = 0
    for (i, poly) in pairs(table.geometry)
        buildedges!(res, offset, poly, i)
        offset += length(poly.points) - length(poly.parts)
    end

    return res
end

function buildedges!(edges, offset, poly, id)
    points = poly.points
    parts = poly.parts

    @inbounds for k in eachindex(parts)
        i1 = parts[k] + 2
        i2 = k == length(parts) ? length(points) : parts[k + 1]
        prev = V(points[i1 - 1].x, points[i1 - 1].y)
        for i in i1:i2
            p = points[i]
            cur = V(p.x, p.y)
            if prev < cur
                edges[i - k + offset] = EData(E(prev, cur), -1, id)
            else
                edges[i - k + offset] = EData(E(cur, prev), id, -1)
            end
            prev = cur
        end
    end

    return edges
end

function buildedges(poly, id)
    points = poly.points
    parts = poly.parts
    res = Vector{EData{typeof(points[1].x)}}(undef, length(points) - length(parts))
    @inbounds for k in eachindex(parts)
        i1 = parts[k] + 2
        i2 = k == length(parts) ? length(points) : parts[k + 1]
        prev = V(points[i1 - 1].x, points[i1 - 1].y)
        for i in i1:i2
            p = points[i]
            cur = V(p.x, p.y)
            if prev < cur
                res[i - k] = EData(E(prev, cur), -1, id)
            else
                res[i - k] = EData(E(cur, prev), id, -1)
            end
            prev = cur
        end
    end

    return res
end

function names2indices(baddies, table)
    m = Dict(map(x -> x[2] => x[1], enumerate(table.tzid)))
    res = Tuple{Int, Int}[]
    for x in baddies
        push!(res, (m[x[1]], m[x[2]]))
    end
    return res
end

function mergeedges(edges, table; precedence = PRECEDENCE)
    precedence = names2indices(precedence, table)
    baddies = Set{Tuple{Int, Int}}()
    edges = sort(edges, by = x -> edge(x))
    res = similar(edges, 0)
    prev = edges[1]
    for i in 2:length(edges)
        cur = edges[i]
        if edge(cur) != edge(prev)
            push!(res, prev)
            prev = cur
        else
            if lhs(cur) != -1 && lhs(prev) == -1
                @set! prev.l = cur.l
            elseif cur.r != -1 && prev.r == -1
                @set! prev.r = cur.r
            else
                caught = false
                for prec in precedence
                    if (cur.l == prec[1] && prev.l == prec[2]) || (cur.l == prec[2] && prev.l == prec[1])
                        @set! prev.l = prec[1]
                        caught = true
                    elseif (cur.r == prec[1] && prev.r == prec[2]) || (cur.r == prec[2] && prev.r == prec[1])
                        @set! prev.r = prec[1]
                        caught = true
                    end
                    caught && break
                end
                if !caught
                    if cur.l != -1 && prev.l != -1
                        m1 = min(cur.l, prev.l)
                        m2 = max(cur.l, prev.l)
                        push!(baddies, (m1, m2))
                    elseif cur.r != -1 && prev.r != -1
                        m1 = min(cur.r, prev.r)
                        m2 = max(cur.r, prev.r)
                        push!(baddies, (m1, m2))
                    end
                end
            end
        end
    end
    push!(res, prev)

    return res, baddies
end

function boundingbox(poly)
    points = poly.points
    p, rest = Iterators.peel(points)
    p = V(p.x, p.y)
    leftp = p
    rightp = p
    topy = p.y
    bottomy = p.y
    for p in rest
        p = V(p.x, p.y)
        if p < leftp
            leftp = p
        end
        if p > rightp
            rightp = p
        end
        if p.y > topy
            topy = p.y
        end
        if p.y < bottomy
            bottomy = p.y
        end
    end
    top = E(V(leftp.x, topy), V(rightp.x, topy))
    bottom = E(V(leftp.x, bottomy), V(rightp.x, bottomy))

    return Trapezoid(leftp, rightp, top, bottom, -1)
end

function boundingbox(v::AbstractVector{T}) where {T <: EData}
    ed, rest = Iterators.peel(v)
    p = edge(ed).p
    leftp = p
    rightp = p
    topy = p.y
    bottomy = p.y
    for ed in rest
        p = edge(ed).p
        if p < leftp
            leftp = p
        end
        if p > rightp
            rightp = p
        end
        if p.y > topy
            topy = p.y
        end
        if p.y < bottomy
            bottomy = p.y
        end
    end
    top = E(V(leftp.x, topy), V(rightp.x, topy))
    bottom = E(V(leftp.x, bottomy), V(rightp.x, bottomy))

    return Trapezoid(leftp, rightp, top, bottom, -1)
end

########################################
# Angle sorting facilities
########################################

function pseudoangle1(v)
    atan(v.y, v.x)
end

function pseudoangle2(v)
    copysign(1. - v.x/(abs(v.x)+abs(v.y)), v.y)
end

# http://steve.hollasch.net/cgindex/math/fowler.html
# https://stackoverflow.com/a/27253590/1494135
function pseudoangle3(v)
    # 1 for above, 0 for below the diagonal/anti-diagonal
    diag  = v.x > v.y
    adiag = v.x > -v.y

    r = !adiag ? 4 : 0

    if (v.y == 0)
        return r
    end

    if (diag ⊻ adiag)
        r += 2 - v.x / v.y
    else
        r += v.y / v.x
    end

    return r
end

function isangleless(v1, v2)
    a1 = abs(v1.x) + abs(v1.y)
    a2 = abs(v2.x) + abs(v2.y)

    a2*copysign(a1 - v1.x, v1.y) < a1*copysign(a2 - v2.x, v2.y)
end


text:
struct TimeSeriesKey <: InfrastructureSystemsType
    time_series_type::Type{<:TimeSeriesMetadata}
    name::String
end

function TimeSeriesKey(; time_series_type::Type{<:TimeSeriesMetadata}, name::String)
    return TimeSeriesKey(time_series_type, name)
end

function TimeSeriesKey(data::TimeSeriesData)
    metadata_type = time_series_data_to_metadata(typeof(data))
    return TimeSeriesKey(metadata_type, get_name(data))
end

const TimeSeriesByType = Dict{TimeSeriesKey, TimeSeriesMetadata}

"""
Time series container for a component.
"""
mutable struct TimeSeriesContainer
    data::TimeSeriesByType
    time_series_storage::Union{Nothing, TimeSeriesStorage}
end

function TimeSeriesContainer()
    return TimeSeriesContainer(TimeSeriesByType(), nothing)
end

Base.length(container::TimeSeriesContainer) = length(container.data)
Base.isempty(container::TimeSeriesContainer) = isempty(container.data)

function set_time_series_storage!(
    container::TimeSeriesContainer,
    storage::Union{Nothing, TimeSeriesStorage},
)
    if !isnothing(container.time_series_storage) && !isnothing(storage)
        throw(
            ArgumentError(
                "The time_series_storage reference is already set. Is this component being " *
                "added to multiple systems?",
            ),
        )
    end

    container.time_series_storage = storage
end

function add_time_series!(
    container::TimeSeriesContainer,
    ts_metadata::T;
    skip_if_present = false,
) where {T <: TimeSeriesMetadata}
    key = TimeSeriesKey(T, get_name(ts_metadata))
    if haskey(container.data, key)
        if skip_if_present
            @warn "time_series $key is already present, skipping overwrite"
        else
            throw(ArgumentError("time_series $key is already stored"))
        end
    else
        container.data[key] = ts_metadata
    end
end

function remove_time_series!(
    container::TimeSeriesContainer,
    ::Type{T},
    name::AbstractString,
) where {T <: TimeSeriesMetadata}
    key = TimeSeriesKey(T, name)
    if !haskey(container.data, key)
        throw(ArgumentError("time_series $key is not stored"))
    end

    pop!(container.data, key)
end

function clear_time_series!(container::TimeSeriesContainer)
    empty!(container.data)
end

function get_time_series(
    ::Type{T},
    container::TimeSeriesContainer,
    name::AbstractString,
) where {T <: TimeSeriesMetadata}
    key = TimeSeriesKey(T, name)
    if !haskey(container.data, key)
        throw(ArgumentError("time_series $key is not stored"))
    end

    return container.data[key]
end

function get_time_series_names(
    ::Type{T},
    container::TimeSeriesContainer,
) where {T <: TimeSeriesMetadata}
    names = Set{String}()
    for key in keys(container.data)
        if key.time_series_type <: T
            push!(names, key.name)
        end
    end

    return Vector{String}(collect(names))
end

function has_time_series_internal(
    container::TimeSeriesContainer,
    ::Type{T},
    name::AbstractString,
) where {T <: TimeSeriesMetadata}
    return haskey(container.data, TimeSeriesKey(T, name))
end

function serialize(container::TimeSeriesContainer)
    # Store a flat array of time series. Deserialization can unwind it.
    return serialize_struct.(values(container.data))
end

function deserialize(::Type{TimeSeriesContainer}, data::Vector)
    container = TimeSeriesContainer()
    for ts_dict in data
        type = get_type_from_serialization_metadata(get_serialization_metadata(ts_dict))
        time_series = deserialize(type, ts_dict)
        add_time_series!(container, time_series)
    end

    return container
end


text:
cd(@__DIR__); include("setups/path512.jl")
pyplot(dpi = 200)

##
@time VM_NGWP = vm_ngwp(𝚽, GP_dual)
@time LP_NGWP = lp_ngwp(𝚽, Gstar_Sig.W, GP_dual; ϵ = 0.3)

dict_VM = Dict()
dict_LP = Dict()
for j in [3, 4, 5]
   for k in [1, 2, 3]
      if GP_dual.rs[k + 1, j] - GP_dual.rs[k, j] < 1
         continue
      end
      WW_VM = sort_wavelets(VM_NGWP[GP_dual.rs[k, j]:(GP_dual.rs[k + 1, j] - 1), j, :]')
      WW_LP = sort_wavelets(LP_NGWP[GP_dual.rs[k, j]:(GP_dual.rs[k + 1, j] - 1), j, :]')
      Y = [-0.5, 0.6]
      l = Int(floor(size(WW_LP, 2) / 2))
      dict_VM[(j-1, k-1, l-1)] = WW_VM[:, l]
      dict_LP[(j-1, k-1, l-1)] = WW_LP[:, l]
      plt = plot(size = (400, 350), layout = Plots.grid(2, 1), ylim = Y)
         plot!(WW_VM[:, l], c = :black, grid = false, frame = :box, lw = 0.8, xtick = false,
         lab = latexstring("\\psi^{($(j-1))}_{$(k-1), $(l-1)} (\\mathrm{VM})"), legendfontsize = 11)
         plot!(WW_LP[:, l], c = :black, grid = false, frame = :box, lw = 0.8,
         lab = latexstring("\\psi^{($(j-1))}_{$(k-1), $(l-1)} (\\mathrm{LP})"), legendfontsize = 11,
         subplot = 2)
         xticks!([1; 64:64:N], vcat(string(1), [string(k) for k in 64:64:N]),
         subplot = 2)
      savefig(plt, "../figs/Path512_VM_LP_j$(j-1)k$(k-1)l$(l-1).png")
   end
end

## Table 8.2
using PrettyTables
header = ["basis vector" "main support width" "sidelobe attenuation"]
res = Matrix{String}(undef, 0, 3)
iter = sort(collect(keys(dict_VM)))
for γ in iter
   global res
   ms = find_mainsupport(dict_VM[γ]; ϵ = 0.01)
   sa = sidelobe_attenuation(dict_VM[γ])
   j, k, l = γ
   res = vcat(res, ["ψ^{$(j)}_{$(k), $(l)} (VM)" ms[2] - ms[1] + 1 round(sa; digits = 4)])
   ms = find_mainsupport(dict_LP[γ]; ϵ = 0.01)
   sa = sidelobe_attenuation(dict_LP[γ])
   res = vcat(res, ["ψ^{$(j)}_{$(k), $(l)} (LP)" ms[2] - ms[1] + 1 round(sa; digits = 4)])
end
pretty_table(res; header = tuple(header), alignment = :c)


text:
########################## conditional distribution ############################
# TODO: add iterative method conditioning
# TODO: specialize conditional for ConditionalKernel GP
# IDEA: generalize gp conditional with arbitrary linear functional
# IDEA: internally normalize data and covariance before inference for stability
# IDEA have separate type for GP to make conditional type dispatch easier?
# nice syntax
Base.:|(G::AbstractGaussian, t::Tuple) = conditional(G, t...)

# if noise covariance is passed separately, first merge it with the observations
function conditional(G::Gaussian, x, y::AbstractVector, Σ::Union{Real, AbstractMatOrFacOrUni})
    conditional(G, x, 𝒩(y, Σ))
end
# converting vector of noise variances to diagonal covariance matrix via constructor
function conditional(G::Gaussian, x::AbstractVector, y::AbstractVector, σ²::AbstractVector)
    conditional(G, x, 𝒩(y, σ²))
end

# conditional with matrix input X for GPs
function conditional(G::Gaussian, X::AbstractMatrix, y::AbstractVector; tol = tol)
	conditional(G, [c for c in eachcol(X)], y, tol = tol)
end
function conditional(G::Gaussian, X::AbstractMatrix, y::Multivariate; tol = tol)
	conditional(G, [c for c in eachcol(X)], y, tol = tol)
end

############################### Univariate #####################################
conditional(G::Univariate, x::Real) = Gaussian(x, zero(0))
conditional(G::Univariate, x::Real, σ²::Real) = conditional(G, 𝒩(x, σ²))
function conditional(G::Univariate, x::Univariate)
    α = var(G) / (var(G) + var(x))
    μ = (1-α) * mean(G) + α * mean(x)
    σ² = (1-α) * var(G)
    Gaussian(μ, σ²)
end
############################### Multivariate ###################################
conditional(G::Multivariate, i::Int, x::Real) = conditional(G, [i], [x])

# pre-processes index array to be bitarray to allow for "not i" indexing
function conditional(G::Multivariate, i::AbstractVector{<:Integer},
                     Gx::Union{AbstractVector, Multivariate})
    length(i) == length(Gx) || error("number of indices and values in x do not match")
    isdata = [any(==(k), i) for k in eachindex(mean(G))] # TODO: compute in quasi-linear time
    conditional(G, isdata, Gx)
end

# conditional distribution, after knowing that G[i] = x + ε where x + ε ∼ Gx
function conditional(G::Multivariate, i::AbstractVector{<:Bool},
                     Gx::Union{AbstractVector, Multivariate})
    conditional!(copy(G), i, Gx)
end

function conditional!(G::Multivariate, i::AbstractVector{<:Bool},
                      Gx::Union{AbstractVector, Multivariate}; tol::Real = tol)
    any(i) || return G # return if there is nothing to condition on
    j = .!i
    x = Gx isa Gaussian ? mean(Gx) : Gx
    Gi, Gj = marginal(G, i), marginal(G, j)
    Gi = factorize(Gi, tol = tol)
    Σii = Gx isa Gaussian ? factorize(Woodbury(cov(Gx), cov(Gi))) : cov(Gi)
    Λi = inverse(Σii) # pseudo-inverse?
    Σi = @view cov(G)[i, :]

    # assumes G.μ, (G.Σ) contains prior mean (covariance)
    G.μ .+= Σi' * (Λi * difference(x, mean(Gi)))
    G.Σ .-= *(Σi', Λi, Σi) # check if we can pre-allocate memory here. woodbury?

    if x isa AbstractVector # collapse of uncertainty of observed noiseless data
        G.μ[i] .= x
        G.Σ[i, :] .= 0; G.Σ[:, i] .= 0
    end
    return G
end

################################ Process #######################################
# conditional distribution of Gaussian Process
# ASSUMPTION: x and y are vectors of data, not data vectors
# FIXME: does not work if G is a conditional kernel (gives rise to Woodbury K, which does not have factorize(K, tol) signature)
function conditional(G::Gaussian, x::AbstractVector,
                     Gy::Union{AbstractVector, Multivariate}; tol::Real = tol)
    length(x) > 0 || return G # return if an empty vector is passed
	if G.Σ isa ConditionalKernel
		throw("iterative conditioning is not yet supported, call conditional with entire data instead.")
		# TODO: if applicable, use cholesky updating formula, need to keep track of all targets y?
	end
    K = gramian(cov(G), x)
    K = length(x) == 1 ? Matrix(K) : factorize(K, tol = tol)
    y = Gy
    if Gy isa Gaussian # if there is no noise variance, gramian has to be full rank
        y, Σ = mean(Gy), cov(Gy)
        K = length(x) == 1 ? (K + Σ) : factorize(Woodbury(Σ, K)) # TODO allow for iterative inference, like in multi.jl: specialize factorize for LazyMatrixSum
    else
        rank(K) == size(K, 1) || throw(RankDeficientException(1))
    end
    μ = ConditionalMean(G, x, y, K)
    Σ = ConditionalKernel(G, x, K)
    return Gaussian(μ, Σ)
end

######################## conditioning on single datapoints #####################
# IDEA: try to consolidate
# process with 1d input and output
function conditional(G::Gaussian, x::Real, y::Real, σ²::Real)
    conditional(G, [x], [y], σ²)
end
conditional(G::Gaussian, x::Real, y::Real) = conditional(G, [x], [y])

# process with nd input and 1d ouput
function conditional(G::Gaussian, x::AbstractVector{<:Real}, y::Real, σ²::Real)
    conditional(G, [x], [y], σ²)
end
function conditional(G::Gaussian, x::AbstractVector{<:Real}, y::Real)
    conditional(G, [x], [y])
end

############################### Conditional Mean ###############################
# here, T would be the output type ..., or should it be field type?
struct ConditionalMean{M, K, U<:AbstractVector, V<:AbstractVector}
    μ::M
	kxu::K
	u::U
    α::V # α = K \ (y - μ(x)) IDEA: sparsity?
end
ConditionalMean(G::Gaussian, u...) = ConditionalMean(G.μ, G.Σ, G.μ, u...)
function ConditionalMean(μx, kxu, μu, u::AbstractVector, y::AbstractVector, K)
    μ = marginal(μu, u)
    α = K \ difference(y, μ)
    ConditionalMean(μx, kxu, u, α)
end
function (M::ConditionalMean)(x)
	value = M.μ(x)
	if length(M.u) < length(M.α) # this means M.kxu will return a vector
		G = gramian(M.kxu, [x], M.u)
		g = vec(Matrix(G)) # because of BlockFactorization
		value += dot(g, M.α)
	else
	    for i in eachindex(M.α)
	        value += M.kxu(x, M.u[i]) * M.α[i]
	    end
	end
	return value
end
# not possible because of BlockFactorization
# IDEA: don't return it if one dimension is one!
# (M::ConditionalMean)(x) = M.μ(x) + dot(vec(gramian(M.kxu, [x], M.u)), M.α)
# function (M::ConditionalMean)(x) # this doesn't work for vector-valued observations
#     value = M.μ(x)
#     for i in eachindex(M.α)
#         value += M.kxu(x, M.u[i]) * M.α[i]
#     end
#     return value
# end
function marginal(M::ConditionalMean, x::VecOfVec)
	μ = marginal(M.μ, x)
	K = gramian(M.kxu, x, M.u)
	return mul!(μ, K, M.α, 1, 1)
end

########################## Conditional Kernel ##################################
# corresponds to the kernel of the conditional GP, with observations at x
# ASSUMPTION: x is vector of data, even if it is just a single element
struct ConditionalKernel{T, K, KX, KY, U<:AbstractVector,
                    	 S<:AbstractMatOrFac{T}} <: MercerKernel{T}
    kxx::K # prior covariance function of predictions
	kxu::KX # covariance function between observations and predictions
	kuy::KY # covariance function between observations and predictions
    u::U # input data on which we are conditioning
    Kuu::S # covariance between u and u
end
# default to having same covariance between observations and predictions
function ConditionalKernel(k, u::AbstractVector, Kuu::AbstractMatOrFac)
	ConditionalKernel(k, k, k, u, Kuu)
end
function ConditionalKernel(k, u::AbstractVector)
	Kuu = factorize(gramian(k, u))
	ConditionalKernel(k, u, Kuu)
end
function ConditionalKernel(G::Gaussian, u::AbstractVector, Kuu::AbstractMatOrFac)
    ConditionalKernel(cov(G), u, Kuu)
end

# might have to force @nospecialize for K, since recursive Conditional usage
# could be stressing the compiler
# TODO: why lazy evaluation here? for multi-output GPs we have VectorValuedGaussian ...
function (C::ConditionalKernel)(x, y)
    Kxu = gramian(C.kxu, [x], C.u) # no need to evaluate lazily because output dimension is one
    Kuy = (x ≡ y) && (C.kxu ≡ C.kuy) ? Kxu' : gramian(C.kuy, C.u, [y])
	Kxu, Kuy = vec(Matrix(Kuy)), vec(Matrix(Kxu)) # vec(Kxu), vec(Kuy) # necessary because BlockFactorization does not allow vec currently
    C.kxx(x, y) - dot(Kxu, inverse(C.Kuu), Kuy) # bottleneck
end

# for efficiency, gramian should dispatch on the type of the factorization of c
# does it make sense to return a lazy matrix? given each element requires O(n^2) cost ...
import CovarianceFunctions: gramian
function gramian(C::ConditionalKernel, x::AbstractVector, y::AbstractVector)
   Kxu = gramian(C.kxu, x, C.u)
   Kuy = (x ≡ y) && (C.kxu ≡ C.kuy) ? Kxu' : gramian(C.kuy, C.u, y)
   Woodbury(gramian(C.kxx, x, y), Kxu, inverse(C.Kuu), Kuy, -1)
end
# IDEA: could specialize ternary * for Woodbury as middle argument to speed up SoR

################### conditioning on linear operator data #######################
# computes the posterior distribution x | A, b given prior on x and b
# IDEA: the posterior precision matrix is simple: (A'A + cov(x)), exploit this?
conditional(x::Multivariate, A::AbstractMatrix, b::Uni) = _conditional(x, A, b)
conditional(x::Multivariate, A::AbstractMatrix, b::Multivariate) = _conditional(x, A, b)
function _conditional(x::Multivariate, A::AbstractMatrix, b::Union{Uni, Multi})
    ACx = A*cov(x)
    y = difference(mean(b), A*mean(x)) # this allocates anyway, so use mul!?
    Cy = Woodbury(cov(b), A, cov(x), A') # covariance matrix of b
    Cy = factorize(Cy)
    μ = mean(x) + ACx' * (Cy \ y)
    Σ = cov(x) - *(ACx', inverse(Cy), ACx) # Woodbury?
    Gaussian(μ, Σ)
end

# conditioning on noiseless observations
# avoiding type-ambiguity with line 19 by creating "hidden" _conditional
conditional(x::Multivariate, A::AbstractMatrix, b::Real) = _conditional(x, A, b)
conditional(x::Multivariate, A::AbstractMatrix, b::AbstractVector) = _conditional(x, A, b)
function _conditional(x::Multivariate, A::AbstractMatrix, b::Union{Real, AbstractVector})
    μ, Σ = copy(b), zero(cov(x))
    if !issquare(A)
        ACx = A*cov(x)
        y = difference(b, A*mean(x))
        Cy = factorize(ACx * A') # scales cubically in observations *(A, cov(x), A')
        μ = mean(x) + ACx' * (Cy \ y)
        Σ = Symmetric(cov(x) - *(ACx', inverse(Cy), ACx))
    end
    Gaussian(μ, Σ)
end

# convenience for linear functionals parameterized by a
# IDEA: try to consolidate
function conditional(x::Multivariate, a::AbstractVector{<:Real}, b::Real)
    conditional(x, a', b)
end
function conditional(x::Multivariate, a::AbstractVector{<:Real}, b::Uni)
    conditional(x, a', b)
end


text:
"A convolutional filter. One would normally fit many such filters as part of a convolutional layer."
immutable ConvFilter{T,N,P<:Params} <: Learnable
    sizein::NTuple{N,Int}
    sizeout::NTuple{N,Int}
    sizefilter::NTuple{N,Int}
    stride::NTuple{N,Int}
    input::SumNode{T,N}
    output::OutputNode{T,N}
    params::P
end
ConvFilter(args...) = ConvFilter(Float64, args...)

function ConvFilter{T,N}(::Type{T}, sizein::NTuple{N,Int}, sizefilter::NTuple{N,Int}, stride::NTuple{N,Int} = (1,1))
    input = InputNode(T, sizein...)
    nr_x,nc_x = sizein
    nr_f,nc_f = sizefilter
    sizeout = (nr_x-nr_f+1, nc_x-nc_f+1)
    output = OutputNode(T, sizeout...)
    nparams = prod(sizefilter) + 1
    params = Params(T, nparams, (sizefilter, (1,)))
    w, b = params.views
    initialize_weights!(w)
    initialize_bias!(b)
    ConvFilter(sizein, sizeout, sizefilter, stride, input, output, params)
end

# non-allocating sum(x .* y)
function summult{T}(x::AbstractMatrix{T}, y::AbstractMatrix{T})
    tot = zero(T)
    nr, nc = size(x)
    @assert size(x) == size(y)
    for c=1:nc, r=1:nr
        tot += x[r,c]*y[r,c]
    end
    tot
end

"allows reading out of bounds (returns `zero(T)`)"
immutable TileView{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    a::A
    pos::NTuple{N,Int}
    sz::NTuple{N,Int}
end
Base.size(tile::TileView) = tile.sz
function Base.getindex{T}(tile::TileView{T,2}, i::Int, j::Int)
    ai = tile.pos[1] + i - 1
    aj = tile.pos[2] + j - 1
    if ai <= 0 || aj <= 0 || ai > size(tile.a,1) || aj > size(tile.a,2)
        zero(T)
    else
        tile.a[ai,aj]
    end
end

# specialized for matrix... TODO: make this a generated function
function transform!{T}(filter::ConvFilter{T,2})
    w, b = filter.params.views
    x = filter.input.val
    y = filter.output.val

    # start with the bias
    fill!(y, b[1])

    # now apply the filter, passing over the input
    nr_y, nc_y = filter.sizeout
    nr_f, nc_f = filter.sizefilter
    for r=1:nr_y, c=1:nc_y
        x_tile = view(x, r:r+nr_f-1, c:c+nc_f-1)
        # x_tile = TileView{T,2,typeof(x)}(x, (r,c), filter.sizefilter)
        y[r, c] = summult(w, x_tile) + b[1]
    end
    y
end

function grad!{T}(filter::ConvFilter{T,2})
    w, b = filter.params.views
    ∇w, ∇b = filter.params.∇_views
    x = filter.input.val
    ∇x = filter.input.∇
    ∇y = filter.output.∇

    # reset the gradients to 0
    fill!(∇x, zero(T))
    fill!(filter.params.∇, zero(T))

    nr_y, nc_y = filter.sizeout
    nr_f, nc_f = filter.sizefilter
    for ry=1:nr_y, cy=1:nc_y, i=1:nr_f, j=1:nc_f
        rx = ry+i-1
        cx = cy+j-1
        ∇x[rx,cx] += w[i,j] * ∇y[ry,cy]
        ∇w[i,j] += x[rx,cx] * ∇y[ry,cy]
    end

    ∇b[1] = sum(∇y)
    return
end

# ----------------------------------------------------------------------------

"Wraps several filters for a convolution layer"
immutable ConvLayer{T,N,P<:Params} <: Learnable
    sizein::NTuple{N,Int}
    sizeout::NTuple{N,Int}
    input::SumNode{T,N}
    output::OutputNode{T,N}
    filters::Vector{ConvFilter{T,N,P}}
    params::P
end

function ConvLayer{T,N}(::Type{T}, sizein::NTuple{N,Int}, sizefilter::NTuple{N,Int}, numfilters::Int)
    filters = [ConvFilter(T, sizein, sizefilter) for i=1:numfilters]
    sizeout = filters[1].sizeout
    input = InputNode(T,sizein)
    output = OutputNode(T,sizeout..., numfilters)
    params = consolidate_params(T, filters)

    # TODO: need a concatenating layer of some sort... hmmm
    # or maybe we need a "link" function that connects a view to part of the array?
    for filter in filters
        link_nodes!(input, filter.input)
        # link_nodes!(output, filter.output)
    end

    ConvLayer(sizein, sizeout, input, output, filters, params)
end

function transform!(conv::ConvLayer)
    foreach(transform!, conv.filters)
    value(output)
end

function grad!(conv::ConvLayer)
    foreach(grad!, conv.filters)
end


text:
using Test
using Overseer

@testset "Indices" begin include("test_indices.jl") end
@testset "Components" begin include("test_components.jl") end
@testset "Ledger" begin include("test_ledger.jl") end