text:
include("struct.jl")#
include("input_data/functions.jl")#
include("wind/functions.jl")#
include("wind/wind_module.jl")#
include("eens/functions.jl")#
include("database/functions.jl")
#include("input_data/test_cases.jl")#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#NOTE EENS is set to zero for flexlan !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

################################################################################
############################ Support Functions #################################
################################################################################
# Equipment Cost Calculations are from Techno-Economic Analysis of HVAC, HVDC and OFAC Offshore Wind Power Connections. S. Hardy, K. Van Brusselen, S. Hendrix
# and suppliment data is from North sea grid Annex

################################### Cables #####################################

#calculates the cost of a given hvac cost
#**
function cost_hvac_cable(cbl,ks)
    #cost of losses in the cable
    cbl.costs.rlc=cost_rlc(cbl,ks)
    #cost and size of cable compensation placed on OSS - divide by 2 for each OSS or PCC
    cbl.costs.qc,cbl.reactors=cost_qc_hvac(cbl,ks)
    #cost of switchgear placed on OSS
    cbl.costs.sg=cost_sg_hvac(cbl)
    #capex of cable
    cbl.costs.cpx_p,cbl.costs.cpx_i=capex_cable(cbl,ks.lac)
    #cost of corrective maintenance
    cbl.costs.cm=cost_cm(cbl.costs.cpx_p,ks.opx_c)
    #cost of expected energy not served
    cbl.costs.eens=0#cost_eens(cbl,ks)
    #totals the cable cost
    cbl.costs.ttl=cost_cbl_sum(cbl)
    cbl.costs.grand_ttl=cbl.costs.ttl
    return cbl
end

##MVAC
#calculates the cost of a given mvac cable
#**
function cost_mvac_cable(cbl,ks)
    #cost of losses in the cable
    cbl.costs.rlc=cost_rlc(cbl,ks)
    #cost of switchgear placed on OSS and OWPP
    cbl.costs.sg=cost_sg_hvac(cbl)
    #capex of cable
    cbl.costs.cpx_p,cbl.costs.cpx_i=capex_cable(cbl,ks.lmac)
    #cost of corrective maintenance
    cbl.costs.cm=cost_cm(cbl.costs.cpx_p,ks.opx_c)
    #cost of expected energy not served
    cbl.costs.eens=0#cost_eens(cbl,ks)
    #totals the cable cost
    cbl.costs.ttl=cost_cbl_sum(cbl)
    cbl.costs.grand_ttl=cbl.costs.ttl
    return cbl
end
#calculates the cost of a given hvdc cable
#**
function cost_hvdc_cable(cbl,ks)
    #cost of losses in the hvdc cable
    cbl.costs.rlc=cost_rlc_hvdc(cbl,ks)
    #capex of cable
    cbl.costs.cpx_p,cbl.costs.cpx_i=capex_cable(cbl,ks.ldc)
    #cost of corrective maintenance
    cbl.costs.cm=cost_cm(cbl.costs.cpx_p,ks.opx_c)
    #cost of expected energy not served
    cbl.costs.eens=0#cost_eens(cbl,ks)
    #totals the cable cost
    cbl.costs.ttl=cost_cbl_sum(cbl)
    cbl.costs.grand_ttl=cbl.costs.ttl
    return cbl
end

#CAPEX of cables
#**
function capex_cable(cbl,lay)
    cpx_p=cbl.length*cbl.num*cbl.costs.perkm_cpx#1.775 is installation cost per core per km ref: North Sea Grid Final Report Annex
    k=0.985^(cbl.num)
    cpx_i=cbl.length*cbl.num*lay*k#1.775=DC, 1.905=HVAC, 1.69=MVAC
    return cpx_p,cpx_i
end

#AC RLC Calculation
#**
function cost_rlc(cbl,ks)
    #AC cable loss cost **
    eta=0.994#assumed efficiency of transformers
    A=cbl.mva*eta#eta is the efficiency of a feeder transformer, s power to transmit
    B=cbl.elec.volt*sqrt(3)
    I=A/B#cable current
    R=(cbl.length*cbl.elec.ohm)/(cbl.num)#cable resistance
    #delta is related to wind profile, T_op lifetime hours and E_op cost of energy
    rlc=I^2*R*npv_hours()*ks.E_op*wind_module.wind_profs[cbl.wnd].delta#I^2R losses times cost factors
    return rlc
end

#AC RLC Calculation
#**
function cost_rlc_hvdc(cbl,ks)
    #AC cable loss cost **
    eta=0.9828#assumed efficiency of transformers
    A=cbl.mva*eta#eta is the efficiency of a feeder transformer, s power to transmit
    B=cbl.elec.volt*2
    I=A/B#cable current
    R=(cbl.length*cbl.elec.ohm)/(cbl.num)#cable resistance
    #delta is related to wind profile, T_op lifetime hours and E_op cost of energy
    rlc=I^2*R*npv_hours()*ks.E_op*wind_module.wind_profs[cbl.wnd].delta#I^2R losses times cost factors
    return rlc
end

#HVAC compensation cost oss/PCC
#**
function cost_qc_hvac(cbl,ks)
    f=cbl.elec.freq
    A=cbl.elec.farrad*cbl.length*cbl.num
    Q=2*pi*f*cbl.elec.volt^2*A
    qc,reactors=get_reactors(Q/2)
    opx=(qc*0.0015)*2*npv_years()
    cbl.costs.cm=cbl.costs.cm+opx
    return 2*qc,reactors
end

#corrective maintenance of all equipment
#North Sea Grid OPEX values used
#**
function cost_cm(cpx,k)
    cm=cpx*k*npv_years()
    return cm
end

#Cost of switch gear is added to the cables directly
#**
function cost_sg_hvac(cbl)
    sg=0
    if (cbl.elec.volt==400.0)
        sg=2*4.545
    elseif (cbl.elec.volt==220.0)
        sg=2*3.09
    elseif (cbl.elec.volt==132.0)
        sg=2*2.4
    elseif (cbl.elec.volt==66.0)
        sg=2*1.8
    else
        println(cbl.size)
        println(cbl.num)
        println(cbl.mva)
        println(cbl.costs.ttl)
        println(cbl.elec.mva)
        println("Problem Placing Switchgear")
    end
    #should add OPEX
    opx=(sg*0.007)*npv_years()
    cbl.costs.cm=cbl.costs.cm+opx
    return sg
end

#sums all cable costs and returns the total**
#**
function cost_cbl_sum(cbl)
    ttl=cbl.costs.rlc+cbl.costs.cpx_i+cbl.costs.cpx_p+cbl.costs.cm+cbl.costs.qc+cbl.costs.eens+cbl.costs.sg
    if (ttl==NaN)
        ttl=Inf
    end
    cbl.costs.grand_ttl=ttl
    return ttl
end

################################# Converters ###################################
#totals the cost of HVDC converter
#**
function cost_conv_sum(conv)
    ttl=conv.costs.tlc+conv.costs.cpx+conv.costs.cm+conv.costs.eens
    return ttl
end
#capex of a converter
#**
function capex_hvdc(conv,ks)
    conv.costs.cpx = ks.conv_c*conv.mva
    return conv.costs.cpx
end

################################ Transformers ##################################
##Transformers
#calculates the cost of a given transformer on an OSS
#**
function cost_xfo_oss(xfo,ks)
    #capex oss
    xfo.costs.cpx_p,xfo.costs.cpx_i=capex_oss(xfo,ks)
    #cost of losses
    xfo.costs.tlc=cost_tlc(xfo,ks)
    #corrective maintenance
    xfo.costs.cm=cost_cm(xfo.costs.cpx_p,ks.opx_x)
    #eens calculation
    xfo.costs.eens=0#cost_eens(xfo,ks)
    #totals the xfo cost
    xfo.costs.ttl=cost_xfo_sum(xfo)
    return xfo
end

#calculates the cost of a given transformer at the PCC
#**
function cost_xfo_pcc(xfo,ks)
    #capex pcc
    xfo.costs.cpx_p,xfo.costs.cpx_i=capex_pcc(xfo)
    #cost of losses
    xfo.costs.tlc=cost_tlc(xfo,ks)
    #corrective maintenance
    xfo.costs.cm=cost_cm(xfo.costs.cpx_p,ks.opx_x)
    #eens calculation
    xfo.costs.eens=0#cost_eens(xfo,ks)
    #totals the xfo cost
    xfo.costs.ttl=cost_xfo_sum(xfo)
    return xfo
end

#Offshore transformer CAPEX
#**
function capex_oss(xfo,ks)
    A=(1+ks.dc*(xfo.num-2))
    B=(ks.f_ct+ks.p_ct)
    cpx=A*B*xfo.num*xfo.elec.mva
    return cpx/2,cpx/2#cost is split 50-50 over procurement and installation
end

#PCC transformer CAPEX
#**
function capex_pcc(xfo)
    cpx=xfo.num*0.0071*(xfo.elec.mva)^1.05
    return cpx*0.8,cpx*0.2#cost is split 80-20 over procurement and installation
end

#xfo Losses Calculation
#power factor of 1 assumed
#**
function cost_tlc(xfo,ks)
    pf=1
    tlc=xfo.mva*pf*(1-xfo.eta)*npv_hours()*ks.E_op*wind_module.wind_profs[xfo.wnd].delta
    return tlc
end

#sums all xfo costs and returns the total**
#**
function cost_xfo_sum(xfo)
    ttl=xfo.costs.tlc+xfo.costs.cpx_p+xfo.costs.cpx_i+xfo.costs.cm+xfo.costs.eens
    return ttl
end

################################### platform ###################################

#depth ref A review of foundations of offshore wind energy convertors: Current status and future perspectives
#AC platform incremental capex
#**
#depth scale factors from Europe's onshore and offshore wind energy potential An assessment of environmental and economic constraints
#EEA Technical report No 6/2009, ISSN 1725-2237
#10–20 20–30 30–40 40–50
#1.000 1.067 1.237 1.396
#multiple is approximated by 0.0136*depth_in_meters + 0.7676
function capex_plat_ac(plat,ks)
    plat.costs.cpx=ks.pac_e*plat.mva#*(plat.depth/15)#15m set as base depth - original with ks.pac_e=16
    #plat.costs.cpx=ks.pac_e*plat.mva^1.3#new to reflect better the base cost
    return plat.costs.cpx
end

#depth ref A review of foundations of offshore wind energy convertors: Current status and future perspectives
#DC platform incremental capex
#**
function capex_plat_dc(plat,ks)
    plat_base=platform()
    plat_base.mva=plat.mva
    plat_base.wnd=plat.wnd
    plat.costs.cpx=(ks.pdc_g*plat.mva)-capex_plat_ac(plat,ks)#*(plat.depth/15)#15m set as base depth
    return plat.costs.cpx
end

############################### Net Present value ##############################
#Changes to a net present value
#**
function npv_years()
    #1+1/(1.04)^1+1/(1.04)^2+...+1/(1.04)^25
    return 16.62
end

#Changes to a net present value
#**
function npv_hours()
    #16.62*365.25*24=145690.91
    return 145691
end


text:
typealias Chars Union(UInt8,UInt16,UInt32,Char)

immutable String{E<:Encoding}
  ptr::Ptr{E} # pointer to code units of encoding E
  len::Int    # number of code units
end

# core access method; unsafe if !(0 < i < s.len)
getcodeunit{E}(s::String{E},i::Integer=1) = unsafe_load(convert(Ptr{codeunit(E)},s.ptr), i)
codeunits(s::String) = s.len
Base.sizeof{E}(s::String{E}) = codeunits(s) * sizeof(E)

# constructors
String{E<:Encoding}(::Type{E}, ptr::Ptr, len::Int) = String(convert(Ptr{E},ptr), len)

# function String{E<:Encoding}(::Type{E}, ptr::Ptr)
#     # look for NUL terminate
# end

emtpy{E<:Encoding}(::Type{E}) = String(convert(Ptr{E},C_NULL), 0)

