text:
"""

    spring_layout(mpr[; C=2.0, maxiter=100, seed=0])

Calculate a 1-skeleton of the Mapper complex `mpr` layout using spring/repulsion model.

This borrows with modifications from [GraphPlot](https://github.com/JuliaGraphs/GraphPlot.jl)

# Arguments
- `mpr::Mapper`: the Mapper complex.
- `C::Real`: the constant to fiddle with density of resulting layout.
- `maxiter::Integer`: the number of iterations we apply the forces.
- `starttemp::Real`: the initial "temperature", controls movement per iteration.
- `seed::Integer`: the seed for pseudorandom generation of locations (default = 0).
"""
function spring_layout(adj::AbstractMatrix, locs_x, locs_y; C=2.0, maxiter=100, starttemp=2.0) where {T<:Integer}

    N = size(adj, 1)

    # The optimal distance bewteen vertices
    K = C * sqrt(4.0 / N)

    # Store forces and apply at end of iteration all at once
    force_x = zeros(N)
    force_y = zeros(N)

    @inbounds for iter = 1:maxiter
        # Calculate forces
        for i = 1:N
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j = 1:N
                i == j && continue
                d_x = locs_x[j] - locs_x[i]
                d_y = locs_y[j] - locs_y[i]
                d   = sqrt(d_x^2 + d_y^2)
                if adj[i,j] != zero(eltype(adj)) || adj[j,i] != zero(eltype(adj))
                    # F = d^2 / K - K^2 / d
                    F_d = d / K - K^2 / d^2
                else
                    # Just repulsive
                    # F = -K^2 / d^
                    F_d = -K^2 / d^2
                end
                # d  /          sin θ = d_y/d = fy/F
                # F /| dy fy    -> fy = F*d_y/d
                #  / |          cos θ = d_x/d = fx/F
                # /---          -> fx = F*d_x/d
                # dx fx
                force_vec_x += F_d*d_x
                force_vec_y += F_d*d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        TEMP = starttemp / iter
        # Now apply them, but limit to temperature
        for i = 1:N
            force_mag  = sqrt(force_x[i]^2 + force_y[i]^2)
            scale      = min(force_mag, TEMP)/force_mag
            locs_x[i] += force_x[i] * scale
            #locs_x[i]  = max(-1.0, min(locs_x[i], +1.0))
            locs_y[i] += force_y[i] * scale
            #locs_y[i]  = max(-1.0, min(locs_y[i], +1.0))
        end
    end

    # Scale to unit square
    min_x, max_x = minimum(locs_x), maximum(locs_x)
    min_y, max_y = minimum(locs_y), maximum(locs_y)
    function scaler(z, a, b)
        2.0*((z - a)/(b - a)) - 1.0
    end
    map!(z -> scaler(z, min_x, max_x), locs_x, locs_x)
    map!(z -> scaler(z, min_y, max_y), locs_y, locs_y)

    return locs_x, locs_y
end

import Random
function spring_layout(adj::AbstractMatrix; seed::Integer=0, xs=nothing, ys=nothing, kwargs...)
    N = size(adj, 1)
    rng = Random.MersenneTwister(seed)
    spring_layout(adj, xs === nothing ? 2 .* rand(rng, N) .- 1.0 : xs,
                        ys === nothing ? 2 .* rand(rng, N) .- 1.0 : ys ; kwargs...)
end

function spring_layout(mpr::Mapper; seed::Integer=0, random::Bool=false, kwargs...)
    xs = random ? nothing : mpr.centers[1,:]
    ys = random ? nothing : mpr.centers[2,:]
    spring_layout(mpr.adj; xs=xs, ys=ys, kwargs...)
end

"""
    circular_layout(mpr)

Arranges verticies of a 1-skeleton of the Mapper complex `mpr` on a circle.
"""
circular_layout(adj::AbstractMatrix) = circlepoints(size(adj, 1)+1, 1.0)
circular_layout(mpr::Mapper) = circular_layout(mpr.adj)

function constant_layout(data::AbstractMatrix)
    inner_layout(cplx::SimplicialComplex) = (data[1,:], data[2,:])
    return inner_layout
end
constant_layout(mpr::Mapper) = (mpr.centers[1,:], mpr.centers[2,:])


text:
using Nihilist

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

##Simple chart make, message send with upper/lowercase letters, numbers, period and special chars
function run_test1()
	chart = Nihilist.chart_maker()
	cleartext="The Nihilist cipher was used extensively by the Soviets during and after the war in the 40s."
	ciphertext = Nihilist.encode(cleartext, chart)
	plaintext = Nihilist.decode(ciphertext, chart)
	strip_text=join(split(lowercase(cleartext)))
	return isequal(strip_text, plaintext[1:length(strip_text)])
end

function run_test2()
	chart = Nihilist.chart_maker(file_loc="~/gitHub/chart.csv")
	cleartext="The Nihilist cipher was used extensively by the Soviets during and after the war in the 40s."
	ciphertext = Nihilist.encode(cleartext, "~/gitHub/chart.csv", spacing_enable=false, num_repeat=4, file_out="output.txt")
	plaintext = Nihilist.decode("output.txt", "~/gitHub/chart.csv", num_repeat=4, file_out="clearoutput.txt")
	strip_text=join(split(lowercase(cleartext)))
	return isequal(strip_text, plaintext[1:length(strip_text)])
end

function run_test3()
	chart = Nihilist.chart_maker()
	cleartext="The Nihilist cipher was used extensively by the Soviets during and after the war in the 40s."
	ciphertext = Nihilist.encode(cleartext, chart, group_length=3, group_per_line=4, num_repeat=3, file_out="output2.txt")
	plaintext = Nihilist.decode("output2.txt", chart, num_repeat=3, file_out="clearoutput2.txt")
	strip_text=join(split(lowercase(cleartext)))
	return isequal(strip_text, plaintext[1:length(strip_text)])
end

function run_test4()
	chart = chart_maker()
	cleartext="The VIC cipher was a pencil and paper cipher used by the Soviet spy Reino Hayhanen codenamed VICTOR.
				If the cipher were to be given a modern technical name it would be known as a straddling bipartite
				monoalphabetic substitution superenciphered by modified double transposition.1 However by general
				classification it is part of the Nihilist family of ciphers. It was arguably the most complex handoperated
				cipher ever seen when it was first discovered. The initial analysis done by the American National Security
				Agency NSA in 1953 did not absolutely conclude that it was a hand cipher but its placement in a hollowed
				out 5c coin implied it could be broken by pencil and paper. The VIC cipher remained unbroken until more
				information about its structure was available. Although certainly not as complex or secure as modern
				computer operated stream ciphers or block ciphers in practice messages protected by it resisted all
				attempts at cryptanalysis by at least the NSA from its discovery in 1953 until Hayhanens defection
				in 1957."
	ciphertext = encode(cleartext, chart, padding=false)
	plaintext = decode(ciphertext, chart)
	strip_text=join(split(lowercase(cleartext)))
	return isequal(strip_text, plaintext[1:length(strip_text)])
end

@testset "Encode/Decode" begin

	@test run_test1()
	@test run_test2()
	@test run_test3()
	@test run_test4()
end

#= output-

3×11 DataFrames.DataFrame
│ Row │ 1 │ 3 │ 7 │ 6 │ 2 │ 4 │ 5 │ 8 │ 0 │ 9 │ CODE │
├─────┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┼──────┤
│ 1   │ X │ X │ t │ n │ r │ i │ a │ o │ e │ s │ X    │
│ 2   │ z │ . │ f │ y │ q │ d │ l │ b │ p │ h │ 3    │
│ 3   │ v │ k │ u │ g │ m │ j │ c │ x │ / │ w │ 1    │

The Nihilist cipher was used extensively by the Soviets during and after the war in the 40s.
 =>
73906 43935 49715 43039 02195 91790 34018 70694
11035 36383 67390 98114 07934 17246 16563 45377
02739 01952 46739 01044 41000 09337 87615 48565
 =>
thenihilistcipherwasusedextensivelybythesovietsduringandafterthewarinthe40s.totncioan
=#

ext: jl


text:
__precompile__(true)

module ClassImbalance

export smote, rose

include("utils.jl")
include("smote_exs.jl")
include("ub_smote.jl")
include("rose.jl")

end # end module ClassImbalance


text:
# This file includes modified code from https://github.com/odow/SDDP.jl as at
# commit a6662eecff6d5de5c499b760773a8569a0979b9b.

#  Copyright 2017,  Oscar Dowson, Eyob Zewdie
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

