text:
include("../src/io/util.jl")
include("../src/io/args.jl")
include("../src/io/netcdf.jl")
include("../src/rasterizer/rasterization.jl")

using Test
using Base.Threads

ATOL=0.05

function norm(vec)
    mag = sqrt(sum(abs2, vec))
    return Tuple(el / mag for el ∈ vec)
end

function make_context(simplices, points;
                      sampling_rate=(100., 100., 100.),
                      mem_limit=4*1024^3,
                      dyn_var_mapping=Dict("W"=>"d"),
                      stat_var_mapping=Dict("b"=>"b"),
                      z_range=Rasterization.z_all,
                      tanioka=false,
                      water_height=0.,
                      t_begin=0,
                      t_end=0,
                      times=[0.],
                      out_filename="temp.nc",
                      custom_domain=((-Inf, Inf), (-Inf, Inf)),
                      max_threads=typemax(Int))

    Rasterization.RasterizationContext(simplices, points, sampling_rate, mem_limit,
                                       dyn_var_mapping, stat_var_mapping, z_range, tanioka,
                                       water_height, t_begin, t_end, times, out_filename, custom_domain, max_threads=max_threads)
end

@testset "Rasterization" begin
    if nthreads() == 1
        error("Binning tests only work in a multi-threaded environment. Please run with at least 3 threads.")
    end

    simplices2D = [[1, 2, 3   ] [2, 3, 4   ]]
    simplices3D = [[1, 2, 3, 4] [2, 3, 4, 5]]

    A = [1., 1., 1.]
    B = [2., 1., 2.]
    C = [1., 2., 1.]
    D = [2., 1., 1.]
    E = [2., 2., 2.]
    points2D = [A B C   E]
    points3D = [A B C D E]

    dimensions = Dict(2=>(simplices=simplices2D, points=points2D),
                      3=>(simplices=simplices3D, points=points3D))

    #============================================#
    # calculate_hnfs!
    #============================================#

    # Each case is [Dimension, Simplex Index, Simplex HNFs]
    hnf_test_cases = [
        (dim=2,
         simp_idx=1,
         hnfs=[Rasterization.HesseNormalForm(norm((-1., -1.,  0.)), -2.12),
               Rasterization.HesseNormalForm(     ( 1.,  0.,  0.) ,  1.00),
               Rasterization.HesseNormalForm(     ( 0.,  1.,  0.) ,  1.00)]),

        (dim=3,
         simp_idx=1,
         hnfs=[Rasterization.HesseNormalForm(norm((-1., -1.,  0.)), -2.12),
               Rasterization.HesseNormalForm(     ( 0.,  0.,  1.) ,  1.00),
               Rasterization.HesseNormalForm(     ( 0.,  1.,  0.) ,  1.00),
               Rasterization.HesseNormalForm(norm(( 1.,  0., -1.)),  0.00)])
    ]

    @testset "calculate_hnfs! ($(test_case.dim)D, simplex $(test_case.simp_idx))" for test_case ∈ hnf_test_cases
        dimension = dimensions[test_case.dim]
        ctx = make_context(dimension.simplices, dimension.points)
        l_simp_points = [@view ctx.points[:, p] for p ∈ ctx.simplices[:, test_case.simp_idx]]
        l_simp_hnfs = Array{Rasterization.HesseNormalForm, 1}(undef, ctx.n_simplex_points)

        Rasterization.calculate_hnfs!(ctx, l_simp_points, l_simp_hnfs)

        for p ∈ 1:ctx.n_simplex_points
            @test test_case.hnfs[p].d ≈ l_simp_hnfs[p].d atol=ATOL

            for dim ∈ 1:3
                @test test_case.hnfs[p].n0[dim] ≈ l_simp_hnfs[p].n0[dim] atol=ATOL
            end
        end
    end

    #============================================#
    # is_in_simplex
    #============================================#

    is_in_test_cases = [
        (dim=2,
        simp_idx=1,
        points=[A, B, C, (A .+ B .+ C) / 3., A .- [1., .5, .2], B .+ [2., 0., 0.], C .+ [.1, .2, 0.]],
        is_in=[true, false, false, true, false, false, false]) # B and C intersect an which has a negative normal. Thus, they are not in the triangle themselves

        (dim=3,
        simp_idx=1,
        points=[A, B, C, D, (A .+ B .+ C .+ D) / 4., (A .+ C .+ D) / 3., A .- [1., .5, .2], B .- [.5, .5, 0.], D .+ [0., .2, 0.]],
        is_in=[true, false, false, false, true, true, false, false, false])
    ]

    @testset "is_in_simplex ($(test_case.dim)D, simplex $(test_case.simp_idx))" for test_case ∈ is_in_test_cases
        dimension = dimensions[test_case.dim]
        ctx = make_context(dimension.simplices, dimension.points)
        l_simp_points = [@view ctx.points[:, p] for p ∈ ctx.simplices[:, test_case.simp_idx]]
        l_simp_hnfs = Array{Rasterization.HesseNormalForm, 1}(undef, ctx.n_simplex_points)

        Rasterization.calculate_hnfs!(ctx, l_simp_points, l_simp_hnfs)

        for (point, result) ∈ zip(test_case.points, test_case.is_in)
            @test result == Rasterization.is_in_simplex(ctx, l_simp_hnfs, point)
        end
    end

    #============================================#
    # is_bathy
    #============================================#

    is_bathy_test_cases = [
        (bathy_heights=[1998., 2001., 2001.], water_height=2000., is_bathy=true),
        (bathy_heights=[-20., -19., -17.], water_height=0., is_bathy=true),
        (bathy_heights=[-20., -20., -20.], water_height=0., is_bathy=true),
        (bathy_heights=[0., 0., 0.], water_height=20., is_bathy=true),
        (bathy_heights=[0., 0., 0.], water_height=0., is_bathy=false),
        (bathy_heights=[2000., 2000., 2000.], water_height=2000., is_bathy=false),
        (bathy_heights=[2001., 2001., 2001.], water_height=2000., is_bathy=true)
    ]

    @testset "is_bathy (wh=$(test_case.water_height), b=$(test_case.bathy_heights))" for test_case ∈ is_bathy_test_cases
        bathys = test_case.bathy_heights
        points = [[1., 1., bathys[1]] [4., 4., bathys[2]] [2., 2., bathys[3]]]
        @test test_case.is_bathy == Rasterization.is_bathy(points, test_case.water_height)
    end

    #============================================#
    # bin
    #============================================#

    nBins = 3
    cX = nBins*3
    cY = 9
    cZ = 1
    points = zeros((3, (cX+1)*(cY+1)*(cZ+1)))
    for i ∈ 0:size(points,2)-1
        points[:, i+1] = [i % (cX+1), i ÷ (cX+1) % (cY+1), i ÷ ((cX+1) * (cY+1)) % (cZ+1)]
    end

    tup2idx = (x, y, z) -> (x + (cX+1) * (y + (cY+1) * z) + 1)

    for z ∈ 0:cZ, y ∈ 0:cY, x ∈ 0:cX
        p = points[:,tup2idx(x,y,z)]
        if tuple(p...) != (x,y,z)
            println("$((x,y,z)) != $p")
        end
    end

    simplices2D = []
    bin_ids = []
    z = 0
    water_height = 1.
    for z ∈ 0:cZ, y ∈ 0:0, x ∈ 0:cX-1
        # For each (x,y,z), which we define as the left upper point of the triangle, add two triangles:
        # (x,y  ,z)  o---o (x+1,y  ,z)
        #            |\ 1|
        #            | \ |
        #            |2 \|
        # (x,y+1,z)  o---o (x+1,y+1,z)
        push!(simplices2D, [tup2idx(x, y, z), tup2idx(x+1, y, z), tup2idx(x+1, y+1, z)])
        push!(simplices2D, [tup2idx(x, y, z), tup2idx(x+1, y+1, z), tup2idx(x, y+1, z)])
        bin_id = z < water_height ? 1 + x ÷ (cX ÷ nBins) : 2 * nBins + 1
        push!(bin_ids, bin_id)
        push!(bin_ids, bin_id)
    end

    for z ∈ 0:cZ, y ∈ 0:0, xBin ∈ 1:nBins-1
        # Triangles all span one bucket border
        x=(cX ÷ nBins)*xBin-1
        push!(simplices2D, [tup2idx(x, y, z), tup2idx(x+2, y, z), tup2idx(x+2, y+1, z)])
        push!(simplices2D, [tup2idx(x, y, z), tup2idx(x+2, y+1, z), tup2idx(x, y+1, z)])
        bin_id = z < water_height ? nBins + xBin : 2 * nBins + 1
        push!(bin_ids, bin_id)
        push!(bin_ids, bin_id)
    end

    for z ∈ 0:cZ, y ∈ 0:0
        # Triangles all span multiple borders
        push!(simplices2D, [tup2idx(1, y, z), tup2idx(cX, y, z), tup2idx(cX, y+1, z)])
        push!(simplices2D, [tup2idx(1, y, z), tup2idx(cX, y+1, z), tup2idx(1, y+1, z)])
        bin_id = z < water_height ? 2 * nBins : 2 * nBins + 1
        push!(bin_ids, bin_id)
        push!(bin_ids, bin_id)
    end

    simplices2D = [simplices2D[i][dim] for dim ∈ 1:3, i ∈ 1:length(simplices2D)]

    @testset "bin" begin
        ctx = make_context(simplices2D, points, sampling_rate=(1., 1., 1.), water_height=water_height, z_range=Rasterization.z_floor, max_threads=3)
        l_bin_ids, l_bin_counts = Rasterization.bin(ctx)

        for i ∈ 1:size(simplices2D,2)
            @test bin_ids[i] == l_bin_ids[i]
        end
    end

    @testset "iterate" begin end
    Rasterization.iterate
    @testset "iterate_simps!" begin end
    Rasterization.iterate_simps!
    @testset "cross3!" begin end
    Rasterization.cross3!
    @testset "dot3" begin end
    Rasterization.dot3

    #============================================#
    # bary_coords_2d
    #============================================#

    triangle = [[-2., 1., 14.] [1., 3., 13.] [-1., -1., 12.]]

    bary_test_cases = [
        (point=triangle[:, 1],                           bary=(  1.,   0.,   0.)),
        (point=triangle[:, 2],                           bary=(  0.,   1.,   0.)),
        (point=triangle[:, 3],                           bary=(  0.,   0.,   1.)),
        (point=sum(triangle,           dims=2)[:,1] / 3, bary=(1/3 , 1/3 , 1/3 )),
        (point=sum(triangle[:, 1:2],   dims=2)[:,1] / 2, bary=(1/2 , 1/2 ,   0.)),
        (point=sum(triangle[:, 2:3],   dims=2)[:,1] / 2, bary=(  0., 1/2 , 1/2 )),
        (point=sum(triangle[:, 1:2:3], dims=2)[:,1] / 2, bary=(1/2 ,   0., 1/2 )),

        (point=triangle[:, 1] + 2 * (triangle[:, 2] - triangle[:, 1]), bary=(-1.,  2., 0.)),
        (point=triangle[:, 1] + 2 * (triangle[:, 3] - triangle[:, 1]), bary=(-1.,  0., 2.)),
        (point=triangle[:, 2] + 2 * (triangle[:, 3] - triangle[:, 2]), bary=( 0., -1., 2.)),
    ]

    @testset "bary_coords_2d" for test_case ∈ bary_test_cases
        calculated_bary = Rasterization.bary_coords_2d(test_case.point, triangle)
        for dim ∈ 1:3
            @test test_case.bary[dim] ≈ calculated_bary[dim] atol=1e-6
        end
    end
end

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
export forw2spot, spot2forw

## capital market -----------------------------------------------
"""
`getyield(τ, stock::Stock)`

calculates the yield of `stock` during year `τ`
"""
function getyield(τ, stock::Stock)
  if τ > 1
    return stock.x[τ] / stock.x[τ - 1] - 1
  elseif τ == 1
    return stock.x[τ] / stock.x_0 - 1
  else
    return stock.yield_0
  end
end

"""
`getyield(τ, rfr::RiskFreeRate)`

calculates the yield of `rfr` during year `τ`
"""
function getyield(τ, rfr::RiskFreeRate)
  if τ >= 1
    return rfr.x[τ]
  else
    return rfr.yield_0
  end
end

"""
`forw2spot(f::Vector{Float64})`

calculates the spot rate `s` rate from the forward rate `f`:

  `(1+f[1])(1+f[2])...(1+f[n]) = (1+s[n])^n`
"""
forw2spot(f::Vector{Float64}) =
  cumprod(1 .+ f) .^ (1 ./ collect(1:length(f))) .- 1

"""
`spot2forw(s::Vector{Float64})`

calculates the forward rate `f` rate from the spot rate `s`:

  `(1 + s[n-1])^(n-1) * (1+f[n]) = (1+s[n])^n`
"""
function spot2forw(s::Vector{Float64})
  f = zeros(Float64, length(s))
  for 𝑛 ∈ length(s):-1:2
    f[𝑛] = (1 + s[𝑛])^𝑛 / (1 + s[𝑛-1])^(𝑛 - 1) -1
  end
  f[1] = s[1]
  return f
end

## investments --------------------------------------------------
"""
`project!(τ::Int, mv_boy::Float64, invest::Invest)`

project `invest::Invest` one year for a given
initial market value.

**Changed**:  `invest`
"""
function project!(τ::Int, mv_boy::Float64, invest::Invest)
  invest.mv[τ] = (1 + getyield(τ, invest.proc)) * mv_boy