function String{E<:Encoding,C<:Chars}(::Type{E}, s::Vector{C})
    # isvalid(E, s) || throw(ArgumentError("vector is not valid $E"))
    ptr, len = storebytes!(s)
    return String{E}(convert(Ptr{E},ptr),div(len,sizeof(E)))
end

# TODO: remove when no longer needed
String(s::ASCIIString) = String(ASCII,s.data)
String(s::UTF8String)  = String(UTF8, s.data)
String(s::UTF16String) = String(UTF16,s.data)
String(s::UTF32String) = String(UTF32,s.data)

# indexing by code point
# DirectIndexedEncoding: code unit == code point
Base.getindex{E<:DirectIndexedEncoding}(s::String{E}, i::Integer=1) = (c = getcodeunit(s, i); isvalid(E, c) ? Char(c) : '\uffd')
Base.endof{E<:DirectIndexedEncoding}(s::String{E}) = codeunits(s)
Base.length{E<:DirectIndexedEncoding}(s::String{E}) = codeunits(s)

# substring
function Base.getindex{E<:DirectIndexedEncoding}(s::String{E}, r::UnitRange{Int})
    n = length(r)
    (0 < first(r) <= endof(s) && last(r) <= endof(s)) || throw(BoundsError())
    isempty(r) && return empty(E)

    if n < div(endof(s),4) || n < 16 # < 25% of original string size or really small
        # just make a copy
        ptr, len = storebytes!(convert(Ptr{UInt8},s.ptr), n, first(r)-1)
        return String(E, ptr, len)
    else
        # share data with original string
        return String(E, s.ptr + UInt(first(r)-1), n)
    end
end

function Base.getindex(s::String,i::Integer=1)
    error("must implement getindex(s::$(typeof(s)), i) --> Char representing the ith code point in s")
end
Base.endof(s::String) = error("must implement endof(s::$(typeof(s))) --> # of code points in s")
function Base.length(s::String)
    i = start(s)
    done(s,i) && return 0
    n = 1
    while true
        c, j = next(s,i)
        done(s,j) && return n
        n += 1
        i = j
    end
end

Base.start(::String) = 1
Base.done(s::String, i::Integer) = i > endof(s)
Base.next(s::String, i::Integer) = (getindex(s, i), i+1)

function =={E}(a::String{E}, b::String{E})
    na = sizeof(a)
    nb = sizeof(b)
    na != nb && return false
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              a.ptr, b.ptr, min(na,nb))
    return c == 0
end

function Base.show(io::IO, x::String)
    if codeunits(x) == 0
        print(io, '"', '"')
    else
        print(io, '"')
        for c in x
            print(io, c)
        end
        print(io, Char('"'))
    end
    return
end

# concatenation
function Base.string{E}(strs::String{E}...)
    length(strs) == 1 && return strs[1]
    @inbounds begin
    N = length(strs)
    n = 0
    for i = 1:N
        n += strs[i].len
    end
    ensureroom!(n)
    # unsafe_copy!(dest::Ptr{T}, src::Ptr{T}, N)
    starting_ptr = pointer(POOL.pool[POOL.ind])+UInt(POOL.pos)
    ptr = starting_ptr
    for i = 1:N
        unsafe_copy!(ptr, convert(Ptr{UInt8},strs[i].ptr), sizeof(E) * strs[i].len)
        ptr += strs[i].len
    end
    end
    return String{E}(starting_ptr,n)
end

ext: jl


text:
using LinearAlgebra
using Dates
using Plots

include("qr_factorization.jl")
include("truncated_svd.jl")

"""
verifies <num_tests> times that ||Q*R - A||_2 <= min_error
"""
function test_qr_A(num_tests, min_error= 1e-13)
    errors = 0

    for i = 1:num_tests

        #generating random matrix
        rows = rand(10:100)
        columns = rand(10:rows)
        r = rand(10:columns)
        #using SVD to obtain a rank r matrix
        A = approx_matrix_svd(rand(rows,columns),r)

        #calculating QR factorization
        QR, v, u = allocate_matrices(A)
        Q = rand(rows, columns)
        qr_factorization!(A,QR,v,u)
        get_Q!(QR,Q)
        R = triu(QR)[1:columns,:]

        #checking ||QR - A||
        if (norm(Q*R - A) > min_error)
            errors += 1
        end

    end
    return errors
end

"""
verifies <num_tests> times that ||R - R_j||_2 <= min_error where R is the R
computed by our QR factorization while R_j is the R computed by the
JULIA's implementation
"""
function test_qr_R(num_tests, min_error = 1e-10)
    errors = 0

    for i = 1:num_tests

        #generating random matrix
        rows = rand(10:100)
        columns = rand(10:rows)
        r = rand(10:columns)
        A = approx_matrix_svd(rand(rows,columns),r)

        #calculating QR factorization
        QR, v, u = allocate_matrices(A)
        qr_factorization!(A,QR,v,u)
        Q_j,R_j = qr(A)
        R = triu(QR)[1:columns,:]

        #checking ||R - R_j||
        if (norm(R - R_j) > min_error)
            errors += 1
        end

    end
    return errors
end

"""
verifies Q_t_times_A() function
"""
function test_qr_multiplication(num_tests, min_error = 1e-10)
    errors = 0

    for i = 1:num_tests
        #generating random matrix
        rows = rand(10:100)
        columns = rand(10:rows)
        r = rand(10:columns)
        A = approx_matrix_svd(rand(rows,columns),r)

        #calculating QR factorization
        QR, v, u = allocate_matrices(A)
        qr_factorization!(A,QR,v,u)
        Q_j,R_j = qr(A)
        W = zeros(rows,columns)
        Q_t_times_A!(QR,A,W)
        W_j = Q_j'*A

        #checking ||W - W_j||
        if (norm(W - W_j) > min_error)
            errors += 1
        end

    end
    return errors
end

function qr_test(A)
    QR, v, u = allocate_matrices(A)
    qr_factorization!(A,QR,v,u)
end

"""
Draws a plot that compares the time (in seconds) required by our implementation and the one built-in in Julia.

The algorithm generates Int((max_dimension - min_dimension) / step) + 1 random matrices of increasing dimension.

At step k, the matrix shape will be (min dimension+((k - 1) * step) x min dimension+((k - 1) * step)). The algorithm terminates when
the matrix shapes reach max dimension.

"""
function test_time_qr(min_dimension, max_dimension; step = 15)
    plotly()
    dim = Int((max_dimension - min_dimension) / step) + 1
    x = zeros(dim)
    y = zeros(dim,2)
    for i = 0:dim-1
        x[i+1] = i * step + min_dimension
        A = rand(Int(x[i+1]),Int(x[i+1]))
        y[i+1,1] = @elapsed qr_test(A)
        y[i+1,2] = @elapsed qr(A)
    end
    plot(x, y, label = ["Our QR" "Julia's QR"], lw = 3)
end


text:
using DataFrames
using CSV

data = readlines(joinpath(@__DIR__, "original.txt"))
data

df = DataFrame(implementation = String[], max_edit = Int[], prefix_length = Int[],
    dict_name = String[], verbosity = String[], type = String[],
    res = Int[], min_val = Float64[], med_val = Float64[], avg_val = Float64[])

for s in data
    if startswith(s, "Precalculation")
        m = match(r"MaxEditDistance=([0-9]+) prefixLength=([0-9]+) dict=(.*)$", s)
        println("$(m[1])\t$(m[2])\t$(m[3])")
        global max_edit = parse(Int, m[1])
        global prefix_length = parse(Int, m[2])
        global dict_name = m[3]
    end

    # Lookup instance 1 results 0.000063ms/op verbosity=Closest query=exact
    if startswith(s, "Lookup")
        m = match(r"Lookup ([\w]*)\s*([\d\,]+)\s*results\s*([\d\.]+)ms/op\s*verbosity=([\w]*)\s*query=([\w\-]*)$", s)
        implementation = "CS$(m[1])"
        res = parse(Int, replace(m[2], ','=>""))
        avg_val = parse(Float64, m[3])
        @show m[3], parse(Float64, m[3])
        verbosity = lowercase(m[4])
        query = m[5]

        # pushing avg_val to other columns, which is not exactly correct
        push!(df, (implementation, max_edit, prefix_length, dict_name, verbosity, query, res,
            avg_val, avg_val, avg_val))
    end
end

CSV.write(joinpath(@__DIR__, "original_bench.csv"), df)

##################
# Sanity check

df2 = CSV.read(joinpath(@__DIR__, "original_bench.csv"))

df == df2   # true


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Define the type for martingale increments
################################################################################

immutable Mart

    ## Members
    mart::Array{Float64} # martingale increments
    method::String # slow, neighbours or discrete distribution
    pnum::Int64 # number of positions
    anum::Int64 # number of actions
    pathnum::Int64 # number of paths
    subsimnum::Int64 # number of subsimulation for each path
    tnum::Int64 ## number of time points
    fullcontrol::Bool

    ## Using only row rearrangement

    function Mart(bellman::Bellman, path::Path, subsim::Array{Float64, 5},
                  subsim_weight::Array{Float64, 1})

        method = "row-rearrange"

        ## Extract sizes from CSS
        pnum = bellman.pnum
        anum = bellman.anum
        pathnum = path.pathnum
        subsimnum = size(subsim)[4]
        tnum = bellman.tnum
        fullcontrol = bellman.fullcontrol

        ## Container
        result = zeros(Array{Float64}(tnum - 1, pnum, pathnum))
        path_state = Array{Float64}(bellman.snum)

        print("Martingales: ")
        ##  Compute martingales
        for t = 1:(tnum - 1)
            print(t, ".")
            for pos = 1:pnum
                for p = 1:pathnum
                    path_state = path.sample[t, p, :]
                    ## Finding the average
                    for h = 1:subsimnum
                        result[t, pos, p] += subsim_weight[h] *
                            maximum(bellman.value[:, :, pos, t + 1] *
                                    subsim[:, :, p, h, t] * path_state)
                    end
                    ## Subtracting the path realization
                    path_state = path.sample[t + 1, p, :]
                    result[t, pos, p] -=
                        maximum(bellman.value[:, :, pos, t + 1] * path_state)
                end
            end
        end

        if fullcontrol
            print("End.\n")
            return new(result, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)
        else # partial control of positions
            ## Container
            result2 = Array{Float64}(tnum - 1, pnum, anum, pathnum)
            prob_weight = Array{Float64}(pnum)
            ## Adjust the martingales for partial control
            for t = 1:(tnum - 1)
                for pos = 1:pnum
                    for a = 1:anum
                        prob_weight = bellman.control[pos, a, :]
                        result2[t, pos, a, :] = transpose(prob_weight) * result[t, :, :]
                    end
                end
            end

            print("End.\n")
            return new(result2, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)

        end

    end

    ############################################################################

    ## Using nearest neighbours

    function Mart(bellman::Bellman, path::Path, subsim::Array{Float64, 5},
                  subsim_weight::Array{Float64, 1}, search::Function,
                  grid::Array{Float64, 2})

        method = "neighbors"

        ## Extract sizes from CSS
        pnum = bellman.pnum
        anum = bellman.anum
        pathnum = path.pathnum
        subsimnum = size(subsim)[4]
        tnum = bellman.tnum
        fullcontrol = bellman.fullcontrol

        ## Container
        result = zeros(Array{Float64}(tnum - 1, pnum, pathnum))
        path_state = Array{Float64}(bellman.snum)
        disturb_state = Array{Float64}(bellman.snum)

        ## Tree
        tgrid = transpose(grid)
        tree = search(tgrid)

        print("Martingales: ")
        ##  Compute martingales
        for t = 1:(tnum - 1)
            print(t, ".")
            for pos = 1:pnum
                for p = 1:pathnum
                    path_state = path.sample[t, p, :]
                    ## Finding the average
                    for h = 1:subsimnum
                        disturb_state = subsim[:, :, p, h, t] * path_state
                        nn = knn(tree, disturb_state, 1, true)[1]
                        result[t, pos, p] += subsim_weight[h] *
                            (bellman.value[nn, :, pos, t + 1] * disturb_state)[1]
                    end
                    ## Subtracting the path realization
                    path_state = path.sample[t + 1, p, :]
                    nn = knn(tree, path_state, 1 ,true)[1]
                    result[t, pos, p] -=
                        (bellman.value[nn, :, pos, t + 1] * path_state)[1]
                end
            end
        end

        if fullcontrol
            print("End.\n")
            return new(result, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)
        else # partial control of positions
            ## Container
            result2 = Array{Float64}(tnum - 1, pnum, anum, pathnum)
            prob_weight = Array{Float64}(pnum)
            ## Adjust the martingales for partial control
            for t = 1:(tnum - 1)
                for pos = 1:pnum
                    for a = 1:anum
                        prob_weight = bellman.control[pos, a, :]
                        result2[t, pos, a, :] = transpose(prob_weight) * result[t, :, :]
                    end
                end
            end

            print("End.\n")
            return new(result2, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)

        end

    end

       ############################################################################

       ## Finite distribution case

       function Mart(bellman::Bellman, path::Path, path_disturb::Array{Float64, 4},
                     search::Function, grid::Array{Float64, 2})

        method = "finite-distribution"

        ## Extract sizes from CSS
        pnum = bellman.pnum
        anum = bellman.anum
        pathnum = path.pathnum
        subsimnum = 0
        tnum = bellman.tnum
        fullcontrol = bellman.fullcontrol

        ## Container
        result = zeros(Array{Float64}(tnum - 1, pnum, pathnum))
        path_state = Array{Float64}(bellman.snum)
        host = Array{Float64}(bellman.snum)
        disturb_host = Array{Float64}(bellman.snum)

        ## Tree
        tgrid = transpose(grid)
        tree = search(tgrid)

        print("Martingales: ")
        ##  Compute martingales
        for t = 1:(tnum - 1)
            print(t, ".")
            for p = 1:pathnum
                path_state = path.sample[t, p, :]
                nn1 = knn(tree, path_state, 1, true)[1]
                host = transpose(grid[nn1, :])
                disturb_host = path_disturb[:, :, t, p] * host
                nn2 = knn(tree, disturb_host, 1, true)[1][1]
                for pos = 1:pnum
                    ## Fining the average
                    result[t, pos, p] = (bellman.evalue[nn1, :, pos, t] * host)[1]
                    ## Subtracting the path realization
                    result[t, pos, p] -=
                        (bellman.value[nn2, :, pos, t + 1] * disturb_host)[1]

                end
            end
        end

        if fullcontrol
            print("End.\n")
            return new(result, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)
        else # partial control of positions
            ## Container
            result2 = Array{Float64}(tnum - 1, pnum, anum, pathnum)
            prob_weight = Array{Float64}(pnum)
            ## Adjust the martingales for partial control
            for t = 1:(tnum - 1)
                for pos = 1:pnum
                    for a = 1:anum
                        prob_weight = bellman.control[pos, a, :]
                        result2[t, pos, a, :] = transpose(prob_weight) * result[t, :, :]
                    end
                end
            end

            print("End.\n")
            return new(result2, method, pnum, anum, pathnum, subsimnum, tnum,
                       fullcontrol)

        end

    end