#=
This problem is a version of the Ambulance dispatch problem. A hospital is
located at 0 on the number line that stretches from 0 to 100. Ambulance bases
are located at points 20, 40, 60, 80, and 100. When not responding to a call,
Ambulances must be located at a base, or the hospital. In this example there are
three ambulances.

Example location:

    H       B       B       B       B       B
    0 ---- 20 ---- 40 ---- 60 ---- 80 ---- 100

Each stage, a call comes in from somewhere on the number line. The agent must
decide which ambulance to dispatch. They pay the cost of twice the driving
distance. If an ambulance is not dispatched in a stage, the ambulance can be
relocated to a different base in preparation for future calls. This incurrs a
cost of the driving distance.
=#
using SDDP, JuMP, GLPKMathProgInterface, Base.Test, SDDiP

function vehiclelocationmodel(nvehicles, baselocations, requestlocations)
    Locations = baselocations       # Points on the number line where emergency bases are located
    Vehicles = 1:nvehicles          # ambulances
    nbases = length(baselocations)
    Bases = 1:nbases                # base 1 = hostpital
    Requests = collect(requestlocations)     # Points on the number line where calls come from

    shiftcost(src, dest) = abs(Locations[src] - Locations[dest])
    dispatchcost(request, base) = 2 * (abs(request - Locations[1]) + abs(request-Locations[base]))

    #Initial State of emergency vehicles at bases
    Q0 = zeros(Int64, (length(Bases), length(Vehicles)))
    Q0[1,:] = 1 # all ambulances start at hospital

    m = SDDPModel(
                 stages = 10,
        objective_bound = 0.0,
                  sense = :Min,
                 solver = GLPKSolverMIP()
                            ) do sp, t

        # Vehicles at which bases?
        @binarystate(sp, 0 <= q[b=Bases, v=Vehicles] <= 1, q0 == Q0[b, v])

        @variables(sp, begin
            # which vehicle is dipatched?
            0 <= dispatch[src=Bases, v=Vehicles] <= 1
            # shifting vehicles between bases
            0 <= shift[src=Bases, v=Vehicles, dest=Bases] <= 1
        end)

        @expression(sp, basebalance[b in Bases, v in Vehicles],
            # flow of vehicles in and out of bases:
            #initial - dispatch      - shifted from      + shifted to
            q0[b, v] - dispatch[b,v] - sum(shift[b,v,:]) + sum(shift[:,v,b])
        )
        @constraints(sp, begin
            # only one vehicle dispatched to call
            sum(dispatch) == 1
            # can only dispatch vehicle from base if vehicle is at that base
            [b in Bases, v in Vehicles], dispatch[b,v] <= q0[b,v]
            # can only shift vehicle if vehicle is at that src base
            # use formulation that gives stronger bounds
            [b in Bases, v in Vehicles, c in Bases], shift[b,v,c] <= q0[b, v]
            # can only shift vehicle if vehicle is not being dispatched
            [b in Bases, v in Vehicles, c in Bases], shift[b,v,c] + dispatch[b,v] <= 1
            # can't shift to same base
            [b in Bases, v in Vehicles], shift[b,v,b] == 0
            # Update states for non-home/non-hospital bases
            [b in Bases[2:end], v in Vehicles], q[b, v] == basebalance[b,v]
            # Update states for home/hospital bases
            [v in Vehicles], q[1, v] == basebalance[1,v] + sum(dispatch[:,v])
        end)

        if t == 1
            @stageobjective(sp, sum(
                    #distance to travel from base to emergency and back to home base
                    dispatch[b,v] * dispatchcost(50, b) +
                    #distance travelled by vehilces relocating bases
                    sum(shiftcost(b, dest) * shift[b, v, dest] for dest in Bases)
                for b in Bases, v in Vehicles)
            )
        else
            @stageobjective(sp, request=Requests, sum(
                    #distance to travel from base to emergency and back to home base
                    dispatch[b,v] * dispatchcost(request, b) +
                    #distance travelled by vehilces relocating bases
                    sum(shiftcost(b, dest) * shift[b, v, dest] for dest in Bases)
                for b in Bases, v in Vehicles)
        )
        end
        setSDDiPsolver!(sp,
            method=KelleyMethod(),
            LPsolver = GLPKSolverLP(),
            pattern=Pattern(benders=4, lagrangian=1, strengthened_benders=1)
            )
    end
end

srand(1234)
ambulancemodel = vehiclelocationmodel(3, [0, 20, 40, 60, 80, 100], 0:10:100)
@assert solve(ambulancemodel, max_iterations=20) == :iteration_limit
@test getbound(ambulancemodel) >= 1206.0

# Symmetry breaking constraints
# @variables(sp, begin
#     # Let s[u, v, b] = 1 if vehicles u and v are both at base b
#     s[u=1:nvehicles-1, v=u+1:nvehicles, Bases], Bin
# end)
# @constraints(sp, begin
#     # Definition of s
#     [u=1:nvehicles-1, v=u+1:nvehicles, b in Bases], s[u, v, b] <= q0[b, u]
#     [u=1:nvehicles-1, v=u+1:nvehicles, b in Bases], s[u, v, b] <= q0[b, v]
#     [u=1:nvehicles-1, v=u+1:nvehicles, b in Bases], q0[b, u] + q0[b, v] <= s[u, v, b] + 1
#     # If multiple vehicles are at the same base, dispatch vehicle v before v+i
#     [b in Bases, u=1:nvehicles-1, v=u+1:nvehicles], dispatch[b, u] - dispatch[b, v] >= s[u,v, b] - 1
#     # From any base, shift u to b and v to c if b has a lower index than c, unless u is being dispatched
#     [b in Bases, u=1:nvehicles-1, v=u+1:nvehicles, c=1:nbases, d=c:nbases],
#         shift[b, u, c] + dispatch[b, u] - shift[b, v, d] >= s[u, v, b] - 1
# end)