end

"""
`project!(τ::Int, mv_bop_total::Float64, ig::InvestGroup)`

project `ig::InvestGroup` one year for a given
initial market value.

**Changed**:  `ig`
"""
function project!(τ::Int, mv_bop_total::Float64, ig::InvestGroup)
  mv_bop = ig.alloc.total[τ] * mv_bop_total
  ig.mv[τ] = 0
  for (𝑖, 𝑖𝑔_𝑖𝑛𝑣𝑒𝑠𝑡) ∈ enumerate(ig.investments)
    project!(τ, ig.alloc.all[τ, 𝑖] * mv_bop, 𝑖𝑔_𝑖𝑛𝑣𝑒𝑠𝑡)
    ig.mv[τ] += 𝑖𝑔_𝑖𝑛𝑣𝑒𝑠𝑡.mv[τ]
  end
  ig.cost.total[τ] =
    ig.cost.abs[τ] * ig.cost.cum_infl_abs[τ] +
    mv_bop * ig.cost.rel[τ] * ig.cost.cum_infl_rel[τ]
end

"""
`alloc!(τ, cap_mkt::CapMkt, invs::InvPort)`

Dynamically re-allocate investments `invs::InvPort` at the
beginning of year `τ`

**Changed**:  `invs`
"""
function alloc!(τ, cap_mkt::CapMkt, invs::InvPort)
  if τ > 1
    invs.igs[:IGStock].alloc.total[τ] =
      0.5 * (1 - exp( -max(0, dynstateavg(τ, cap_mkt))))
    invs.igs[:IGCash].alloc.total[τ] =
      1 - invs.igs[:IGStock].alloc.total[τ]
    ## we leave the allocations within each group unchanged:
    for 𝑠𝑦𝑚𝑏 ∈ [:IGCash, :IGStock]
      for 𝑖 = 1:size(invs.igs[𝑠𝑦𝑚𝑏].alloc.all, 2)
        invs.igs[𝑠𝑦𝑚𝑏].alloc.all[τ, 𝑖] =
          invs.igs[𝑠𝑦𝑚𝑏].alloc.all[τ-1, 𝑖]
      end
    end
  end
end

"""
`project!(τ::Int, mv_boy::Float64, invs::InvPort)`

Project `invs::InvPort`  one year for a given
initial market value.

**Changed**:  `invs`
"""
function project!(τ::Int, mv_boy::Float64, invs::InvPort)
  invs.mv[τ] = 0.0
  invs.cost[τ] = 0.0
  invs.mv_boy[τ] = mv_boy
  for 𝑖𝑔 ∈ values(invs.igs)
    project!(τ, mv_boy, 𝑖𝑔)
    invs.mv[τ] += 𝑖𝑔.mv[τ]
    invs.cost[τ] += 𝑖𝑔.cost.total[τ]
  end
  invs.yield[τ] = invs.mv[τ] / mv_boy - 1
end

## insurance liabilities  ---------------------------------------
"""
`premium(ins_sum, rfr, prob, β, λ)`

Calculate the premium of a product
"""
function premium(ins_sum, rfr, prob, β, λ)
  lx_boy = [1; cumprod(prob[:px])[1:end-1]]
  v_eoy = 1 ./ cumprod(1 .+ rfr)
  v_boy = [1; v_eoy[1:end-1]]
  num =
    sum(lx_boy .* ins_sum .*
        (v_boy .* λ[:boy] .* λ[:cum_infl] ./ (1 .+ λ[:infl]) .+
           v_eoy .* (λ[:eoy] .* λ[:cum_infl] .+
                       prob[:px] .* β[:px] .+
                       prob[:qx] .* β[:qx])
         ))
  denom =
    sum(lx_boy .* β[:prem] .*
        (v_boy - v_eoy .* prob[:sx] .* β[:sx]))
  return num / denom
end

"""
`tpgrec(τ, tpg, rfr, prob, β, λ)`

One year backward recursion formula for technical provisions
of guaranteed benefits
"""
function tpgrec(τ, tpg, rfr, prob, β, λ)
  disc = 1/(1 + rfr[τ + 1])
  -β[τ + 1, :prem] +
    λ[τ + 1, :boy] *
    λ[τ + 1, :cum_infl] / (1 + λ[τ + 1, :infl]) +
    disc * (λ[τ + 1, :eoy] * λ[τ + 1, :cum_infl] +
              prob[τ + 1, :qx] * β[τ + 1, :qx] +
              prob[τ + 1, :sx] * β[τ + 1, :sx] +
              prob[τ + 1, :px] * (β[τ + 1, :px] + tpg))
end

"""
`tpg(τ, rfr, prob, β, λ)`

Technical provisions of guaranteed benefits
at the end of year `τ`
"""
function tpg(τ, rfr, prob, β, λ)
  dur = nrow(β)
  res = 0.0
  if τ >= dur
    return 0.0
  else
    for 𝑠 ∈ (dur-1):-1:τ
      res = tpgrec(𝑠, res, rfr, prob, β, λ)
    end
    return res
  end
end

"""
`tpg(τ, rfr, mp)`

Best estimate guaranteed technical provisions
of a model point `mp` at the end of year `τ`
"""
tpg(τ, rfr, mp) = tpg(τ, rfr, mp.prob, mp.β, mp.λ)

"""
`tpgfixed(τ, rfr, fixed_cost_gc::Vector)`

Technical provisions at the end of year `τ`
for future fixed (going concern) costs
"""
function tpgfixed(τ, rfr, fixed_cost_gc::Vector)
  dur = length(fixed_cost_gc)
  if dur < τ + 1
    return 0.0
  else
    disc = cumprod(1 ./ (1 .+ rfr[(τ + 1):dur]))
    return fixed_cost_gc[(τ + 1):dur] ⋅ disc
  end
end

## other liabilities --------------------------------------------
"""
`pv(τ::Int, cap_mkt::CapMkt, debt::Debt)`

Present value of `debt` at the end of year `τ`.
Any servicing of this debt during year `τ` has occured beforehand
"""
function pv(τ::Int, cap_mkt::CapMkt, debt::Debt)
  ## calculate pv at the end of the year after servicing debt
  if (debt.τ_init > τ) | (debt.τ_mat <= τ)
    return 0.0
  else
    p_v = debt.nominal
    for 𝑠 ∈ (debt.τ_mat - 1) : -1 : τ
      p_v = (debt.coupon + p_v) / (1 + cap_mkt.rfr.x[𝑠 + 1])
    end
    return p_v
  end
end

"""
`pv(τ::Int, cap_mkt::CapMkt, l_other::LiabOther)`

present value of other liabilities `l_other` at the end of year
`τ`. Any servicing of debt within this portfolio during year `τ`
has occured beforehand
"""
function pv(τ::Int, cap_mkt::CapMkt, l_other::LiabOther)
  p_v = 0.0
  for 𝑑𝑒𝑏𝑡 ∈ l_other.subord
    p_v += pv(τ, cap_mkt, 𝑑𝑒𝑏𝑡)
  end
  return p_v
end

"""
`paycoupon(τ::Int, debt::Debt)`

Coupon payment for `debt` at the end of year `τ`
"""
paycoupon(τ::Int, debt::Debt) =
  (debt.τ_init <= τ <= debt.τ_mat ? debt.coupon : 0.0)

"""
`paycoupon(τ::Int, l_other::LiabOther)`

Coupon payment for all debts within the portfolio of
other liabilities `l_other` at the end of year `τ`
"""
function paycoupon(τ::Int, l_other::LiabOther)
  pay = 0.0
  for 𝑑𝑒𝑏𝑡 ∈ l_other.subord
    pay += paycoupon(τ, 𝑑𝑒𝑏𝑡)
  end
  return pay
end

"""
`payprincipal(τ::Int, debt::Debt)`

Payment of the principal of `debt` at the end of year `τ`,
if the debt matures at this point in time
"""
payprincipal(τ::Int, debt::Debt) =
  (τ == debt.τ_mat ? debt.nominal : 0.0)

"""
`payprincipal(τ::Int, l_other::LiabOther)`

Payment of the total principal of all debts within the
portfolio of other liabilities `l_other`, which mature at
the end of year `τ`
"""
function payprincipal(τ::Int, l_other::LiabOther)
  pay = 0.0
  for 𝑑𝑒𝑏𝑡 ∈ l_other.subord
    pay += payprincipal(τ, 𝑑𝑒𝑏𝑡)
  end
  return pay
end

"""
`getloan(τ::Int, debt::Debt)`

Get the nominal of `debt` at time `τ`, if it has been
taken out at this point in time, otherwise return 0.0
"""
getloan(τ::Int, debt::Debt) =
  (τ == debt.τ_init ? debt.nominal : 0.0)

"""
`getloan(τ::Int, l_other::LiabOther)`

Get the total nominal of all debts within the
portfolio of other liabilities `l_other`,
which are taken out at the beginning of year `τ`
"""
function getloan(τ::Int, l_other::LiabOther)
  nominal = 0.0
  for 𝑑𝑒𝑏𝑡 ∈ l_other.subord
    nominal += getloan(τ, 𝑑𝑒𝑏𝑡)
  end
  return nominal
end

"""
`goingconcern(debts::Vector{Debt}, Δgc::Vector{Float64})`

Calculates a vector of debts from an existing vector of debts
according to the going concern assumption. The total initial debt
is the same, but the new debt vectors mature earlier so that the
total nominal decreases according to the going concern factors.
We do not input the factors `gc` directly but their year on year
differences `Δgc`.
"""
function goingconcern(debts::Vector{Debt}, Δgc::Vector{Float64})
  new_debt_vec = Array{Debt}(undef, 0)
  for 𝑑𝑒𝑏𝑡 ∈ debts
    if 𝑑𝑒𝑏𝑡.nominal > 0.0
      τ_init = max(1, 𝑑𝑒𝑏𝑡.τ_init)
      diff_nom = -Δgc * 𝑑𝑒𝑏𝑡.nominal
      for 𝜏 ∈ τ_init:𝑑𝑒𝑏𝑡.τ_mat
        t = 𝑑𝑒𝑏𝑡.t_init + 𝜏 - 𝑑𝑒𝑏𝑡.τ_init
        push!(new_debt_vec,
              Debt(𝑑𝑒𝑏𝑡.name,
                   𝑑𝑒𝑏𝑡.t_init,
                   t,
                   𝑑𝑒𝑏𝑡.τ_init,
                   𝜏,
                   diff_nom[𝜏],
                   𝑑𝑒𝑏𝑡.coupon * diff_nom[𝜏] / 𝑑𝑒𝑏𝑡.nominal))
      end
    end
  end
  return(new_debt_vec)
end

"""
`goingconcern!(l_other::LiabOther, Δgc::Vector{Float64})`

Transforms a portfolio of other liabilities `l_other::LiabOther`
according to the going concern assumption. We do not input the
factors `gc` directly but their year on year differences `Δgc`.

**Changed**: `l_other`
"""
function goingconcern!(l_other::LiabOther, Δgc::Vector{Float64})
  l_other.subord = goingconcern(l_other.subord, Δgc)
end

## dynamics -----------------------------------------------------
"""
`dynstate(τ, cap_mkt::CapMkt)`

indicator for the state of the economy at the end of year `τ`
"""
dynstate(τ, cap_mkt::CapMkt) =
  getyield(τ, cap_mkt.stock) / max(getyield(τ, cap_mkt.rfr), eps()) - 1

"""
`dynstateavg(τ, cap_mkt::CapMkt)`

Two year average of the indicator for the state of the economy
at the end of year `τ`.
"""
dynstateavg(τ, cap_mkt::CapMkt) =
  0.5 * (getyield(τ - 1, cap_mkt.stock) /
           max(getyield(τ-1, cap_mkt.rfr), eps())
         + getyield(τ, cap_mkt.stock) /
           max(getyield(τ, cap_mkt.rfr), eps())) - 1

"""
`bonusrate(yield_eoy, rfr_price, bonus_factor)`

Helper function for dynamic bonus rate declaration
"""
bonusrate(yield_eoy, rfr_price, bonus_factor) =
  max(bonus_factor * (yield_eoy - rfr_price), 0.0)

"""
`bonusrate(τ, yield_eoy, mp::ModelPoint, dyn)`

Dynamic bonus rate declaration for a model point at the end
of year `τ`
"""
bonusrate(τ, yield_eoy, mp::ModelPoint, dyn) =
  bonusrate(yield_eoy,
            ( τ == 0 ? mp.rfr_price_0 : mp.rfr_price[τ]),
            dyn.bonus_factor)

"""
`biquotient(τ, yield_eoy, cap_mkt, invs, mp, dyn)`

Indicator for bonus rate expectation at the end of year `τ`
"""
function biquotient(τ, yield_eoy, cap_mkt, invs, mp, dyn)
  if τ ≤ mp.dur
    ind_bonus =
      getyield(τ, cap_mkt.stock) /
      max(0.0,
          bonusrate(τ - 1, yield_eoy, mp, dyn) + mp.rfr_price[τ])
    ind_bonus_hypo =
      getyield(0, cap_mkt.stock) /
      max(eps(), mp.bonus_rate_hypo + mp.rfr_price_0)
    return ind_bonus / ind_bonus_hypo
  else
    return 0.0
  end
end