end


text:
function _end(obj::FieldVectorPreservingVisitor)
    return jcall(obj, "end", FieldElement, ())
end

function start(obj::FieldVectorPreservingVisitor, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "start", void, (jint, jint, jint), arg0, arg1, arg2)
end

function visit(obj::FieldVectorPreservingVisitor, arg0::jint, arg1::FieldElement)
    return jcall(obj, "visit", void, (jint, FieldElement), arg0, arg1)
end


text:
using Test
using IncrementalInference
using Random

##

@testset "test basic multihypo" begin

## A simple multihypo example

Random.seed!(42) # The answer to reproducable noise

fg = LocalDFG(solverParams=SolverParams(graphinit=false, gibbsIters=5, spreadNH=5.0))

pRight = 0.99
pWrong = 0.01
pr_noise = 0.01
od_noise = 0.1
lm_noise = 0.01

# true positions
# x0 at 0
# x1 at 1
# l1 at 1
# l2 at 2

#x0 prior
addVariable!(fg, :x0, ContinuousScalar)
prpo = Normal(0.0, pr_noise)
addFactor!(fg, [:x0], Prior(Normal(rand(prpo), pr_noise)))

#l1 and l2
addVariable!(fg, :l1, ContinuousScalar, tags=[:LANDMARK])
addVariable!(fg, :l2, ContinuousScalar, tags=[:LANDMARK])

#x0 to l1 or l2
p2ln = Normal(1.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x0; :l1; :l2], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x0; :l1], p2p) #this one used for sanity check

#x0 to x1
addVariable!(fg, :x1, ContinuousScalar)
pp = Normal(1.0, od_noise)
addFactor!(fg, [:x0,:x1], LinearRelative(Normal(rand(pp), od_noise)))

#x1 to l1 or l2
p2ln = Normal(0.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l1; :l2], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x1; :l1], p2p) #this one used for sanity check

#x1 to l2 or l1
p2ln = Normal(1.0, lm_noise)
p2p = LinearRelative(Normal(rand(p2ln), lm_noise))
addFactor!(fg, [:x1; :l2; :l1], p2p, multihypo = [1, pRight, pWrong])
# addFactor!(fg, [:x1; :l2], p2p) #this one used for sanity check

##

# prescribe an elimination order to get a single clique
eo = [:l2,:x1,:x0,:l1]
# fg.solverParams.graphinit=true
smtasks = Task[]
tree = solveTree!(fg, eliminationOrder=eo) #, smtasks=smtasks, recordcliqs=ls(fg));

# hists = fetchCliqHistoryAll!(smtasks)

# plotKDE(fg, ls(fg))

##

@test isapprox(DFG.getPPESuggested(fg, :x0)[], 0, atol = 0.2)
@test isapprox(DFG.getPPESuggested(fg, :x1)[], 1, atol = 0.2)
@test isapprox(DFG.getPPESuggested(fg, :l1)[], 1, atol = 0.2)

L2 = getBelief(fg, :l2)
npts = length(getPoints(L2))
pts = [2.0.+0.1*randn(1) for _ in 1:npts]
L2_ = manikde!(ContinuousScalar, pts)

# test that there is at least a mode present
@test mmd(L2_, L2, ContinuousScalar) < 1e-3
# @test isapprox(DFG.getPPESuggested(fg, :l2)[], 2, atol = 0.2)

##

end

@testset "test multihypo chain example (see #462)..." begin

##

l1 = -10.0
l2 = +10.0
lnoise = 1.0
x1 = 0
x2 = 0
Onoise = 0.1

fg = initfg()

addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :l1, ContinuousScalar)
addVariable!(fg, :l1_0, ContinuousScalar)
addVariable!(fg, :l2, ContinuousScalar)
addVariable!(fg, :l2_0, ContinuousScalar)

# priors on two landmarks only
addFactor!(fg, [:l1], Prior(Normal(l1, lnoise)))
addFactor!(fg, [:l2], Prior(Normal(l2, lnoise)))

# relative constraints
addFactor!(fg, [:x1;:l1;:l1_0], LinearRelative(Normal(l1-x1, lnoise)), multihypo=[1;1/2;1/2])
addFactor!(fg, [:x2;:l2;:l2_0], LinearRelative(Normal(l2-x2, lnoise)), multihypo=[1;1/2;1/2])
addFactor!(fg, [:x1;:x2], LinearRelative(Normal(0, Onoise)))

tree = solveTree!(fg)

# drawTree(tree, show=true)

# expect x1 x2 to have at least one mode at 0

@test getPPE(fg, :x1).suggested[1] - x1 |> abs < 1.2
@test getPPE(fg, :x2).suggested[1] - x2 |> abs < 1.2

@test getPPE(fg, :l1).suggested[1] - l1 |> abs < 1.2
@test getPPE(fg, :l2).suggested[1] - l2 |> abs < 1.2

# l1_0, l2_0 should be nearby around l1 and l2, but cannot confirm 100%
@test getPPE(fg, :l1_0).suggested[1] - l1 |> abs < 10
@test getPPE(fg, :l2_0).suggested[1] - l2 |> abs < 10

##

end

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)
#
# plotKDE(fg, [:l1;:l2])
# plotKDE(fg, [:l1_0;:l2_0])
# plotKDE(fg, [:x1;:x2])

#


text:
"""
    flatten([eltype=Float64], x)

Returns a "flattened" representation of `x` as a vector of real numbers, and a function
`unflatten` that takes a vector of reals of the same length and returns an object of the
same type as `x`.

`unflatten` is the inverse of `flatten`, so
```julia
julia> x = (randn(5), 5.0, (a=5.0, b=randn(2, 3)));

julia> v, unflatten = flatten(x);

julia> x == unflatten(v)
true
```
"""
function flatten end

flatten(x) = flatten(Float64, x)

function flatten(::Type{T}, x::Integer) where T<:Real
    v = T[]
    unflatten_to_Integer(v::Vector{T}) = x
    return v, unflatten_to_Integer
end

function flatten(::Type{T}, x::R) where {T<:Real, R<:Real}
    v = T[x]
    unflatten_to_Real(v::Vector{T}) = convert(R, only(v))
    return v, unflatten_to_Real
end

flatten(::Type{T}, x::Vector{R}) where {T<:Real, R<:Real} = (Vector{T}(x), Vector{R})