text:
"""
    LiftingLine(submodels)

Construct a lifting line model using the aerodynamic models in `models`. State variables,
inputs, and parameters correspond to the state variables, inputs, and parameters of each of
the cross sections concatenated
"""
function LiftingLine(submodels)

    # residual
    f = (resid, dx, x, y, p, t) -> liftingline_residual!(resid, dx, x, y, p, t; submodels = submodels)

    # dimensions
    nx = sum(number_of_states.(submodels))
    ny = sum(number_of_inputs.(submodels))
    np = sum(number_of_parameters.(submodels))

    # rate jacobian
    if all(isempty.(rate_jacobian_type.(submodels)))
        ratejac = Empty()
    elseif all(iszero.(rate_jacobian_type.(submodels)))
        ratejac = Zeros()
    elseif all(isidentity.(rate_jacobian_type.(submodels)))
        ratejac = Identity()
    elseif all(isinvariant.(rate_jacobian_type.(submodels)))
        ratejac = Invariant(liftingline_rate_jacobian(; submodels))
    elseif all(isconstant.(rate_jacobian_type.(submodels)))
        ratejac = Constant((J, p) -> liftingline_rate_jacobian!(J, p; submodels))
    elseif all(islinear.(rate_jacobian_type.(submodels)))
        ratejac = Linear(
            (J, dx, x, y, p, t) -> liftingline_rate_jacobian!(J, dx, x, y, p, t; submodels)
        )
    else
        ratejac = Nonlinear(
            (J, dx, x, y, p, t) -> liftingline_rate_jacobian!(J, dx, x, y, p, t; submodels)
        )
    end

    # state jacobian
    if all(isempty.(state_jacobian_type.(models)))
        statejac = Empty()
    elseif all(iszero.(state_jacobian_type.(models)))
        statejac = Zeros()
    elseif all(isidentity.(state_jacobian_type.(models)))
        statejac = Identity()
    elseif all(isinvariant.(state_jacobian_type.(models)))
        statejac = Invariant(liftingline_state_jacobian(; submodels))
    elseif all(isconstant.(state_jacobian_type.(models)))
        statejac = Constant((J, p) -> liftingline_state_jacobian!(J, p; submodels))
    elseif all(islinear.(state_jacobian_type.(submodels)))
        statejac = Linear(
            (J, dx, x, y, p, t) -> liftingline_state_jacobian!(J, dx, x, y, p, t; submodels)
        )
    else
        statejac = Nonlinear(
            (J, dx, x, y, p, t) -> liftingline_state_jacobian!(J, dx, x, y, p, t; submodels)
        )
    end

    # input jacobian
    if all(isempty.(input_jacobian_type.(models)))
        inputjac = Empty()
    elseif all(iszero.(input_jacobian_type.(models)))
        inputjac = Zeros()
    elseif all(isidentity.(input_jacobian_type.(models)))
        inputjac = Identity()
    elseif all(isinvariant.(input_jacobian_type.(models)))
        inputjac = Invariant(liftingline_input_jacobian(; submodels))
    elseif all(isconstant.(input_jacobian_type.(models)))
        inputjac = Constant((J, p) -> liftingline_input_jacobian!(J, p; submodels))
    elseif all(islinear.(input_jacobian_type.(submodels)))
        inputjac = Linear(
            (J, dx, x, y, p, t) -> liftingline_input_jacobian!(J, dx, x, y, p, t; submodels)
        )
    else
        inputjac = Nonlinear(
            (J, dx, x, y, p, t) -> liftingline_input_jacobian!(J, dx, x, y, p, t; submodels)
        )
    end

    # parameter jacobian
    if all(isempty.(parameter_jacobian_type.(models)))
        paramjac = Empty()
    elseif all(iszero.(parameter_jacobian_type.(models)))
        paramjac = Zeros()
    elseif all(isidentity.(parameter_jacobian_type.(models)))
        paramjac = Identity()
    elseif all(isinvariant.(parameter_jacobian_type.(models)))
        paramjac = Invariant(liftingline_parameter_jacobian(; submodels))
    elseif all(isconstant.(parameter_jacobian_type.(models)))
        paramjac = Constant((J, p) -> liftingline_parameter_jacobian!(J, p; submodels))
    elseif all(islinear.(parameter_jacobian_type.(submodels)))
        paramjac = Linear(
            (J, dx, x, y, p, t) -> liftingline_parameter_jacobian!(J, dx, x, y, p, t; submodels)
        )
    else
        paramjac = Nonlinear(
            (J, dx, x, y, p, t) -> liftingline_parameter_jacobian!(J, dx, x, y, p, t; submodels)
        )
    end

    # time gradient
    if all(isempty.(time_gradient_type.(models)))
        tgrad = Empty()
    elseif all(iszero.(time_gradient_type.(models)))
        tgrad = Zeros()
    elseif all(isinvariant.(time_gradient_type.(models)))
        tgrad = Invariant(liftingline_time_gradient(; submodels))
    elseif all(isconstant.(time_gradient_type.(models)))
        tgrad = Constant((dT, p) -> liftingline_time_gradient!(dT, p; submodels))
    elseif all(islinear.(time_gradient_type.(models)))
        paramjac = Linear(
            (dT, dx, x, y, p, t) -> liftingline_time_gradient!(dT, dx, x, y, p, t; submodels)
        )
    else
        paramjac = Nonlinear(
            (dT, dx, x, y, p, t) -> liftingline_time_gradient!(dT, dx, x, y, p, t; submodels)
        )
    end

    # input/output functions
    setstate = (x; kwargs...) -> liftingline_setstate!(x; submodels, kwargs...)
    setinput = (y; kwargs...) -> liftingline_setinput!(y; submodels, kwargs...)
    setparam = (p; kwargs...) -> liftingline_setparam!(p; submodels, kwargs...)
    sepstate = (x) -> liftingline_sepstate(x; submodels)
    sepinput = (y) -> liftingline_sepinput(y; submodels)
    sepparam = (p) -> liftingline_sepparam(p; submodels)

    return Model{true}(fresid, nx, ny, np;
        ratejac = ratejac,
        statejac = statejac,
        inputjac = inputjac,
        paramjac = paramjac,
        tgrad = tgrad,
        setstate = setstate,
        setinput = setinput,
        setparam = setparam,
        sepstate = sepstate,
        sepinput = sepinput,
        sepparam = sepparam,
    )
end

function liftingline_residual!(resid, dx, x, y, p, t; submodels)

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    resids = view.(Ref(resid), ix)
    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_residual!.(resids, submodels, dxs, xs, ys, ps, t)

    return resid
end

function liftingline_rate_jacobian(; submodels)
    Jii = get_rate_jacobian.(submodels)
    return block_diagonal(Jii)
end

function liftingline_rate_jacobian!(J, p; submodels)

    J .= 0

    ix = state_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ix)
    ps = view.(Ref(p), ip)

    get_rate_jacobian!.(Js, submodels, ps)

    return J
end

function liftingline_rate_jacobian!(J, dx, x, y, p, t; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ix)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_rate_jacobian!.(Js, submodels, dxs, xs, ys, ps, t)

    return J
end

function liftingline_state_jacobian(; submodels)
    Jii = get_state_jacobian.(submodels)
    return block_diagonal(Jii)
end

function liftingline_state_jacobian!(J, p; submodels)

    J .= 0

    ix = state_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ix)
    ps = view.(Ref(p), ip)

    get_state_jacobian!.(Js, submodels, ps)

    return J
end

function liftingline_state_jacobian!(J, dx, x, y, p, t; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ix)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_state_jacobian!.(Js, submodels, dxs, xs, ys, ps, t)

    return J
end

function liftingline_input_jacobian(; submodels)
    Jyii = get_input_jacobian.(submodels)
    return block_diagonal(Jyii)
end

function liftingline_input_jacobian!(J, p; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, iy)
    ps = view.(Ref(p), ip)

    get_input_jacobian!.(Js, submodels, ps)

    return J
end

function liftingline_input_jacobian!(J, dx, x, y, p, t; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, iy)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_input_jacobian!.(Js, submodels, dxs, xs, ys, ps, t)

    return J
end

function liftingline_parameter_jacobian(; submodels)
    Jii = get_parameter_jacobian.(submodels)
    return block_diagonal(Jii)
end

function liftingline_parameter_jacobian!(J, p; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ip)
    ps = view.(Ref(p), ip)

    get_parameter_jacobian!.(Js, submodels, ps)

    return J
end

function liftingline_parameter_jacobian!(J, dx, x, y, p, t; submodels)

    J .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    Js = view.(Ref(J), ix, ip)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_parameter_jacobian!.(Js, submodels, dxs, xs, ys, ps, t)

    return J
end

liftingline_time_gradient(; submodels) = vcat(get_time_gradient.(submodels)...)

function liftingline_time_gradient!(dT, p; submodels)

    J .= 0

    ix = state_indices(submodels)
    ip = parameter_indices(submodels)

    dTs = view.(Ref(dT), ix)
    ps = view.(Ref(p), ip)

    get_time_gradient!.(dTs, submodels, ps)

    return dT
end

function liftingline_time_gradient!(dT, dx, x, y, p, t; submodels)

    dT .= 0

    ix = state_indices(submodels)
    iy = input_indices(submodels)
    ip = parameter_indices(submodels)

    dTs = view.(Ref(dT), ix)

    dxs = view.(Ref(dx), ix)
    xs = view.(Ref(x), ix)
    ys = view.(Ref(y), iy)
    ps = view.(Ref(p), ip)

    get_time_gradient!.(dTs, submodels, dxs, xs, ys, ps, t)

    return J
end

function liftingline_setstate!(x; section_states, submodels)
    section_indices = state_indices(submodels)
    vxs = view.(Ref(x), section_indices)
    bfn! = (x, model, states) -> set_states!(x, model; states...)
    bfn!.(vxs, submodels, section_states)
    return x
end

function liftingline_setinput!(y; section_inputs, submodels)
    section_indices = input_indices(submodels)
    vys = view.(Ref(y), section_indices)
    bfn! = (y, model, inputs) -> set_inputs!(y, model; inputs...)
    bfn!.(vys, submodels, section_inputs)
    return y
end

function liftingline_setparam!(p; section_parameters, submodels)
    section_indices = parameter_indices(submodels)
    vps = view.(Ref(p), section_indices)
    bfn! = (p, model, parameters) -> set_parameters!(p, model; parameters...)
    bfn!.(vps, submodels, section_parameters)
    return p
end

function liftingline_sepstate(x; submodels)
    section_indices = state_indices(submodels)
    vxs = view.(Ref(x), section_indices)
    return (section_states = separate_states.(submodels, vxs),)
end

function liftingline_sepinput(y; submodels)
    section_indices = input_indices(submodels)
    vys = view.(Ref(y), section_indices)
    return (section_inputs = separate_inputs.(submodels, vys),)
end

function liftingline_sepparam(p; submodels)
    section_indices = parameter_indices(submodels)
    vps = view.(Ref(p), section_indices)
    return (section_parameters = separate_parameters.(submodels, vps),)
end


text:
module PATHSolver

using ForwardDiff