"""
`δsx(τ, cap_mkt, invs, mp, dyn)`

Dynamic lapse probability factor to adjust the initial estimate
"""
function δsx(τ, cap_mkt, invs, mp, dyn)
  yield_eoy =
    invs.igs[:IGStock].alloc.total[τ] * getyield(τ, cap_mkt.stock) +
    invs.igs[:IGCash].alloc.total[τ] * getyield(τ, cap_mkt.rfr)
  if τ - 1 > mp.t_start
    bi_quot =  biquotient(τ, yield_eoy, cap_mkt, invs, mp, dyn)
  else
    bi_quot = 1
  end
  state_quot = dynstate(τ, cap_mkt) / dynstate(0, cap_mkt)
  δ_SX = 1.0
  if state_quot < 0.5
    δ_SX += 0.15
  elseif  state_quot > 2.0
    δ_SX -= 0.15
  end
  δ_SX += 0.25 * min(4.0, max(0.0, bi_quot - 1.2))
  return δ_SX
end

"""
`freesurp(dyn, invest_pre, liab)`

Helper function for the free surplus calculation
"""
freesurp(dyn, invest_pre, liab) =
 max(0, invest_pre - (1 + dyn.quota_surp) * liab)

"""
`freesurp(τ, proj::Projection, dyn)`

Free surplus for the dynamic dividend declaration
"""
freesurp(τ, proj::Projection, dyn) =
  if τ == 1
      freesurp(dyn,
               proj.val_0[1, :invest],
               proj.val_0[1, :tpg] + proj.val_0[1, :l_other])
  else
      freesurp(dyn,
               proj.val[τ-1, :invest],
               proj.val[τ-1, :tpg] + proj.val[τ-1, :l_other])
  end

"""
`update!(τ, proj::Projection, dyn::Dynamic)`

Update dynamic parameters
"""
function update!(τ, proj::Projection, dyn::Dynamic)
  if τ == 1
    dyn.free_surp_boy[τ] =
      freesurp(dyn,
               proj.val_0[1, :invest],
               proj.val_0[1, :tpg] + proj.val_0[1, :l_other])
  else
    dyn.free_surp_boy[τ] =
      freesurp(dyn,
               proj.val[τ-1, :invest],
               proj.val[τ-1, :tpg] + proj.val[τ-1, :l_other])
  end
end

## cashflow projection  -----------------------------------------
"""
`val0!(cap_mkt::CapMkt, invs::InvPort, liabs::LiabIns,
  l_other::LiabOther, proj::Projection)`

Valuation at time `t_0`

**Changed:** `proj::Projection`
"""
function val0!(cap_mkt::CapMkt,
               invs::InvPort,
               liabs::LiabIns,
               l_other::LiabOther,
               proj::Projection)
  proj.val_0[1, :invest] = invs.mv_0
  for 𝑚𝑝 ∈ liabs.mps
    if 0 <= 𝑚𝑝.dur
      proj.val_0[1, :tpg] += tpg(0, cap_mkt.rfr.x, 𝑚𝑝)
    end
  end
  proj.val_0[1, :l_other] = pv(0, cap_mkt, l_other)
  proj.val_0[1, :surplus] =
    proj.val_0[1, :invest] -
    proj.val_0[1, :tpg] -
    proj.val_0[1, :l_other]
end

"""
`projectboy!(τ, proj::Projection, liabs::LiabIns)`

Project one year, update values at the beginning of the year `τ`

**Changed:** `proj::Projection`, `liabs::LiabIns  (liabs.mp)`
"""
function projectboy!(τ, proj::Projection, liabs::LiabIns)
  proj.cf[τ, :prem] = 0.0
  proj.cf[τ, :λ_boy] = 0.0
  for 𝑚𝑝 ∈ liabs.mps
    if τ <= 𝑚𝑝.dur
      𝑚𝑝.lx_boy[τ] = (τ == 1 ? 1 : 𝑚𝑝.lx_boy_next)
      proj.cf[τ, :prem] += 𝑚𝑝.lx_boy[τ] * 𝑚𝑝.β[τ, :prem]
      proj.cf[τ, :λ_boy] -=
        𝑚𝑝.lx_boy[τ] * 𝑚𝑝.λ[τ, :boy] *
        𝑚𝑝.λ[τ, :cum_infl] / (1 + 𝑚𝑝.λ[τ, :infl])
    end
  end
end

"""
`projecteoy!(τ, cap_mkt::CapMkt, invs::InvPort, liabs::LiabIns,
  dyn::Dynamic, proj::Projection)`

Project one year, update values at the end of the year `τ`

**Changed:** `proj::Projection`, `liabs::LiabIns  (liabs.mp)`
"""
function projecteoy!(τ,
                     cap_mkt::CapMkt,
                     invs::InvPort,
                     liabs::LiabIns,
                     dyn::Dynamic,
                     proj::Projection)
  tpg_price_positive = 0.0
  for 𝑚𝑝 ∈ liabs.mps
    if τ <= 𝑚𝑝.dur
      tpg_price_positive +=
        𝑚𝑝.lx_boy[τ] *
        max(0, (τ == 1 ? 𝑚𝑝.tpg_price_0 : 𝑚𝑝.tpg_price[τ-1]))
    end
  end
  for 𝑚𝑝 ∈ liabs.mps
    if τ <= 𝑚𝑝.dur
      prob = deepcopy(𝑚𝑝.prob)
      prob[:,:sx] *=
        δsx(τ, cap_mkt, invs, 𝑚𝑝, dyn)
      prob[:,:px] = 1 .- prob[:,:qx] - prob[:,:sx]
      𝑚𝑝.lx_boy_next = 𝑚𝑝.lx_boy[τ] * prob[τ, :px]
      for 𝑤𝑥 ∈ [:qx, :sx, :px]
        proj.cf[τ, 𝑤𝑥] -=
          𝑚𝑝.lx_boy[τ] * prob[τ, 𝑤𝑥] * 𝑚𝑝.β[τ, 𝑤𝑥]
      end
      proj.cf[τ, :λ_eoy] -=
        𝑚𝑝.lx_boy[τ] * 𝑚𝑝.λ[τ, :eoy] * 𝑚𝑝.λ[τ, :cum_infl]
      proj.val[τ, :tpg] +=
        𝑚𝑝.lx_boy[τ] * prob[τ, :px] *
        tpg(τ, cap_mkt.rfr.x, prob, 𝑚𝑝.β, 𝑚𝑝.λ)
    end
  end
  proj.cf[τ, :Δtpg] =
    -(proj.val[τ, :tpg] -
        (τ == 1 ? proj.val_0[1, :tpg] : proj.val[τ - 1, :tpg]))

end

"""
`project!(τ, cap_mkt::CapMkt, invs::InvPort, dyn::Dynamic,
  proj::Projection)`

Project one year, investment results from the year `τ`

**Changed:** `proj::Projection`, `invs::InvPort`
"""
function project!(τ,
                  cap_mkt::CapMkt,
                  invs::InvPort,
                  dyn::Dynamic,
                  proj::Projection)
  mv_boy =
    (τ == 1 ? proj.val_0[1, :invest] : proj.val[τ - 1, :invest])
  mv_boy += proj.cf[τ, :prem] + proj.cf[τ, :λ_boy]
  alloc!(τ, cap_mkt, invs)
  project!(τ, mv_boy, invs)
  proj.cf[τ, :invest] = invs.mv[τ] - mv_boy
end

"""
`bonus!(τ, invs::InvPort, liabs::LiabIns, dyn::Dynamic,
  proj, surp_pre_profit_tax_bonus)`

Bonus at the end of year `τ`

**Changed:** `proj::Projection`
"""
function bonus!(τ,
                invs::InvPort,
                liabs::LiabIns,
                dyn::Dynamic,
                proj,
                surp_pre_profit_tax_bonus)
  for 𝑚𝑝 ∈ liabs.mps
    if τ <= 𝑚𝑝.dur
      proj.cf[τ, :bonus] -=
        min(surp_pre_profit_tax_bonus,
            𝑚𝑝.lx_boy[τ] *
              bonusrate(τ, invs.yield[τ], 𝑚𝑝, dyn) *
              max(0, (τ == 1 ?
                        𝑚𝑝.tpg_price_0 :
                        𝑚𝑝.tpg_price[τ-1])))
    end
  end
end

"""
`investpredivid(τ, invs::InvPort, proj::Projection)`

Market value of assets before payment of dividends
"""
function investpredivid(τ, invs::InvPort, proj::Projection)
  invs.mv_boy[τ] +
    sum(convert(Array,
                proj.cf[τ,
                        [:invest, :qx, :sx, :px, :λ_eoy, :bonus,
                         :l_other, :tax, :gc]]))
end

"""
`project!(τ, cap_mkt::CapMkt, invs::InvPort, liabs::LiabIns,
  liab_other::LiabOther, dyn::Dynamic, proj::Projection)`

Project one year

**Changed:** `proj::Projection`, `invs::InvPort`, `dyn::Dynamic`
"""
function project!(τ,
                  cap_mkt::CapMkt,
                  invs::InvPort,
                  liabs::LiabIns,
                  liab_other::LiabOther,
                  dyn::Dynamic,
                  proj::Projection)
  projectboy!(τ, proj, liabs)
  proj.cf[τ, :new_debt] = -getloan(τ, liab_other)
  project!(τ, cap_mkt, invs, dyn, proj)
  update!(τ, proj, dyn)
  proj.cf[τ, :λ_eoy] = -invs.cost[τ]
  proj.cf[τ, :l_other] = -paycoupon(τ, liab_other)
  proj.cf[τ, :l_other] -= payprincipal(τ, liab_other)
  proj.val[τ, :l_other] = pv(τ, cap_mkt, liab_other)
  projecteoy!(τ, cap_mkt, invs, liabs, dyn, proj)

  proj.cf[τ,:tax] = 0.0
  proj.cf[τ,:bonus] = 0.0
  surp_pre_profit_tax_bonus =
    max(0,
        investpredivid(τ, invs, proj) -
          proj.val[τ, :tpg] -
          proj.val[τ, :l_other]  )
  bonus!(τ, invs, liabs, dyn, proj, surp_pre_profit_tax_bonus)
  proj.cf[τ, :profit] =
    sum(convert(Array, proj.cf[τ, [:prem, :invest,
                                   :qx, :sx, :px, :λ_boy, :λ_eoy,
                                   :Δtpg, :bonus, :l_other]]))
  tax = proj.tax_rate * proj.cf[τ, :profit] ## could be negative
  tax_credit_pre =
    (τ == 1 ? proj.tax_credit_0 : proj.tax_credit[τ - 1])
  proj.cf[τ, :tax] = -max(0, tax - tax_credit_pre)
  proj.tax_credit[τ] =  ## no new tax credit generated
  tax_credit_pre - tax - proj.cf[τ, :tax]

  proj.val[τ, :invest] = investpredivid(τ, invs, proj)
  proj.cf[τ, :divid] =
    -freesurp(dyn,
              proj.val[τ, :invest],
              proj.val[τ, :tpg] + proj.val[τ, :l_other])
  proj.val[τ, :invest] +=  proj.cf[τ, :divid]

  proj.val[τ, :surplus] =
    proj.val[τ, :invest] -
    proj.val[τ, :tpg] -
    proj.val[τ, :l_other]
end

"""
`pvprev(rfr, cf, pv) = (cf + pv) /  (1 + rfr)`

Recursive step in generic present value calculation
"""
pvprev(rfr, cf, pv) = (cf + pv) /  (1 + rfr)

"""
`pvvec(rfr::Vector{Float64}, cf)``

Generic present value calculation, where `rfr` denotes the risk
free rate for each year and `cf`. The length of the vector `rfr`
may not be smaller than the length of the cashflow.
Each payment occurs at the end of the corresponding year.

Output: A vector of the same length as `cf` contaning the present
value at the end of each year. The last component is zero.
"""
function pvvec(rfr::Vector{Float64}, cf)
  T = length(cf)
  val = zeros(Float64, T)
  for 𝑡 ∈ reverse(collect(1:(T-1))) # [T-1:-1:1] #r
    val[𝑡] = pvprev(rfr[𝑡 + 1], val[𝑡 + 1], cf[𝑡 + 1])
  end
  return val
end

"""
`valbonus!(rfr::Vector{Float64}, proj::Projection) `

Provisions for future bonus payments (at each time `τ`).
Needs to be called after the projection is completed

**Changed:** proj::Projection
"""
function valbonus!(rfr::Vector{Float64},
                   proj::Projection)
  proj.val[:bonus] = pvvec(rfr,  -proj.cf[:bonus])
  proj.val_0[:bonus] =
    pvprev(rfr[1], -proj.cf[1, :bonus], proj.val[1, :bonus])
end

"""
`valcostprov!(rfr::Vector{Float64}, invs::InvPort,
  proj::Projection)`

Provisions for future costs (at each time `τ`). Absolute costs
(including absolute investment costs from *all* investments) and
relative investment costs for provisions are considered.
It is assumed that provisions are backed by cash investments
Needs to be called after the projection is completed

**Changed:** proj::Projection
"""
function valcostprov!(rfr::Vector{Float64},
                      invs::InvPort,
                      proj::Projection)
  cash_cost = deepcopy(invs.igs[:IGCash].cost)
  proj.cf[1, :cost_prov] =
    proj.fixed_cost_gc[1] +
    sum(convert(Array,
                proj.val_0[1, [:tpg, :bonus, :l_other]])) *
    cash_cost.cum_infl_rel[1] * cash_cost.rel[1]
  for 𝑡 ∈ 2:proj.dur
    proj.cf[𝑡, :cost_prov] =
      proj.fixed_cost_gc[𝑡] +
      sum(convert(Array,
                  proj.val[𝑡 - 1, [:tpg, :bonus, :l_other]])) *
      cash_cost.cum_infl_rel[𝑡] * cash_cost.rel[𝑡]
  end
  proj.val[:cost_prov] =
    pvvec(rfr - cash_cost.rel,  proj.cf[:cost_prov])
  proj.val_0[:cost_prov] = pvprev(rfr[1] - cash_cost.rel[1],
                                  proj.cf[1, :cost_prov],
                                  proj.val[1, :cost_prov])