function flatten(::Type{T}, x::AbstractVector) where T<:Real
    x_vecs_and_backs = map(val -> flatten(T, val), x)
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = _cumsum(map(length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return oftype(x, x_Vec)
    end
    return reduce(vcat, x_vecs), Vector_from_vec
end

function flatten(::Type{T}, x::AbstractArray) where T<:Real
    x_vec, from_vec = flatten(T, vec(x))
    Array_from_vec(x_vec) = oftype(x, reshape(from_vec(x_vec), size(x)))
    return x_vec, Array_from_vec
end

function flatten(::Type{T}, x::Tuple) where T<:Real
    x_vecs_and_backs = map(val -> flatten(T, val), x)
    x_vecs, x_backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    lengths = map(length, x_vecs)
    sz = _cumsum(lengths)
    function unflatten_to_Tuple(v::Vector{T})
        map(x_backs, lengths, sz) do x_back, l, s
            return x_back(v[s - l + 1:s])
        end
    end
    return reduce(vcat, x_vecs), unflatten_to_Tuple
end

function flatten(::Type{T}, x::NamedTuple) where T<:Real
    x_vec, unflatten = flatten(T, values(x))
    function unflatten_to_NamedTuple(v::Vector{T})
        v_vec_vec = unflatten(v)
        return typeof(x)(v_vec_vec)
    end
    return x_vec, unflatten_to_NamedTuple
end

function flatten(::Type{T}, d::Dict) where T<:Real
    d_vec, unflatten = flatten(T, collect(values(d)))
    function unflatten_to_Dict(v::Vector{T})
        v_vec_vec = unflatten(v)
        return Dict(key => v_vec_vec[n] for (n, key) in enumerate(keys(d)))
    end
    return d_vec, unflatten_to_Dict
end

_cumsum(x) = cumsum(x)
if VERSION < v"1.5"
    _cumsum(x::Tuple) = (_cumsum(collect(x))..., )
end


text:
# building fixed rate bonds
using QuantLib
using Dates
using Printf

function build_bonds(bond_mats::Vector{Date}, bond_rates::Vector{Float64}, tenor::QuantLib.Time.TenorPeriod, conv::QuantLib.Time.BusinessDayConvention,
                    rule::QuantLib.Time.DateGenerationRule, calendar::QuantLib.Time.BusinessCalendar, dc::QuantLib.Time.DayCount, freq::QuantLib.Time.Frequency, issue_date::Date)
  bonds = Vector{FixedRateBondHelper}(undef, length(bond_mats))
  pricing_engine = DiscountingBondEngine()

  for i =1:length(bond_mats)
    term_date = bond_mats[i]
    # rate = bond_rates[i] / 100
    rate = bond_rates[i]
    sched = QuantLib.Time.Schedule(issue_date, term_date, tenor, conv, conv, rule, true)
    bond_help = FixedRateBondHelper(Quote(100.0), FixedRateBond(0, 100.0, sched, rate, dc, conv, 100.0, issue_date, calendar, pricing_engine))
    bonds[i] = bond_help
  end

  return bonds
end

function build_depos(depo_quotes::Vector{Float64}, depo_tenors::Vector{P}, dc::QuantLib.Time.DayCount, conv::QuantLib.Time.BusinessDayConvention, calendar::QuantLib.Time.BusinessCalendar, fixing_days::Int) where {P <: Dates.Period}
  depos = Vector{DepositRateHelper}(undef, length(depo_quotes))
  for i = 1:length(depo_quotes)
    depo_quote = Quote(depo_quotes[i])
    depo_tenor = QuantLib.Time.TenorPeriod(depo_tenors[i])
    depo = DepositRateHelper(depo_quote, depo_tenor, fixing_days, calendar, conv, true, dc)
    depos[i] = depo
  end

  return depos
end

function build_swaps(swap_quotes::Vector{Float64}, swap_tenors::Vector{P}, fixed_dc::QuantLib.Time.DayCount, fixed_conv::QuantLib.Time.BusinessDayConvention, calendar::QuantLib.Time.BusinessCalendar, fixed_freq::QuantLib.Time.Frequency, float_index::IborIndex, forward_start::Int) where {P <: Dates.Period}
  forward_start_period = Dates.Day(forward_start)
  swaps = Vector{SwapRateHelper}(undef, length(swap_quotes))
  pricing_engine = DiscountingSwapEngine()

  for i = 1:length(swap_quotes)
    swaps[i] = SwapRateHelper(swap_quotes[i], swap_tenors[i], calendar, fixed_freq, fixed_conv, fixed_dc, float_index, 0.0, forward_start_period, pricing_engine)
  end

  return swaps
end

function get_npvs(bonds, issue_date, calendar, dc, freq)
  pricing_engine = DiscountingBondEngine()
  rate_quote = Quote(0.05)
  compounding = CompoundedCompounding()
  yts = FlatForwardTermStructure(0, issue_date, calendar, rate_quote, dc, compounding, freq)
  npvs = zeros(length(bonds))
  for i=1:length(bonds)
    npvs[i] = calculate(pricing_engine, yts, bonds[i])
  end

  return npvs
end

function par_rate(yts::YieldTermStructure, dates::Vector{Date}, dc::QuantLib.Time.DayCount)
  sum = 0.0
  for i = 2:length(dates)
    dt = QuantLib.Time.year_fraction(dc, dates[i - 1], dates[i])
    sum += discount(yts, dates[i]) * dt
  end

  result = discount(yts, dates[1]) - discount(yts, dates[end])
  return result/sum
end

function setup()
  today = now()
  issue_date = Date(Dates.Year(today), Dates.Month(today), Dates.Day(today))
  bond_mats = [issue_date + Dates.Year(i) for i in range(2, step=2, length=15)]
  # bond_rates = [5.75, 6.0, 6.25, 6.5, 6.75, 6.80, 7.00, 7.1, 7.15, 7.2]
  bond_rates = [0.0200, 0.0225, 0.0250, 0.0275, 0.0300, 0.0325, 0.0350, 0.0375, 0.0400, 0.0425, 0.0450, 0.0475, 0.0500, 0.0525, 0.0550]
  set_eval_date!(settings, issue_date)

  freq = QuantLib.Time.Annual()
  tenor = QuantLib.Time.TenorPeriod(freq)
  conv = QuantLib.Time.Unadjusted()
  rule = QuantLib.Time.DateGenerationBackwards()
  calendar = QuantLib.Time.USGovernmentBondCalendar()
  dc = QuantLib.Time.SimpleDayCount()
  bonds = build_bonds(bond_mats, bond_rates, tenor, conv, rule, calendar, dc, freq, issue_date)

  return issue_date, bonds, dc, calendar
end

function piecewise_yld_curve()
  issue_date, bonds, dc, calendar = setup()

  interp = QuantLib.Math.LogLinear()
  trait = Discount()
  bootstrap = IterativeBootstrap()

  yts = PiecewiseYieldCurve(issue_date, bonds, dc, interp, trait, 0.00000000001, bootstrap)

  solver = QuantLib.Math.BrentSolver()
  solver2 = QuantLib.Math.FiniteDifferenceNewtonSafe()
  calculate!(IterativeBootstrap(), yts, solver2, solver)

  # println(yts.data)

  for bond in bonds
    date_vec = Vector{Date}(undef, length(bond.cashflows.coupons) + 1)
    date_vec[1] = issue_date
    for (i, cf) in enumerate(bond.cashflows.coupons)
      date_vec[i+1] = date(cf)
    end
    par = par_rate(yts, date_vec, dc)
    println(100.0 * par)
  end
end

  # npvs = get_npvs(bonds, issue_date, calendar, dc, freq)
function fitted_bond_curve_exp()
  issue_date, bonds, dc, calendar = setup()

  esf = ExponentialSplinesFitting(true, length(bonds))
  fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, esf, 1e-10, 5000, 1.0)
  initialize!(fitted_curve)
  calculate!(fitted_curve)

  println(fitted_curve.fittingMethod.commons.minimumCostValue)
  println(fitted_curve.fittingMethod.commons.numberOfIterations)
  println(fitted_curve.fittingMethod.commons.guessSolution)

  disc = discount(fitted_curve, 1.0)
  println(disc)
end

function fitted_bond_curve_simp()
  issue_date, bonds, dc, calendar = setup()

  spf = SimplePolynomialFitting(true, 3, length(bonds))
  fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, spf, 1e-10, 5000, 1.0)
  initialize!(fitted_curve)
  calculate!(fitted_curve)

  # println(fitted_curve.fittingMethod.minimumCostValue)
  println(fitted_curve.fittingMethod.commons.numberOfIterations)
  # println(fitted_curve.fittingMethod.guessSolution)

  disc = discount(fitted_curve, 1.0)
  println(disc)
end

function fitted_bond_curve_ns()
  issue_date, bonds, dc, calendar = setup()

  nsf = NelsonSiegelFitting(length(bonds))
  fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, nsf, 1e-10, 5000, 1.0)
  initialize!(fitted_curve)
  calculate!(fitted_curve)

  # println(fitted_curve.fittingMethod.minimumCostValue)
  println(fitted_curve.fittingMethod.commons.numberOfIterations)
  # println(fitted_curve.fittingMethod.guessSolution)

  disc = discount(fitted_curve, 1.0)
  println(disc)
end

function fitted_bond_curve_sven()
  issue_date, bonds, dc, calendar = setup()

  sf = SvenssonFitting(length(bonds))
  fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, sf, 1e-10, 5000, 1.0)
  initialize!(fitted_curve)
  calculate!(fitted_curve)

  # println(fitted_curve.fittingMethod.minimumCostValue)
  println(fitted_curve.fittingMethod.commons.numberOfIterations)
  # println(fitted_curve.fittingMethod.guessSolution)

  disc = discount(fitted_curve, 1.0)
  println(disc)
end

function fitted_bond_curve_cbspline()
  issue_date, bonds, dc, calendar = setup()

  knots = [-30.0, -20.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0]
  cbsf = CubicBSplinesFitting(true, knots, length(bonds))
  fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, cbsf, 1e-10, 5000, 1.0)
  initialize!(fitted_curve)
  calculate!(fitted_curve)

  println(fitted_curve.fittingMethod.commons.numberOfIterations)
  disc = discount(fitted_curve, 1.0)
  println(disc)
end

function fitted_bond_curve_all()
  tic()
  issue_date, bonds, dc, calendar = setup()

  esf = ExponentialSplinesFitting(true, length(bonds))
  esf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, esf, 1e-10, 5000, 1.0)
  initialize!(esf_fitted_curve)
  calculate!(esf_fitted_curve)

  spf = SimplePolynomialFitting(true, 3, length(bonds))
  spf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, spf, 1e-10, 5000, 1.0)
  initialize!(spf_fitted_curve)
  calculate!(spf_fitted_curve)

  nsf = NelsonSiegelFitting(length(bonds))
  nsf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, nsf, 1e-10, 5000, 1.0)
  initialize!(nsf_fitted_curve)
  calculate!(nsf_fitted_curve)

  knots = [-30.0, -20.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0]
  cbsf = CubicBSplinesFitting(true, knots, length(bonds))
  cbsf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, cbsf, 1e-10, 5000, 1.0)
  initialize!(cbsf_fitted_curve)
  calculate!(cbsf_fitted_curve)

  sf = SvenssonFitting(length(bonds))
  sf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, sf, 1e-10, 5000, 1.0)
  initialize!(sf_fitted_curve)
  calculate!(sf_fitted_curve)
  toc()

  println("Exponential Splines: $(esf_fitted_curve.fittingMethod.commons.numberOfIterations)")
  println("Simple Polynomial: $(spf_fitted_curve.fittingMethod.commons.numberOfIterations)")
  println("Nelson Siegel: $(nsf_fitted_curve.fittingMethod.commons.numberOfIterations)")
  println("Cubic B-Splines: $(cbsf_fitted_curve.fittingMethod.commons.numberOfIterations)")
  println("Svensson Fitting: $(sf_fitted_curve.fittingMethod.commons.numberOfIterations)")
end

function generate_floatingrate_bond(yts::YieldTermStructure, disc_yts::YieldTermStructure, cap_vol::OptionletVolatilityStructure)
  # Floating Rate bond
  settlement_days = 3
  face_amount = 100.0
  fb_issue_date = Date(2005, 10, 21)
  bond_engine = DiscountingBondEngine(disc_yts)
  fb_dc = QuantLib.Time.Actual360()
  conv = QuantLib.Time.ModifiedFollowing()
  fb_schedule = QuantLib.Time.Schedule(Date(2005, 10, 21), Date(2010, 10, 21), QuantLib.Time.TenorPeriod(QuantLib.Time.Quarterly()),
                                        QuantLib.Time.Unadjusted(), QuantLib.Time.Unadjusted(), QuantLib.Time.DateGenerationBackwards(), false,
                                        QuantLib.Time.USNYSECalendar())
  fixing_days = 2
  in_arrears = true
  gearings = ones(length(fb_schedule.dates) - 1)
  spreads = fill(0.001, length(fb_schedule.dates) - 1)
  libor_3m = usd_libor_index(QuantLib.Time.TenorPeriod(Dates.Month(3)), yts)
  floating_bond = FloatingRateBond(settlement_days, face_amount, fb_schedule, libor_3m, fb_dc, conv, fixing_days, fb_issue_date, bond_engine, in_arrears, 100.0, gearings, spreads, cap_vol=cap_vol)

  return floating_bond
end

function generate_fixedrate_bond(yts::YieldTermStructure)
  settlement_days = 3
  face_amount = 100.0

  fx_schedule = QuantLib.Time.Schedule(Date(2007, 5, 15), Date(2017, 5, 15), QuantLib.Time.TenorPeriod(QuantLib.Time.Semiannual()),
                                        QuantLib.Time.Unadjusted(), QuantLib.Time.Unadjusted(), QuantLib.Time.DateGenerationBackwards(), false,
                                        QuantLib.Time.USGovernmentBondCalendar())

  pe = DiscountingBondEngine(yts)

  fixedrate_bond = FixedRateBond(settlement_days, face_amount, fx_schedule, 0.045, QuantLib.Time.ISMAActualActual(), QuantLib.Time.ModifiedFollowing(),
                                100.0, Date(2007, 5, 15), fx_schedule.cal, pe)

  return fixedrate_bond
end