if isfile(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
    include(joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl"))
else
    error("PATHSolver not properly installed. Please run Pkg.build(\"PATHSolver\")")
end

export solveMCP, path_options

function solveMCP(f_eval::Function, lb::Vector, ub::Vector)
    j_eval = x -> ForwardDiff.jacobian(f_eval, x)
    return solveMCP(f_eval, j_eval, lb, ub)
end

function solveMCP(f_eval::Function, j_eval::Function, lb::Vector, ub::Vector)

    global user_f = f_eval
    global user_j = j_eval

    f_user_cb = cfunction(f_user_wrap, Cint, (Cint, Ptr{Cdouble}, Ptr{Cdouble}))
    j_user_cb = cfunction(j_user_wrap, Cint, (Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}))

    n = length(lb)
    z = copy(lb)
    f = zeros(n)

    J0 = j_eval(z)
    s_col, s_len, s_row, s_data = sparse_matrix(J0)
    nnz = length(s_data)

    t = ccall( (:path_solver, "libpath47julia"), Cint,
                (Cint, Cint,
                 Ptr{Cdouble}, Ptr{Cdouble},
                 Ptr{Cdouble}, Ptr{Cdouble},
                 Ptr{Void},Ptr{Void}),
                 n, nnz, z, f, lb, ub, f_user_cb, j_user_cb)

    status =
        [ :Solved,              # 1 - solved
          :StationaryPoint,     # 2 - stationary point found
          :MajorIterLimit,      # 3 - major iteration limit
          :MinorIterLimit,      # 4 - cumulative minor iteration limit
          :TimeLimit,           # 5 - time limit
          :Interrupt,           # 6 - user interrupt
          :BoundError,          # 7 - bound error (lb is not less than ub)
          :DominaError,         # 8 - domain error (could not find a starting point)
          :InternalError        # 9 - internal error
         ]

    remove_option_file()

    return status[t], z, f
end

function remove_option_file()
    if isfile("path.opt")
        rm("path.opt")
    end
end

function path_options(options...)
    opt_file = open("path.opt", "w")
    println(opt_file, "* Generated by PATHSolver.jl. Do not edit.")
    for opt in options
        println(opt_file, opt)
    end
    close(opt_file)
end

###############################################################################
# wrappers for callback functions
###############################################################################
# static int (*f_eval)(int n, double *z, double *f);
# static int (*j_eval)(int n, int nnz, double *z, int *col_start, int *col_len,
            # int *row, double *data);

function f_user_wrap(n::Cint, z::Ptr{Cdouble}, f::Ptr{Cdouble})
    F = user_f(pointer_to_array(z, Int(n)))
    unsafe_store_vector!(f, F)
    return Cint(0)
end

function j_user_wrap(n::Cint, nnz::Cint, z::Ptr{Cdouble},
    col_start::Ptr{Cint}, col_len::Ptr{Cint}, row::Ptr{Cint}, data::Ptr{Cdouble})

    J = user_j(pointer_to_array(z, Int(n)) )

    s_col, s_len, s_row, s_data = sparse_matrix(J)

    unsafe_store_vector!(col_start, s_col)
    unsafe_store_vector!(col_len, s_len)
    unsafe_store_vector!(row, s_row)
    unsafe_store_vector!(data, s_data)

    return Cint(0)
end
###############################################################################

###############################################################################
# Converting the Jacobian matrix to the sparse matrix format of the PATH Solver
###############################################################################
function sparse_matrix(A::AbstractSparseArray)
    m, n = size(A)
    @assert m==n

    col_start = Array{Int}(n)
    col_len = Array{Int}(n)
    row = Array{Int}(0)
    data = Array{Float64}(0)
    for j in 1:n
        if j==1
            col_start[j] = 1
        else
            col_start[j] = col_start[j-1] + col_len[j-1]
        end

        col_len[j] = 0
        for i in 1:n
            if A[i,j] != 0.0
                col_len[j] += 1
                push!(row, i)
                push!(data, A[i,j])
            end
        end
    end

    return col_start, col_len, row, data
end

function sparse_matrix(A::Matrix)
    return sparse_matrix(sparse(A))
    # m, n = size(A)
    # @assert m==n
    #
    # col_start = Array{Int}(n)
    # col_len = Array{Int}(n)
    # row = Array{Int}(0)
    # data = Array{Float64}(0)
    # for j in 1:n
    #     if j==1
    #         col_start[j] = 1
    #     else
    #         col_start[j] = col_start[j-1] + col_len[j-1]
    #     end
    #
    #     col_len[j] = 0
    #     for i in 1:n
    #         if A[i,j] != 0.0
    #             col_len[j] += 1
    #             push!(row, i)
    #             push!(data, A[i,j])
    #         end
    #     end
    # end
    #
    # return col_start, col_len, row, data
end
###############################################################################

function unsafe_store_vector!(x_ptr::Ptr{Cint}, x_val::Vector)
    for i in 1:length(x_val)
        unsafe_store!(x_ptr, x_val[i], i)
    end
    return
end

function unsafe_store_vector!(x_ptr::Ptr{Cdouble}, x_val::Vector)
    for i in 1:length(x_val)
        unsafe_store!(x_ptr, x_val[i], i)
    end
    return
end

end # Module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LoadBalancerProbeListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{Probe} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function LoadBalancerProbeListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(LoadBalancerProbeListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(LoadBalancerProbeListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type LoadBalancerProbeListResult

const _property_map_LoadBalancerProbeListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_LoadBalancerProbeListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{Probe}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ LoadBalancerProbeListResult }) = collect(keys(_property_map_LoadBalancerProbeListResult))
Swagger.property_type(::Type{ LoadBalancerProbeListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerProbeListResult[name]))}
Swagger.field_name(::Type{ LoadBalancerProbeListResult }, property_name::Symbol) =  _property_map_LoadBalancerProbeListResult[property_name]

function check_required(o::LoadBalancerProbeListResult)
    true
end

function validate_property(::Type{ LoadBalancerProbeListResult }, name::Symbol, val)
end


text:
"""
    MattssonSvärdShoeybi2008()

Coefficients of the SBP operators given in
- Mattsson, Svärd, Shoeybi (2008)
  Stable and accurate schemes for the compressible Navier-Stokes equations.
  Journal of Computational Physics 227, pp. 2293-2316.
"""
struct MattssonSvärdShoeybi2008 <: SourceOfCoefficients end

function Base.show(io::IO, source::MattssonSvärdShoeybi2008)
    if get(io, :compact, false)
      summary(io, source)
    else
        print(io,
            "Mattsson, Svärd, Shoeybi (2008) \n",
            "  Stable and accurate schemes for the compressible Navier-Stokes equations. \n",
            "  Journal of Computational Physics 227, pp. 2293-2316. \n",
            "See also \n",
            "  Mattsson, Nordström (2004) \n",
            "  Summation by parts operators for finite difference approximations of second \n",
            "    derivatives. \n",
            "  Journal of Computational Physics 199, pp. 503-540.")
    end
end