end

ext: jl


text:
# Note: Terminology overall taken from:
# https://en.wikipedia.org/wiki/Interval_(mathematics)#Terminology

"""
    Bound <: Any

Abstract type representing all possible endpoint classifications (e.g. open, closed,
unbounded).
"""
:Bound

"""
    Bounded <: Bound

Abstract type indicating that the endpoint of an interval is not unbounded (e.g. open or
closed).
"""
:Bounded

"""
    Closed <: Bounded <: Bound

Type indicating that the endpoint of an interval is closed (the endpoint value is *included*
in the interval).
"""
:Closed

"""
    Open <: Bounded <: Bound

Type indicating that the endpoint of an interval is open (the endpoint value is *not
included* in the interval).
"""
:Open

"""
    Unbounded <: Bound

Type indicating that the endpoint of an interval is unbounded (the endpoint value is
effectively infinite).
"""
:Unbounded

"""
    first(interval::AbstractInterval{T}) -> Union{T,Nothing}

The value of the lower endpoint. When the lower endpoint is unbounded `nothing` will be
returned.
"""
Base.first(::AbstractInterval)

"""
    last(interval::AbstractInterval{T}) -> Union{T,Nothing}

The value of the upper endpoint. When the upper endpoint is unbounded `nothing` will be
returned.
"""
Base.last(::AbstractInterval)

"""
    span(interval::AbstractInterval) -> Any

The delta between the upper and lower endpoints. For bounded intervals returns a
non-negative value while intervals with any unbounded endpoints will throw an
`ArgumentError`.

To avoid having to capture the exception use the pattern:
```julia
!isbounded(interval) ? span(interval) : infinity`
```
Where `infinity` is a variable representing the value you wish to use to represent an
unbounded, or infinite, span.
"""
span(::AbstractInterval)

"""
    isclosed(interval) -> Bool

Is a closed-interval: includes both of its endpoints.
"""
isclosed(::AbstractInterval)

"""
    isopen(interval) -> Bool

Is an open-interval: excludes both of its endpoints.
"""
Base.isopen(::AbstractInterval)

"""
    isunbounded(interval) -> Bool

Is an unbounded-interval: unbounded at both ends.
"""
isunbounded(::AbstractInterval)

"""
    isbounded(interval) -> Bool

Is a bounded-interval: either open, closed, left-closed/right-open, or
left-open/right-closed.

Note using `!isbounded` is commonly used to determine if any end of the interval is
unbounded.
"""
isbounded(::AbstractInterval)


text:
# Assorted utility functions.

"""
    diag(A::Array, dim1::Int, dim2::Int)

Collapse `dim1` and `dim2` of `A` into a single dimension containing only the
elements at equal indices.
"""
function diag(A::Array, dim1::Int, dim2::Int)
    s = [size(A)...]
    s[dim1] == s[dim2] || throw(DomainError(s, "Dimensions must have the same size."))

    deleteat!(s, dim2)

    result = similar(A, s...)
    for idx1 in CartesianIndices(A)
        idx1[dim1] == idx1[dim2] || continue
        idx2 = [idx1.I...]
        deleteat!(idx2, dim2)
        result[idx2...] = A[idx1]
    end

    result
end

"""
    @cat A B

Append `B` along the last dimension of `A`, modifying the binding of `A`.
"""
macro cat(A::Symbol, B)
    quote
        $(esc(A)) = cat($(esc(A)), $(esc(B)); dims=length(size($(esc(A)))))
    end
end

Maybe{T} = Union{T,Nothing} where {T}

function show_value(io::IO, x::Float64)
    if x == 0.0
        print(io, "    .          ")
    else
        @printf(io, "% 15.10f", x)
    end
    nothing
end

function show_vector(io::IO, xs::Vector{Float64})
    for i in 1:length(xs)
        show_value(io, xs[i])
    end
    println(io)
    nothing
end

function show_matrix(io::IO, xs::Matrix{Float64})
    for i in 1:size(xs, 2)
        show_vector(io, xs[:, i])
    end
    nothing
end


text:
const γ_exact = 7 // 5
const μ_exact = 1 // 100
@noinline ρ_g(t, x, y, z, ::Val{2}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) + 3
@noinline U_g(t, x, y, z, ::Val{2}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y)
@noinline V_g(t, x, y, z, ::Val{2}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y)
@noinline W_g(t, x, y, z, ::Val{2}) = 0
@noinline E_g(t, x, y, z, ::Val{2}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) + 100
@noinline Sρ_g(t, x, y, z, ::Val{2}) =
    pi * (
        2 * sin(2 * pi * x) - 2 * sin(2 * pi * y) - sin(pi * (2 * t - 2 * x)) +
        sin(pi * (2 * t + 2 * x)) +
        sin(pi * (2 * t - 2 * y)) - sin(pi * (2 * t + 2 * y)) +
        2 * sin(pi * (2 * x + 2 * y)) +
        sin(pi * (-2 * t + 2 * x + 2 * y)) +
        sin(pi * (2 * t + 2 * x + 2 * y)) +
        10 * cos(pi * (-t + x + y)) - 2 * cos(pi * (t - x + y)) +
        2 * cos(pi * (t + x - y)) +
        14 * cos(pi * (t + x + y))
    ) / 8
@noinline SU_g(t, x, y, z, ::Val{2}) =
    pi * (
        -600 * sin(pi * t) * sin(pi * x)^2 * cos(pi * t) * cos(pi * y)^2 -
        900 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        900 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        1800 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) +
        540 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 +
        1080 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 +
        7 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        pi * sin(pi * y) * cos(pi * t) * cos(pi * x) +
        120 * cos(pi * t) * cos(pi * x) * cos(pi * y)
    ) / 300
@noinline SV_g(t, x, y, z, ::Val{2}) =
    pi * (
        -600 * sin(pi * t) * sin(pi * x)^2 * cos(pi * t) * cos(pi * y)^2 -
        900 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        540 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        1080 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) +
        900 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 -
        120 * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        1800 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 +
        7 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        pi * sin(pi * y) * cos(pi * t) * cos(pi * x)
    ) / 300
@noinline SW_g(t, x, y, z, ::Val{2}) = 0
@noinline SE_g(t, x, y, z, ::Val{2}) =
    pi * (
        -300 * sin(pi * t) * sin(pi * x) * cos(pi * y) +
        480 * sin(pi * x)^4 * sin(pi * y) * cos(pi * t)^4 * cos(pi * y)^3 +
        1080 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        480 * sin(pi * x)^3 * cos(pi * t)^4 * cos(pi * x) * cos(pi * y)^4 -
        7 * pi * sin(pi * x)^2 * sin(pi * y)^2 * cos(pi * t)^2 -
        840 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) -
        1080 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 +
        14 * pi * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        4 *
        pi *
        sin(pi * x) *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y) - 42000 * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        840 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 -
        7 * pi * cos(pi * t)^2 * cos(pi * x)^2 * cos(pi * y)^2 +
        42000 * cos(pi * t) * cos(pi * x) * cos(pi * y)
    ) / 300
@noinline ρ_g(t, x, y, z, ::Val{3}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3
@noinline U_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y) *
    cos(pi * z)
@noinline V_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y) *
    cos(pi * z)
@noinline W_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    sin(pi * z) *
    cos(pi * t) *
    cos(pi * y)
@noinline E_g(t, x, y, z, ::Val{3}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 100
@noinline Sρ_g(t, x, y, z, ::Val{3}) =
    pi * (
        -sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) -
        2 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        3 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        2 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        3 * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) +
        3 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    )
@noinline SU_g(t, x, y, z, ::Val{3}) =
    pi * (
        -600 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        900 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) -
        900 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        600 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) +
        300 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 -
        1800 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        180 *
        sin(pi * x)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z) -
        900 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        540 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        900 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        360 *
        sin(pi * x) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 +
        1080 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        10 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) +
        pi * sin(pi * y) * cos(pi * t) * cos(pi * x) * cos(pi * z) -
        pi * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z) +
        120 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    ) / 300
@noinline SV_g(t, x, y, z, ::Val{3}) =
    pi * (
        -600 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        900 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) +
        180 *
        sin(pi * x)^3 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z) -
        540 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        600 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) +
        300 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 +
        360 *
        sin(pi * x)^2 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * y) -
        1080 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        900 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        900 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        900 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        120 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        pi * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        1800 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        10 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) +
        pi * sin(pi * y) * cos(pi * t) * cos(pi * x) * cos(pi * z)
    ) / 300
@noinline SW_g(t, x, y, z, ::Val{3}) =
    pi *
    (
        -600 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z) - 900 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        900 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        240 * sin(pi * x)^3 * sin(pi * z)^2 * cos(pi * t)^3 * cos(pi * y)^3 +
        840 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^2 -
        1800 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z) +
        900 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^2 +
        2160 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z) -
        pi * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        1800 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z) - 120 * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        10 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        pi * cos(pi * t) * cos(pi * x) * cos(pi * y)
    ) *
    sin(pi * z) / 300
@noinline SE_g(t, x, y, z, ::Val{3}) =
    pi * (
        45 *
        (1 - cos(2 * pi * x))^2 *
        (1 - cos(4 * pi * z)) *
        (cos(2 * pi * t) + 1)^2 *
        (cos(2 * pi * y) + 1)^2 / 128 -
        300 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) +
        240 *
        sin(pi * x)^4 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^4 *
        cos(pi * y)^3 *
        cos(pi * z)^2 +
        480 *
        sin(pi * x)^4 *
        sin(pi * y) *
        cos(pi * t)^4 *
        cos(pi * y)^3 *
        cos(pi * z)^4 +
        60 * sin(pi * x)^4 * sin(pi * z)^4 * cos(pi * t)^4 * cos(pi * y)^4 -
        120 * sin(pi * x)^4 * cos(pi * t)^4 * cos(pi * y)^4 * cos(pi * z)^4 +
        540 *
        sin(pi * x)^3 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z) +
        1080 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        240 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^4 *
        cos(pi * x) *
        cos(pi * y)^4 *
        cos(pi * z)^2 +
        180 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) -
        480 *
        sin(pi * x)^3 *
        cos(pi * t)^4 *
        cos(pi * x) *
        cos(pi * y)^4 *
        cos(pi * z)^4 -
        360 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 -
        3 * pi * sin(pi * x)^2 * sin(pi * y)^2 * sin(pi * z)^2 * cos(pi * t)^2 -
        7 * pi * sin(pi * x)^2 * sin(pi * y)^2 * cos(pi * t)^2 * cos(pi * z)^2 -
        7 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * y) -
        840 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        3 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        540 *
        sin(pi * x)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z) -
        420 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        4 * pi * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 -
        1080 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        16 * pi * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 +
        420 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 +
        4 *
        pi *
        sin(pi * x) *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y) *
        cos(pi * z)^2 -
        42000 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        7 *
        pi *
        sin(pi * x) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 +
        3 *
        pi *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        840 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        42000 * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) -
        3 * pi * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * x)^2 * cos(pi * y)^2 -
        7 * pi * cos(pi * t)^2 * cos(pi * x)^2 * cos(pi * y)^2 * cos(pi * z)^2 +
        42000 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    ) / 300


text:
using CuArrays

# Test GPU movement inside the call to `gradient`
@testset "GPU movement" begin
  r = rand(Float32, 3,3)
  @test gradient(x -> sum(cu(x)), r)[1] isa AbstractArray
end

@testset "basic bcasting" begin
  a = cu(collect(1:9))
  v(x, n) = x .^ n
  pow_grada = cu(Float32[7.0, 448.0, 5103.0, 28672.0, 109375.0, 326592.0, 823543.0, 1.835008e6, 3.720087e6])
  @test gradient(x -> v(x, 7) |> sum, a) == (pow_grada,)
  w(x) = broadcast(log, x)
  log_grada = cu(Float32[1.0, 0.5, 0.33333334, 0.25, 0.2, 0.16666667, 0.14285715, 0.125, 0.11111111])
  @test gradient(x -> w(x) |> sum, a) == (log_grada,)
end


text:
importall Base

abstract ImmutableArray{T,N} <: AbstractArray{T,N}
typealias ImmutableVector{T} ImmutableArray{T,1}
typealias ImmutableMatrix{T} ImmutableArray{T,2}

export unit, row, column