function generate_discounting_ts(sett::Date)
  settlement_date = sett
  freq = QuantLib.Time.Semiannual()
  tenor = QuantLib.Time.TenorPeriod(freq)
  conv = QuantLib.Time.Unadjusted()
  conv_depo = QuantLib.Time.ModifiedFollowing()
  rule = QuantLib.Time.DateGenerationBackwards()
  calendar = QuantLib.Time.USGovernmentBondCalendar()
  dc_depo = QuantLib.Time.Actual365()
  dc = QuantLib.Time.ISDAActualActual()
  dc_bond = QuantLib.Time.ISMAActualActual()
  fixing_days = 3

  # build depos
  depo_rates = [0.0096, 0.0145, 0.0194]
  depo_tens = [Dates.Month(3), Dates.Month(6), Dates.Month(12)]

  # build bonds
  issue_dates = [Date(2005, 3, 15), Date(2005, 6, 15), Date(2006, 6, 30), Date(2002, 11, 15), Date(1987, 5, 15)]
  mat_dates = [Date(2010, 8, 31), Date(2011, 8, 31), Date(2013, 8, 31), Date(2018, 8, 15), Date(2038, 5, 15)]

  coupon_rates = [0.02375, 0.04625, 0.03125, 0.04000, 0.04500]
  market_quotes = [100.390625, 106.21875, 100.59375, 101.6875, 102.140625]

  insts = Vector{BootstrapHelper}(undef, length(depo_rates) + length(issue_dates))
  for i = 1:length(depo_rates)
    depo_quote = Quote(depo_rates[i])
    depo_tenor = QuantLib.Time.TenorPeriod(depo_tens[i])
    depo = DepositRateHelper(depo_quote, depo_tenor, fixing_days, calendar, conv_depo, true, dc_depo)
    insts[i] = depo
  end

  pricing_engine = DiscountingBondEngine()

  for i =1:length(coupon_rates)
    term_date = mat_dates[i]
    # rate = bond_rates[i] / 100
    rate = coupon_rates[i]
    issue_date = issue_dates[i]
    market_quote = market_quotes[i]
    sched = QuantLib.Time.Schedule(issue_date, term_date, tenor, conv, conv, rule, true)
    bond = FixedRateBondHelper(Quote(market_quote), FixedRateBond(3, 100.0, sched, rate, dc_bond, conv, 100.0, issue_date, calendar, pricing_engine))
    insts[i + length(depo_rates)] = bond
  end

  interp = QuantLib.Math.LogLinear()
  trait = Discount()
  bootstrap = IterativeBootstrap()

  yts = PiecewiseYieldCurve(settlement_date, insts, dc, interp, trait, 0.00000000001, bootstrap)

  # solver = QuantLib.Math.BrentSolver()
  # solver2 = QuantLib.Math.FiniteDifferenceNewtonSafe()
  calculate!(yts)

  return yts
end

function main()
  issue_date, bonds, dc, calendar = setup()

  println("Today's date: $issue_date")
  println("Calculating fit for 15 bonds....")

  interp = QuantLib.Math.LogLinear()
  trait = Discount()
  bootstrap = IterativeBootstrap()

  yts = PiecewiseYieldCurve(issue_date, bonds, dc, interp, trait, 0.00000000001, bootstrap)

  # solver = QuantLib.Math.BrentSolver()
  # solver2 = QuantLib.Math.FiniteDifferenceNewtonSafe()
  calculate!(yts)

  esf = ExponentialSplinesFitting(true, length(bonds))
  esf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, esf, 1e-10, 5000, 1.0)
  # @code_warntype initialize!(esf_fitted_curve)
  initialize!(esf_fitted_curve)
  calculate!(esf_fitted_curve)

  println("(a) exponential splines")
  println("reference date : ", esf_fitted_curve.referenceDate)
  println("number of iterations : ", esf_fitted_curve.fittingMethod.commons.numberOfIterations)

  spf = SimplePolynomialFitting(true, 3, length(bonds))
  spf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, spf, 1e-10, 5000, 1.0)
  initialize!(spf_fitted_curve)
  calculate!(spf_fitted_curve)

  println("(b) simple polynomial")
  println("reference date : ", spf_fitted_curve.referenceDate)
  println("number of iterations : ", spf_fitted_curve.fittingMethod.commons.numberOfIterations)

  nsf = NelsonSiegelFitting(length(bonds))
  nsf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, nsf, 1e-10, 5000, 1.0)
  initialize!(nsf_fitted_curve)
  calculate!(nsf_fitted_curve)

  println("(c) Nelson-Siegel")
  println("reference date : ", nsf_fitted_curve.referenceDate)
  println("number of iterations : ", nsf_fitted_curve.fittingMethod.commons.numberOfIterations)

  knots = [-30.0, -20.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0]
  cbsf = CubicBSplinesFitting(true, knots, length(bonds))
  cbsf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, cbsf, 1e-10, 5000, 1.0)
  initialize!(cbsf_fitted_curve)
  calculate!(cbsf_fitted_curve)

  println("(d) cubic B-splines")
  println("reference date : ", cbsf_fitted_curve.referenceDate)
  println("number of iterations : ", cbsf_fitted_curve.fittingMethod.commons.numberOfIterations)

  sf = SvenssonFitting(length(bonds))
  sf_fitted_curve = FittedBondDiscountCurve(0, issue_date, calendar, bonds, dc, sf, 1e-10, 5000, 1.0)
  initialize!(sf_fitted_curve)
  calculate!(sf_fitted_curve)

  println("(e) Svensson")
  println("reference date : ", sf_fitted_curve.referenceDate)
  println("number of iterations : ", sf_fitted_curve.fittingMethod.commons.numberOfIterations)

  println("Output par rates for each curve.  In this case, par rates should equal coupons for these par bonds")

  println(" tenor | coupon | bstrap |  (a)  |  (b)  |  (c)  |  (d)  |  (e)  ")

  for bh in bonds
    bond = bh.bond
    date_vec = Vector{Date}(undef, length(bond.cashflows.coupons) + 1)
    date_vec[1] = issue_date
    for (i, cf) in enumerate(bond.cashflows.coupons)
      date_vec[i+1] = date(cf)
    end

    if bh.bond.cashflows.redemption != nothing
      push!(date_vec, date(bh.bond.cashflows.redemption))
    end
    tenor = QuantLib.Time.year_fraction(dc, issue_date, date(bond.cashflows.coupons[end]))
    println(@sprintf(" %.2f  | %.3f | %.3f | %.3f | %.3f | %.3f | %.3f | %.3f ",
            tenor,
            100.0 * bond.cashflows.coupons[end].rate.rate,
            100.0 * par_rate(yts, date_vec, dc),
            100.0 * par_rate(esf_fitted_curve, date_vec, dc),
            100.0 * par_rate(spf_fitted_curve, date_vec, dc),
            100.0 * par_rate(nsf_fitted_curve, date_vec, dc),
            100.0 * par_rate(cbsf_fitted_curve, date_vec, dc),
            100.0 * par_rate(sf_fitted_curve, date_vec, dc)
            )
      )
  end
end

# function main2()
#   settlement_date = Date(2008, 9, 18)
#   calendar = QuantLib.Time.USGovernmentBondCalendar()
#   set_eval_date!(settings, settlement_date - Base.Dates.Day(3))
#
#   yts = generate_discounting_ts(settlement_date)
#   pe = DiscountingBondEngine(yts)
#
#   # build zero coupon bond
#   zcb = ZeroCouponBond(3, calendar, 100.0, Date(2013, 8, 15), QuantLib.Time.Following(), 116.92, Date(2003, 8, 15), pe)
#
#   # println(npv(zcb, pricing_engine, yts))
#   # println(clean_price(zcb))
#   # println(dirty_price(zcb))
#   return npv(zcb, zcb.pricingEngine), clean_price(zcb), dirty_price(zcb)
# end

function main3()
  settlement_date = Date(2008, 9, 18)
  set_eval_date!(settings, settlement_date - Dates.Day(3))
  cal = QuantLib.Time.TargetCalendar()
  dc = QuantLib.Time.ISDAActualActual()

  # Build deposits
  depo_quotes = [0.043375, 0.031875, 0.0320375, 0.03385, 0.0338125, 0.0335125]
  depo_tenors = [Dates.Week(1), Dates.Month(1), Dates.Month(3), Dates.Month(6), Dates.Month(9), Dates.Year(1)]
  deposit_dc = QuantLib.Time.Actual360()
  fixing_days = 3
  biz_conv = QuantLib.Time.ModifiedFollowing()

  depos = build_depos(depo_quotes, depo_tenors, deposit_dc, biz_conv, cal, fixing_days)

  # build swaps
  fixedLegFreq = QuantLib.Time.Annual()
  fixedLegConv = QuantLib.Time.Unadjusted()
  fixedLegDC = QuantLib.Time.EuroThirty360()
  floatingLegIndex = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(6)))
  forwardStart = 1

  swap_quotes = [0.0295, 0.0323, 0.0359, 0.0412, 0.0433]
  swap_tenors = [Dates.Year(2), Dates.Year(3), Dates.Year(5), Dates.Year(10), Dates.Year(15)]

  swaps = build_swaps(swap_quotes, swap_tenors, fixedLegDC, fixedLegConv, cal, fixedLegFreq, floatingLegIndex, forwardStart)

  insts = Vector{BootstrapHelper}(undef, length(swap_quotes) + length(depo_quotes))

  insts[1:length(depo_quotes)] = depos
  insts[length(depo_quotes) + 1: end] = swaps

  interp = QuantLib.Math.LogLinear()
  trait = Discount()
  bootstrap = IterativeBootstrap()

  yts = PiecewiseYieldCurve(settlement_date, insts, dc, interp, trait, 1e-15, bootstrap)

  # solver = QuantLib.Math.BrentSolver()
  # solver2 = QuantLib.Math.FiniteDifferenceNewtonSafe()
  calculate!(yts)

  disc_yts = generate_discounting_ts(settlement_date)

  cap_vol = ConstantOptionVolatility(3, cal, QuantLib.Time.ModifiedFollowing(), 0.0, QuantLib.Time.Actual365())
  # Floating Rate Bond
  fb = generate_floatingrate_bond(yts, disc_yts, cap_vol)

  # fb.iborIndex.ts = yts
  # fb.pricingEngine.yts = disc_yts

  # Zero Coupon Bond
  zcb_pe = DiscountingBondEngine(disc_yts)
  zcb_cal = QuantLib.Time.USGovernmentBondCalendar()

  # build zero coupon bond
  zcb = ZeroCouponBond(3, zcb_cal, 100.0, Date(2013, 8, 15), QuantLib.Time.Following(), 116.92, Date(2003, 8, 15), zcb_pe)

  # fixed rate bond
  fxb = generate_fixedrate_bond(disc_yts)
  # fxb.pricingEngine.yts = disc_yts

  println("Today's date: ", settings.evaluation_date)
  println("Settlement date: ", settlement_date)
  println("")
  println("              ZC      Fixed    Floating  ")
  println("-----------------------------------------")
  println(@sprintf("  NPV       %.2f   %.2f   %.2f", npv(zcb), npv(fxb), npv(fb)))
  println(@sprintf(" Clean      %.2f   %.2f   %.2f", clean_price(zcb), clean_price(fxb), clean_price(fb)))
  println(@sprintf(" Dirty      %.2f   %.2f   %.2f", dirty_price(zcb), dirty_price(fxb), dirty_price(fb)))
  println(@sprintf("accrued      %.2f     %.2f     %.2f", accrued_amount(zcb, settlement_date), accrued_amount(fxb, settlement_date), accrued_amount(fb, settlement_date)))
  println(@sprintf(" Next C      N/A      %.2f%%    %.2f%%", next_coupon_rate(fxb.cashflows, settlement_date) * 100.0, next_coupon_rate(fb.cashflows, settlement_date) * 100.0))
  println(@sprintf(" Yield      %.2f%%    %.2f%%    %.2f%%",
          QuantLib.yield(zcb, clean_price(zcb), QuantLib.Time.Actual360(), CompoundedCompounding(), QuantLib.Time.Annual(), settlement_date) * 100.0,
          QuantLib.yield(fxb, clean_price(fxb), QuantLib.Time.Actual360(), CompoundedCompounding(), QuantLib.Time.Annual(), settlement_date) * 100.0,
          QuantLib.yield(fb, clean_price(fb), QuantLib.Time.Actual360(), CompoundedCompounding(), QuantLib.Time.Annual(), settlement_date) * 100.0))

  # return npv(fb, fb.pricingEngine), clean_price(fb), dirty_price(fb), QuantLib.yield(fb, clean_price(fb), QuantLib.Time.Actual360(), CompoundedCompounding(), QuantLib.Time.Annual(), settlement_date), next_coupon_rate(fb.cashflows, settlement_date)