function first_derivative_coefficients(source::MattssonSvärdShoeybi2008, order::Int,
                                       T=Float64, mode=FastMode())
    if order == 2
        left_boundary = (
            DerivativeCoefficientRow{T,1,2}(SVector(-one(T), one(T))),
        )
        right_boundary = (
            -left_boundary[1],
        )
        upper_coef = SVector(T(1//2))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 4
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,4}(SVector(T(-24//17),
                                                    T(59//34),
                                                    T(-4//17),
                                                    T(-3//34) )),
            # q2
            DerivativeCoefficientRow{T,1,3}(SVector(T(-1//2),
                                                    T(0),
                                                    T(1//2))),
            # q3
            DerivativeCoefficientRow{T,1,5}(SVector(T(4//43),
                                                    T(-59//86),
                                                    T(0),
                                                    T(59//86),
                                                    T(-4//43))),
            # q4
            DerivativeCoefficientRow{T,1,6}(SVector(T(3//98),
                                                    T(0),
                                                    T(-59//98),
                                                    T(0),
                                                    T(32//49),
                                                    T(-4//49))),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(2//3), T(-1//12))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(17//48),
                                T(59//48),
                                T(43//48),
                                T(49//48) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 6
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,6}(SVector(T(-21600//13649),
                                                    T(104009//54596),
                                                    T(30443//81894),
                                                    T(-33311//27298),
                                                    T(16863//27298),
                                                    T(-15025//163788) )),
            # q2
            DerivativeCoefficientRow{T,1,6}(SVector(T(-104009//240260),
                                                    T(0),
                                                    T(-311//72078),
                                                    T(20229//24026),
                                                    T(-24337//48052),
                                                    T(36661//360390) )),
            # q3
            DerivativeCoefficientRow{T,1,6}(SVector(T(-30443//162660),
                                                    T(311//32532),
                                                    T(0),
                                                    T(-11155//16266),
                                                    T(41287//32532),
                                                    T(-21999//54220) )),
            # q4
            DerivativeCoefficientRow{T,1,7}(SVector(T(33311//107180),
                                                    T(-20229//21436),
                                                    T(485//1398),
                                                    T(0),
                                                    T(4147//21436),
                                                    T(25427//321540),
                                                    T(72//5359) )),
            # q5
            DerivativeCoefficientRow{T,1,8}(SVector(T(-16863//78770),
                                                    T(24337//31508),
                                                    T(-41287//47262),
                                                    T(-4147//15754),
                                                    T(0),
                                                    T(342523//472620),
                                                    T(-1296//7877),
                                                    T(144//7877) )),
            # q6
            DerivativeCoefficientRow{T,1,9}(SVector(T(15025//525612),
                                                    T(-36661//262806),
                                                    T(21999//87602),
                                                    T(-25427//262806),
                                                    T(-342523//525612),
                                                    T(0),
                                                    T(32400//43801),
                                                    T(-6480//43801),
                                                    T(720//43801) ))
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(13649//43200),
                                T(12013//8640),
                                T(2711//4320),
                                T(5359//4320),
                                T(7877//8640),
                                T(43801//43200) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 8
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,8}(SVector(T(-2540160//1498139),
                                                    T(515174618457408334671//81327545493318772580),
                                                    T(-319653445960068136431//16265509098663754516),
                                                    T(1756838164798071830975//48796527295991263548),
                                                    T(-132855587210457907869//4066377274665938629),
                                                    T(205031990442925032261//16265509098663754516),
                                                    T(-147365687//1707878460),
                                                    T(-13621478277275954493//16265509098663754516) )),
            # q2
            DerivativeCoefficientRow{T,1,8}(SVector(T(-515174618457408334671//420776658856749814780),
                                                    T(0),
                                                    T(335158710375082986831//30055475632624986770),
                                                    T(-81614753349680015895//3005547563262498677),
                                                    T(996306452293020120805//36066570759149984124),
                                                    T(-338011221795904141851//30055475632624986770),
                                                    T(1109310//21038833),
                                                    T(503120655918727839389//631164988285124722170) )),
            # q3
            DerivativeCoefficientRow{T,1,8}(SVector(T(106551148653356045477//4733507987139058764),
                                                    T(-111719570125027662277//1690538566835378130),
                                                    T(0),
                                                    T(66611314247587471205//338107713367075626),
                                                    T(-177186717090009550283//676215426734151252),
                                                    T(13307564278128858831//112702571122358542),
                                                    T(-354462//1972295),
                                                    T(-7342453837193297117//788917997856509794) )),
            # q4
            DerivativeCoefficientRow{T,1,8}(SVector(T(-1756838164798071830975//297539771191584439476),
                                                    T(81614753349680015895//3542140133233148089),
                                                    T(-199833942742762413615//7084280266466296178),
                                                    T(0),
                                                    T(334165099234925485449//14168560532932592356),
                                                    T(-97726720448760690009//7084280266466296178),
                                                    T(-2934266//74384943),
                                                    T(65502786250599341781//49589961865264073246) )),
            # q5
            DerivativeCoefficientRow{T,1,9}(SVector(T(132855587210457907869//5691012984537517679),
                                                    T(-996306452293020120805//9756022259207173164),
                                                    T(531560151270028650849//3252007419735724388),
                                                    T(-334165099234925485449//3252007419735724388),
                                                    T(0),
                                                    T(68552420813742056079//3252007419735724388),
                                                    T(-2343561//22764052),
                                                    T(-197397151320811747355//68292155814450212148),
                                                    T(-2592//299527) )),
            # q6
            DerivativeCoefficientRow{T,1,10}(SVector(T(-9763428116329763441//3358016562304840404),
                                                     T(112670407265301380617//8395041405762101010),
                                                     T(-13307564278128858831//559669427050806734),
                                                     T(32575573482920230003//1679008281152420202),
                                                     T(-22850806937914018693//3358016562304840404),
                                                     T(0),
                                                     T(5346432//9794215),
                                                     T(1388463154122240//14728142817126493),
                                                     T(3072//103097),
                                                     T(-288//103097) )),
             # q7
             DerivativeCoefficientRow{T,1,11}(SVector(T(21052241//763903740),
                                                      T(-1109310//12731729),
                                                      T(3190158//63658645),
                                                      T(2934266//38195187),
                                                      T(2343561//50926916),
                                                      T(-48117888//63658645),
                                                      T(0),
                                                      T(52309152//63658645),
                                                      T(-145152//670091),
                                                      T(27648//670091),
                                                      T(-2592//670091) )),
              # q8
              DerivativeCoefficientRow{T,1,12}(SVector(T(13621478277275954493//55672594705880416916),
                                                       T(-503120655918727839389//417544460294103126870),
                                                       T(66082084534739674053//27836297352940208458),
                                                       T(-65502786250599341781//27836297352940208458),
                                                       T(197397151320811747355//167017784117641250748),
                                                       T(-87473178709701120//732534140866847591),
                                                       T(-366164064//487135205),
                                                       T(0),
                                                       T(4064256//5127739),
                                                       T(-1016064//5127739),
                                                       T(193536//5127739),
                                                       T(-18144//5127739) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(1498139//5080320),
                                T(1107307//725760),
                                T(20761//80640),
                                T(1304999//725760),
                                T(299527//725760),
                                T(103097//80640),
                                T(670091//725760),
                                T(5127739//5080320) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

function second_derivative_coefficients(source::MattssonSvärdShoeybi2008, order::Int,
                                        T=Float64, mode=FastMode())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(1))
        central_coef = T(-2)
        lower_coef = upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
                                                    T(2),
                                                    T(-1//2) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,4}(SVector(T(2),
                                                    T(-5),
                                                    T(4),
                                                    T(-1) )),
            # d2
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
            # d3
            DerivativeCoefficientRow{T,1,5}(SVector(T(-4//43),
                                                    T(59//43),
                                                    T(-110//43),
                                                    T(59//43),
                                                    T(-4//43) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(-1//49),
                                                    T(0),
                                                    T(59//49),
                                                    T(-118//49),
                                                    T(64//49),
                                                    T(-4//49) ))
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(4//3), T(-1//12))
        central_coef = T(-5//2)
        lower_coef = upper_coef
        left_weights = SVector(T(17//48), T(59//48), T(43//48), T(49//48))
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6),
                                                    T(3),
                                                    T(-3//2),
                                                    T(1//3) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(114170//40947),
                                                    T(-438107//54596),
                                                    T(336409//40947),
                                                    T(-276997//81894),
                                                    T(3747//13649),
                                                    T(21035//163788) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(6173//5860),
                                                    T(-2066//879),
                                                    T(3283//1758),
                                                    T(-303//293),
                                                    T(2111//3516),
                                                    T(-601//4395) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(-52391//81330),
                                                    T(134603//32532),
                                                    T(-21982//2711),
                                                    T(112915//16266),
                                                    T(-46969//16266),
                                                    T(30409//54220) )),
            # d4
            DerivativeCoefficientRow{T,1,7}(SVector(T(68603//321540),
                                                    T(-12423//10718),
                                                    T(112915//32154),
                                                    T(-75934//16077),
                                                    T(53369//21436),
                                                    T(-54899//160770),
                                                    T(48//5359) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(-7053//39385),
                                                    T(86551//94524),
                                                    T(-46969//23631),
                                                    T(53369//15754),
                                                    T(-87904//23631),
                                                    T(820271//472620),
                                                    T(-1296//7877),
                                                    T(96//7877) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(21035//525612),
                                                    T(-24641//131403),
                                                    T(30409//87602),
                                                    T(-54899//131403),
                                                    T(820271//525612),
                                                    T(-117600//43801),
                                                    T(64800//43801),
                                                    T(-6480//43801),
                                                    T(480//43801) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(3//2), T(-3//20), T(1//90))
        central_coef = T(-49//18)
        lower_coef = upper_coef
        left_weights = SVector( T(13649//43200),
                                T(12013//8640),
                                T(2711//4320),
                                T(5359//4320),
                                T(7877//8640),
                                T(43801//43200) )
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12),
                                                    T(4),
                                                    T(-3),
                                                    T(4//3),
                                                    T(-1//4) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 8
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,7}(SVector(T(4870382994799//1358976868290),
                                                    T(-893640087518//75498714905),
                                                    T(926594825119//60398971924),
                                                    T(-1315109406200//135897686829),
                                                    T(39126983272//15099742981),
                                                    T(12344491342//75498714905),
                                                    T(-451560522577//2717953736580) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(333806012194//390619153855),
                                                    T(-154646272029//111605472530),
                                                    T(1168338040//33481641759),
                                                    T(82699112501//133926567036),
                                                    T(-171562838//11160547253),
                                                    T(-28244698346//167408208795),
                                                    T(11904122576//167408208795),
                                                    T(-2598164715//312495323084) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(7838984095//52731029988),
                                                    T(1168338040//5649753213),
                                                    T(-88747895//144865467),
                                                    T(423587231//627750357),
                                                    T(-43205598281//22599012852),
                                                    T(4876378562//1883251071),
                                                    T(-5124426509//3766502142),
                                                    T(10496900965//39548272491) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(-94978241528//828644350023),
                                                    T(82699112501//157837019052),
                                                    T(1270761693//13153084921),
                                                    T(-167389605005//118377764289),
                                                    T(48242560214//39459254763),
                                                    T(-31673996013//52612339684),
                                                    T(43556319241//118377764289),
                                                    T(-44430275135//552429566682) )),
            # d5
            DerivativeCoefficientRow{T,1,9}(SVector(T(1455067816//21132528431),
                                                    T(-171562838//3018932633),
                                                    T(-43205598281//36227191596),
                                                    T(48242560214//9056797899),
                                                    T(-52276055645//6037865266),
                                                    T(57521587238//9056797899),
                                                    T(-80321706377//36227191596),
                                                    T(8078087158//21132528431),
                                                    T(-1296//299527) )),
            # d6
            DerivativeCoefficientRow{T,1,10}(SVector(T(10881504334//327321118845),
                                                     T(-28244698346//140280479505),
                                                     T(4876378562//9352031967),
                                                     T(-10557998671//12469375956),
                                                     T(57521587238//28056095901),
                                                     T(-278531401019//93520319670),
                                                     T(73790130002//46760159835),
                                                     T(-137529995233//785570685228),
                                                     T(2048//103097),
                                                     T(-144//103097) )),
             # d7
             DerivativeCoefficientRow{T,1,11}(SVector(T(-135555328849//8509847458140),
                                                      T(11904122576//101307707835),
                                                      T(-5124426509//13507694378),
                                                      T(43556319241//60784624701),
                                                      T(-80321706377//81046166268),
                                                      T(73790130002//33769235945),
                                                      T(-950494905688//303923123505),
                                                      T(239073018673//141830790969),
                                                      T(-145152//670091),
                                                      T(18432//670091),
                                                      T(-1296//670091) )),
              # d8
              DerivativeCoefficientRow{T,1,12}(SVector(T(0),
                                                       T(-2598164715//206729925524),
                                                       T(10496900965//155047444143),
                                                       T(-44430275135//310094888286),
                                                       T(425162482//2720130599),
                                                       T(-137529995233//620189776572),
                                                       T(239073018673//155047444143),
                                                       T(-144648000000//51682481381),
                                                       T(8128512//5127739),
                                                       T(-1016064//5127739),
                                                       T(129024//5127739),
                                                       T(-9072//5127739) )),

        )
        right_boundary = left_boundary
        upper_coef = SVector(T(8//5), T(-1//5), T(8//315), T(-1//560))
        central_coef = T(-205//72)
        lower_coef = upper_coef
        left_weights = SVector( T(1498139//5080320),
                                T(1107307//725760),
                                T(20761//80640),
                                T(1304999//725760),
                                T(299527//725760),
                                T(103097//80640),
                                T(670091//725760),
                                T(5127739//5080320) )
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,7}(SVector(T(-4723//2100),
                                                    T(839//175),
                                                    T(-157//35),
                                                    T(278//105),
                                                    T(-103//140),
                                                    T(-1//175),
                                                    T(6//175) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end


text:
# Retrieve number of variables from equation instance
@inline nvariables(::AbstractEquations{NDIMS, NVARS}) where {NDIMS, NVARS} = NVARS

# TODO: Taal performance, 1:NVARS vs. Base.OneTo(NVARS) vs. SOneTo(NVARS)
@inline eachvariable(equations::AbstractEquations) = Base.OneTo(nvariables(equations))

# Add methods to show some information on systems of equations.
function Base.show(io::IO, equations::AbstractEquations)
  # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.
  @nospecialize equations # reduce precompilation time

  print(io, get_name(equations), " with ")
  if nvariables(equations) == 1
    print(io, "one variable")
  else
    print(io, nvariables(equations), " variables")
  end
end

function Base.show(io::IO, ::MIME"text/plain", equations::AbstractEquations)
  # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.
  @nospecialize equations # reduce precompilation time

  if get(io, :compact, false)
    show(io, equations)
  else
    summary_header(io, get_name(equations))
    summary_line(io, "#variables", nvariables(equations))
    for variable in eachvariable(equations)
      summary_line(increment_indent(io),
                   "variable " * string(variable),
                   varnames(cons2cons, equations)[variable])
    end
    summary_footer(io)
  end
end

@inline Base.ndims(::AbstractEquations{NDIMS}) where NDIMS = NDIMS

"""
    flux(u, orientation, equations)

Given the conservative variables `u`, calculate the (physical) flux in spatial
direction `orientation` for the coressponding set of governing `equations`
`orientation` is `1`, `2`, and `3` for the x-, y-, and z-directions, respectively.
"""
function flux(u, orientation, equations) end

# set sensible default values that may be overwritten by specific equations
have_nonconservative_terms(::AbstractEquations) = Val(false)
have_constant_speed(::AbstractEquations) = Val(false)

default_analysis_errors(::AbstractEquations)     = (:l2_error, :linf_error)
default_analysis_integrals(::AbstractEquations)  = (entropy_timederivative,)

"""
    flux_central(u_ll, u_rr, orientation, equations::AbstractEquations)

The classical central numerical flux `f((u_ll) + f(u_rr)) / 2`. When this flux is
used as volume flux, the discretization is equivalent to the classical weak form
DG method (except floating point errors).
"""
@inline function flux_central(u_ll, u_rr, orientation, equations::AbstractEquations)
  # Calculate regular 1D fluxes
  f_ll = flux(u_ll, orientation, equations)
  f_rr = flux(u_rr, orientation, equations)

  # Average regular fluxes
  return 0.5 * (f_ll + f_rr)
end

@inline cons2cons(u, ::AbstractEquations) = u
function cons2prim(u, ::AbstractEquations) end
@inline Base.first(u, ::AbstractEquations) = first(u)

# FIXME: Deprecations introduced in v0.3
@deprecate varnames_cons(equations) varnames(cons2cons, equations)
@deprecate varnames_prim(equations) varnames(cons2prim, equations)
@deprecate flux_upwind(u_ll, u_rr, orientation, equations) flux_godunov(u_ll, u_rr, orientation, equations)
@deprecate calcflux(u, orientation, equations) flux(u, orientation, equations)

####################################################################################################
# Include files with actual implementations for different systems of equations.

# Linear scalar advection
abstract type AbstractLinearScalarAdvectionEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("linear_scalar_advection_1d.jl")
include("linear_scalar_advection_2d.jl")
include("linear_scalar_advection_3d.jl")

# Inviscid Burgers
abstract type AbstractInviscidBurgersEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("inviscid_burgers_1d.jl")

# CompressibleEulerEquations
abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("compressible_euler_1d.jl")
include("compressible_euler_2d.jl")
include("compressible_euler_3d.jl")

# CompressibleEulerMulticomponentEquations
abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: AbstractEquations{NDIMS, NVARS} end
include("compressible_euler_multicomponent_2d.jl")

# Retrieve number of components from equation instance for the multicomponent case
@inline ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP}) where {NDIMS, NVARS, NCOMP} = NCOMP
@inline eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) = Base.OneTo(ncomponents(equations))

# Ideal MHD
abstract type AbstractIdealGlmMhdEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("ideal_glm_mhd_1d.jl")
include("ideal_glm_mhd_2d.jl")
include("ideal_glm_mhd_3d.jl")

# Diffusion equation: first order hyperbolic system
abstract type AbstractHyperbolicDiffusionEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("hyperbolic_diffusion_1d.jl")
include("hyperbolic_diffusion_2d.jl")
include("hyperbolic_diffusion_3d.jl")

# Lattice-Boltzmann equation (advection part only)
abstract type AbstractLatticeBoltzmannEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("lattice_boltzmann_2d.jl")
include("lattice_boltzmann_3d.jl")

# Acoustic perturbation equations
abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
include("acoustic_perturbation_2d.jl")


text:
using DynamicDecisionNetworks
using Test

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


text:
# compare FD to a traditional rank-k approximation
# Danny Perry (dperry@cs.utah.edu)
# Oct 2014
#

using PyPlot
using Debug
using RDatasets

using Sketch

X = []
Labels = []

#name = "iris"
name = "baseball"

if name == "iris"
	iris = dataset("datasets","iris")

	X = [ convert(Array{Float64,1}, iris[1]) convert(Array{Float64,1}, iris[2]) convert(Array{Float64,1}, iris[3]) convert(Array{Float64,1}, iris[4]) ]
	n = size(X,1)
	println("X: ", size(X))
	Labels = ASCIIString[]
	for i=1:n
		push!(Labels, iris[5][i])
	end
else

	baseball = dataset("plyr","baseball")

	numerical_cols = [1,4,8,9,10,11,12,13,17]
	X = zeros(size(baseball,1),length(numerical_cols))

	for i=1:length(numerical_cols)
		println(i,": ",numerical_cols[i])
		X[:,i] = convert(Array{Float64,1}, baseball[numerical_cols[i]])
	end

end

n = size(X,1)
X_centered = X .- (ones(n,1)*mean(X,1))

@debug function compare(X)

	n = size(X,1)
	d = size(X,2)

	errs = Float64[]
	minerrs = Float64[]
	maxerrs = Float64[]

	figure()
	for l=2:20
		B = FD(X,l)
		BB = B'*B
		U,S,V = svd(BB)

		projerr = sum((X .- (X*V*V')).^2,2)
		err = mean(projerr)
		minerr = minimum(projerr)
		maxerr = maximum(projerr)

		push!(errs,err)
		push!(minerrs,minerr)
		push!(maxerrs,maxerr)

		clf()
		plot(1:length(errs),errs,"r-")
		hold(:on)
		plot(1:length(errs),minerrs,"r--")
		plot(1:length(errs),maxerrs,"r--")
		title("FD projection error")
		xlabel("size of sketch")
		ylabel("proj err squared")
		sleep(0.0001)
	end

	@bp
end

compare(X_centered)


text:
import OnlineStats: OnlineStat, Extrema, fit!, value, HistogramStat, Ash
import ...Cubes.Axes: CategoricalAxis, RangeAxis
import IterTools
using WeightedOnlineStats
import ProgressMeter: next!, Progress, ProgressUnknown

import WeightedOnlineStats: WeightedOnlineStat
abstract type TableAggregator end
struct OnlineAggregator{O,S} <: TableAggregator
    o::O
end
function OnlineAggregator(O::OnlineStat, s::Symbol)
    OnlineAggregator{typeof(O),s}(copy(O))
end
cubeeltype(t::OnlineAggregator) = Float64
function fitrow!(o::OnlineAggregator{T,S}, r) where {T<:OnlineStat,S}
    v = getproperty(r, S)
    !ismissing(v) && fit!(o.o, v)
end
value(o::OnlineAggregator) = value(o.o)
struct WeightOnlineAggregator{O,S,W} <: TableAggregator
    o::O
    w::W
end
function WeightOnlineAggregator(O::WeightedOnlineStat, s::Symbol, w)
    WeightOnlineAggregator{typeof(O),s,typeof(w)}(copy(O), w)
end
cubeeltype(t::WeightOnlineAggregator{T}) where {T} = cubeeltype(T)
value(o::WeightOnlineAggregator) = value(o.o)
function fitrow!(o::WeightOnlineAggregator{T,S}, r) where {T<:OnlineStat,S}
    v = getproperty(r, S)
    w = o.w(r)
    if !checkmiss(v) && !ismissing(w)
        fit!(o.o, v, w)
    end
end
checkmiss(v) = ismissing(v)
checkmiss(v::AbstractVector) = any(ismissing, v)
struct GroupedOnlineAggregator{O,S,BY,W,C} <: TableAggregator
    d::O
    w::W
    by::BY
    cloneobj::C
end
value(o::GroupedOnlineAggregator) = Dict(zip(keys(o.d), map(value, (values(o.d)))))
struct SymType{S} end
SymType(s::Symbol) = SymType{s}()
(f::SymType{S})(x) where {S} = getproperty(x, S)
getbytypes(et, by) =
    Tuple{map(i -> Base.nonmissingtype(Base.return_types(i, Tuple{et})[1]), by)...}
cubeeltype(t::GroupedOnlineAggregator{T}) where {T} = cubeeltype(T)
cubeeltype(t::Type{<:Dict{<:Any,T}}) where {T} = cubeeltype(T)
cubeeltype(t::Type{<:WeightedOnlineStat{T}}) where {T} = T
cubeeltype(t::Type{<:OnlineStat{<:Number}}) = Float64
cubeeltype(t::Type{<:WeightedCovMatrix{T}}) where {T} = T
cubeeltype(t::Type{<:Extrema{T}}) where {T} = T

function GroupedOnlineAggregator(O::OnlineStat, s::Symbol, by, w, iter)
    ost = typeof(O)
    et = eltype(iter)
    bytypes =
        Tuple{map(i -> Base.nonmissingtype(Base.return_types(i, Tuple{et})[1]), by)...}
    d = Dict{bytypes,ost}()
    GroupedOnlineAggregator{typeof(d),s,typeof(by),typeof(w),ost}(d, w, by, O)
end

dicteltype(::Type{<:Dict{K,V}}) where {K,V} = V
dictktype(::Type{<:Dict{K,V}}) where {K,V} = K
actval(v) = v
actval(v::SentinelMissings.SentinelMissing) = v[]
function fitrow!(o::GroupedOnlineAggregator{T,S,BY,W}, r) where {T,S,BY,W}
    v = getproperty(r, S)
    if !ismissing(v)
        w = o.w(r)
        if w===nothing
            bykey = map(i->actval(i(r)),o.by)
            if !any(ismissing,bykey)
                if haskey(o.d,bykey)
                    fit!(o.d[bykey],actval(v))
                else
                    o.d[bykey] = copy(o.cloneobj)
                    fit!(o.d[bykey],actval(v))
                end
            end
        else
            if !ismissing(w)
                bykey = map(i -> actval(i(r)), o.by)
                if !any(ismissing, bykey)
                    if haskey(o.d, bykey)
                        fit!(o.d[bykey], v, w)
                    else
                        o.d[bykey] = copy(o.cloneobj)
                        fit!(o.d[bykey], v, w)
                    end
                end
            end
        end
    end
end
export TableAggregator, fittable, cubefittable
function TableAggregator(iter, O, fitsym; by = (), weight = nothing)
    !isa(by, Tuple) && (by = (by,))
    if !isempty(by)
        weight === nothing && (weight = (i -> nothing))
        by = map(i -> isa(i, Symbol) ? (SymType(i)) : i, by)
        GroupedOnlineAggregator(O, fitsym, by, weight, iter)
    else
        if weight === nothing
            OnlineAggregator(O, fitsym)
        else
            WeightOnlineAggregator(O, fitsym, weight)
        end
    end
end

function tooutaxis(
    ::SymType{s},
    iter::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,S},
    k,
    ibc,
) where {s,S}
    ichosen = findfirst(i -> i === s, fieldnames(S))
    if ichosen <= length(iter.inars)
        bycube = iter.dc.incubes[ichosen].cube
        if haskey(bycube.properties, "labels")
            idict = bycube.properties["labels"]
            axname = get(bycube.properties, "name", "Label")
            outAxis = CategoricalAxis(axname, collect(String, values(idict)))
            convertdict = Dict(k => i for (i, k) in enumerate(keys(idict)))
        else
            sort!(k)
            outAxis = CategoricalAxis("Label$(ibc)", k)
            convertdict = Dict(k => i for (i, k) in enumerate(k))
        end
    else
        iax = findAxis(string(s), iter.dc.LoopAxes)
        outAxis = iter.dc.LoopAxes[iax]
        convertdict = Dict(k => i for (i, k) in enumerate(outAxis.values))
    end
    outAxis, convertdict
end
function tooutaxis(f, iter::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,S}, k, ibc) where {S}
    sort!(k)
    outAxis = CategoricalAxis("Category$(ibc)", k)
    convertdict = Dict(k => i for (i, k) in enumerate(k))
    outAxis, convertdict
end

varsym(::WeightOnlineAggregator{<:Any,S}) where {S} = S
varsym(::OnlineAggregator{<:Any,S}) where {S} = S
varsym(::GroupedOnlineAggregator{<:Any,S}) where {S} = S
axt(::CategoricalAxis) = CategoricalAxis
axt(::RangeAxis) = RangeAxis
getStatType(::WeightOnlineAggregator{T}) where {T} = T
getStatType(::OnlineAggregator{T}) where {T} = T
getStatType(t::GroupedOnlineAggregator{T}) where {T} = getStatType(T)
getStatType(t::Type{<:Dict{<:Any,T}}) where {T} = T

getStatOutAxes(tab, agg) = getStatOutAxes(tab, agg, getStatType(agg))
getStatOutAxes(tab, agg, ::Type{<:OnlineStat}) = ()
function getStatOutAxes(tab, agg, ::Type{<:Extrema})
    (CategoricalAxis(:Extrema, ["min", "max"]),)
end
function getStatOutAxes(tab, agg, ::Type{<:WeightedCovMatrix})
    varn = fieldnames(eltype(tab))
    s = varsym(agg)
    icube = findfirst(isequal(s), varn)
    ax = tab.dc.incubes[icube].axesSmall[1]
    oldname = YAXArrays.Cubes.Axes.axname(ax)
    coname = string("Co", oldname)
    v = ax.values
    axtype = axt(ax)
    a1 = axtype(oldname, copy(v))
    a2 = axtype(coname, copy(v))
    (a1, a2)
end
function getStatOutAxes(tab,agg,::Type{<:Union{Ash,HistogramStat, WeightedAdaptiveHist}})
    nbin = getnbins(agg)
    a1 = RangeAxis("Bin", 1:nbin)
    a2 = CategoricalAxis("Hist", ["MidPoints", "Frequency"])
    (a1, a2)
end
function getByAxes(iter, agg::GroupedOnlineAggregator)
    by = agg.by
    ntuple(length(by)) do ibc
        bc = agg.by[ibc]
        tooutaxis(bc, iter, unique(map(i -> i[ibc], collect(keys(agg.d)))), ibc)
    end
end
getByAxes(iter, agg) = ()

function tooutcube(agg, iter, post)
    outaxby = getByAxes(iter, agg)
    axby = map(i -> i[1], outaxby)
    convdictall = map(i -> i[2], outaxby)

    outaxstat = getStatOutAxes(iter, agg)
    outax = (outaxstat..., axby...)
    snew = map(length, outax)
    aout = fill!(zeros(Union{cubeeltype(agg),Missing}, snew), missing)
    filloutar(aout, convdictall, agg, map(i -> 1:length(i), outaxstat), post)

    YAXArray(collect(outax), aout)
end
function filloutar(aout, convdictall, agg::GroupedOnlineAggregator, s, post)
    for (k, v) in agg.d
        i = CartesianIndices((
            s...,
            map(
                (i, d) -> d[convert(keytype(d), i)]:d[convert(keytype(d), i)],
                k,
                convdictall,
            )...,
        ))
        aout[i.indices...] .= post(v)
    end
end
function filloutar(aout, convdictall, agg, g, post)
    copyto!(aout, post(agg.o))
end
"""
    fittable(tab,o,fitsym;by=(),weight=nothing)

Loops through an iterable table `tab` and thereby fitting an OnlineStat `o` with the values
specified through `fitsym`. Optionally one can specify a field (or tuple) to group by.
Any groupby specifier can either be a symbol denoting the entry to group by or an anynymous
function calculating the group from a table row.

For example the following would caluclate a weighted mean over a cube weighted by grid cell
area and grouped by country and month:

````julia
fittable(iter,WeightedMean,:tair,weight=(i->abs(cosd(i.lat))),by=(i->month(i.time),:country))
````
"""
function fittable(tab, o, fitsym; by = (), weight = nothing, showprog = false)
    agg = TableAggregator(tab, o, fitsym, by = by, weight = weight)
    if showprog
        runfitrows_progress(agg, tab)
    else
        foreach(i -> fitrow!(agg, i), tab)
    end
    agg
end
fittable(tab, o::Type{<:OnlineStat}, fitsym; kwargs...) =
    fittable(tab, o(), fitsym; kwargs...)

getmeter(tab) = getmeter(Base.IteratorSize(tab), tab)
getmeter(::Union{Base.HasLength,Base.HasShape}, tab) = Progress(length(tab))
getmeter(::Base.SizeUnknown, tab) = ProgressUnknown("Rows processed: ")

@noinline function runfitrows_progress(agg, tab)
    p = getmeter(tab)
    every = 0
    for row in tab
        fitrow!(agg, row)
        every += 1
        if every == 100
            next!(p, step = 100)
            every = 0
        end
    end
end

struct collectedValue{V,S,SY}
    value::V
    laststruct::S
end

function Base.getproperty(s::collectedValue{<:Any,<:Any,SY}, sy::Symbol) where {SY}
    if sy == SY
        getfield(s, :value)
    else
        getproperty(getfield(s, :laststruct), sy)
    end
end

function collectval(row::Union{Tuple,Vector}, ::Val{SY}) where {SY}
    nvars = length(row)
    v = ntuple(i -> getfield(row[i], SY), nvars) |> collect
    val = collectedValue{typeof(v),typeof(row[end]),SY}(v, row[end])
end

getpostfunction(s::OnlineStat) = getpostfunction(typeof(s))
getpostfunction(::Type{<:OnlineStat}) = value
function getpostfunction(hist::Union{Ash, HistogramStat, WeightedAdaptiveHist})
    nb = getnbins(hist)
    i->begin
        r = hcat(value(i)...)
        if size(r,1)<nb
            r = vcat(r,zeros(eltype(r),nb-size(r,1),2))
        end
        r
    end
end
getnbins(f::GroupedOnlineAggregator) = getnbins(f.cloneobj)
getnbins(f::TableAggregator) = getnbins(f.o)
getnbins(histogram::HistogramStat) = histogram.k
getnbins(whist::WeightedAdaptiveHist) = whist.alg.b
getnbins(a::Ash) = length(a.density)

fitfun(o) = fitfun(typeof(o))
fitfun(::Type{<:Any}) = fittable

"""
    cubefittable(tab,o,fitsym;post=getpostfunction(o),kwargs...)

Executes [`fittable`](@ref) on the [`@CubeTable`](@ref) `tab` with the
(Weighted-)OnlineStat `o`, looping through the values specified by `fitsym`.
Finally, writes the results from the `TableAggregator` to an output data cube.
"""
function cubefittable(tab, o, fitsym; post = getpostfunction(o), kwargs...)
    agg = fitfun(o)(tab, o, fitsym; showprog = true, kwargs...)
    tooutcube(agg, tab, post)
end

function tupleeltypebyname(::Type{NamedTuple{names,tt}}, s::Symbol) where {names,tt}
    i = findfirst(isequal(s), names)
    fieldtype(tt, i)
end


text:
using Empirikos
using Distributions
using Test
using Roots

biases = 0:0.1:2.0
standard_errors = 0.1:0.1:3.0
αs = 0.01:0.01:0.2

function gaussian_ci2(se; maxbias=0.0, α=0.05)
    level = 1 - α
    maxbias = abs(maxbias)
    rel_bias = maxbias/se
    zz = fzero( z-> cdf(Normal(), rel_bias-z) + cdf(Normal(), -rel_bias-z) +  level -1,
        0, rel_bias - quantile(Normal(),(1- level)/2.1))
    zz*se
end

for bias in biases, se in standard_errors, α in αs
    ci1 = Empirikos.gaussian_ci(se; maxbias=bias, α=α)
    ci2 = gaussian_ci2(se; maxbias=bias, α=α)
    @test ci1 ≈ ci2
end

for se in standard_errors, α in αs
    ci1 = Empirikos.gaussian_ci(se; α=α)
    ci2 = Empirikos.gaussian_ci(se; maxbias =0.0, α=α)
    ci3 = quantile(Normal(), 1-α/2)*se
    @test ci1 ≈ ci2
    @test ci1 ≈ ci3
end

ci = Empirikos.BiasVarianceConfidenceInterval(estimate=1.0, se=1.0)
@test ci.lower ≈ 1 - 1.96 atol = 0.001
@test ci.upper ≈ 1 + 1.96 atol = 0.001

ci_no_se = Empirikos.BiasVarianceConfidenceInterval(estimate=1.0, se=0.0, maxbias=1.0)
@test ci_no_se.lower ≈ 1.0 - 1.0
@test ci_no_se.upper ≈ 1.0 + 1.0

@test Empirikos.gaussian_ci(0.0; maxbias=2.0) == 2.0


text:
tolerance(expected, observed) = maxpcterr(expected, observed) < ε

maxpcterr(expected::Array, observed::Array) = maximum(vec(pcterr(expected, observed)))

maxpcterr{T<:Number}(expected::T, observed::T) = maximum(pcterr(expected, observed))

pcterr(expected, observed) = abs(expected - observed) ./ expected

function hcaulk(a, b)
    z = zeros(size(a, 1))
    display(round([a z b], 4))
    println("")
end

function vcaulk(a, b)
    z = zeros(size(a, 2))'
    display(round([a; z; b], 4))
    println("")
end