function generate_arrays(maxSz::Integer)

    # operations
    const unaryOps = (:-, :~, :conj, :abs,
                      :sin, :cos, :tan, :sinh, :cosh, :tanh,
                      :asin, :acos, :atan, :asinh, :acosh, :atanh,
                      :sec, :csc, :cot, :asec, :acsc, :acot,
                      :sech, :csch, :coth, :asech, :acsch, :acoth,
                      :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
                      :sind, :tand, :acosd, :acotd, :acscd, :asecd,
                      :asind, :atand, :radians2degrees, :degrees2radians,
                      :log, :log2, :log10, :log1p, :exponent, :exp,
                      :exp2, :expm1, :cbrt, :sqrt, :square, :erf,
                      :erfc, :erfcx, :erfi, :dawson, :ceil, :floor,
                      :trunc, :round, :significand, :lgamma, :hypot,
                      :gamma, :lfact, :frexp, :modf, :airy, :airyai,
                      :airyprime, :airyaiprime, :airybi, :airybiprime,
                      :besselj0, :besselj1, :bessely0, :bessely1,
                      :eta, :zeta, :digamma)

    # vec-vec and vec-scalar
    const binaryOps = (:.+, :.-,:.*, :./, :.\, :.^,
                       :.==, :.!=, :.<, :.<=, :.>, :.>=,
                       :min, :max,
                       :div, :fld, :rem, :mod, :mod1, :cmp,
                       :atan2, :besselj, :bessely, :hankelh1, :hankelh2,
                       :besseli, :besselk, :beta, :lbeta)

    # vec-vec only
    const binaryOps2 = (:+,:-)

    const reductions = ((:sum,:+),(:prod,:*),(:minimum,:min),(:maximum,:max))

    # expression functions
    vecTyp(n) = symbol(string("Vector",n))
    vecTypT(n) = Expr(:curly, vecTyp(n), :T)
    matTyp(r,c) = symbol(string("Matrix",r,"x",c))
    matTypT(r,c) = Expr(:curly, matTyp(r,c,), :T)
    elt(i) = symbol(string("e",i))
    col(i) = symbol(string("c",i))
    mem(s,e) = Expr(:.,s,Expr(:quote,e))
    velt(v,i) = mem(v,elt(i))
    melt(m,i,j) = mem(mem(m,col(j)),elt(i))

    # vector types
    for sz = 1:maxSz
        local Typ = vecTyp(sz)
        local TypT = vecTypT(sz)

        # the body of the type definition
        local defn = :(immutable $TypT <: ImmutableVector{T} end)

        # the members of the type
        for i = 1:sz
            local e = elt(i)
            push!(defn.args[3].args, :($e::T))
        end
        # instantiate the type definition
        eval(defn)

        # unary and n-ary constructors
        ctorn = :($TypT($([:($(symbol("a$i"))::T) for i=1:sz]...)) = $TypT($([:($(symbol("a$i"))) for i=1:sz]...)))
        ctor1 = :($TypT(a::T) = $TypT($([:(a) for i=1:sz]...)))

        eval(ctorn)
        eval(ctor1)

        # construct or convert from other vector types
        @eval $Typ(a::AbstractVector) = $Typ(ntuple($sz, x-> a[x])...)
        @eval convert{T}(::Type{$TypT}, x::AbstractVector) = $Typ(x)

        # convert to Array
        @eval begin
            function convert{T}(::Type{Vector{T}}, v::$TypT)
                a = Array(T,$sz)
                for i = 1:$sz
                    a[i] = v[i]
                end
                a
            end
        end

        # getindex
        local getix = :(error(BoundsError))
        for i = sz:-1:1
            local val = mem(:v,elt(i))
            getix = :(ix == $i ? $val : $getix)
        end
        getix = :(getindex{T}(v::$TypT, ix::Integer) = $getix)
        eval(getix)

        # helper for defining maps
        mapBody(f,j) = begin
            mp = :($Typ())
            for i = 1:sz
                local ff = copy(f)
                ff.args[j] = mem(:v,elt(i))
                push!(mp.args, ff)
            end
            mp
        end

        for op = unaryOps
            local bdy = mapBody(:($op(x)),2)
            @eval $op(v::$Typ) = $bdy
        end

        for op = binaryOps
            local bdy = :($Typ())
            for i = 1:sz
                local mem1 = mem(:v1,elt(i))
                local mem2 = mem(:v2,elt(i))
                push!(bdy.args, Expr(:call,op,mem1,mem2))
            end
            @eval $op(v1::$Typ,v2::$Typ) = $bdy

            bdy = mapBody(:($op(s,x)),3)
            if op == :.^ # special version for MathConst{:e}
                @eval $op(s::MathConst{:e},m::$Typ) = $bdy
            end
            if op == :min || op == :max
                @eval $op{T2<:Real}(s::T2,v::$Typ) = $bdy
            else
                @eval $op(s::Number,v::$Typ) = $bdy
            end

            bdy = mapBody(:($op(x,s)),2)
            if op == :min || op == :max
                @eval $op{T2<:Real}(v::$Typ,s::T2) = $bdy
            else
                @eval $op(v::$Typ,s::Number) = $bdy
            end
        end

        for op = binaryOps2
            local bdy = :($Typ())
            for i = 1:sz
                local mem1 = mem(:v1,elt(i))
                local mem2 = mem(:v2,elt(i))
                push!(bdy.args, Expr(:call,op,mem1,mem2))
            end
            @eval $op(v1::$Typ,v2::$Typ) = $bdy
        end

        for pr = reductions
            local bdy = Expr(:call,pr[2])
            for i = 1:sz
                push!(bdy.args, mem(:v,elt(i)))
            end
            local meth = pr[1]
            @eval $meth(v::$Typ) = $bdy
        end

        # convert to column matrix
        local colMatT = matTypT(sz,1)
        @eval column{T}(v::$TypT) = $colMatT(v)

        # convert to row matrix
        local rowMat = Expr(:call,matTyp(1,sz))
        for i = 1:sz
            local val = mem(:v,elt(i))
            push!(rowMat.args, :(Vector1{T}($val)))
        end
        @eval row{T}(v::$TypT) = $rowMat

        # vector norms
        @eval norm{T}(v::$TypT) = sqrt(dot(v,v))
        @eval norm{T}(v::$TypT,p::Number) = begin
            if p == 1
                sum(abs(v))
            elseif p == 2
                norm(v)
            elseif p == Inf
                max(abs(v))
            else
                norm(copy(v),p)
            end
        end

        # standard basis vectors
        local bdy = :($TypT())
        for j = 1:sz
            push!(bdy.args, :(i==$j?one(T):zero(T)))
        end
        @eval unit{T}(::Type{$TypT}, i::Integer) = $bdy

        # diagonal matrix
        sqMatTypT = matTypT(sz,sz)
        bdy = :($sqMatTypT())
        for i = 1:sz
            push!(bdy.args, :(v[$i].*unit($TypT,$i)))
        end
        @eval diagm{T}(v::$TypT) = $bdy

        # elementwise type conversion
        bdy = mapBody(:(convert(T,x)),3)
        @eval convert{T}(::Type{$TypT}, v::$Typ) = $bdy

        # some one-liners
        @eval similar{T}(::$TypT, t::DataType, dims::Dims) = Array(t, dims)
        @eval size(::$Typ) = ($sz,)
        @eval zero{T}(::Type{$TypT}) = $Typ(zero(T))
        @eval dot{T}(v1::$TypT,v2::$TypT) = sum(v1.*conj(v2))
        @eval unit{T}(v::$TypT) = v/norm(v)
    end

    # matrix types
    for rSz = 1:maxSz, cSz = 1:maxSz
        local Typ = matTyp(rSz,cSz)
        local TypT = matTypT(rSz,cSz)
        local ColTyp = vecTyp(rSz)
        local ColTypT = vecTypT(rSz)
        local RowTyp = vecTyp(cSz)
        local RowTypT = vecTypT(cSz)

        # the body of the type definition
        local defn = :(immutable $TypT <: ImmutableMatrix{T} end)

        # the members of the type
        for i = 1:cSz
            local c = col(i)
            push!(defn.args[3].args, :($c::$ColTypT))
        end

        # instantiate the type definition
        eval(defn)

        # unary and n-ary constructors
        ctorn = :($TypT($([:($(symbol("a$i"))::$ColTypT) for i=1:cSz]...)) = $TypT($([:($(symbol("a$i"))) for i=1:cSz]...)))
        ctor1 = :($TypT(a::$ColTypT) = $TypT($([:(a) for i=1:cSz]...)))

        eval(ctorn)
        eval(ctor1)

        # construction from a scalar
        @eval $TypT(a::T) = $Typ($ColTyp(a))

        # construct or convert from other matrix types
        @eval $Typ(a::AbstractMatrix) = $Typ(ntuple($cSz, c->
            $ColTyp(ntuple($rSz, r-> a[r,c])...))...)
        @eval convert{T}(::Type{$TypT}, x::AbstractMatrix) = $Typ(x)

        # convert to Array
        @eval begin
            function convert{T}(::Type{Matrix{T}}, m::$TypT)
                a = Array(T,$rSz,$cSz)
                for i = 1:$rSz, j = 1:$cSz
                    a[i,j] = m[i,j]
                end
                a
            end
        end

        # column access
        local cl = :(error(BoundsError))
        for j = cSz:-1:1
            local val = mem(:m,col(j))
            cl = :(ix == $j ? $val : $cl)
        end
        @eval column{T}(m::$TypT, ix::Integer) = $cl

        # row access
        local rw = :(error(BoundsError))
        for i = rSz:-1:1
            local rowexp = :($RowTypT())
            for j = 1:cSz
                push!(rowexp.args, mem(mem(:m,col(j)),elt(i)))
            end
            rw = :(ix == $i ? $rowexp : $rw)
        end
        @eval row{T}(m::$TypT, ix::Integer) = $rw

        # getindex
        @eval getindex{T}(m::$TypT, i::Integer, j::Integer) = column(m,j)[i]
        @eval getindex{T}(m::$TypT, ix::Integer) =
            getindex(m,mod(ix-1,$rSz)+1,div(ix-1,$rSz)+1)

        # ctranspose
        local bdy = Expr(:call, matTypT(cSz,rSz))
        for i = 1:rSz
            local rw = :($RowTypT())
            for j = 1:cSz
                local val = mem(mem(:m,col(j)),elt(i))
                push!(rw.args, :(conj($val)))
            end
            push!(bdy.args, rw)
        end
        @eval ctranspose{T}(m::$TypT) = $bdy

        # helper for defining maps
        mapBody(f,k) = begin
            local bdy = :($Typ())
            for j = 1:cSz
                local cl = :($ColTyp())
                for i = 1:rSz
                    local ff = copy(f)
                    ff.args[k] = mem(mem(:m,col(j)),elt(i))
                    push!(cl.args, ff)
                end
                push!(bdy.args, cl)
            end
            bdy
        end

        for op = unaryOps
            local bdy = mapBody(:($op(x)),2)
            @eval $op(m::$Typ) = $bdy
        end

        for op = binaryOps
            local bdy = :($Typ())
            for j = 1:cSz
                local cl = :($ColTyp())
                for i = 1:rSz
                    push!(cl.args,
                          Expr(:call,op,
                               mem(mem(:m1,col(j)),elt(i)),
                               mem(mem(:m2,col(j)),elt(i))))
                end
                push!(bdy.args, cl)
            end
            @eval $op(m1::$Typ,m2::$Typ) = $bdy

            bdy = mapBody(:($op(s,x)),3)
            if op == :.^ # special version for MathConst{:e}
                @eval $op(s::MathConst{:e},m::$Typ) = $bdy
            end
            if op == :min || op == :max
                @eval $op{T2<:Real}(s::T2,m::$Typ) = $bdy
            else
                @eval $op(s::Number,m::$Typ) = $bdy
            end

            bdy = mapBody(:($op(x,s)),2)
            if op == :min || op == :max
                @eval $op{T2<:Real}(m::$Typ,s::T2) = $bdy
            else
                @eval $op(m::$Typ,s::Number) = $bdy
            end
        end

        for op = binaryOps2
            local bdy = :($Typ())
            for j = 1:cSz
                local cl = :($ColTyp())
                for i = 1:rSz
                    push!(cl.args,
                          Expr(:call,op,
                               mem(mem(:m1,col(j)),elt(i)),
                               mem(mem(:m2,col(j)),elt(i))))
                end
                push!(bdy.args, cl)
            end
            @eval $op(m1::$Typ,m2::$Typ) = $bdy
        end

        for pr = reductions
            local bdy = Expr(:call,pr[2])
            for i = 1:rSz, j = 1:cSz
                push!(bdy.args, mem(mem(:m,col(j)),elt(i)))
            end
            local meth = pr[1]
            @eval $meth(m::$Typ) = $bdy
        end

        # vector-matrix multiplication
        bdy = :($RowTypT())
        for j = 1:cSz
            local e = :(+())
            for i = 1:rSz
                push!(e.args,
                      Expr(:call, :*,
                           mem(:v,elt(i)),
                           mem(mem(:m,col(j)),elt(i))))
            end
            push!(bdy.args, e)
        end
        @eval *{T}(v::$ColTypT,m::$TypT) = $bdy

        # matrix-vector multiplication
        bdy = :($ColTypT())
        for i = 1:rSz
            local e = :(+())
            for j = 1:cSz
                push!(e.args,
                      Expr(:call, :*,
                           melt(:m,i,j),
                           velt(:v,j)))
            end
            push!(bdy.args, e)
        end
        @eval *{T}(m::$TypT,v::$RowTypT) = $bdy

        # vector-matrix-vector multiplication
        bdy = :(+())
        for i = 1:rSz, j = 1:cSz
            push!(bdy.args,
                  Expr(:call, :*,
                       mem(:vl,elt(i)),
                       mem(mem(:m,col(j)),elt(i)),
                       mem(:vr,elt(j))))
        end
        @eval *{T}(vl::$ColTypT,m::$TypT,vr::$RowTypT) = $bdy

        # identity
        bdy = :($TypT())
        for j = 1:cSz
            push!(bdy.args, :(unit($ColTypT,$j)))
        end
        @eval eye{T}(::Type{$TypT}) = $bdy

        # matrix diagonal
        diagSz = min(rSz,cSz)
        diagTypT = vecTypT(diagSz)
        bdy = :($diagTypT())
        for i = 1:diagSz
            push!(bdy.args, mem(mem(:m,col(i)),elt(i)))
        end
        @eval diag{T}(m::$TypT) = $bdy

        # elementwise type conversion
        bdy = mapBody(:(convert(T,x)),3)
        @eval convert{T}(::Type{$TypT}, m::$Typ) = $bdy

        # some one-liners
        @eval similar{T}(::$TypT, t::DataType, dims::Dims) = Array(t, dims)
        @eval size(::$Typ) = ($rSz,$cSz)
        @eval zero{T}(::Type{$TypT}) = $Typ($ColTyp(zero(T)))
    end

    # matrix-matrix multiplication
    for n = 1:maxSz, p = 1:maxSz, m = 1:maxSz
        local bdy = Expr(:call, matTypT(n,m))
        for j = 1:m
            local c = Expr(:call, vecTypT(n))
            for i = 1:n
                local e = :(+())
                for k = 1:p
                    push!(e.args,
                          Expr(:call, :*,
                               melt(:m1,i,k),
                               melt(:m2,k,j)))
                end
                push!(c.args, e)
            end
            push!(bdy.args, c)
        end
        local m1T = matTypT(n,p)
        local m2T = matTypT(p,m)
        @eval *{T}(m1::$m1T,m2::$m2T) = $bdy
    end

    # matrix determinant and inverse
    for sz in 1:maxSz
        local Typ = matTyp(sz,sz)
        local TypT = matTypT(sz,sz)

        @eval det{T}(a::$TypT) = det(convert(Array{T,2},a))
        @eval inv{T}(a::$TypT) = $Typ(inv(convert(Array{T,2},a)))
    end

    # cross products
    if maxSz >= 2
        @eval cross(a::Vector2,b::Vector2) = a.e1*b.e2-a.e2*b.e1
    end

    if maxSz >= 3
        @eval cross(a::Vector3,b::Vector3) = Vector3(a.e2*b.e3-a.e3*b.e2,
                                                     a.e3*b.e1-a.e1*b.e3,
                                                     a.e1*b.e2-a.e2*b.e1)
    end