end


text:
module Sch

using Distributed
import MemPool: DRef

import ..Dagger
import ..Dagger: Context, Processor, Thunk, ThunkFuture, ThunkFailedException, Chunk, OSProc, order, free!, dependents, noffspring, istask, inputs, affinity, tochunk, @dbg, @logmsg, timespan_start, timespan_end, unrelease, procs, move, capacity, chunktype, choose_processor, default_enabled, get_processors, execute!, rmprocs!, addprocs!, thunk_processor

const OneToMany = Dict{Thunk, Set{Thunk}}

include("util.jl")
include("fault-handler.jl")
include("dynamic.jl")

"""
    ComputeState

The internal state-holding struct of the scheduler.

Fields:
- uid::UInt64 - Unique identifier for this scheduler instance
- dependents::OneToMany - The result of calling `dependents` on the DAG
- finished::Set{Thunk} - The set of completed `Thunk`s
- waiting::OneToMany - Map from downstream `Thunk` to upstream `Thunk`s that still need to execute
- waiting_data::OneToMany - Map from upstream `Thunk` to all downstream `Thunk`s, accumulating over time
- ready::Vector{Thunk} - The list of `Thunk`s that are ready to execute
- cache::Dict{Thunk, Any} - Maps from a finished `Thunk` to it's cached result, often a DRef
- running::Set{Thunk} - The set of currently-running `Thunk`s
- thunk_dict::Dict{Int, Any} - Maps from thunk IDs to a `Thunk`
- node_order::Any - Function that returns the order of a thunk
- worker_procs::Dict{Int,OSProc} - Maps from worker ID to root processor
- worker_pressure::Dict{Int,Dict{Type,Float64}} - Cache of worker pressure
- worker_capacity::Dict{Int,Dict{Type,Float64}} - Maps from worker ID to capacity
- worker_chans::Dict{Int, Tuple{RemoteChannel,RemoteChannel}} - Communication channels between the scheduler and each worker
- halt::Ref{Bool} - Flag indicating, when set, that the scheduler should halt immediately
- lock::ReentrantLock() - Lock around operations which modify the state
- futures::Dict{Thunk, Vector{ThunkFuture}} - Futures registered for waiting on the result of a thunk.
- errored::Set{Thunk} - Thunks that threw an error
- chan::Channel - Channel for receiving completed thunks
"""
struct ComputeState
    uid::UInt64
    dependents::OneToMany
    finished::Set{Thunk}
    waiting::OneToMany
    waiting_data::OneToMany
    ready::Vector{Thunk}
    cache::Dict{Thunk, Any}
    running::Set{Thunk}
    thunk_dict::Dict{Int, Any}
    node_order::Any
    worker_procs::Dict{Int,OSProc}
    worker_pressure::Dict{Int,Dict{Type,Float64}}
    worker_capacity::Dict{Int,Dict{Type,Float64}}
    worker_chans::Dict{Int, Tuple{RemoteChannel,RemoteChannel}}
    halt::Ref{Bool}
    lock::ReentrantLock
    futures::Dict{Thunk, Vector{ThunkFuture}}
    errored::Set{Thunk}
    chan::Channel{Any}
end

"""
    SchedulerOptions

Stores DAG-global options to be passed to the Dagger.Sch scheduler.

# Arguments
- `single::Int=0`: Force all work onto worker with specified id. `0` disables
this option.
- `proclist=nothing`: Force scheduler to use one or more processors that are
instances/subtypes of a contained type. Alternatively, a function can be
supplied, and the function will be called with a processor as the sole
argument and should return a `Bool` result to indicate whether or not to use
the given processor. `nothing` enables all default processors.
- `allow_errors::Bool=true`: Allow thunks to error without affecting
non-dependent thunks.
- `checkpoint=nothing`: If not `nothing`, Saves the final result of
the current scheduler invocation to persistent storage, for later retrieval by
`restore`.
- `restore=nothing`: If not `nothing`, returns the (cached) final result of
the current scheduler invocation, were it to execute. If this returns a
result, all thunks will be skipped. If this throws an error, restoring will be
skipped, the error will be displayed, and the scheduler will execute as usual.
"""
Base.@kwdef struct SchedulerOptions
    single::Int = 0
    proclist = nothing
    allow_errors::Bool = false
    checkpoint = nothing
    restore = nothing
end

"""
    ThunkOptions

Stores Thunk-local options to be passed to the Dagger.Sch scheduler.

# Arguments
- `single::Int=0`: Force thunk onto worker with specified id. `0` disables this
option.
- `proclist=nothing`: Force thunk to use one or more processors that are
instances/subtypes of a contained type. Alternatively, a function can be
supplied, and the function will be called with a processor as the sole
argument and should return a `Bool` result to indicate whether or not to use
the given processor. `nothing` enables all default processors.
- `procutil::Dict{Type,Any}=Dict{Type,Any}()`: Indicates the maximum expected
processor utilization for this thunk. Each keypair maps a processor type to
the utilization, where the value can be a real (approx. the number of processors
of this type utilized), or `MaxUtilization()` (utilizes all processors of this
type). By default, the scheduler assumes that this thunk only uses one
processor.
- `allow_errors::Bool=true`: Allow this thunk to error without affecting
non-dependent thunks.
- `checkpoint=nothing`: If not `nothing`, saves the final result of the thunk
to persistent storage, for later retrieval by `restore`.
- `restore=nothing`: If not `nothing`, returns the (cached) result of this
thunk, were it to execute. If this returns a result, this thunk will be
skipped. If this throws an error, restoring will be skipped, the error will be
displayed, and the scheduler will execute this thunk as usual.
"""
Base.@kwdef struct ThunkOptions
    single::Int = 0
    proclist = nothing
    procutil::Dict{Type,Any} = Dict{Type,Any}()
    allow_errors::Bool = false
    checkpoint = nothing
    restore = nothing
end

# Eager scheduling
include("eager.jl")

"Combine `SchedulerOptions` and `ThunkOptions` into a new `ThunkOptions`."
function merge(sopts::SchedulerOptions, topts::ThunkOptions)
    single = topts.single != 0 ? topts.single : sopts.single
    allow_errors = sopts.allow_errors || topts.allow_errors
    proclist = topts.proclist !== nothing ? topts.proclist : sopts.proclist
    ThunkOptions(single, proclist, topts.procutil, allow_errors, topts.checkpoint, topts.restore)
end

function isrestricted(task::Thunk, proc::OSProc)
    if (task.options !== nothing) && (task.options.single != 0) &&
       (task.options.single != proc.pid)
        return true
    end
    return false
end
merge(sopts::SchedulerOptions, ::Nothing) =
    ThunkOptions(sopts.single, sopts.proclist, Dict{Type,Any}())

function cleanup(ctx)
end

function init_proc(state, p)
    # Initialize pressure and capacity
    proc = OSProc(p.pid)
    lock(state.lock) do
        state.worker_procs[p.pid] = proc
        state.worker_pressure[p.pid] = Dict{Type,Float64}()
        state.worker_capacity[p.pid] = Dict{Type,Float64}()
        for T in unique(typeof.(get_processors(proc)))
            state.worker_pressure[p.pid][T] = 0
            state.worker_capacity[p.pid][T] = capacity(proc, T)
        end
        state.worker_pressure[p.pid][OSProc] = 0
        state.worker_capacity[p.pid][OSProc] = 0
    end
    cap = remotecall_fetch(capacity, p.pid)
    @async lock(state.lock) do
        state.worker_capacity[p.pid] = cap
    end

    # Setup worker-to-scheduler channels
    inp_chan = RemoteChannel(p.pid)
    out_chan = RemoteChannel(p.pid)
    lock(state.lock) do
        state.worker_chans[p.pid] = (inp_chan, out_chan)
    end
end

"Process-local count of actively-executing Dagger tasks per processor type."
const ACTIVE_TASKS = Dict{UInt64,Dict{Type,Ref{Float64}}}()
const ACTIVE_TASKS_LOCK = ReentrantLock()

"Process-local condition variable indicating task completion."
const TASK_SYNC = Condition()

"Indicates that a thunk uses all processors of a given type."
struct MaxUtilization end

function compute_dag(ctx, d::Thunk; options=SchedulerOptions())
    if options === nothing
        options = SchedulerOptions()
    end
    ctx.options = options
    if options.restore !== nothing
        try
            return options.restore()
        catch err
            @error "Scheduler restore failed" exception=(err,catch_backtrace())
        end
    end
    master = OSProc(myid())
    @dbg timespan_start(ctx, :scheduler_init, 0, master)

    chan = Channel{Any}(32)
    deps = dependents(d)
    ord = order(d, noffspring(deps))

    node_order = x -> -get(ord, x, 0)
    state = start_state(deps, node_order, chan)

    # setup thunk_dict mappings
    for node in keys(deps)
        state.thunk_dict[node.id] = node
        for dep in deps[node]
            state.thunk_dict[dep.id] = dep
        end
    end

    # Initialize procs, pressure, and capacity
    @sync for p in procs_to_use(ctx)
        @async init_proc(state, p)
    end

    # setup dynamic listeners
    dynamic_listener!(ctx, state)

    @dbg timespan_end(ctx, :scheduler_init, 0, master)

    # start off some tasks
    # Note: procs_state may be different things for different contexts. Don't touch it out here!
    procs_state = assign_new_procs!(ctx, state, procs_to_use(ctx))

    safepoint(state)

    # Loop while we still have thunks to execute
    while !isempty(state.ready) || !isempty(state.running)
        if !isempty(state.ready)
            # Nothing running, so schedule up to N thunks, 1 per N workers
            schedule!(ctx, state, procs_state)
        end

        # This is a bit redundant as the @async task below does basically the
        # same job Without it though, testing of process modification becomes
        # non-deterministic (due to sleep in CI environment) which is why it is
        # still here.
        procs_state = assign_new_procs!(ctx, state, procs_state)

        check_integrity(ctx)

        # Check periodically for new workers in a parallel task so that we
        # don't accidentally end up having to wait for `take!(chan)` on some
        # large task before new workers are put to work. Locking is used to
        # stop this task as soon as something pops out from the channel to
        # minimize risk that the task schedules thunks simultaneously as the
        # main task (after future refactoring).
        newtasks_lock = ReentrantLock()
        @async while !isempty(state.ready) || !isempty(state.running)
            sleep(1)
            islocked(newtasks_lock) && return
            procs_state = lock(newtasks_lock) do
                assign_new_procs!(ctx, state, procs_state)
            end
        end

        isempty(state.running) && continue
        pid, proc, thunk_id, (res, metadata) = take!(chan) # get result of completed thunk
        gproc = OSProc(pid)
        lock(newtasks_lock) # This waits for any assign_new_procs! above to complete and then shuts down the task
        safepoint(state)
        thunk_failed = false
        if res isa Exception
            if unwrap_nested_exception(res) isa Union{ProcessExitedException, Base.IOError}
                @warn "Worker $(pid) died on thunk $thunk_id, rescheduling work"

                # Remove dead worker from procs list
                remove_dead_proc!(ctx, state, gproc)

                lock(state.lock) do
                    handle_fault(ctx, state, state.thunk_dict[thunk_id], gproc)
                end
                continue
            else
                if ctx.options.allow_errors || state.thunk_dict[thunk_id].options.allow_errors
                    thunk_failed = true
                else
                    throw(res)
                end
            end
        end
        state.worker_pressure[pid][typeof(proc)] = metadata.pressure
        node = state.thunk_dict[thunk_id]
        state.cache[node] = res
        if node.options !== nothing && node.options.checkpoint !== nothing
            try
                node.options.checkpoint(node, res)
            catch err
                @error "Thunk checkpoint failed" exception=(err,catch_backtrace())
            end
        end

        # FIXME: Move log start and lock to before error check
        @dbg timespan_start(ctx, :finish, thunk_id, master)
        lock(state.lock) do
            finish_task!(state, node, thunk_failed)
        end
        @dbg timespan_end(ctx, :finish, thunk_id, master)

        safepoint(state)
    end
    value = state.cache[d] # TODO: move(OSProc(), state.cache[d])
    if d in state.errored
        throw(value)
    end
    if options.checkpoint !== nothing
        try
            options.checkpoint(value)
        catch err
            @error "Scheduler checkpoint failed" exception=(err,catch_backtrace())
        end
    end
    value