end


text:
@doc raw"""
`DAE`: Differential Algebraic Equation

Defines a differential algebraic initial value problem
```math
\begin{align*}
\dot{q} (t) &= v(t, q(t)) + u(t, q(t), \lambda(t)) , & q(t_{0}) &= q_{0} , \\
0 &= \phi (t, q(t), \lambda(t)) , & \lambda(t_{0}) &= \lambda_{0} ,
\end{align*}
```
with vector field ``v``, projection ``u``, algebraic constraint ``\phi=0``,
initial conditions ``q_{0}`` and ``\lambda_{0}``, the dynamical variable ``q``
taking values in ``\mathbb{R}^{m}`` and the algebraic variable ``\lambda``
taking values in ``\mathbb{R}^{n}``.

### Fields

* `d`: dimension of dynamical variable ``q`` and the vector field ``v``
* `m`: dimension of algebraic variable ``\lambda`` and the constraint ``\phi``
* `n`: number of initial conditions
* `v`: function computing the vector field
* `u`: function computing the projection
* `ϕ`: algebraic constraint
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `λ₀`: initial condition for algebraic variable ``\lambda``

The function `v`, providing the vector field, takes three arguments,
`v(t, q, v)`, the functions `u` and `ϕ`, providing the projection and the
algebraic constraint take four arguments, `u(t, q, λ, u)` and `ϕ(t, q, λ, ϕ)`,
where `t` is the current time, `q` and `λ` are the current solution vectors,
and `v`, `u` and `ϕ` are the vectors which hold the result of evaluating the
vector field ``v``, the projection ``u`` and the algebraic constraint ``\phi``
on `t`, `q` and `λ`.

### Example

```julia
    function v(t, q, v)
        v[1] = q[1]
        v[2] = q[2]
    end

    function u(t, q, λ, u)
        u[1] = +λ[1]
        u[2] = -λ[1]
    end

    function ϕ(t, q, λ, ϕ)
        ϕ[1] = q[2] - q[1]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ₀ = [0.]

    dae = DAE(v, u, ϕ, t₀, q₀, λ₀)

```
"""
struct DAE{dType <: Number, tType <: Number, vType <: Function, uType <: Function,
           ϕType <: Function, pType <: Union{Tuple,Nothing}, N} <: AbstractEquationDAE{dType, tType}

    d::Int
    m::Int
    n::Int
    v::vType
    u::uType
    ϕ::ϕType
    t₀::tType
    q₀::Array{dType, N}
    λ₀::Array{dType, N}
    parameters::pType
    periodicity::Vector{dType}

    function DAE(DT::DataType, N::Int, d::Int, m::Int, n::Int,
                 v::vType, u::uType, ϕ::ϕType, t₀::tType,
                 q₀::AbstractArray{dType}, λ₀::AbstractArray{dType};
                 parameters=nothing, periodicity=zeros(DT,d)) where {
                        dType <: Number, tType <: Number, vType <: Function,
                        uType <: Function, ϕType <: Function}

        @assert d == size(q₀,1)
        @assert m == size(λ₀,1)
        @assert n == size(q₀,2) == size(λ₀,2)
        @assert d ≥ m
        @assert ndims(q₀) == ndims(λ₀) == N ∈ (1,2)

        new{DT, tType, vType, uType, ϕType, typeof(parameters), N}(d, m, n, v, u, ϕ, t₀,
                convert(Array{DT}, q₀), convert(Array{DT}, λ₀), parameters, periodicity)
    end
end

function DAE(v, u, ϕ, t₀, q₀::AbstractArray{DT}, λ₀::AbstractArray{DT}; kwargs...) where {DT}
    DAE(DT, ndims(q₀), size(q₀,1), size(λ₀,1), size(q₀,2), v, u, ϕ, t₀, q₀, λ₀; kwargs...)
end

function DAE(v, u, ϕ, q₀, λ₀; kwargs...)
    DAE(v, u, ϕ, zero(eltype(q₀)), q₀, λ₀; kwargs...)
end

Base.hash(dae::DAE, h::UInt) = hash(dae.d, hash(dae.m, hash(dae.n, hash(dae.v,
        hash(dae.u, hash(dae.t₀, hash(dae.q₀, hash(dae.λ₀,
        hash(dae.periodicity, hash(dae.parameters, h))))))))))

Base.:(==)(dae1::DAE, dae2::DAE) = (
                                dae1.d == dae2.d
                             && dae1.m == dae2.m
                             && dae1.n == dae2.n
                             && dae1.v == dae2.v
                             && dae1.u == dae2.u
                             && dae1.t₀ == dae2.t₀
                             && dae1.q₀ == dae2.q₀
                             && dae1.λ₀ == dae2.λ₀
                             && dae1.parameters == dae1.parameters
                             && dae1.periodicity == dae1.periodicity)

function Base.similar(dae::DAE, q₀, λ₀=get_λ₀(q₀, dae.λ₀); kwargs...)
    similar(dae, dae.t₀, q₀, λ₀; kwargs...)
end

function Base.similar(dae::DAE, t₀::TT, q₀::AbstractArray{DT}, λ₀::AbstractArray{DT}=get_λ₀(q₀, dae.λ₀);
                      parameters=dae.parameters, periodicity=dae.periodicity) where {DT  <: Number, TT <: Number}
    @assert dae.d == size(q₀,1)
    @assert dae.m == size(λ₀,1)
    DAE(dae.v, dae.u, dae.ϕ, t₀, q₀, λ₀; parameters=parameters, periodicity=periodicity)
end

@inline Base.ndims(dae::DAE) = dae.d

@inline periodicity(equation::DAE) = equation.periodicity


text:
using CompEcon
using Dierckx
using Statistics

linspace(a,b,n) = range(a, stop=b, length=n)

# Create functions to test against
@. f1(x, y) = sin(x) * cos(y)
@. f2(x, y) = exp(x) * cos(y)
@. f3(x, y) = 1.5 + log(x * y)
fs = [f1, f2, f3]
nf = length(fs)

x = collect(linspace(.5, 2, 35))
y = collect(linspace(.2, 4, 35))

function compare_D_CE_levels(f::Function, x::Array{Float64, 1}, y::Array{Float64, 1}, order::Int)

    # Create Dierckx Spline (and data)
    n = length(x)
    z = f(x, y')
    dspl = Spline2D(x, y, z; s=0.0, kx=order, ky=order)

    # Create CompEcon Spline (and data)
    cebasis_x = Basis(SplineParams(x, 0, order))
    cebasis_y = Basis(SplineParams(y, 0, order))
    cebasis = Basis(cebasis_x, cebasis_y)
    xx = nodes(cebasis)[1]
    yy = f(xx[:, 1], xx[:, 2])
    c, bs = funfitxy(cebasis, xx, yy)

    # Evaluate Splines on finer grid
    xfine = collect(linspace(x[1], x[end], 2*n + 1))
    yfine = collect(linspace(y[1], y[end], 2*n + 1))
    zfine = f(xfine, yfine')

    deval = Dierckx.evalgrid(dspl, xfine, yfine)
    derr = abs.(zfine .- deval)

    ceeval = CompEcon.funeval(c, cebasis, gridmake(xfine, yfine))
    ceerr = abs.(zfine .- reshape(ceeval, 2*n+1, 2*n+1))

    return derr, ceerr
end

for i=1:nf

    curr_f = fs[i]

    derr1_levels, ceerr1_levels = compare_D_CE_levels(curr_f, x, y, 1)
    derr2_levels, ceerr2_levels = compare_D_CE_levels(curr_f, x, y, 2)
    derr3_levels, ceerr3_levels = compare_D_CE_levels(curr_f, x, y, 3)

    println("Linear Approximation to Function")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr1_levels)) \n\t Min: $(minimum(derr1_levels)) \n\t Mean: $(mean(derr1_levels))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr1_levels)) \n\t Min: $(minimum(ceerr1_levels)) \n\t Mean: $(mean(ceerr1_levels))")

    println("Quadratic Approximation to Function")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr2_levels)) \n\t Min: $(minimum(derr2_levels)) \n\t Mean: $(mean(derr2_levels))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr2_levels)) \n\t Min: $(minimum(ceerr2_levels)) \n\t Mean: $(mean(ceerr2_levels))")

    println("Cubic Approximation to Function")
    println("Dierckx Level Errors \n\t Max: $(maximum(derr3_levels)) \n\t Min: $(minimum(derr3_levels)) \n\t Mean: $(mean(derr3_levels))")
    println("CompEcon Level Errors \n\t Max: $(maximum(ceerr3_levels)) \n\t Min: $(minimum(ceerr3_levels)) \n\t Mean: $(mean(ceerr3_levels))")

end


text:
# ------------ Lightweight visitor ------------ #

const Visitor = DynamicTarget

function visit!(vs::Visitor, addr)
    haskey(vs.tree, addr) && error("VisitorError (visit!): already visited address $(addr).")
    push!(vs, addr)
end

function compare(am::AddressMap, vs::Visitor)
    missed = []
    for (k, v) in shallow_iterator(am)
        k in vs || push!(missed, k)
    end
    isempty(missed), missed
end


text:
module testGrouping

using Base.Test
using DataFrames
using TimeData
using Datetime
include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running grouping tests\n")

## create core data
nObs = 10
varNames = [:a, :b, :c, :d, :e, :f]
simVals = rand(nObs, 6)
dats = Date{ISOCalendar}[date(2013, 7, ii) for ii=1:nObs]

## create groups
groups = {:group1 => [:a, :c, :e], :group2 => [:b, :d],
          :group3 => [:f]}

## test groupVars
tm = Timematr(simVals, varNames, dats)
tmGroup = AssetMgmt.groupVars(tm, groups, x -> sum(x, 2))

## test groupAssets
invs = AssetMgmt.equWgtInvestments(tm)
groupedInvs = AssetMgmt.groupAssets(invs, groups)
groupedInvsSorted = groupedInvs[[:group1, :group2, :group3]]

## calculate expected result
expDf = AssetMgmt.composeDataFrame(repmat([0.5 1/3 1/6], nObs, 1),
                                   [:group1 :group2 :group3])
expInvs = AssetMgmt.Investments(expDf, dats)

@test isequal(expInvs, groupedInvsSorted)

groupedTm = AssetMgmt.groupAssets(invs, groups, x -> sum(x, 2))
end


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

# AEDAT 1.0
v1_data = [0x7f, 0x22, 0x01, 0x39, 0xa4, 0xd0]

## Test reading bytes (correct byte order etc):
evt = read(IOBuffer(v1_data), AEDATEvent{UInt16})
@assert evt == AEDATEvent{UInt16}(0x7f22, 20554960)

## Determine file version:
dat = AEDATFile(IOBuffer(v1_data))
@assert isa(dat, AEDATFile{1,0})

## Test parsing an event:
dat = AEDATFile(IOBuffer(v1_data))
evt = read_event(dat)
@assert AER.isevent(DVS128, evt)

evt2 = convert(AEDATEvent{DVS128}, evt)
@assert evt2 == AEDATEvent{DVS128}(DVS128(false, 17, 127, false), 20554960)

dat = AEDATFile(IOBuffer(v1_data))
evt3 = read_event(dat, DVS128)
@assert evt2 == evt3

# AEDAT 2.0
v2_data = [0xa9, 0xcc, 0xe2, 0x00, 0x10, 0x21, 0x07, 0xb2]
v2_header = "#!AER-DAT2.0\r\n# Foo\r\n"

## Test reading bytes (correct byte order etc):
evt = read(IOBuffer(v2_data), AEDATEvent{UInt32})
@assert evt == AEDATEvent{UInt32}(0xa9cce200, 270600114)

## Test reading file version and header:
function v2_io()
    io = IOBuffer()
    print(io, v2_header)
    write(io, v2_data)
    seekstart(io)
    io
end
dat = AEDATFile(v2_io())
@assert isa(dat, AEDATFile{2,0}) dat
@assert dat.header == [" Foo\r\n"]

## Test default event types:
mark(dat.io)
evt1 = read_event(dat)
reset(dat.io)
mark(dat.io)
evt2 = read_event(dat, UInt32)
@assert evt1 == evt2

## Test parsing an event:
@assert AER.isevent(DAVIS240_APS, evt1)
reset(dat.io)
mark(dat.io)
evt = read_event(dat, DAVIS240_APS)
@assert evt == AEDATEvent{DAVIS240_APS}(DAVIS240_APS(167, 206, 0, 512), 270600114)


text:
# Persistent landscape

scale(p::Pair) = (p.first + p.second)/2 => (p.second - p.first)/2
scale(ps::Vector{<:Pair}) = map(scale, sort(ps))

struct Landscape
    levels::Vector{Vector{<:Pair}}
end
Landscape(ints::Vector{<:Pair}...) = Landscape([ints...])
getindex(ls::Landscape, i::Integer) = ls.levels[i]
length(ls::Landscape) = length(ls.levels)
show(io::IO, ls::Landscape) = print(io, "λ[k = $(length(ls.levels))]")
function Base.dump(io::IOContext, ls::Landscape; kwargs...)
    print(io, "Landscape:")
    for (i, l) in enumerate(ls.levels)
        print(io, "\n\tk = $i: ")
        for itr in l
            print(io, "$itr ")
        end
    end
end
function iterate(ls::Landscape, state=nothing)
    state === nothing && return iterate(ls, 1)
    state > length(ls) && return nothing
    return ls[state], state+1
end

function value(i1::Pair, i2::Pair, x)
    a = (i2.second - i1.second)/(i2.first - i1.first)
    b = i1.second - (a == 0 ? 0 : a*i1.first)
    return a*x+b
end

"""Construction of a persistent landscape from a barcode

Algorithm is adopted form "A Persistence Landscapes Toolbox for Topological Statistics" by  Bubenik & Dlotko, p. 22
"""
function landscape(bar::Vector{<:Pair})
    A = sort(bar)
    k = 0
    L = Vector{Vector{Pair}}()
    while length(A) > 0
        k += 1
        push!(L, Vector{Pair}())
        b, d  = A[1]
        deleteat!(A, 1)
        p = 1
        if (b, d) == (-Inf, Inf)
            push!(L[k],  -Inf => Inf)
        else
            if d == Inf
                push!(L[k], -Inf  => 0.0)
                push!(L[k], b => 0.0)
                push!(L[k], Inf  => Inf)
            else
                if b == -Inf
                    push!(L[k], -Inf => Inf)
                else
                    push!(L[k], -Inf => 0.0)
                    push!(L[k], b => 0.0)
                    push!(L[k], scale(b => d))
                end
            end
        end
        while !(L[k][end] == (Inf => 0.0) || L[k][end] == (Inf => Inf))
            i = findfirst(x->x.second > d, A[p:end])
            if i === nothing
                push!(L[k], d => 0.0)
                push!(L[k], Inf => 0.0)
            else
                b′, d′  = A[p+i-1]
                deleteat!(A, p+i-1)
                b′ > d && push!(L[k], d => 0.0)
                if b′ >= d
                    push!(L[k], b′ => 0.0)
                else
                    push!(L[k], scale(b′ => d))
                    insert!(A, p+i-1, b′ => d)
                    p += 1
                end
                if d′ == Inf
                    push!(L[k], Inf => Inf)
                else
                    push!(L[k], scale(b′ => d′))
                end
                b, d  = b′, d′
            end
        end
    end
    return Landscape(L)
end

"""
    reduce(op, l1::Landscape, x::Number)

Perform a reduction `op` between landscape and number.
"""
function reduce(op, l::Landscape, x::Number)
    n = length(l)
    L = Vector{Vector{Pair}}(undef, n)
    for i in 1:n
        λᵢ = Pair[]
        for j in 1:length(l[i])
            push!(λᵢ, l[i][j].first => op(l[i][j].second, x))
        end
        L[i] = λᵢ
    end
    return Landscape(L);
end

"""
    reduce(op, l1::Landscape, l2::Landscape)

Perform a reduction `op` of two persistence landscapes.
"""
function reduce(op, l1::Landscape, l2::Landscape)
    n = min(length(l1), length(l2))
    m = max(length(l1), length(l2))
    L = Vector{Vector{Pair}}(undef, m)
    for i in 1:n
        λᵢ = Pair[]
        p = q = 1
        while p+1 <= length(l1[i]) && q+1 <= length(l2[i])
            if l1[i][p].first < l2[i][q].first
                v = value(l2[i][q-1], l2[i][q], l1[i][p].first)
                push!(λᵢ, l1[i][p].first => op(l1[i][p].second, v))
                @debug "L" p=p q=q λᵢ=λᵢ[end] val=v
                p+=1
            elseif l1[i][p].first > l2[i][q].first
                v = value(l1[i][p-1], l1[i][p], l2[i][q].first)
                push!(λᵢ, l2[i][q].first => op(v, l2[i][q].second))
                @debug "G" p=p q=q λᵢ=λᵢ[end] val=v
                q+=1
            else # l1[i][p].first == l2[i][q].first
                push!(λᵢ, l2[i][q].first => op(l1[i][p].second, l2[i][q].second))
                @debug "E" p=p q=q λᵢ=λᵢ[end]
                p+=1
                q+=1
            end
        end
        while p+1 < length(l1[i]) && q+1 >= length(l2[i])
            push!(λᵢ, l1[i][p].first => op(l1[i][p].second, 0))
            @debug "L1" p=p q=q λᵢ=λᵢ[end]
            p+=1
        end
        while p+1 >= length(l1[i]) && q+1 < length(l2[i])
            push!(λᵢ, l2[i][q].first => op(0, l2[i][q].second))
            @debug "L2" p=p q=q λᵢ=λᵢ[end]
            q+=1
        end
        push!(λᵢ, Inf => 0.0)
        L[i] = λᵢ
    end
    if length(l1) > n
        for i in n+1:m
            λᵢ = copy(l1[i])
            for j in 1:length(λᵢ)
                λᵢ[j] = l1[i][j].first => op(l1[i][j].second, 0)
            end
            L[i] = λᵢ
            @debug "L1N:" i=i n=n m=m λᵢ=λᵢ
        end
    end
    if length(l2) > n
        for i in n+1:m
            λᵢ = copy(l2[i])
            for j in 1:length(λᵢ)
                λᵢ[j] = l2[i][j].first => op(0, l2[i][j].second)
            end
            L[i] = λᵢ
            @debug "L2N:" i=i n=n m=m λᵢ=λᵢ
        end
    end
    return Landscape(L)
end

"""
    mean(ls::Vector{Landscape})

Calculate avarage landscape.
"""
function mean(ls::Vector{Landscape})
    n = length(ls)
    m = foldl((l,v) -> reduce(+, l, v), ls)
    return reduce(/, m, n)
end


text:
include("testdata.jl")

@testset "expandtable(Single)" begin
    dat = expandtable(ACTmath.scale, ACTmath.xcount)
    ft = freqtab(dat)
    @test ft.table.freq == ACTmath.xcount[begin+1:end]
end


text:
using PartialRejectionSampling
using Distributions
using Random
using Plots

# Pedagogy
function hard_core_spatial_pedagogy(
    pp::PRS.HardCorePointProcess{T};
    win::Union{Nothing,PRS.AbstractWindow}=nothing,
    rng=-1
)::Vector{T} where {T}

    path(i) = joinpath(@__DIR__,
                       "output",
                       "hard_core_spatial",
                       join([lpad(i, 3, "0"), ".pdf"]))

    window_ = isnothing(win) ? PRS.window(pp) : win

    n = rand(rng, Distributions.Poisson(pp.β * PRS.volume(window_)))
    points = Matrix{Float64}(undef, PRS.dimension(pp), n)
    for x in eachcol(points)
        x .= rand(rng, window_)
    end

    i = 0
    while true

        p = pedagogy_plot(points, true, 0, "white", pp.window)
        Plots.savefig(p, path(i))

        p = pedagogy_plot(points, true, pp.r/2, "white", pp.window)
        i+=1
        Plots.savefig(p, path(i))

        bad = vec(any(PRS.pairwise_distances(points) .< pp.r, dims=2))
        !any(bad) && break

        p = pedagogy_plot(points[:, .!bad], true, 0, "white", pp.window)

        pedagogy_plot!(p, points[:, bad], true, pp.r/2, "red", pp.window)
        i+=1
        Plots.savefig(p, path(i))

        p = pedagogy_plot(points[:, .!bad], true, 0, "white", pp.window)
        pedagogy_plot!(p, points[:, bad], false, pp.r, "orange", pp.window)
        i+=1
        Plots.savefig(p, path(i))

        resampled = PRS.generate_sample_poisson_union_balls(rng, pp.β, points[:, bad], pp.r, window_)
        pedagogy_plot!(p, resampled, true, 0, "blue", pp.window)
        i+=1
        Plots.savefig(p, path(i))

        points = hcat(points[:, .!bad], resampled)
        i += 1
    end

    i += 1
    p = pedagogy_plot(points, true, 0, "white", pp.window)
    Plots.savefig(p, path(i))

    return [points[:, i] for i in 1:size(points, 2)]
end

function pedagogy_plot!(
    p,
    points,
    show_center=true,
    radius=0,
    color="white",
    window=SquareWindow(zeros(2), 1)
)
    for x in (points isa Matrix ? eachcol(points) : points)
        if radius > 0
            plot_disk!(p, x, radius, color)
        end
        if show_center
            Plots.scatter!(p, [x[1]], [x[2]], markersize=2, color=color, grid=false)
        end
    end

    Plots.xlims!(window.c[1], window.c[1] + window.w)
    Plots.ylims!(window.c[2], window.c[2] + window.w)

    return p
end

function pedagogy_plot(
    points,
    show_center=true,
    radius=0,
    color="white",
    window=SquareWindow(zeros(2), 1)
)

    p = Plots.plot([0], [0],
            label="", legend=false,
            color="white",
            linewidth=0,
            aspect_ratio=:equal,
            grid=false,
            title="")

    pedagogy_plot!(p, points, show_center, radius, color, window)

    return p
end

function plot_disk!(p, center, radius, color=:white)
    θ = LinRange(0, 2π, 15)
    x, y = center[1] .+ radius .* cos.(θ), center[2] .+ radius .* sin.(θ)

    Plots.plot!(p, x, y, seriestype=[:shape], c=color, linewidth=0.5, legend=false, fillapha=0.2)

    return p
end

function plot(
    pp::PRS.AbstractSpatialPointProcess,
    points;
    title=""
)
    p = Plots.plot([0], [0],
            label="", legend=false,
            color="white",
            linewidth=0,
            aspect_ratio=:equal,
            grid=false,
            title=title)

    θ = collect(range(0, 2π, length=15))
    rad = pp.r / 2  # radius = interaction range / 2
    circ_x, circ_y = rad .* cos.(θ), rad .* sin.(θ)

    for x in points
        Plots.plot!(x[1] .+ circ_x,
                    x[2] .+ circ_y,
                    color="black",
                    linewidth=0.5,
                    grid=false)
    end

    win = pp.window
    Plots.xlims!(win.c[1], win.c[1] + win.w[1])
    Plots.ylims!(win.c[2], win.c[2] + (win.w isa Number ? win.w[1] : win.w[2]))

    return p
end


text:
function generatecolours(N::Int)
    # oversample since the HSV starts from red, and ends at red.
    M = round(Int, 1.3*N)

    colours = range(Colors.HSV(0,1,1), stop = Colors.HSV(-360,1,1), length = M)
    colours = convert(Vector{Colors.RGB}, colours)
    colours = collect( [colours[n].r; colours[n].g; colours[n].b] for n = 1:N )

    return colours
end