end

function procs_to_use(ctx, options=ctx.options)
    return if options.single !== 0
        @assert options.single in vcat(1, workers()) "Sch option `single` must specify an active worker ID."
        OSProc[OSProc(options.single)]
    else
        procs(ctx)
    end
end

check_integrity(ctx) = @assert !isempty(procs_to_use(ctx)) "No suitable workers available in context."

function schedule!(ctx, state, procs=procs_to_use(ctx))
    lock(state.lock) do
        safepoint(state)
        @assert length(procs) > 0
        proc_keys = map(x->x.pid, procs)
        proc_set = Set{Any}()
        for p in proc_keys
            for proc in get_processors(OSProc(p))
                push!(proc_set, p=>proc)
            end
        end
        failed_scheduling = Thunk[]
        while !isempty(state.ready)
            task = pop!(state.ready)
            opts = merge(ctx.options, task.options)
            proclist = opts.proclist
            proc_set_useable = if proclist === nothing
                filter(x->default_enabled(x[2]), proc_set)
            elseif proclist isa Function
                filter(x->proclist(x[2]), proc_set)
            else
                filter(x->typeof(x[2]) in proclist, proc_set)
            end
            if opts.single != 0
                proc_set_useable = filter(x->x[1]==opts.single, proc_set_useable)
            end
            @assert !isempty(proc_set_useable) "No processors available, try making proclist more liberal"
            procutil = opts.procutil
            gproc = nothing
            proc = nothing
            extra_util = nothing
            cap = nothing
            # FIXME: Sort by lowest utilization
            for (gp,p) in proc_set_useable
                T = typeof(p)
                extra_util = get(procutil, T, 1)
                real_util = state.worker_pressure[gp][T]
                cap = capacity(OSProc(gp), T)
                if ((extra_util isa MaxUtilization) && (real_util > 0)) ||
                   ((extra_util isa Real) && (extra_util + real_util > cap))
                    continue
                else
                    gproc = OSProc(gp)
                    proc = p
                    break
                end
            end
            if proc !== nothing
                extra_util = extra_util isa MaxUtilization ? cap : extra_util
                fire_task!(ctx, task, (gproc, proc), state; util=extra_util)
                continue
            else
                push!(failed_scheduling, task)
            end
        end
        append!(state.ready, failed_scheduling)
    end
end

# Main responsibility of this function is to check if new procs have been pushed to the context
function assign_new_procs!(ctx, state, assignedprocs)
    ps = procs_to_use(ctx)
    # Must track individual procs to handle the case when procs are removed
    diffps = setdiff(ps, assignedprocs)
    if !isempty(diffps)
        for p in diffps
            init_proc(state, p)
        end
        schedule!(ctx, state, diffps)
    end
    return ps
end

shall_remove_proc(ctx, proc) = proc ∉ procs_to_use(ctx)

function remove_dead_proc!(ctx, state, proc, options=ctx.options)
    @assert options.single !== proc.pid "Single worker failed, cannot continue."
    rmprocs!(ctx, [proc])
    delete!(state.worker_procs, proc.pid)
    delete!(state.worker_pressure, proc.pid)
    delete!(state.worker_capacity, proc.pid)
    delete!(state.worker_chans, proc.pid)
end

function pop_with_affinity!(ctx, tasks, proc)
    # TODO: use the size
    parent_affinity_procs = Vector(undef, length(tasks))
    # parent_affinity_sizes = Vector(undef, length(tasks))
    for i=length(tasks):-1:1
        t = tasks[i]
        aff = affinity(t)
        aff_procs = first.(aff)
        if proc in aff_procs
            if !isrestricted(t,proc)
                deleteat!(tasks, i)
                return t
            end
        end
        parent_affinity_procs[i] = aff_procs
    end
    for i=length(tasks):-1:1
        # use up tasks without affinities
        # let the procs with the respective affinities pick up
        # other tasks
        aff_procs = parent_affinity_procs[i]
        if isempty(aff_procs)
            t = tasks[i]
            if !isrestricted(t,proc)
                deleteat!(tasks, i)
                return t
            end
        end
        if all(!(p in aff_procs) for p in procs(ctx))
            # no proc is ever going to ask for it
            t = tasks[i]
            if !isrestricted(t,proc)
                deleteat!(tasks, i)
                return t
            end
        end
    end
    return nothing
end

function fire_task!(ctx, thunk, (gproc, proc), state; util=1.0)
    push!(state.running, thunk)
    if thunk.cache && thunk.cache_ref !== nothing
        # the result might be already cached
        data = unrelease(thunk.cache_ref) # ask worker to keep the data around
                                          # till this compute cycle frees it
        if data !== nothing
            # cache hit
            state.cache[thunk] = data
            thunk_failed = thunk in state.errored
            finish_task!(state, thunk, thunk_failed; free=false)
            return
        else
            # cache miss
            thunk.cache_ref = nothing
        end
    end
    if thunk.options !== nothing && thunk.options.restore !== nothing
        try
            result = thunk.options.restore(thunk)
            state.cache[thunk] = result
            finish_task!(state, thunk, false; free=false)
            return
        catch err
            @error "Thunk restore failed" exception=(err,catch_backtrace())
        end
    end

    ids = map(enumerate(thunk.inputs)) do (idx,x)
        istask(x) ? x.id : -idx
    end

    data = map(thunk.inputs) do x
        istask(x) ? state.cache[x] : x
    end
    toptions = thunk.options !== nothing ? thunk.options : ThunkOptions()
    options = merge(ctx.options, toptions)
    @assert (options.single == 0) || (gproc.pid == options.single)
    sch_handle = SchedulerHandle(ThunkID(thunk.id), state.worker_chans[gproc.pid]...)
    state.worker_pressure[gproc.pid][typeof(proc)] += util
    async_apply((gproc, proc), thunk.id, thunk.f, data, state.chan,
                thunk.get_result, thunk.persist, thunk.cache, thunk.meta,
                options, ids, ctx.log_sink, sch_handle, state.uid)
end

function finish_task!(state, node, thunk_failed; free=true)
    pop!(state.running, node)
    if !thunk_failed
        push!(state.finished, node)
    else
        set_failed!(state, node)
    end
    if istask(node) && node.cache
        node.cache_ref = state.cache[node]
    end
    if !thunk_failed
        for dep in sort!(collect(state.dependents[node]), by=state.node_order)
            set = state.waiting[dep]
            node in set && pop!(set, node)
            if isempty(set)
                pop!(state.waiting, dep)
                push!(state.ready, dep)
            end
            # todo: free data
        end
        if haskey(state.futures, node)
            # Notify any listening thunks
            for future in state.futures[node]
                if istask(node) && haskey(state.cache, node)
                    put!(future, state.cache[node])
                else
                    put!(future, nothing)
                end
            end
            delete!(state.futures, node)
        end
    end
    # Internal clean-up
    for inp in filter(istask, inputs(node))
        if inp in keys(state.waiting_data)
            s = state.waiting_data[inp]
            if node in s
                pop!(s, node)
            end
            if free && isempty(s)
                if haskey(state.cache, inp)
                    _node = state.cache[inp]
                    free!(_node, force=false, cache=(istask(inp) && inp.cache))
                    pop!(state.cache, inp)
                end
            end
        end
    end
end

function start_state(deps::Dict, node_order, chan)
    state = ComputeState(rand(UInt64),
                         deps,
                         Set{Thunk}(),
                         OneToMany(),
                         OneToMany(),
                         Vector{Thunk}(undef, 0),
                         Dict{Thunk, Any}(),
                         Set{Thunk}(),
                         Dict{Int, Thunk}(),
                         node_order,
                         Dict{Int,OSProc}(),
                         Dict{Int,Dict{Type,Float64}}(),
                         Dict{Int,Dict{Type,Float64}}(),
                         Dict{Int, Tuple{RemoteChannel,RemoteChannel}}(),
                         Ref{Bool}(false),
                         ReentrantLock(),
                         Dict{Thunk, Vector{ThunkFuture}}(),
                         Set{Thunk}(),
                         chan)

    nodes = sort(collect(keys(deps)), by=node_order)
    # N.B. Using merge! here instead would modify deps
    for (key,val) in deps
        state.waiting_data[key] = copy(val)
    end
    for k in nodes
        if istask(k)
            waiting = Set{Thunk}(Iterators.filter(istask, inputs(k)))
            if isempty(waiting)
                push!(state.ready, k)
            else
                state.waiting[k] = waiting
            end
        end
    end
    state
end

@noinline function do_task(to_proc, thunk_id, f, data, send_result, persist, cache, meta, options, ids, log_sink, sch_handle, uid)
    ctx = Context(Processor[]; log_sink=log_sink)

    from_proc = OSProc()
    # TODO: Time choose_processor
    Tdata = map(x->x isa Chunk ? chunktype(x) : x, data)

    # Fetch inputs
    fetched = if meta
        data
    else
        fetch.(map(Iterators.zip(data,ids)) do (x, id)
            @async begin
                @dbg timespan_start(ctx, :move, (thunk_id, id), (f, id))
                x = move(to_proc, x)
                @dbg timespan_end(ctx, :move, (thunk_id, id), (f, id))
                return x
            end
        end)
    end

    # Check if we'll go over capacity from running this thunk
    extra_util = get(options.procutil, typeof(to_proc), 1)
    real_util = lock(ACTIVE_TASKS_LOCK) do
        AT = get!(()->Dict{Type,Ref{Float64}}(), ACTIVE_TASKS, uid)
        get!(()->Ref{Float64}(0.0), AT, typeof(to_proc))
    end
    cap = Float64(capacity(OSProc(), typeof(to_proc)))
    while true
        lock(ACTIVE_TASKS_LOCK)
        if ((extra_util isa MaxUtilization) && (real_util[] > 0)) ||
           ((extra_util isa Real) && (extra_util + real_util[] > cap))
            # Fully subscribed, wait and re-check
            @debug "($(myid())) $f Waiting for free $(typeof(to_proc)): $extra_util | $(real_util[])/$cap"
            unlock(ACTIVE_TASKS_LOCK)
            wait(TASK_SYNC)
        else
            # Under-subscribed, calculate extra utilization and execute thunk
            @debug "($(myid())) Using available $to_proc: $extra_util | $(real_util[])/$cap"
            extra_util = if extra_util isa MaxUtilization
                count(c->typeof(c)===typeof(to_proc), from_proc.children)
            else
                extra_util
            end
            real_util[] += Float64(extra_util)
            unlock(ACTIVE_TASKS_LOCK)
            break
        end
    end

    @dbg timespan_start(ctx, :compute, thunk_id, f)
    res = nothing
    result_meta = try
        # Set TLS variables
        task_local_storage(:processor, to_proc)
        task_local_storage(:sch_handle, sch_handle)

        # Execute
        res = execute!(to_proc, f, fetched...)

        # Construct result
        send_result || meta ? res : tochunk(res, to_proc; persist=persist, cache=persist ? true : cache)
    catch ex
        bt = catch_backtrace()
        RemoteException(myid(), CapturedException(ex, bt))
    end
    @dbg timespan_end(ctx, :compute, thunk_id, (f, to_proc, typeof(res), sizeof(res)))
    lock(ACTIVE_TASKS_LOCK) do
        real_util[] -= Float64(extra_util)
    end
    @debug "($(myid())) Releasing $(typeof(to_proc)): $extra_util | $(real_util[])/$cap"
    notify(TASK_SYNC)
    metadata = (pressure=real_util[],)
    (result_meta, metadata)
end

@noinline function async_apply((gp,p), thunk_id, f, data, chan, send_res, persist, cache, meta, options, ids, log_sink, sch_handle, uid)
    @async begin
        try
            put!(chan, (gp.pid, p, thunk_id, remotecall_fetch(do_task, gp.pid, p, thunk_id, f, data,
                                                          send_res, persist, cache, meta, options, ids,
                                                          log_sink, sch_handle, uid)))
        catch ex
            bt = catch_backtrace()
            put!(chan, (gp.pid, p, thunk_id, (CapturedException(ex, bt), nothing)))
        end
        nothing
    end
end

end # module Sch


text:
include("planet.jl")
include("subgrid_scale.jl")
include("edmf.jl")
include("microphysics.jl")
include("atmos_subgrid_scale.jl")
include("surface_fluxes.jl")
include("override_defaults.jl")


text:
using MosaicViews
using Test
using ImageCore, ColorVectorSpace
using OffsetArrays

@testset "MosaicView" begin
    @test_throws ArgumentError MosaicView(rand(2,2,2,2,2))

    @testset "1D input" begin
        A = [1,2,3]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (3, 1) # 1D vector is lifted to 2D matrix
        @test axes(mv) == (Base.OneTo(3), Base.OneTo(1))
        @test @inferred(getindex(mv, 1, 1)) === 1
        @test @inferred(getindex(mv, 2, 1)) === 2
    end

    @testset "2D input" begin
        A = [1 2;3 4]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (2, 2)
        @test axes(mv) == (Base.OneTo(2), Base.OneTo(2))
        @test @inferred(getindex(mv, 1, 1)) === 1
        @test @inferred(getindex(mv, 2, 1)) === 3
    end

    @testset "3D input" begin
        A = zeros(Int,2,2,2)
        A[:,:,1] = [1 2; 3 4]
        A[:,:,2] = [5 6; 7 8]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 2)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test_throws BoundsError mv[0,1]
        @test_throws BoundsError mv[1,0]
        @test_throws BoundsError mv[1,3]
        @test_throws BoundsError mv[5,1]
        @test all(mv .== vcat(A[:,:,1],A[:,:,2]))
        # singleton dimension doesn't change anything
        @test mv == MosaicView(reshape(A,2,2,2,1))
    end

    @testset "4D input" begin
        A = zeros(Int,2,2,1,2)
        A[:,:,1,1] = [1 2; 3 4]
        A[:,:,1,2] = [5 6; 7 8]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (2, 4)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test_throws BoundsError mv[0,1]
        @test_throws BoundsError mv[1,0]
        @test_throws BoundsError mv[3,1]
        @test_throws BoundsError mv[1,5]
        @test all(mv .== hcat(A[:,:,1,1],A[:,:,1,2]))
        A = zeros(Int,2,2,2,3)
        A[:,:,1,1] = [1 2; 3 4]
        A[:,:,1,2] = [5 6; 7 8]
        A[:,:,1,3] = [9 10; 11 12]
        A[:,:,2,1] = [13 14; 15 16]
        A[:,:,2,2] = [17 18; 19 20]
        A[:,:,2,3] = [21 22; 23 24]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 6)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test all(mv .== vcat(hcat(A[:,:,1,1],A[:,:,1,2],A[:,:,1,3]), hcat(A[:,:,2,1],A[:,:,2,2],A[:,:,2,3])))
    end
end

@testset "mosaicview" begin
    @testset "1D input" begin
        A1 = [1, 2, 3]
        A2 = [4, 5, 6]
        mv = mosaicview(A1, A2)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A1)
        @test size(mv) == (6, 1)
        @test @inferred(getindex(mv, 1, 1)) == 1

        @test mosaicview(A1, A2; nrow=1) == [1 4; 2 5; 3 6]
        @test mosaicview(A1, A2; nrow=2) == reshape([1, 2, 3, 4, 5, 6], (6, 1))
        @test mosaicview(A1, A2) == mosaicview([A1, A2]) == mosaicview((A1, A2))
    end

    @testset "2D input" begin
        A = [i*ones(Int, 2, 3) for i in 1:4]
        Ao = [i*ones(Int, 0i:1, 0:2) for i in 1:4]

        for B in (A, tuple(A...), Ao)
            @test_throws ArgumentError mosaicview(B, nrow=0)
            @test_throws ArgumentError mosaicview(B, ncol=0)
            @test_throws ArgumentError mosaicview(B, nrow=1, ncol=1)

            mv = mosaicview(B)
            @test mosaicview(B...) == mv
            @test typeof(mv) <: MosaicView
            @test eltype(mv) == eltype(eltype(B))
            @test size(mv) == (8, 3)
            @test @inferred(getindex(mv,3,1)) === 2
            @test collect(mv) == [
                1  1  1
                1  1  1
                2  2  2
                2  2  2
                3  3  3
                3  3  3
                4  4  4
                4  4  4
            ]

            mv = mosaicview(B, nrow=2)
            @test typeof(mv) <: MosaicView
            @test eltype(mv) == eltype(eltype(B))
            @test size(mv) == (4, 6)
            @test collect(mv) == [
             1  1  1  3  3  3
             1  1  1  3  3  3
             2  2  2  4  4  4
             2  2  2  4  4  4
            ]
        end

        @test mosaicview(A...) == mosaicview(A)
        @test mosaicview(A..., nrow=2) == mosaicview(A, nrow=2)
        @test mosaicview(A..., nrow=2, rowmajor=true) == mosaicview(A, nrow=2, rowmajor=true)

        A1 = reshape([1 2 3], (1, 3))
        A2 = reshape([4;5;6], (3, 1))
        @test mosaicview([A1, A2]; center=false) == [
         1 2 3;
         0 0 0;
         0 0 0;
         4 0 0;
         5 0 0;
         6 0 0
        ]
        @test mosaicview([A1, A2]; center=true) == [
         0 0 0;
         1 2 3;
         0 0 0;
         0 4 0;
         0 5 0;
         0 6 0
        ]
        @test mosaicview([A1, A2]) == mosaicview([A1, A2]; center=true)

        # same size but different axes
        A1 = fill(1, 1:2, 1:2)
        A2 = fill(2, 2:3, 2:3)
        @test collect(mosaicview(A1, A2; center=true)) == [
            1 1 0;
            1 1 0;
            0 0 0;
            2 2 0;
            2 2 0;
            0 0 0;
        ]
        @test collect(mosaicview(A1, A2; center=false)) == [
            1 1 0;
            1 1 0;
            0 0 0;
            0 0 0;
            0 2 2;
            0 2 2;
        ]

        @test mosaicview(A1) == A1 # a trivial case
    end

    @testset "3D input" begin
        A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
        B = reshape(A, 2, 3, :)
        @test_throws ArgumentError mosaicview(B, nrow=0)
        @test_throws ArgumentError mosaicview(B, ncol=0)
        @test_throws ArgumentError mosaicview(B, nrow=1, ncol=1)
        @test mosaicview(B, center=2) == mosaicview(B) # no op

        mv = mosaicview(B)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(B)
        @test size(mv) == (8, 3)
        @test @inferred(getindex(mv,3,1)) === 2
        @test mv == [
            1  1  1
            1  1  1
            2  2  2
            2  2  2
            3  3  3
            3  3  3
            5  5  5
            5  5  5
        ]
        mv = mosaicview(B, nrow=2)
        @test mv == MosaicView(A)
        @test typeof(mv) != typeof(MosaicView(A))
        @test parent(parent(mv)).data == B
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(B)
        @test size(mv) == (4, 6)

        @test mosaicview(B, B) == mosaicview(cat(B, B; dims=4))
        @test mosaicview(B, B, nrow=2) == mosaicview(cat(B, B; dims=4), nrow=2)
        @test mosaicview(B, B, nrow=2, rowmajor=true) == mosaicview(cat(B, B; dims=4), nrow=2, rowmajor=true)
    end

    @testset "4D input" begin
        A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
        @test_throws ArgumentError mosaicview(A, nrow=0)
        @test_throws ArgumentError mosaicview(A, ncol=0)
        @test_throws ArgumentError mosaicview(A, nrow=1, ncol=1)
        @test mosaicview(A, center=2) == mosaicview(A) # no op

        mv = mosaicview(A)
        @test mv == MosaicView(A)
        @test typeof(mv) != typeof(MosaicView(A))
        @test parent(parent(mv)).data == reshape(A, 2, 3, :)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 6)
        mv = mosaicview(A, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 7)
        @test mv == mosaicview(A, nrow=2, npad=1)
        @test mv == mosaicview(A, ncol=2, npad=1)
        @test @inferred(getindex(mv,3,1)) === 0
        @test @inferred(getindex(mv,2,5)) === 3
        @test mv == [
            1  1  1  0  3  3  3
            1  1  1  0  3  3  3
            0  0  0  0  0  0  0
            2  2  2  0  5  5  5
            2  2  2  0  5  5  5
        ]
        mv = mosaicview(A, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test mv == [
            1  1  1  0  3  3  3  0  0  0  0
            1  1  1  0  3  3  3  0  0  0  0
            0  0  0  0  0  0  0  0  0  0  0
            2  2  2  0  5  5  5  0  0  0  0
            2  2  2  0  5  5  5  0  0  0  0
        ]
        mv = mosaicview(A, rowmajor=true, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test @inferred(getindex(mv,3,1)) === 0
        @test @inferred(getindex(mv,2,5)) === 2
        @test mv == [
            1  1  1  0  2  2  2  0  3  3  3
            1  1  1  0  2  2  2  0  3  3  3
            0  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0
        ]
        mv = mosaicview(A, rowmajor=true, ncol=3, npad=2)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (6, 13)
        @test mv == [
            1  1  1  0  0  2  2  2  0  0  3  3  3
            1  1  1  0  0  2  2  2  0  0  3  3  3
            0  0  0  0  0  0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0  0  0
        ]
        mv = mosaicview(A, fillvalue=-1.0, rowmajor=true, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test mv == [
             1   1   1  -1   2   2   2  -1   3   3   3
             1   1   1  -1   2   2   2  -1   3   3   3
            -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
             5   5   5  -1  -1  -1  -1  -1  -1  -1  -1
             5   5   5  -1  -1  -1  -1  -1  -1  -1  -1
        ]

        @test mosaicview(A, A) == mosaicview(cat(A, A; dims=5))
        @test mosaicview(A, A, nrow=2) == mosaicview(cat(A, A; dims=4), nrow=2)
        @test mosaicview(A, A, nrow=2, rowmajor=true) == mosaicview(cat(A, A; dims=4), nrow=2, rowmajor=true)
    end

    @testset "Colorant Array" begin
        A = rand(RGB{Float32}, 2, 3, 2, 2)
        mv = mosaicview(A)
        @test eltype(mv) == eltype(A)
        @test mv == @inferred(MosaicView(A))
        mv = mosaicview(A, rowmajor=true, ncol=3)
        @test eltype(mv) == eltype(A)
        @test @inferred(getindex(mv, 3, 4)) == RGB(0,0,0)
        mv = mosaicview(A, fillvalue=colorant"white", rowmajor=true, ncol=3)
        @test eltype(mv) == eltype(A)
        @test @inferred(getindex(mv, 3, 4)) == RGB(1,1,1)

        # this should work regardless they're of different size and color
        @test_nowarn mosaicview(rand(RGB{Float32}, 4, 4),
                                rand(Gray{N0f8}, 5, 5))
    end

    # all arrays should have the same dimension
    @test_throws ArgumentError mosaicview(ones(2), ones(1, 2))
    @test_throws ArgumentError mosaicview((ones(2), ones(1, 2)))
    @test_throws ArgumentError mosaicview([ones(2), ones(1, 2)])
end

@testset "deprecations" begin
    @info "deprecations are expected"
    A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
    mv_old = mosaicview(A, -1.0, rowmajor=true, ncol=3, npad=1)
    mv_new = mosaicview(A, fillvalue=-1.0, rowmajor=true, ncol=3, npad=1)
    @test mv_old == mv_new
end