function visualize2Dpartition(X_parts::Vector{Vector{Vector{T}}},
    y_set,
    t_set,
    fig_num::Int,
    title_string::String;
    display_range = 1:length(X_parts),
    new_fig_flag = true) where T <: Real

    colours_for_pts = generatecolours(length(X_parts))

    if new_fig_flag
        PyPlot.figure(fig_num)
        fig_num += 1
    end

    ax = PyPlot.axes()

    # partitions of points.
    for i in display_range

        # points.
        X = X_parts[i]

        x1 = collect( X[n][1] for n = 1:length(X))
        x2 = collect( X[n][2] for n = 1:length(X))

        PyPlot.scatter(x1, x2, label = "$(i)", color = colours_for_pts[i])

        # annotate centroid.
        z1, z2 = Statistics.mean(X)
        # PyPlot.annotate("$(i)",
        # xy=[z1;z2],
        # xycoords="data",
        # fontsize=20.0,
        # #bbox = ["boxstyle"=>"rarrow,pad=0.3", "fc"=>"cyan", "ec"=>"b", "lw"=>2])
        # bbox = (z1,z2, z1+0.5, z2+0.5) )

        #PyPlot.text(z1, z2, "$(i)", backgroundcolor = "white")
        bbox = Dict("boxstyle"=>"round, pad=0.3", "fc"=>"cyan", "ec"=>"b", "lw"=>2)
        PyPlot.text(z1, z2, "$(i)", bbox = bbox)
    end

    ### boundaries of partitions.
    @assert length(t_set) == length(y_set)
    @assert !isempty(X_parts[1][1])

    #colours_for_boundaries = generatecolours(length(t_set))

    for i = 1:length(t_set)

       # PyPlot.plot(t_set[i], y_set[i], color = colours_for_boundaries[i])
       PyPlot.plot(t_set[i], y_set[i], color = "black")
    end

    PyPlot.axis("scaled")
    PyPlot.title(title_string)
    PyPlot.legend()

    return fig_num, ax
end

function visualizesingleregions(X_set, y_set, t_set, fig_num;
    display_region_range = 1:length(X_set))

    for i in display_region_range
        fig_num, ax = visualize2Dpartition(X_set, y_set, t_set, fig_num,
        "X_set, levels = $(levels)"; display_range = i:i)

        # ## I am here. debug the visualization of the complement. then, mixGP, use X_set as training set.
        # # display the rest of the points.
        # Xc_set = collect( X_set[n] for n = 1:length(X_set) if n != i )
        # Xc_flatten = combinevectors(combinevectors(Xc_set[]))

        # Xc_display = similar(X_set)
        # resize!(Xc_display, 1)
        # Xc_display[1] = Xc_flatten

        # fig_num, ax = visualize2Dpartition(Xc_display, y_set, t_set, fig_num,
        # "Xc_set, levels = $(levels)"; display_range = i:i, new_fig_flag = false)

    end

    return fig_num
end

function visualizemeshgridpcolorx1horizontal( x_ranges::Vector{LinRange{T,L}},
    Y::Matrix{T},
    marker_locations::Vector,
    marker_symbol::String,
    fig_num::Int,
    title_string::String;
    x1_title_string::String = "Dimension 1",
    x2_title_string::String = "Dimension 2",
    cmap = "Greens_r") where {T <: Real,L}

    #
    @assert length(x_ranges) == 2
    x_coords = collect( collect(x_ranges[d]) for d = 1:2 )

    PyPlot.figure(fig_num)
    fig_num += 1
    PyPlot.pcolormesh(x_coords[1], x_coords[2], Y, cmap = cmap, shading = "auto")
    PyPlot.xlabel(x1_title_string)
    PyPlot.ylabel(x2_title_string)
    PyPlot.title(title_string)

    for i = 1:length(marker_locations)
        #pt = reverse(marker_locations[i])
        pt = marker_locations[i]
        PyPlot.annotate(marker_symbol, xy=pt, xycoords="data")
    end

    PyPlot.plt.colorbar()
    PyPlot.axis("scaled")

    return fig_num
end

##### debug.
## debug.
function parsedebugvarsset(N_regions::Int, debug_vars_set, dummy_val::T) where T

    Nq = length(debug_vars_set)

    # allocate.
    Ws_tilde = Vector{Vector{T}}(undef, Nq)
    Us = Vector{Vector{T}}(undef, Nq)
    Rs = Vector{Vector{Int}}(undef, Nq)
    Ps = Vector{Int}(undef, Nq)

    for j = 1:Nq

        ### initialize.
        Ws_tilde[j] = Vector{T}(undef, N_regions)
        Us[j] = Vector{T}(undef, N_regions)

        # debug: initialize to known value.
        fill!(Ws_tilde[j], 0.0)
        fill!(Us[j], Inf)

        ### parse.
        w_tilde = debug_vars_set[j].w_tilde_set[1]

        u = debug_vars_set[j].u_set[1]
        # v = debug_vars_set[j].v_set[1]

        region_inds = debug_vars_set[j].region_inds_set[1]
        p_region_ind = debug_vars_set[j].p_region_ind_set[1]

        @assert length(u) == length(w_tilde) == length(region_inds) + 1

        for i = 1:length(region_inds)

            r = region_inds[i]
            Ws_tilde[j][r] = w_tilde[i]
            Us[j][r] = u[i]
        end
        Ws_tilde[j][p_region_ind] = one(T)
        Us[j][p_region_ind] = u[end]

        Ps[j] = p_region_ind
        Rs[j] = region_inds

        #hps_keep_flags_set = debug_vars_set[j].hps_keep_flags_set_set[1]
        #zs_set = debug_vars_set[j].zs_set_set[1]
        #ts_set = debug_vars_set[j].ts_set_set[1]
    end

    return Ws_tilde, Us, Ps, Rs
end

function findneighbourpartitions(p::Vector{T},
    radius::T,
    root,
    levels,
    hps,
    p_region_ind::Int;
    δ::T = 1e-10) where T

    region_inds = Vector{Int}(undef, length(hps))
    j = 0

    keep_flags = falses(length(hps)) # for debug.
    zs = Vector{Vector{T}}(undef, length(hps)) # for debug.
    ts = Vector{T}(undef, length(hps)) # for debug.

    for i = 1:length(hps)

        # parse.
        u = hps[i].v
        c = hps[i].c

        # intersection with hyperplane.
        t = -dot(u, p) + c
        z = p + t .* u

        zs[i] = z # debug.
        ts[i] = t

        if norm(z-p) < radius

            # get two points along ray normal to the plane, on either side of plane.
            z1 = p + (t+δ) .* u
            z2 = p + (t-δ) .* u

            # get z1's, z2's region.
            z1_region_ind = PatchMixtureKriging.findpartition(z1, root, levels)
            z2_region_ind = PatchMixtureKriging.findpartition(z2, root, levels)

            # println("i = ", i)
            # println("z1 = ", z1)
            # println("z2 = ", z2)
            # println("z1_region_ind = ", z1_region_ind)
            # println("z2_region_ind = ", z2_region_ind)
            # println()

            # do not keep if both are in p's region.
            # keep if either z1 or z2 is in p's region.
            #if !(z2_region_ind == p_region_ind && z1_region_ind == p_region_ind) &&
            #    (z2_region_ind == p_region_ind || z1_region_ind == p_region_ind)
            if xor(z2_region_ind == p_region_ind, z1_region_ind == p_region_ind)

                keep_flags[i] = true

                j += 1
                region_inds[j] = z1_region_ind
                if z1_region_ind == p_region_ind
                    region_inds[j] = z2_region_ind
                end

            end
        end
    end

    resize!(region_inds, j)

    return region_inds, ts, zs, keep_flags
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
todo"add ability to show a Cfunction's signature, so perhaps storing arg names is needed"

struct Cfunction{RetT, ArgsT<:Tuple, ConvT<:Cconvention}
	let constructor = false end
end

Base.Ptr{FuncT}(args...) where {FuncT<:Cfunction} = FuncT(args...)

Cfunction{RetT, ArgsT, ConvT}(ptr::Ptr{Cvoid}) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention} = reinterpret(Ptr{Cfunction{RetT, ArgsT, ConvT}}, ptr)
Cfunction{RetT, ArgsT}(ptr::Ptr{Cvoid}) where {RetT, ArgsT<:Tuple} = reinterpret(Ptr{Cfunction{RetT, ArgsT}}, ptr)

Cfunction{RetT, ArgsT, ConvT}(lib::Clibrary, sym::Symbol) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention} = Cfunction{RetT, ArgsT, ConvT}(sym, lib)
Cfunction{RetT, ArgsT}(lib::Clibrary, sym::Symbol) where {RetT, ArgsT<:Tuple} = Cfunction{RetT, ArgsT}(sym, lib)

Cfunction{RetT, ArgsT, ConvT}(sym::Symbol, lib::Clibrary, libs::Clibrary...) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention} = Cfunction{RetT, ArgsT, ConvT}(_dlsym(sym, lib, libs...))
Cfunction{RetT, ArgsT}(sym::Symbol, lib::Clibrary, libs::Clibrary...) where {RetT, ArgsT<:Tuple} = Cfunction{RetT, ArgsT}(_dlsym(sym, lib, libs...))

# NOTE:  this returns a tuple (since the user must retain a reference to func for the function pointer to remain usable)
Cfunction{RetT, ArgsT}(func::Base.CFunction) where {RetT, ArgsT<:Tuple} = (Cfunction{RetT, ArgsT, default_convention(ArgsT)}(Base.unsafe_convert(Ptr{Cvoid}, func)), func)

@generated function Cfunction{RetT, ArgsT}(func::Function) where {RetT, ArgsT<:Tuple}
	_tuplize(::Type{Tuple{}}) where {T<:Tuple} = ()
	_tuplize(::Type{T}) where {T<:Tuple} = (Base.tuple_type_head(T), _tuplize(Base.tuple_type_tail(T))...,)

	types = _tuplize(ArgsT)
	return quote
		$(Expr(:meta, :inline))
		return Cfunction{RetT, ArgsT}(Base.@cfunction($(Expr(:$, :func)), RetT, ($(types...),)))
	end
end

convention(::Type{Cfunction{RetT, ArgsT, ConvT}}) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention} = convention(ConvT)
convention(::Type{Cconvention{SymT}}) where {SymT} = SymT
convention(conv::Cconvention) = convention(typeof(conv))

todo"determine the correct default conventions for each platform"
todo"some compilers use different calling convention for variadic functions"
todo"identify what is Julia's calling convention"
default_convention(::Type{ArgsT}) where {ArgsT<:Tuple} = ifelse(Sys.iswindows(), STDCALL, CDECL)

# https://www.gnu.org/software/libc/manual/html_node/How-Variadic.html
# it appears that chars and ints are promoted to 32-bit ints and floats are promoted to doubles when used as varargs
cconvert_vararg(::Type{CDECL}, ::Type{T}) where {T} = T
cconvert_vararg(::Type{CDECL}, ::Type{T}) where {T<:Signed} = ifelse(sizeof(T) < sizeof(Cint), Cint, T)
cconvert_vararg(::Type{CDECL}, ::Type{T}) where {T<:Unsigned} = ifelse(sizeof(T) < sizeof(Cuint), Cuint, T)
cconvert_vararg(::Type{CDECL}, ::Type{T}) where {T<:Cchar} = ifelse(sizeof(T) < sizeof(Cint), Cint, T)
cconvert_vararg(::Type{CDECL}, ::Type{T}) where {T<:AbstractFloat} = ifelse(sizeof(T) < sizeof(Cdouble), Cdouble, T)

# these define a type to add to the type tuple when ccalling variadic functions
cconvert_default(::Type{T}) where {T} = error("Sorry, `$(T)` is not yet usable for variadic functions, provide `CBinding.cconvert_default(::Type{$(T)}) = CConvertType` to define how to pass it to a variadic function call")
cconvert_default(::Type{T}) where {T<:AbstractString} = Cstring
cconvert_default(::Type{T}) where {T<:Char} = Cchar
cconvert_default(::Type{T}) where {T<:Union{Integer, AbstractFloat, Ref}} = T   # this covers Ptr too, Ptr <: Ref
cconvert_default(::Type{T}) where {E, T<:AbstractArray{E}} = Ptr{E}

(f::Ptr{Cfunction{RetT, ArgsT}})(args...) where {RetT, ArgsT<:Tuple} = Cfunction{RetT, ArgsT, default_convention(ArgsT)}(reinterpret(Ptr{Cvoid}, f))(args...)
(f::Ptr{Cfunction{RetT, ArgsT, ConvT}})(args...) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention} = invoke(reinterpret(Ptr{Cfunction{concrete(RetT), ArgsT, ConvT}}, f), args...)
@generated function invoke(f::Ptr{Cfunction{RetT, ArgsT, ConvT}}, args...) where {RetT, ArgsT<:Tuple, ConvT<:Cconvention}
	error = nothing
	_tuplize(::Type{Tuple{}}) = ()
	_tuplize(::Type{Tuple{}}, argT, argsT...) = (error = :(throw(MethodError(f, args))) ; ())
	_tuplize(::Type{Tuple{Vararg}}) = (:(Ptr{Cvoid}...),)  # HACK: extra `Ptr{Cvoid}...` is being added to trigger vararg ccall behavior rather than regular behavior (if there is a difference in the backend)
	_tuplize(::Type{Tuple{Vararg}}, argT, argsT...) = (cconvert_vararg(ConvT, cconvert_default(argT)), _tuplize(Tuple{Vararg}, argsT...)...,)
	_tuplize(::Type{T}) where {T<:Tuple} = (error = :(throw(MethodError(f, args))) ; ())
	_tuplize(::Type{T}, argT, argsT...) where {T<:Tuple} = (Base.tuple_type_head(T), _tuplize(Base.tuple_type_tail(T), argsT...)...,)

	types = _tuplize(ArgsT, args...)
	return !isnothing(error) ? error : quote
		$(Expr(:meta, :inline))
		ccall(reinterpret(Ptr{Cvoid}, f), $(convention(ConvT)), RetT, ($(types...),), $(map(i -> :(args[$(i)]), eachindex(args))...),)
	end
end