text:
import LinearAlgebra.BLAS: gemv, gemv!, gemm!, trsm!, axpy!, ger!
import LinearAlgebra: cholesky, Cholesky
import Base: getproperty

Base.@deprecate chol(X) cholesky(X).U

#=
See [1] for implementation details: pages 5-9 in particular. The derivations presented in
[1] assume column-major layout, whereas Julia primarily uses row-major. We therefore
implement both the derivations in [1] and their transpose, which is more appropriate to
Julia.

[1] - "Differentiation of the Cholesky decomposition", Murray 2016
=#

const AM = AbstractMatrix
const UT = UpperTriangular

@explicit_intercepts cholesky Tuple{AbstractMatrix{<:∇Scalar}}
∇(::typeof(cholesky), ::Type{Arg{1}}, p, U::Cholesky, Ū::AM{T}, Σ::AM{T}) where T<:∇Scalar =
    chol_blocked_rev(Matrix(Ū), Matrix(U.U), 25, true)

@explicit_intercepts getproperty Tuple{Cholesky, Symbol} [true, false]
function ∇(::typeof(getproperty), ::Type{Arg{1}}, p, y, ȳ, C::Cholesky, x::Symbol)
    if x === :U
        C.uplo === 'U' ? UpperTriangular(ȳ) : LowerTriangular(ȳ')
    elseif x === :L
        C.uplo === 'L' ? LowerTriangular(ȳ) : UpperTriangular(ȳ')
    else
        throw(ArgumentError("unrecognized field $x; use U or L"))
    end
end

@explicit_intercepts(
    Cholesky,
    Tuple{AbstractMatrix{<:∇Scalar}, Union{Char, Symbol}, Integer},
    [true, false, false],
)
function ∇(
    ::Type{Cholesky},
    ::Type{Arg{1}},
    p,
    C::Cholesky,
    X̄::Union{UpperTriangular, LowerTriangular},
    X::Union{UpperTriangular, LowerTriangular},
    uplo::Union{Char, Symbol},
    info::Integer,
)
    # We aren't doing any actual computation if we've constructed a Cholesky object
    # directly, so just pass through this call and return the sensitivies
    return X̄
end

"""
    level2partition(A::AbstractMatrix, j::Int, upper::Bool)

Returns views to various bits of the lower triangle of `A` according to the
`level2partition` procedure defined in [1] if `upper` is `false`. If `upper` is `true` then
the transposed views are returned from the upper triangle of `A`.
"""
function level2partition(A::AM, j::Int, upper::Bool)

    # Check that A is square and j is a valid index.
    M, N = size(A)
    (0 >= j || j > M) && throw(ArgumentError("j is out of range."))
    M != N && throw(ArgumentError("A is not square."))

    if upper
        r = view(A, 1:j-1, j)
        d = view(A, j, j)
        B = view(A, 1:j-1, j+1:N)
        c = view(A, j, j+1:N)
    else
        r = view(A, j, 1:j-1)
        d = view(A, j, j)
        B = view(A, j+1:N, 1:j-1)
        c = view(A, j+1:N, j)
    end
    return r, d, B, c
end

"""
    level3partition(A::AbstractMatrix, j::Int, k::Int, upper::Bool)

Returns views to various bits of the lower triangle of `A` according to the
`level3partition` procedure defined in [1] if `upper` is `false`. If `upper` is `true` then
the transposed views are returned from the upper triangle of `A`.
"""
function level3partition(A::AM, j::Int, k::Int, upper::Bool)

    # Check that A is square and j is a valid index.
    M, N = size(A)
    (0 >= j || j > M) && throw(ArgumentError("j is out of range."))
    M != N && throw(ArgumentError("A is not square."))

    # Get views into bits of A.
    if upper
        R = view(A, 1:j-1, j:k)
        D = view(A, j:k, j:k)
        B = view(A, 1:j-1, k+1:N)
        C = view(A, j:k, k+1:N)
    else
        R = view(A, j:k, 1:j-1)
        D = view(A, j:k, j:k)
        B = view(A, k+1:N, 1:j-1)
        C = view(A, k+1:N, j:k)
    end
    return R, D, B, C
end

"""
    chol_unblocked_rev!(
        Ā::AbstractMatrix{T},
        L::AbstractMatrix{T},
        upper::Bool
    ) where T<:Real

Compute the reverse-mode sensitivities of the Cholesky factorisation in an unblocked manner.
If `upper` is `false`, then the sensitivites computed from and stored in the lower triangle
of `Ā` and `L` respectively. If `upper` is `true` then they are computed and stored in the
upper triangles. If at input `upper` is `false` and `tril(Ā) = L̄`, at output
`tril(Ā) = tril(Σ̄)`, where `Σ = LLᵀ`. Analogously, if at input `upper` is `true` and
`triu(Ā) = triu(Ū)`, at output `triu(Ā) = triu(Σ̄)` where `Σ = UᵀU`.
"""
function chol_unblocked_rev!(Σ̄::AM{T}, L::AM{T}, upper::Bool) where T<:Real

    # Check that L is square, that Σ̄ is square and that they are the same size.
    M, N = size(Σ̄)
    M != N && throw(ArgumentError("Σ̄ is not square."))

    # Compute the reverse-mode diff.
    j = N
    for ĵ in 1:N
        r, d, B, c = level2partition(L, j, upper)
        r̄, d̄, B̄, c̄ = level2partition(Σ̄, j, upper)

        # d̄ <- d̄ - c'c̄ / d.
        d̄[1] -= dot(c, c̄) / d[1]

        # [d̄ c̄'] <- [d̄ c̄'] / d.
        d̄ ./= d
        c̄ ./= d

        # r̄ <- r̄ - [d̄ c̄'] [r' B']'.
        r̄ = axpy!(-Σ̄[j, j], r, r̄)
        r̄ = gemv!(upper ? 'N' : 'T', -one(T), B, c̄, one(T), r̄)

        # B̄ <- B̄ - c̄ r.
        B̄ = upper ? ger!(-one(T), r, c̄, B̄) : ger!(-one(T), c̄, r, B̄)
        d̄ ./= 2
        j -= 1
    end
    return (upper ? triu! : tril!)(Σ̄)
end
chol_unblocked_rev(Σ̄::AM, L::AM, upper::Bool) = chol_unblocked_rev!(copy(Σ̄), L, upper)

"""
    chol_blocked_rev!(
        Σ̄::AbstractMatrix{T},
        L::AbstractMatrix{T},
        Nb::Int,
        upper::Bool
    ) where T<:∇Scalar

Compute the sensitivities of the Cholesky factorisation using a blocked, cache-friendly
procedure. `Σ̄` are the sensitivities of `L`, and will be transformed into the sensitivities
of `Σ`, where `Σ = LLᵀ`. `Nb` is the block-size to use. If the upper triangle has been used
to represent the factorization, that is `Σ = UᵀU` where `U := Lᵀ`, then this should be
indicated by passing `upper = true`.
"""
function chol_blocked_rev!(Σ̄::AM{T}, L::AM{T}, Nb::Int, upper::Bool) where T<:∇Scalar

    # Check that L is square, that Σ̄ is square and that they are the same size.
    M, N = size(Σ̄)
    M != N && throw(ArgumentError("Σ̄ is not square."))

    tmp = Matrix{T}(undef, Nb, Nb)

    # Compute the reverse-mode diff.
    k = N
    if upper
        for k̂ in 1:Nb:N
            j = max(1, k - Nb + 1)
            R, D, B, C = level3partition(L, j, k, true)
            R̄, D̄, B̄, C̄ = level3partition(Σ̄, j, k, true)

            C̄ = trsm!('L', 'U', 'N', 'N', one(T), D, C̄)
            gemm!('N', 'N', -one(T), R, C̄, one(T), B̄)
            gemm!('N', 'T', -one(T), C, C̄, one(T), D̄)
            chol_unblocked_rev!(D̄, D, true)
            gemm!('N', 'T', -one(T), B, C̄, one(T), R̄)
            if size(D̄, 1) == Nb
                tmp = axpy!(one(T), D̄, transpose!(tmp, D̄))
                gemm!('N', 'N', -one(T), R, tmp, one(T), R̄)
            else
                gemm!('N', 'N', -one(T), R, D̄ + D̄', one(T), R̄)
            end

            k -= Nb
        end
        return triu!(Σ̄)
    else
        for k̂ in 1:Nb:N
            j = max(1, k - Nb + 1)
            R, D, B, C = level3partition(L, j, k, false)
            R̄, D̄, B̄, C̄ = level3partition(Σ̄, j, k, false)

            C̄ = trsm!('R', 'L', 'N', 'N', one(T), D, C̄)
            gemm!('N', 'N', -one(T), C̄, R, one(T), B̄)
            gemm!('T', 'N', -one(T), C̄, C, one(T), D̄)
            chol_unblocked_rev!(D̄, D, false)
            gemm!('T', 'N', -one(T), C̄, B, one(T), R̄)
            if size(D̄, 1) == Nb
                tmp = axpy!(one(T), D̄, transpose!(tmp, D̄))
                gemm!('N', 'N', -one(T), tmp, R, one(T), R̄)
            else
                gemm!('N', 'N', -one(T), D̄ + D̄', R, one(T), R̄)
            end

            k -= Nb
        end
        return tril!(Σ̄)
    end
end
chol_blocked_rev(Σ̄::AbstractMatrix, L::AbstractMatrix, Nb::Int, upper::Bool) =
    chol_blocked_rev!(copy(Σ̄), L, Nb, upper)


text:
using JuMP, GLPK, Ipopt
using DualDecomposition

const DD = DualDecomposition
const parallel = DD.parallel

"""
a: interest rate
π: unit stock price
ρ: unit dividend price

K: number of stages
L: number of stock types
2^L scenarios in each stage
2^L^(K-1)=16 scenarios in total
ρ = 0.05 * π
bank: interest rate 0.01
stock1: 1.03 or 0.97
stock2: 1.06 or 0.94
...

b_k: initial asset (if k=1) and income (else)
B_k: money in bank
x_{k,l}: number of stocks to buy/sell (integer)
y_{k,l}: total stocks

deterministic model:

    max     B_K+∑_{l=1}^{L}π_{K,l}y_{K,l}

    s.t.    B_1+∑_{l=1}^{L}π_{1,l}x_{1,l} = b_1

            b_k+(1+a)B_{k-1}+∑_{l=1}^{L}ρ_{k,l}y_{k-1,l} = B_k+∑_{l=1}^{L}π_{k,l}x_{k,l}, ∀ k=2,…,K

            y_{1,l} = x_{1,l}, ∀ l=1,…,L

            y_{k-1,l}+x_{k,l} = y_{k,l}, ∀ k=2,…,K, l=1,…,L

            x_{k,l} ∈ ℤ , ∀ k=1,…,K, l=1,…,L

            y_{k,l} ≥ 0, ∀ k=1,…,K, l=1,…,L

            B_k ≥ 0, ∀ k=1,…,K.
"""
const K = 3
const L = 2
const a = 0.01
const b_init = 100  # initial capital
const b_in = 30   # income

# iteratively add nodes
# root node
function create_nodes()::DD.Tree
    ξ = Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}(:π => ones(L))
    tree = DD.Tree(ξ)

    #subproblem formulation
    function subproblem_builder(tree::DD.Tree, subtree::DD.SubTree, node::DD.SubTreeNode)
        mdl = subtree.model
        x = @variable(mdl, x[l=1:L], Int, base_name="n1_x")
        #x = @variable(mdl, x[l=1:L], base_name="n1_x")

        y = @variable(mdl, y[l=1:L] >= 0, base_name="n1_y")
        DD.set_output_variable!(node, :y, y)

        B = @variable(mdl, B >= 0, base_name="n1_B")
        DD.set_output_variable!(node, :B, B)

        π = DD.get_scenario(node)[:π]
        @constraints(mdl,
            begin
                B + sum( π[l] * x[l] for l in 1:L) == b_init
                [l=1:L], y[l] - x[l] == 0
            end
        )
        DD.set_stage_objective(node, 0.0)

        JuMP.unregister(mdl, :x)
        JuMP.unregister(mdl, :y)
        JuMP.unregister(mdl, :B)
    end

    DD.set_stage_builder!(tree, 1, subproblem_builder)

    create_nodes!(tree, 1)
    return tree
end

# child nodes
function create_nodes!(tree::DD.Tree, pt::Int)
    for scenario = 1:2^L
        prob = 1/2^L
        π = get_realization(DD.get_scenario(tree, pt)[:π], scenario)
        ξ = Dict{Symbol, Union{Float64,<:AbstractArray{Float64}}}(:π => π)
        id = DD.add_child!(tree, pt, ξ, prob)

        #subproblem formulation
        function subproblem_builder(tree::DD.Tree, subtree::DD.SubTree, node::DD.SubTreeNode)
            mdl = subtree.model
            id = DD.get_id(node)
            #x = @variable(mdl, x[l=1:L], Int, base_name = "n$(id)_x")
            x = @variable(mdl, x[l=1:L], base_name = "n$(id)_x")

            y = @variable(mdl, y[l=1:L] >= 0, base_name = "n$(id)_y")
            DD.set_output_variable!(node, :y, y)

            B = @variable(mdl, B >= 0, base_name = "n$(id)_B")
            DD.set_output_variable!(node, :B, B)

            y_ = @variable(mdl, y_[l=1:L] >= 0, base_name = "n$(id)_y_")
            DD.set_input_variable!(node, :y, y_)

            B_ = @variable(mdl, B_ >= 0, base_name = "n$(id)_B_")
            DD.set_input_variable!(node, :B, B_)

            π = DD.get_scenario(node)[:π]
            pt = DD.get_parent(node)
            ρ = DD.get_scenario(tree, pt)[:π] * 0.05
            @constraint(mdl, B + sum( π[l] * x[l] - ρ[l] * y_[l] for l in 1:L) - (1+a) * B_ == b_in)
            @constraint(mdl, [l=1:L], y[l] - x[l] - y_[l] == 0)

            #dummy bound for input variables to avoid subproblem unboundedness
            @constraint(mdl, [l=1:L], y_[l] <= 500)
            @constraint(mdl, B_ <= 500)
            if DD.get_stage(node) < K
                DD.set_stage_objective(node, 0.0)
            else
                DD.set_stage_objective(node, -(B + sum( π[l] * y[l] for l in 1:L )))
            end
            JuMP.unregister(mdl, :x)
            JuMP.unregister(mdl, :y)
            JuMP.unregister(mdl, :B)
            JuMP.unregister(mdl, :y_)
            JuMP.unregister(mdl, :B_)
        end

        DD.set_stage_builder!(tree, id, subproblem_builder)
        if DD.get_stage(tree, id) < K
            create_nodes!(tree, id)
        end
    end
end

# construct realization event
function get_realization(ξ::Array{Float64,1}, scenario::Int)::Array{Float64,1}
    ret = ones(L)
    multipliers = digits(scenario - 1, base=2, pad=L)*2 - ones(L)
    for l = 1:L
        ret[l] = ξ[l] * (1 + multipliers[l] * l * 0.03)
    end
    return ret
end

tree = create_nodes()
#node_cluster = DD.decomposition_not(tree)
node_cluster = DD.decomposition_scenario(tree)
#node_cluster = DD.decomposition_temporal(tree) #There is a DUAL_INFEASIBLE issue

# Number of block components
NS = length(node_cluster)

# Initialize MPI
parallel.init()

# Create DualDecomposition instance.
algo = DD.LagrangeDual()

# partition scenarios into processes
parallel.partition(NS)

coupling_variables = Vector{DD.CouplingVariableRef}()
models = Dict{Int,JuMP.Model}()

for block_id in parallel.getpartition()
    nodes = node_cluster[block_id]
    subtree = DD.create_subtree!(tree, block_id, coupling_variables, nodes)
    set_optimizer(subtree.model, GLPK.Optimizer)
    DD.add_block_model!(algo, block_id, subtree.model)
    models[block_id] = subtree.model
end

# Set nonanticipativity variables as an array of symbols.
DD.set_coupling_variables!(algo, coupling_variables)

# Lagrange master method
LM = DD.BundleMaster(BM.ProximalMethod, optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
#LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

# Solve the problem with the solver; this solver is for the underlying bundle method.
DD.run!(algo, LM)

# Write timing outputs to files
DD.write_all(algo)

# Finalize MPI
parallel.finalize()

ext: jl


text:
using Base: OneTo

using Base.Threads: @spawn
using Base: @sync

export Panel, to_band!

struct Panel
    w::Int
    h::Int
end

raw"""
Computes the amount of τs we have to store for all
reflectors (I - τvv') in a tsqr. If we have 2 stacked panels
of width w, we need 3 small qr's using w reflectors each, so
3w τs.
x   x
 \ /
  x

Another example is 5 stacked panels, will go like this:
x   x   x   x   x
 \ /     \ /   /
  x       x   /
   \     /   /
    \   /   /
     \ /   /
      x   /
       \ /
        x
It's effectively 9 small qr's with w reflectors, so 9w τs.
Generaly the formula is simply panel_width * (2 * num_panels_vertically - 1)
"""
function reflectors_per_tsqr(matrix_height, p::Panel)
    num_panels_vertically = ÷(matrix_height, p.h, RoundUp)
    reflectors_per_panel = p.w #assumes p.h ≥ p.w
    return reflectors_per_panel * (2 * num_panels_vertically - 1)
end

function to_band!(A::AbstractMatrix{T}, p = Panel(16, 32)) where {T}
    @assert size(A, 1) == size(A, 2)
    @assert p.w ≤ p.h

    # on top we have a section of size p.w x p.w
    m = size(A, 1)

    τs = Vector{T}(undef, reflectors_per_tsqr(m - p.w, p))

    # loop over all column blocks
    for b = OneTo(÷(m - p.w, p.w, RoundUp))

        # divide into the bit we're QR'ing
        # the bit to which we apply it from the left
        # and the bit to which we apply it from the right.
        curr_block = (b - 1) * p.w + 1
        next_block = curr_block + p.w

        A_qr = view(A, next_block:m, curr_block:next_block-1)
        A_left = view(A, next_block:m, next_block:m)
        A_right = view(A, :, next_block:m)

        @timeit TIMER "tsqr!" tsqr!(A_qr, τs, p)
        @timeit TIMER "apply_tsqr_left!" apply_tsqr_left!(A_qr, A_left, τs, p)
        @timeit TIMER "apply_tsqr_right!" apply_tsqr_right!(A_qr, A_right, τs, p)
    end

    return A
end

function exponent_of_2_greater_than_or_equal_to(m)
    k = 0
    two_to_the_k = 1
    while two_to_the_k < m
        two_to_the_k *= 2
        k += 1
    end

    return k
end

function tsqr!(A, τs, p::Panel)
    m, n = size(A)

    @assert n == p.w

    num_panels_vertically = ÷(m, p.h, RoundUp)

    # Keep the τ offset conveniently here, cause
    # we need that when walking up the tree later.
    τ_offset = 1

    # First standard QR on every panel.
    #=@sync=# for b = OneTo(num_panels_vertically)
        b_start = (b - 1) * p.h + 1
        b_end = min(b_start + p.h - 1, m)
        τs′ = view(τs, τ_offset:τ_offset+p.w-1)
        A′ = view(A, b_start:b_end, :)

        #=@spawn=# block_qr!(A′, τs′)

        τ_offset += p.w
    end

    # Then we walk up the tree and do QR on upper triangular blocks.
    tree_height = exponent_of_2_greater_than_or_equal_to(num_panels_vertically)

    # we're combining two panels of height `merged_block_size`
    merged_block_size = p.h

    # we've already done the leaf nodes, so - 1.
    for level = OneTo(tree_height)
        let num_panels_vertically = ÷(m, merged_block_size, RoundUp), merged_block_size=merged_block_size

            # combine all pairs of panels.
            #=@sync=# for b = OneTo(num_panels_vertically ÷ 2)
                b_start = (b - 1) * 2merged_block_size + 1
                b_end = min(b_start + 2merged_block_size - 1, m)

                A′ = view(A, b_start:b_end, :)
                τs′ = view(τs, τ_offset:τ_offset+p.w-1)

                #=@spawn=# merge_block_qr!(A′, τs′, merged_block_size)

                τ_offset += p.w
            end
        end

        merged_block_size *= 2
    end
end

function apply_tsqr_left!(qr, A, τs, p::Panel)
    m, n = size(qr)

    @assert n == p.w

    num_panels_vertically = ÷(m, p.h, RoundUp)

    # Keep the τ offset conveniently here, cause
    # we need that when walking up the tree flater.
    τ_offset = 1

    # First standard QR on every panel.
    #=@sync=# for b = OneTo(num_panels_vertically)
        b_start = (b - 1) * p.h + 1
        b_end = min(b_start + p.h - 1, m)
        τs′ = view(τs, τ_offset:τ_offset+p.w-1)

        qr′ = view(qr, b_start:b_end, :)
        A′ = view(A, b_start:b_end, :)

        #=@spawn=# apply_tiny_Q_left!(qr′, A′, τs′)

        τ_offset += p.w
    end

    # Then we walk up the tree and do QR on upper triangular blocks.
    tree_height = exponent_of_2_greater_than_or_equal_to(num_panels_vertically)

    # we're combining two panels of height `merged_block_size`
    merged_block_size = p.h

    # we've already done the leaf nodes, so - 1.
    for level = OneTo(tree_height)
        let num_panels_vertically = ÷(m, merged_block_size, RoundUp), merged_block_size=merged_block_size

            # combine all pairs of panels.
            #=@sync=# for b = OneTo(num_panels_vertically ÷ 2)
                b_start = (b - 1) * 2merged_block_size + 1
                b_end = min(b_start + 2merged_block_size - 1, m)

                qr′ = view(qr, b_start:b_end, :)
                A′ = view(A, b_start:b_end, :)
                τs′ = view(τs, τ_offset:τ_offset+p.w-1)

                #=@spawn=# apply_tiny_Q_with_gap_left!(qr′, A′, τs′, merged_block_size)

                τ_offset += p.w
            end
        end

        merged_block_size *= 2
    end
end

function apply_tsqr_right!(qr, A, τs, p::Panel)
    m, n = size(qr)

    @assert n == p.w

    num_panels_vertically = ÷(m, p.h, RoundUp)

    # Keep the τ offset conveniently here, cause
    # we need that when walking up the tree flater.
    τ_offset = 1

    # First standard QR on every panel.
    #=@sync=# for b = OneTo(num_panels_vertically)
        b_start = (b - 1) * p.h + 1
        b_end = min(b_start + p.h - 1, m)
        τs′ = view(τs, τ_offset:τ_offset+p.w-1)

        qr′ = view(qr, b_start:b_end, :)
        A′ = view(A, :, b_start:b_end)

        #=@spawn=# apply_tiny_Q_right!(qr′, A′, τs′)

        τ_offset += p.w
    end

    # Then we walk up the tree and do QR on upper triangular blocks.
    tree_height = exponent_of_2_greater_than_or_equal_to(num_panels_vertically)

    # we're combining two panels of height `merged_block_size`
    merged_block_size = p.h

    # we've already done the leaf nodes, so - 1.
    for level = OneTo(tree_height)
        let num_panels_vertically = ÷(m, merged_block_size, RoundUp), merged_block_size=merged_block_size

            # combine all pairs of panels.
            #=@sync=# for b = OneTo(num_panels_vertically ÷ 2)
                b_start = (b - 1) * 2merged_block_size + 1
                b_end = min(b_start + 2merged_block_size - 1, m)

                qr′ = view(qr, b_start:b_end, :)
                A′ = view(A, :, b_start:b_end)
                τs′ = view(τs, τ_offset:τ_offset+p.w-1)

                #=@spawn=# apply_tiny_Q_with_gap_right!(qr′, A′, τs′, merged_block_size)

                τ_offset += p.w
            end
        end

        merged_block_size *= 2
    end
end

### just a simple dense QR.

function block_qr!(A, τs)
    m, n = size(A)

    steps = min(m - 1, n)

    # Create a bunch of reflectors and save the norm value (τ)
    @inbounds for k = OneTo(steps)
        x = view(A, k:m, k)
        τs[k] = reflector!(x)
        reflectorApplyLeft!(view(A, k:m, k+1:n), τs[k], x)
    end

    return A
end

### QR two stacked, upper triangular matrices.

function merge_block_qr!(A, τs, v_block_size)
    m, n = size(A)

    # Create a bunch of reflectors and save the norm value (τ)
    @inbounds for k = OneTo(n)
        gap = v_block_size - k + 1
        to = min(m, k + v_block_size)
        x = view(A, k:to, k)
        τs[k] = reflector!(x, gap)
        reflectorApplyLeft!(view(A, k:to, k+1:n), τs[k], x, gap)
    end

    return A
end

### Applying Q from a tiny QR decomp.

function apply_tiny_Q_left!(qr, A, τs)
    m, n = size(qr)

    steps = min(m - 1, n)

    @inbounds for k = OneTo(steps)
        reflectorApplyLeft!(view(A, k:size(A, 1), :), τs[k], view(qr, k:m, k))
    end

    return A
end

function apply_tiny_Q_right!(qr, A, τs)
    m, n = size(qr)

    steps = min(m - 1, n)

    @inbounds for k = OneTo(steps)
        reflectorApplyRight!(view(A, :, k:size(A, 2)), τs[k], view(qr, k:m, k))
    end

    return A
end

function apply_tiny_Q_with_gap_left!(qr, A, τs, v_block_size)
    m, n = size(qr)

    # Create a bunch of reflectors and save the norm value (τ)
    @inbounds for k = OneTo(n)
        gap = v_block_size - k + 1
        to = min(m, k + v_block_size)
        x = view(qr, k:to, k)
        reflectorApplyLeft!(view(A, k:to, :), τs[k], x, gap)
    end

    return A
end

function apply_tiny_Q_with_gap_right!(qr, A, τs, v_block_size)
    m, n = size(qr)

    # Create a bunch of reflectors and save the norm value (τ)
    @inbounds for k = OneTo(n)
        gap = v_block_size - k + 1
        to = min(m, k + v_block_size)
        x = view(qr, k:to, k)
        reflectorApplyRight!(view(A, :, k:to), τs[k], x, gap)
    end

    return A
end

#### Computing and applying a single reflector

@inline function reflector!(x::AbstractVector, gap = 1)
    @inbounds begin
        n = length(x)
        n == 0 && return zero(eltype(x))
        ξ1 = x[1]
        normu = abs2(ξ1)
        for i = 1+gap:n
            normu += abs2(x[i])
        end
        iszero(normu) && return zero(eltype(x))
        normu = sqrt(normu)
        ν = copysign(normu, real(ξ1))
        ξ1 += ν
        x[1] = -ν
        for i = 1+gap:n
            x[i] /= ξ1
        end
        ξ1/ν
    end
end

@inline function reflectorApplyLeft!(A::AbstractMatrix, τ::Number, x::AbstractVector, gap = 1)
    m, n = size(A)
    m == 0 && return A

    for j = 1:n
        Aj, xj = view(A, 1+gap:m, j), view(x, 1+gap:m)
        vAj = conj(τ)*(A[1, j] + dot(xj, Aj))
        A[1, j] -= vAj
        axpy!(-vAj, xj, Aj)
    end

    return A
end

@inline function reflectorApplyRight!(A::AbstractMatrix, τ::Number, x::AbstractVector, gap = 1)
    m, n = size(A)
    n == 0 && return A
    for j = 1:m
        Aj, xj = view(A, j, 1+gap:n), view(x, 1+gap:n)
        vAj = τ * (A[j, 1] + dot(Aj, xj))
        A[j, 1] -= vAj
        axpy!(-vAj, xj, Aj)
    end
    return A
end


text:
#######################################################################
# Schreier Vectors
#######################################################################

@inline Base.in(x, schr::Schreier) = x in schr.orb
@inline Base.push!(schr::Schreier, x) = push!(schr.orb, x)

function Base.getindex(schr::Schreier, pt)
    gr_word = Word(schr.gens_inv, schr.orb, pt, schr.op)
    return inv(gr_word(schr.gens_inv))
end

@inline Base.iterate(schr::Schreier) = iterate(schr.orb)
@inline Base.iterate(schr::Schreier, s) = iterate(schr.orb, s)
@inline Base.length(schr::Schreier) = length(schr.orb)
@inline Base.eltype(::Schreier{GEl,I}) where {GEl, I} = I
@inline Base.first(schr::Schreier) = first(schr.orb)
@inline Base.last(schr::Schreier) = last(schr.orb)
@inline islast(schr::Schreier, pt) = islast(schr.orb,pt)

function Base.:(==)(s1::Schreier,s2::Schreier)
    s1.gens_inv == s2.gens_inv || return false
    s1.orb == s2.orb || return false
    s1.op == s2.op || return  false
    return true
end

Base.hash(schr::S, h::UInt) where S<:Schreier = hash(S, hash(schr.gens_inv, hash(schr.orb, hash(schr.op, h))))

@doc doc"""
    Schreier(gens, pt[, op=^])
Compute the Schreier Tree of `pt` under group generated by `gens`.
"""
Schreier(gens::Vector{<:GroupElem}, pt, op=^) = Schreier(Orbit, gens, pt, op)

@doc doc"""
    Schreier(OrbT::Type{<:AbstractOrbit}, gens, pt[, op=^])
Compute the Schreier Tree of `pt` under group generated by `gens`.

The type `OrbT` will be used as a base storage for the orbit. `OrbT` must
implement _Orbit Protocol_`.
"""
function Schreier(OrbT::Type{<:AbstractOrbit}, gens::Vector{<:GroupElem}, pt, op=^)
    schr = Schreier(inv.(gens), OrbT(pt, 0), op)
    for o in schr
        for (idx, g) in enumerate(gens)
            γ = op(o, g)
            if γ ∉ schr
                push!(schr,(γ, idx))
            end
        end
    end
    return schr
end

@doc doc"""
    getinv(schr::Schreier{GEl, T}, pt::T)::GEl
While simple indexing into Schreier Trees returns coset representatives `uᵦ`,
this method returns `inv(uᵦ)`, i.e. the group element which takes `pt` to
`first(schr)`, the initial point of `schr`.
"""
function getinv(schr::Schreier, pt)
    gr_word = Word(schr.gens_inv, schr.orb, pt, schr.op)
    g = gr_word(schr.gens_inv)
    return g
end

function orbit_stabilizer(gens::Vector{GEl}, pt, op=^) where GEl<:GroupElem
    schr = Schreier(Orbit1, inv.(gens), pt)
    stab = Vector{GEl}()
    for o in schr
        for (idx,g) in enumerate(gens)
            γ = op(o, g)
            if γ ∉ schr
                push!(schr, (γ, g))
            else
                push!(stab, schr[o]*g*getinv(schr, γ))
            end
        end
    end
    return schr, PermGroup(unique!(stab))
end


text:
module Day03

using AdventOfCode2019

function day03(input::String = readInput(joinpath(@__DIR__, "..", "data", "day03.txt")))
    wires = split.(split(input), ",")
    paths = [Array{Tuple{UnitRange{Int},UnitRange{Int}},1}(undef, length(wire)) for wire in wires]
    corners = [Array{Tuple{Int,Int},1}(undef, length(wire) + 1) for wire in wires]
    corners[1][1], corners[2][1] = (0, 0), (0, 0)
    x, y = 0, 0
    for (iw, wire) in enumerate(wires)
        dirs = [a[1] for a in wire]
        lengths = [parse(Int, a[2:end]) for a in wire]
        for (j, (d, l)) in enumerate(zip(dirs, lengths))
            if d == 'R'
                paths[iw][j] = (x:x+l, y:y)
                corners[iw][j+1] = (x+l, y)
                x += l
            elseif d == 'L'
                paths[iw][j] = (x-l:x, y:y)
                corners[iw][j+1] = (x-l, y)
                x -= l
            elseif d == 'U'
                paths[iw][j] = (x:x, y:y+l)
                corners[iw][j+1] = (x, y+l)
                y += l
            elseif d == 'D'
                paths[iw][j] = (x:x, y-l:y)
                corners[iw][j+1] = (x, y-l)
                y -= l
            else
                throw(AssertionError("Invalid direction"))
            end
        end
        x, y = 0, 0
    end

    # Calculate the intersection points
    common = Array{Tuple{Int,Int},1}(undef,0)
    for seg1 in paths[1]
        for seg2 in paths[2]
            inter = intersect.(seg1, seg2)
            for x in inter[1]
                for y in inter[2]
                    push!(common, (x, y))
                end
            end
        end
    end

    # Evaluate the nearest intersection point
    dists = [abs(x) + abs(y) for (x, y) in common] |> sort
    dists[1] == 0 || throw(AssertionError("Shortest distance is not 0"))

    # Calculate the total distances to each intersection point for both wires
    pathDistances = [zeros(Int, length(common)) for _ in 1:2]
    for (i, path) in enumerate(paths)
        d = 0
        for (k, seg) in enumerate(path)
            for (j, p) in enumerate(common)
                if p[1] in seg[1] && p[2] in seg[2] && pathDistances[i][j] == 0
                    pathDistances[i][j] = d + abs(corners[i][k][1] - p[1]) + abs(corners[i][k][2] - p[2])
                end
            end
            d += sum(length.(seg)) - 2
        end
    end

    totalDistance = pathDistances[1] + pathDistances[2] |> sort
    totalDistance[1] == 0 || throw(AssertionError("Shortest distance is not 0"))

    return [dists[2], totalDistance[2]]
end

end # module


text:
module MarketData

import ..Currencies
using Dates

"""
# Interface

* `MarketData.get_value(provider, serie_id, at, asof::Date; locf::Bool=false) :: Union{Missing, T}`

* `MarketData.get_serie_currency(provider, serie_id) :: Union{Nothing, Currecies.Currency}`

* `MarketData.has_serie(provider, serie_id) :: Bool`

# Arguments

* `at` : a key to the value of the time series.

* `asof` : the closing date at which the time series is being assessed. This is usually set to the pricing date.

* `locf` : last observation carried forward. Use `locf=true` to repeat the latest observation prior to `at`.

# Provided methods

* `MarketData.has_value(provider, serie_id, at, asof::Date; locf::Bool=false) :: Bool`

* `MarketData.get_cash(provider, serie_id, at, asof::Date; locf::Bool=false) :: Union{Missing, Currencies.Cash}`

* `MarketData.assert_has_serie(provider, serie_id)`
"""
abstract type AbstractMarketDataProvider end

"""
    get_value(provider, serie_id, at, asof::Date; locf::Bool=false) :: Union{Missing, T}

Returns a value for `serie_id` at state (or date) `at`.

# Arguments

* `at` : a key to the value of the time series.

* `asof` : the date at which the time series is being assessed.

* `locf` : last observation carried forward. Use `locf=true` to repeat the latest observation prior to `at`.

This method should error if `provider` does not know about `serie_id`. Use [`MarketData.assert_has_serie`](@ref) for that.
"""
function get_value end

"""
    get_serie_currency(provider, serie_id) :: Union{Nothing, Currecies.Currency}

Returns the currency in which `serie_id` is expressed. Returns `nothing` if `serie_id` is not based on a currency.

This method should error if `provider` does not know about `serie_id`. Use [`MarketData.assert_has_serie`](@ref) for that.
"""
function get_serie_currency end
function has_serie end
has_value(provider::AbstractMarketDataProvider, serie_id, at, asof::Date; locf::Bool=false) = !ismissing(get_value(provider, serie_id, at, asof; locf=locf))

function get_cash(provider::AbstractMarketDataProvider, serie_id, at, asof::Date; locf::Bool=false) :: Union{Missing, Currencies.Cash}

    currency = get_serie_currency(provider, serie_id)
    if currency == nothing
        error("Serie $serie_id has no currency")
    end

    serie_value = get_value(provider, serie_id, at, asof; locf=locf)
    if ismissing(serie_value)
        return missing
    else
        return serie_value * currency
    end
end

"""
    assert_has_serie(provider::AbstractMarketDataProvider, serie_id)

Throws `AssertionError` if `has_serie(provider, serie_id)` returns `false`.
"""
function assert_has_serie(provider::AbstractMarketDataProvider, serie_id)
    @assert has_serie(provider, serie_id) "Provider does not know about serie $serie_id"
    nothing
end

"A Market Data provided that has no data. Useful when calculation does not depend on Market Data."
struct EmptyMarketDataProvider <: AbstractMarketDataProvider
end

get_value(::EmptyMarketDataProvider, serie_id, at, asof::Date; locf::Bool=false) = error("EmptyMarketDataProvider has no series")
get_serie_currency(::EmptyMarketDataProvider, serie_id) = error("EmptyMarketDataProvider has no series")
has_serie(::EmptyMarketDataProvider, serie_id) = false

end


text:
include("keyboard.jl")

ansi_moveup(n::Int) = string("\e[", n, "A")
const ansi_movecol1 = "\e[1G"

"""
    play(framestack::AbstractVector{T}; kwargs...) where {T<:AbstractArray}
    play(arr::T, dim=3; kwargs...)

(Experimental) Play a video of a framestack of image arrays, or 3D array along dimension `dim`.

!!! compat "ImageShow 0.3"
    The `play` function requires at least ImageShow 0.3.

# Control keys

- `p` or `space-bar`: pause/resume
- `b`, `←`(left arrow), or `↑`(up arrow): step backward
- `f`, `→`(right arrow), or `↓`(down arrow): step forward
- `ctrl-c` or `q`: exit

# Parameters

- `fps`: frame per second.

# Examples

```julia
using ImageCore, TestImages, ImageShow

img3d = RGB.(testimage("mri-stack"))
ImageShow.play(img3d) # or ImageShow.play(img3d, 3)

framestack = [img3d[:, :, i] for i in axes(img3d, 3)];
ImageShow.play(framestack)
```

See also [`explore`](@ref ImageShow.explore) for a similar version that pauses at the start
and the end.
"""
function play(framestack::AbstractVector{<:AbstractMatrix}; fps::Real=15)
    # NOTE: the default fps is chosen purely by experience and may be changed in the future
    _play(framestack; fps=fps, paused=false, quit_after_play=true)
end
play(img::AbstractArray{<:Colorant, 3}, dim=3; kwargs...) = play(map(i->selectdim(img, dim, i), axes(img, dim)); kwargs...)

"""
    play(f, Xs; kwargs...)
    play(f, Xs, Ys...; kwargs...)

(Experimental) A lazy version of `play([f(X) for X in Xs]; kwargs...)` that allocates memory only when needed.

!!! compat "ImageShow 0.3"
    The `play` function requires at least ImageShow 0.3.

# Parameters

- `fps::Int`: frame per second.

# Examples

Rotate the image and see how things going:

```julia
using TestImages, ImageShow, ImageTransformations
img = testimage("cameraman")
ImageShow.play(-π/4:π/16:π/4]; fps=3) do θ
    imrotate(img, θ, axes(img))
end
```

The following example is less meaningful, but it shows how multiple arguments are passed:

```julia
sizes = 16:4:64
values = range(0, stop=1, length=length(sizes))
ImageShow.play(values, sizes; fps=3) do v, x
    fill(RGB(v, v, v), ntuple(_->x, 2)...)
end
```
"""
play(f, arg1, args...; kwargs...) = play(mappedarray(f, arg1, args...); kwargs...)
# MappedArrays are not efficient here https://github.com/JuliaArrays/MappedArrays.jl/issues/46
play(frames::AbstractMappedArray; kwargs...) = play(collect(frames); kwargs...)

"""
    explore(framestack::AbstractVector{T}; kwargs...) where {T<:AbstractArray}
    explore(arr::T, dim=3; kwargs...)

(Experimental) Play a video of a framestack of image arrays, or 3D array along dimension `dim`.

!!! compat "ImageShow 0.3"
    The `play` function requires at least ImageShow 0.3.

Same as [`play`](@ref ImageShow.play), but will pause at the start and the end of the play. For the detailed
usage, please see the [`play` documentation](@ref ImageShow.play).
"""
function explore(framestack::AbstractVector{<:AbstractMatrix}; fps::Real=15)
    _play(framestack; fps=fps, paused=true, quit_after_play=false)
end
explore(img::AbstractArray{<:Colorant, 3}, dim=3; kwargs...) = explore(map(i->selectdim(img, dim, i), axes(img, dim)); kwargs...)

"""
    explore(f, Xs; kwargs...)
    explore(f, Xs, Ys...; kwargs...)

(Experimental)  A lazy version of `explore([f(X) for X in Xs]; kwargs...)` that allocates memory only when needed.
"""
explore(f, arg1, args...; kwargs...) = explore(mappedarray(f, arg1, args...); kwargs...)
# MappedArrays are not efficient here https://github.com/JuliaArrays/MappedArrays.jl/issues/46
explore(frames::AbstractMappedArray; kwargs...) = explore(collect(frames); kwargs...)

function _play(
        framestack::AbstractVector{<:AbstractMatrix};
        fps, paused, quit_after_play,
        # The following keywords are for advanced usages(e.g., test), common users are not expected
        # to use them directly.
        display_io::Union{Nothing, IO}=nothing,
        summary_io::IO=stdout,
        keyboard_io::IO=stdin
)
    nframes = length(framestack)

    # vars
    frame_idx = 1
    actual_fps = 0
    should_exit = false

    function render_frame(frame_idx, actual_fps; first_frame)
        frame = framestack[frame_idx]
        cols, rows = size(frame)

        if !first_frame
            print(summary_io, ansi_moveup(2), ansi_movecol1)
        end
        println(summary_io, "Frame: $frame_idx/$nframes FPS: $(round(actual_fps, digits=1))", " "^5)
        println(summary_io, "exit: \"q\" play/pause: \"space-bar\" seek: \"arrow keys\"")

        # When calling `display(MIME"image/png"(), img)`, VSCode/IJulia/Atom will eventually
        # create an `IOBuffer` to get the Base64+PNG encoded data, and send the encoded data to
        # the outside display pipeline, e.g., as JSON message.
        # For test purpose, we could directly show it to our manually created IO.
        isnothing(display_io) ? display(frame) : show(display_io, MIME"image/png"(), frame)
    end
    # These codes live in ImageShow and thus MIME"image/png" is always showable
    @assert showable(MIME"image/png"(), framestack[frame_idx])
    render_frame(frame_idx, actual_fps; first_frame=true)

    keytask = @async begin
        while !should_exit
            control_value = read_key(keyboard_io)

            if control_value == :CONTROL_BACKWARD
                frame_idx = max(frame_idx-1, 1)
            elseif control_value == :CONTROL_FORWARD
                frame_idx = min(frame_idx+1, nframes)
            elseif control_value == :CONTROL_PAUSE
                paused = !paused
            elseif control_value == :CONTROL_EXIT
                should_exit = true
            elseif control_value == :CONTROL_VOID
                nothing
            else
                error("Control value $control_value not recognized.")
            end
            sleep(1e-2) # 10ms should be enough for most keyboard event
            # @show control_value frame_idx paused should_exit
        end
    end
    keytask_channel = make_channel(keytask)

    try
        last_frame_idx = frame_idx

        while !should_exit && 1<= frame_idx <= nframes
            # when paused, only update the frame when last_frame_idx changes, i.e., only when
            # user hit arrow keys.
            # Otherwise the same frame will be rendered again and again and again and increases
            # the plot count in the plotpane endlessly.
            if frame_idx != last_frame_idx
                fps_value = paused ? 0 : actual_fps
                actual_fps = fixed_fps(fps) do
                    render_frame(frame_idx, actual_fps; first_frame=false)
                end
                last_frame_idx = frame_idx
            end
            if !quit_after_play && frame_idx == nframes
                # don't immediately quit the play
                paused = true
            end
            paused || (frame_idx += 1)
            # Wait for keyboard event update
            # This does not mean fps is at most 100
            paused && sleep(1e-2)
        end
    catch e
        e isa InterruptException || rethrow(e)
    finally
        # Stop the running read_key task so that REPL/stdin is not blocked.
        # If it's an IOBuffer then there's no need to do so becaused it will exit eventually
        # at `should_exit`.
        if !isa(keyboard_io, IOBuffer)
            safe_kill(keytask, keytask_channel)
        end
    end
    return nothing
end

"""
    fixed_fps(f::Function, fps::Real)

Run function f() at a fixed fps rate if possible.

Example:

The following codes render the frames at a given fps 60.

```julia
while true
    actual_fps = fixed_fps(60) do
        render_frame(...)
    end
end
```
"""
function fixed_fps(f, fps::Real)
    tim = Timer(1/fps)
    t = @elapsed f()
    wait(tim)
    close(tim)
    return 1/t
end


text:
#using Revise
using Plots, CSV, DataFrames, Serialization
include("draftpolicy.jl")
path = "/scratch/users/munro/draft-policy/"

function simulate_season(year, mult=10)
    team_id = CSV.read(string(path,"data/cleaned/teams_",
                            string(year),".csv"))
    season = CSV.read(string(path,"data/cleaned/season_",year,".csv"))
    matches = convert(Matrix, season[[:hID,:aID]])
    wins = convert(Array, season[:win])
    teams = convert(Array,team_id[:id])
    s = Season(teams,matches)
    record,_=calculate_draft_rule(s,wins,mult,16)
    stoptime=record.stop
    println(year)
    println(stoptime)
    end_wins = record.cwins[:,size(record.cwins,2)]
    end_draft = record.draft_prob[:,size(record.draft_prob,2)]
    print("got here")
    df = DataFrame(id = copy(team_id[:id]),team = copy(team_id[:team]),wins = copy(end_wins),score = copy(end_draft))
    sort!(df, (:wins), rev=(false))
    bottom_id = df[:id][1:3]
    bottom_name = df[:team][1:3]
    prob = record.draft_prob[bottom_id,:]
    wins = record.cwins[bottom_id,:]
    serialize(string(path,"object",string(year),".jb"),record)
    CSV.write(string(path,"record",string(year),".csv"),df)
    df
end

for year in [1985,1986,1987,1988,1989]
    simulate_season(year)
end

max_stars_repo_name: evanmunro/draft-policy

max_issues_repo_name: evanmunro/draft-policy

max_forks_repo_name: evanmunro/draft-policy


text:
function registry_url_list()::Vector{String}
    registry_url_list_raw::Vector{String} = String[
        "https://github.com/JuliaRegistries/General.git",
        "https://github.com/bcbi/BCBIRegistry.git",
        ]
    registry_url_list::Vector{String} = strip.(registry_url_list_raw)
    return registry_url_list
end


text:
abstract type RPRMaterial end

struct Material{Typ} <: RPRMaterial
    node::MaterialNode
    parameters::Dict{Symbol,Any}
    matsys::MaterialSystem
end

function set!(shape::Shape, material::Material)
    return set!(shape, material.node)
end

function set!(m::MaterialNode, input::RPR.rpr_material_node_input, m2::Material)
    return set!(m, input, m2.node)
end

const ALL_MATERIALS = [:Diffuse => RPR.RPR_MATERIAL_NODE_DIFFUSE,
                       :Microfacet => RPR.RPR_MATERIAL_NODE_MICROFACET,
                       :Reflection => RPR.RPR_MATERIAL_NODE_REFLECTION,
                       :Refraction => RPR.RPR_MATERIAL_NODE_REFRACTION,
                       :MicrofacetRefraction => RPR.RPR_MATERIAL_NODE_MICROFACET_REFRACTION,
                       :Transparent => RPR.RPR_MATERIAL_NODE_TRANSPARENT,
                       :Emissive => RPR.RPR_MATERIAL_NODE_EMISSIVE, :Ward => RPR.RPR_MATERIAL_NODE_WARD,
                       :Add => RPR.RPR_MATERIAL_NODE_ADD, :Blend => RPR.RPR_MATERIAL_NODE_BLEND,
                       :Arithmetic => RPR.RPR_MATERIAL_NODE_ARITHMETIC,
                       :Fresnel => RPR.RPR_MATERIAL_NODE_FRESNEL,
                       :NormalMap => RPR.RPR_MATERIAL_NODE_NORMAL_MAP,
                       :ImageTexture => RPR.RPR_MATERIAL_NODE_IMAGE_TEXTURE,
                       :Noise2dTexture => RPR.RPR_MATERIAL_NODE_NOISE2D_TEXTURE,
                       :DotTexture => RPR.RPR_MATERIAL_NODE_DOT_TEXTURE,
                       :GradientTexture => RPR.RPR_MATERIAL_NODE_GRADIENT_TEXTURE,
                       :CheckerTexture => RPR.RPR_MATERIAL_NODE_CHECKER_TEXTURE,
                       :ConstantTexture => RPR.RPR_MATERIAL_NODE_CONSTANT_TEXTURE,
                       :InputLookup => RPR.RPR_MATERIAL_NODE_INPUT_LOOKUP,
                       :BlendValue => RPR.RPR_MATERIAL_NODE_BLEND_VALUE,
                       :Passthrough => RPR.RPR_MATERIAL_NODE_PASSTHROUGH,
                       :Orennayar => RPR.RPR_MATERIAL_NODE_ORENNAYAR,
                       :FresnelSchlick => RPR.RPR_MATERIAL_NODE_FRESNEL_SCHLICK,
                       :DiffuseRefraction => RPR.RPR_MATERIAL_NODE_DIFFUSE_REFRACTION,
                       :BumpMap => RPR.RPR_MATERIAL_NODE_BUMP_MAP, :Volume => RPR.RPR_MATERIAL_NODE_VOLUME,
                       :MicrofacetAnisotropicReflection => RPR.RPR_MATERIAL_NODE_MICROFACET_ANISOTROPIC_REFLECTION,
                       :MicrofacetAnisotropicRefraction => RPR.RPR_MATERIAL_NODE_MICROFACET_ANISOTROPIC_REFRACTION,
                       :Twosided => RPR.RPR_MATERIAL_NODE_TWOSIDED,
                       :UvProcedural => RPR.RPR_MATERIAL_NODE_UV_PROCEDURAL,
                       :MicrofacetBeckmann => RPR.RPR_MATERIAL_NODE_MICROFACET_BECKMANN,
                       :Phong => RPR.RPR_MATERIAL_NODE_PHONG,
                       :BufferSampler => RPR.RPR_MATERIAL_NODE_BUFFER_SAMPLER,
                       :UvTriplanar => RPR.RPR_MATERIAL_NODE_UV_TRIPLANAR,
                       :AoMap => RPR.RPR_MATERIAL_NODE_AO_MAP,
                       :UserTexture_0 => RPR.RPR_MATERIAL_NODE_USER_TEXTURE_0,
                       :UserTexture_1 => RPR.RPR_MATERIAL_NODE_USER_TEXTURE_1,
                       :UserTexture_2 => RPR.RPR_MATERIAL_NODE_USER_TEXTURE_2,
                       :UserTexture_3 => RPR.RPR_MATERIAL_NODE_USER_TEXTURE_3,
                       :Uber => RPR.RPR_MATERIAL_NODE_UBERV2, :Transform => RPR.RPR_MATERIAL_NODE_TRANSFORM,
                       :RgbToHsv => RPR.RPR_MATERIAL_NODE_RGB_TO_HSV,
                       :HsvToRgb => RPR.RPR_MATERIAL_NODE_HSV_TO_RGB,
                       :UserTexture => RPR.RPR_MATERIAL_NODE_USER_TEXTURE,
                       :ToonClosure => RPR.RPR_MATERIAL_NODE_TOON_CLOSURE,
                       :ToonRamp => RPR.RPR_MATERIAL_NODE_TOON_RAMP,
                       :VoronoiTexture => RPR.RPR_MATERIAL_NODE_VORONOI_TEXTURE,
                       :GridSampler => RPR.RPR_MATERIAL_NODE_GRID_SAMPLER,
                       :Blackbody => RPR.RPR_MATERIAL_NODE_BLACKBODY,
                       :MatxDiffuseBrdf => RPR.RPR_MATERIAL_NODE_MATX_DIFFUSE_BRDF,
                       :MatxDielectricBrdf => RPR.RPR_MATERIAL_NODE_MATX_DIELECTRIC_BRDF,
                       :MatxGeneralizedSchlickBrdf => RPR.RPR_MATERIAL_NODE_MATX_GENERALIZED_SCHLICK_BRDF,
                       :MatxNoise3d => RPR.RPR_MATERIAL_NODE_MATX_NOISE3D,
                       :MatxTangent => RPR.RPR_MATERIAL_NODE_MATX_TANGENT,
                       :MatxNormal => RPR.RPR_MATERIAL_NODE_MATX_NORMAL,
                       :MatxPosition => RPR.RPR_MATERIAL_NODE_MATX_POSITION,
                       :MatxRoughnessAnisotropy => RPR.RPR_MATERIAL_NODE_MATX_ROUGHNESS_ANISOTROPY,
                       :MatxRotate3d => RPR.RPR_MATERIAL_NODE_MATX_ROTATE3D,
                       :MatxNormalize => RPR.RPR_MATERIAL_NODE_MATX_NORMALIZE,
                       :MatxIfgreater => RPR.RPR_MATERIAL_NODE_MATX_IFGREATER,
                       :MatxSheenBrdf => RPR.RPR_MATERIAL_NODE_MATX_SHEEN_BRDF,
                       :MatxDiffuseBtdf => RPR.RPR_MATERIAL_NODE_MATX_DIFFUSE_BTDF,
                       :MatxConvert => RPR.RPR_MATERIAL_NODE_MATX_CONVERT,
                       :MatxSubsurfaceBrdf => RPR.RPR_MATERIAL_NODE_MATX_SUBSURFACE_BRDF,
                       :MatxDielectricBtdf => RPR.RPR_MATERIAL_NODE_MATX_DIELECTRIC_BTDF,
                       :MatxConductorBrdf => RPR.RPR_MATERIAL_NODE_MATX_CONDUCTOR_BRDF,
                       :MatxFresnel => RPR.RPR_MATERIAL_NODE_MATX_FRESNEL,
                       :MatxLuminance => RPR.RPR_MATERIAL_NODE_MATX_LUMINANCE,
                       :MatxFractal3d => RPR.RPR_MATERIAL_NODE_MATX_FRACTAL3D,
                       :MatxMix => RPR.RPR_MATERIAL_NODE_MATX_MIX, :Matx => RPR.RPR_MATERIAL_NODE_MATX]

const ALL_MATERIALS_DICT = Dict(ALL_MATERIALS)

function Material(matsys::MaterialSystem, @nospecialize(attributes))
    type = if haskey(attributes, :type)
        attributes[:type]
    else
        :Uber # default, since it supports most attributes & is the only one supported by all backends
    end
    material_enum = get(ALL_MATERIALS_DICT, type) do
        error("Material type: $(type) not supported. Choose from RadeonProRender.ALL_MATERIALS")
    end
    delete!(attributes, :type)
    return Material{material_enum}(matsys; attributes...)
end

for (name, enum) in ALL_MATERIALS
    @eval const $(Symbol(name, :Material)) = Material{$(enum)}
end

function snake2camel(str)
    str = replace(lowercase(str), r"_(.)" => r -> uppercase(r[2]))
    return uppercase(str[1]) * str[2:end]
end

function Base.show(io::IO, material::Material{Typ}) where {Typ}
    typ_str = snake2camel(replace(string(Typ), "RPR_MATERIAL_NODE_" => ""))
    println(io, typ_str, "Material:")
    params = getfield(material, :parameters)
    maxlen = 1
    if !isempty(keys(params))
        maxlen = maximum(x -> length(string(x)), keys(params)) + 1
    end
    fstr = Printf.Format("    %-$(maxlen)s %s\n")
    for (k, v) in params
        Printf.format(io, fstr, string(k, ":"), v)
    end
end

function (::Type{Material{Typ}})(matsys::MaterialSystem; kw...) where {Typ}
    mat = Material{Typ}(MaterialNode(matsys, Typ), Dict{Symbol,Any}(kw), matsys)
    for (key, value) in kw
        setproperty!(mat, key, value)
    end
    return mat
end

function Base.setproperty!(material::T, field::Symbol, value::Vec4) where {T<:Material}
    set!(material.node, field2enum(T, field), value...)
    getfield(material, :parameters)[field] = value
    return
end

function Base.setproperty!(material::T, field::Symbol, value::Vec3) where {T<:Material}
    setproperty!(material, field, Vec4f(value..., 0.0))
    return
end

convert_to_type(matsys::MaterialSystem, ::Type{<:Number}, value) = Vec4f(value)
convert_to_type(matsys::MaterialSystem, ::Type{<:Number}, value::Material) = value
convert_to_type(matsys::MaterialSystem, ::Type{<:Colorant}, value) = to_color(value)
convert_to_type(matsys::MaterialSystem, ::Type{<:Colorant}, value::Material) = value

function convert_to_type(matsys::MaterialSystem, ::Type{<: Number}, value::AbstractMatrix)
    convert_to_type(matsys, Material, Float32.(value))
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Number}, value::AbstractMatrix{<: Colorant})
    convert_to_type(matsys, RGB, value)
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Colorant}, value::AbstractMatrix{<:Number})
    convert_to_type(matsys, Float32, value)
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Colorant}, value::AbstractMatrix{<:Color3})
    convert_to_type(matsys, Material, convert(Matrix{RGB{Float32}}, value))
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Colorant}, value::AbstractMatrix{<:TransparentColor})
    convert_to_type(matsys, Material, convert(Matrix{RGBA{Float32}}, value))
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Material}, value::AbstractMatrix)
    return Texture(matsys, value)
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Image}, value::Image)
    return value
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Image}, value::AbstractMatrix)
    return Texture(matsys, value)
end

function convert_to_type(matsys::MaterialSystem, ::Type{<: Union{Colorant, Image}}, value::Material)
    return value
end

function convert_to_type(matsys::MaterialSystem, typ, value)
    return value
end

function Base.setproperty!(material::T, field::Symbol, value) where {T<:Material}
    info = material_field_info()
    field_enum = field2enum(T, field)
    field_type = info[field_enum]
    if field_type isa Tuple
        type, range = field_type
        field_type = type # ignore range for now
    end
    matsys = material.matsys
    set!(material.node, field_enum, convert_to_type(matsys, field_type, value))
    getfield(material, :parameters)[field] = value
    return
end

function Base.setproperty!(material::T, field::Symbol, c::Color3) where {T<:Material}
    set!(material.node, field2enum(T, field), red(c), green(c), blue(c), 0.0)
    getfield(material, :parameters)[field] = c
    return
end

function Base.setproperty!(material::T, field::Symbol, value::Nothing) where {T<:Material}
    # TODO, can you actually unset a material property?
end

function field2enum(::Type{T}, field::Symbol) where {T}
    info = material_info(T)
    if haskey(info, field)
        return getfield(info, field)
    else
        error("Material $(T) doesn't have the property $(field)")
    end
end

function Base.propertynames(m::T) where {T<:Material}
    info = material_info(T)
    return propertynames(info)
end

function material_info(::Type{UberMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_COLOR,
            diffuse_weight=RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_WEIGHT,
            diffuse_roughness=RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_ROUGHNESS,
            diffuse_normal=RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_NORMAL,
            reflection_color=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_COLOR,
            reflection_weight=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_WEIGHT,
            reflection_roughness=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ROUGHNESS,
            reflection_anisotropy=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ANISOTROPY,
            reflection_anisotropy_rotation=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ANISOTROPY_ROTATION,
            reflection_mode=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_MODE,
            reflection_ior=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_IOR,
            reflection_metalness=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_METALNESS,
            reflection_normal=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_NORMAL,
            reflection_dielectric_reflectance=RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_DIELECTRIC_REFLECTANCE,
            refraction_color=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_COLOR,
            refraction_weight=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_WEIGHT,
            refraction_roughness=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ROUGHNESS,
            refraction_ior=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_IOR,
            refraction_normal=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_NORMAL,
            refraction_thin_surface=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_THIN_SURFACE,
            refraction_absorption_color=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_COLOR,
            refraction_absorption_distance=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_DISTANCE,
            refraction_caustics=RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_CAUSTICS,
            coating_color=RPR.RPR_MATERIAL_INPUT_UBER_COATING_COLOR,
            coating_weight=RPR.RPR_MATERIAL_INPUT_UBER_COATING_WEIGHT,
            coating_roughness=RPR.RPR_MATERIAL_INPUT_UBER_COATING_ROUGHNESS,
            coating_mode=RPR.RPR_MATERIAL_INPUT_UBER_COATING_MODE,
            coating_ior=RPR.RPR_MATERIAL_INPUT_UBER_COATING_IOR,
            coating_metalness=RPR.RPR_MATERIAL_INPUT_UBER_COATING_METALNESS,
            coating_normal=RPR.RPR_MATERIAL_INPUT_UBER_COATING_NORMAL,
            coating_transmission_color=RPR.RPR_MATERIAL_INPUT_UBER_COATING_TRANSMISSION_COLOR,
            coating_thickness=RPR.RPR_MATERIAL_INPUT_UBER_COATING_THICKNESS,
            sheen=RPR.RPR_MATERIAL_INPUT_UBER_SHEEN, sheen_tint=RPR.RPR_MATERIAL_INPUT_UBER_SHEEN_TINT,
            sheen_weight=RPR.RPR_MATERIAL_INPUT_UBER_SHEEN_WEIGHT,
            emission_color=RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_COLOR,
            emission_weight=RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_WEIGHT,
            emission_mode=RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_MODE,
            transparency=RPR.RPR_MATERIAL_INPUT_UBER_TRANSPARENCY,
            sss_scatter_color=RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_COLOR,
            sss_scatter_distance=RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DISTANCE,
            sss_scatter_direction=RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DIRECTION,
            sss_weight=RPR.RPR_MATERIAL_INPUT_UBER_SSS_WEIGHT,
            sss_multiscatter=RPR.RPR_MATERIAL_INPUT_UBER_SSS_MULTISCATTER,
            backscatter_weight=RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_WEIGHT,
            backscatter_color=RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_COLOR,
            fresnel_schlick_approximation=RPR.RPR_MATERIAL_INPUT_UBER_FRESNEL_SCHLICK_APPROXIMATION)
end

function material_info(::Type{DiffuseMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL,
            roughness=RPR.RPR_MATERIAL_INPUT_ROUGHNESS)
end

function material_info(::Type{DiffuseRefractionMaterial})
    return material_info(DiffuseMaterial)
end

function material_info(::Type{BumpMapMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, scale=RPR.RPR_MATERIAL_INPUT_SCALE)
end

function material_info(::Type{VolumeMaterial})
    return (scattering=RPR.RPR_MATERIAL_INPUT_SCATTERING, absorption=RPR.RPR_MATERIAL_INPUT_ABSORBTION,
            emission=RPR.RPR_MATERIAL_INPUT_EMISSION, scatter_direction=RPR.RPR_MATERIAL_INPUT_G,
            multiscatter=RPR.RPR_MATERIAL_INPUT_MULTISCATTER)
end

function material_info(::Type{MicrofacetAnisotropicReflectionMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL,
            ior=RPR.RPR_MATERIAL_INPUT_IOR, roughness=RPR.RPR_MATERIAL_INPUT_ROUGHNESS,
            anisotropic=RPR.RPR_MATERIAL_INPUT_ANISOTROPIC, rotation=RPR.RPR_MATERIAL_INPUT_ROTATION)
end

function material_info(::Type{PhongMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL,
            ior=RPR.RPR_MATERIAL_INPUT_IOR, roughness=RPR.RPR_MATERIAL_INPUT_ROUGHNESS)
end

function material_info(::Type{BlendMaterial})
    return (color1=RPR.RPR_MATERIAL_INPUT_COLOR0, color2=RPR.RPR_MATERIAL_INPUT_COLOR1,
            weight=RPR.RPR_MATERIAL_INPUT_WEIGHT,
            transmission_color=RPR.RPR_MATERIAL_INPUT_TRANSMISSION_COLOR,
            thickness=RPR.RPR_MATERIAL_INPUT_THICKNESS)
end

function material_info(::Type{MicrofacetAnisotropicRefractionMaterial})
    return material_info(MicrofacetAnisotropicReflectionMaterial)
end

function material_info(::Type{MicrofacetMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL,
            ior=RPR.RPR_MATERIAL_INPUT_IOR, roughness=RPR.RPR_MATERIAL_INPUT_ROUGHNESS)
end

function material_info(::Type{TwosidedMaterial})
    return (frontface=RPR.RPR_MATERIAL_INPUT_FRONTFACE, backface=RPR.RPR_MATERIAL_INPUT_BACKFACE)
end

function material_info(::Type{EmissiveMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR,)
end

function material_info(::Type{RefractionMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL,
            ior=RPR.RPR_MATERIAL_INPUT_IOR, caustics=RPR.RPR_MATERIAL_INPUT_CAUSTICS)
end

function material_info(::Type{ReflectionMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR, normal=RPR.RPR_MATERIAL_INPUT_NORMAL)
end

function material_info(::Type{TransparentMaterial})
    return (color=RPR.RPR_MATERIAL_INPUT_COLOR,)
end

function material_info(::Type{ImageTextureMaterial})
    return (data=RPR.RPR_MATERIAL_INPUT_DATA, uv=RPR.RPR_MATERIAL_INPUT_UV)
end

function material_info(::Type{InputLookupMaterial})
    return (value=RPR.RPR_MATERIAL_INPUT_VALUE,)
end

function material_info(::Type{ArithmeticMaterial})
    return (color1=RPR.RPR_MATERIAL_INPUT_COLOR0, color2=RPR.RPR_MATERIAL_INPUT_COLOR1,
            color3=RPR.RPR_MATERIAL_INPUT_COLOR2, color4=RPR.RPR_MATERIAL_INPUT_COLOR3,
            op=RPR.RPR_MATERIAL_INPUT_OP)
end

function material_field_info()
    f01 = (Float32, (0, 1))
    return Dict(RPR.RPR_MATERIAL_INPUT_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_COLOR0 => Material,
                RPR.RPR_MATERIAL_INPUT_COLOR1 => Material,
                RPR.RPR_MATERIAL_INPUT_NORMAL => Vec3,
                RPR.RPR_MATERIAL_INPUT_UV => Vec2,
                RPR.RPR_MATERIAL_INPUT_DATA => Image,
                RPR.RPR_MATERIAL_INPUT_ROUGHNESS => f01,
                RPR.RPR_MATERIAL_INPUT_IOR => (Float32, (1, 5)),
                RPR.RPR_MATERIAL_INPUT_ROUGHNESS_X => f01,
                RPR.RPR_MATERIAL_INPUT_ROUGHNESS_Y => f01,
                RPR.RPR_MATERIAL_INPUT_ROTATION => f01,
                RPR.RPR_MATERIAL_INPUT_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_OP => RPR.rpr_material_node_arithmetic_operation,
                # RPR.RPR_MATERIAL_INPUT_INVEC => ,
                # RPR.RPR_MATERIAL_INPUT_UV_SCALE => ,
                RPR.RPR_MATERIAL_INPUT_VALUE => RPR.rpr_material_node_lookup_value,
                RPR.RPR_MATERIAL_INPUT_REFLECTANCE => f01,
                RPR.RPR_MATERIAL_INPUT_SCALE => f01,
                RPR.RPR_MATERIAL_INPUT_SCATTERING => RGBA,
                RPR.RPR_MATERIAL_INPUT_ABSORBTION => RGBA,
                RPR.RPR_MATERIAL_INPUT_EMISSION => RGBA,
                 RPR.RPR_MATERIAL_INPUT_G => (Float32, (-1, 1)), #	Forward or back scattering.
                RPR.RPR_MATERIAL_INPUT_MULTISCATTER => Bool,
                RPR.RPR_MATERIAL_INPUT_COLOR2 => RGBA,
                RPR.RPR_MATERIAL_INPUT_COLOR3 => RGBA,
                RPR.RPR_MATERIAL_INPUT_ANISOTROPIC => (Float32, (-1, 1)), # amount forwards/backward
                RPR.RPR_MATERIAL_INPUT_FRONTFACE => Material,
                RPR.RPR_MATERIAL_INPUT_BACKFACE => Material,
                # RPR.RPR_MATERIAL_INPUT_ORIGIN => ,
                # RPR.RPR_MATERIAL_INPUT_ZAXIS => ,
                # RPR.RPR_MATERIAL_INPUT_XAXIS => ,
                # RPR.RPR_MATERIAL_INPUT_THRESHOLD => ,
                # RPR.RPR_MATERIAL_INPUT_OFFSET => ,
                # RPR.RPR_MATERIAL_INPUT_UV_TYPE => ,
                # RPR.RPR_MATERIAL_INPUT_RADIUS => ,
                # RPR.RPR_MATERIAL_INPUT_SIDE => ,
                RPR.RPR_MATERIAL_INPUT_CAUSTICS => Bool,
                RPR.RPR_MATERIAL_INPUT_TRANSMISSION_COLOR => RGBA,
                RPR.RPR_MATERIAL_INPUT_THICKNESS => f01,
                # RPR.RPR_MATERIAL_INPUT_0 => ,
                # RPR.RPR_MATERIAL_INPUT_1 => ,
                # RPR.RPR_MATERIAL_INPUT_2 => ,
                # RPR.RPR_MATERIAL_INPUT_3 => ,
                # RPR.RPR_MATERIAL_INPUT_4 => ,
                RPR.RPR_MATERIAL_INPUT_SCHLICK_APPROXIMATION => Float32,
                # RPR.RPR_MATERIAL_INPUT_APPLYSURFACE => ,
                # RPR.RPR_MATERIAL_INPUT_TANGENT => ,
                # RPR.RPR_MATERIAL_INPUT_DISTRIBUTION => ,
                # RPR.RPR_MATERIAL_INPUT_BASE => ,
                # RPR.RPR_MATERIAL_INPUT_TINT => ,
                # RPR.RPR_MATERIAL_INPUT_EXPONENT => ,
                # RPR.RPR_MATERIAL_INPUT_AMPLITUDE => ,
                # RPR.RPR_MATERIAL_INPUT_PIVOT => ,
                # RPR.RPR_MATERIAL_INPUT_POSITION => ,
                # RPR.RPR_MATERIAL_INPUT_AMOUNT => ,
                # RPR.RPR_MATERIAL_INPUT_AXIS => ,
                # RPR.RPR_MATERIAL_INPUT_LUMACOEFF => ,
                # RPR.RPR_MATERIAL_INPUT_REFLECTIVITY => ,
                # RPR.RPR_MATERIAL_INPUT_EDGE_COLOR => ,
                # RPR.RPR_MATERIAL_INPUT_VIEW_DIRECTION => ,
                # RPR.RPR_MATERIAL_INPUT_INTERIOR => ,
                # RPR.RPR_MATERIAL_INPUT_OCTAVES => ,
                # RPR.RPR_MATERIAL_INPUT_LACUNARITY => ,
                # RPR.RPR_MATERIAL_INPUT_DIMINISH => ,
                # RPR.RPR_MATERIAL_INPUT_WRAP_U => ,
                # RPR.RPR_MATERIAL_INPUT_WRAP_V => ,
                # RPR.RPR_MATERIAL_INPUT_WRAP_W => ,
                # RPR.RPR_MATERIAL_INPUT_5 => ,
                # RPR.RPR_MATERIAL_INPUT_6 => ,
                # RPR.RPR_MATERIAL_INPUT_7 => ,
                # RPR.RPR_MATERIAL_INPUT_8 => ,
                # RPR.RPR_MATERIAL_INPUT_9 => ,
                # RPR.RPR_MATERIAL_INPUT_10 => ,
                # RPR.RPR_MATERIAL_INPUT_11 => ,
                # RPR.RPR_MATERIAL_INPUT_12 => ,
                # RPR.RPR_MATERIAL_INPUT_13 => ,
                # RPR.RPR_MATERIAL_INPUT_14 => ,
                # RPR.RPR_MATERIAL_INPUT_15 => ,
                # RPR.RPR_MATERIAL_INPUT_DIFFUSE_RAMP => ,
                # RPR.RPR_MATERIAL_INPUT_SHADOW => ,
                # RPR.RPR_MATERIAL_INPUT_MID => ,
                # RPR.RPR_MATERIAL_INPUT_HIGHLIGHT => ,
                # RPR.RPR_MATERIAL_INPUT_POSITION1 => ,
                # RPR.RPR_MATERIAL_INPUT_POSITION2 => ,
                # RPR.RPR_MATERIAL_INPUT_RANGE1 => ,
                # RPR.RPR_MATERIAL_INPUT_RANGE2 => ,
                # RPR.RPR_MATERIAL_INPUT_INTERPOLATION => ,
                # RPR.RPR_MATERIAL_INPUT_RANDOMNESS => ,
                # RPR.RPR_MATERIAL_INPUT_DIMENSION => ,
                # RPR.RPR_MATERIAL_INPUT_OUTTYPE => ,
                # RPR.RPR_MATERIAL_INPUT_DENSITY => ,
                # RPR.RPR_MATERIAL_INPUT_DENSITYGRID => ,
                # RPR.RPR_MATERIAL_INPUT_DISPLACEMENT => ,
                # RPR.RPR_MATERIAL_INPUT_TEMPERATURE => ,
                # RPR.RPR_MATERIAL_INPUT_KELVIN => ,
                RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_ROUGHNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_DIFFUSE_NORMAL => Vec3f,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ROUGHNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ANISOTROPY => (Float32, (-1, 1)),
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_ANISOTROPY_ROTATION => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_MODE => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_IOR => (Float32, (1, 5)),
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_METALNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_NORMAL => Vec3f,
                RPR.RPR_MATERIAL_INPUT_UBER_REFLECTION_DIELECTRIC_REFLECTANCE => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ROUGHNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_IOR => (Float32, (1, 5)),
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_NORMAL => Vec3f,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_THIN_SURFACE => Bool,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_ABSORPTION_DISTANCE => Float32, # soft 0-10
                RPR.RPR_MATERIAL_INPUT_UBER_REFRACTION_CAUSTICS => Bool,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_ROUGHNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_MODE => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_IOR => (Float32, (1, 5)),
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_METALNESS => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_NORMAL => Vec3f,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_TRANSMISSION_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_COATING_THICKNESS => Float32, # soft 0-10
                RPR.RPR_MATERIAL_INPUT_UBER_SHEEN => RGB, RPR.RPR_MATERIAL_INPUT_UBER_SHEEN_TINT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_SHEEN_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_EMISSION_MODE => Bool,
                RPR.RPR_MATERIAL_INPUT_UBER_TRANSPARENCY => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_COLOR => RGB,
                RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DISTANCE => Vec3f, # soft 0-10
                RPR.RPR_MATERIAL_INPUT_UBER_SSS_SCATTER_DIRECTION => (Float32, (-1, 1)),
                RPR.RPR_MATERIAL_INPUT_UBER_SSS_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_SSS_MULTISCATTER => Bool,
                RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_WEIGHT => f01,
                RPR.RPR_MATERIAL_INPUT_UBER_BACKSCATTER_COLOR => RGB
                # RPR.RPR_MATERIAL_INPUT_ADDRESS => ,
                # RPR.RPR_MATERIAL_INPUT_TYPE => ,
                # RPR.RPR_MATERIAL_INPUT_UBER_FRESNEL_SCHLICK_APPROXIMATION => ,
                )
end

to_color(c::Colorant) = convert(RGBA{Float32}, c)
to_color(c::Union{Symbol,String}) = parse(RGBA{Float32}, string(c))

function LayerMaterial(a, b; weight=Vec3f(0.5), transmission_color=nothing, thickness=Vec3f(0.1))
    return BlendMaterial(a.matsys; color1=a.node, color2=b.node, weight=weight,
                         transmission_color=transmission_color, thickness=thickness)
end

function ShinyMetal(matsys; kw...)
    return MicrofacetMaterial(matsys; color=to_color(:white), roughness=Vec3f(0.0), ior=Vec3f(1.5), # Aluminium
                              kw...)
end

function Chrome(matsys; kw...)
    return MicrofacetMaterial(matsys; color=to_color(:gray), roughness=Vec3f(0.01), ior=Vec3f(1.390), # Aluminium
                              kw...)
end

function Plastic(matsys; kw...)
    return UberMaterial(matsys; color=to_color(:yellow), reflection_color=Vec4f(1),
                        reflection_weight=Vec4f(1), reflection_roughness=Vec4f(0),
                        reflection_anisotropy=Vec4f(0), reflection_anisotropy_rotation=Vec4f(0),
                        reflection_metalness=Vec4f(0), reflection_ior=Vec4f(1.5), refraction_weight=Vec4f(0),
                        coating_weight=Vec4f(0), sheen_weight=Vec4f(0), emission_weight=Vec3f(0),
                        transparency=Vec4f(0), reflection_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR),
                        emission_mode=UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_SINGLESIDED),
                        coating_mode=UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR), sss_multiscatter=false,
                        refraction_thin_surface=false)
end

function Glass(matsys; kw...)
    return UberMaterial(matsys; reflection_color=Vec4f(0.9), reflection_weight=Vec4f(1.0),
                        reflection_roughness=Vec4f(0.0), reflection_anisotropy=Vec4f(0.0),
                        reflection_anisotropy_rotation=Vec4f(0.0), reflection_mode=1,
                        reflection_ior=Vec4f(1.5), refraction_color=Vec4f(0.9), refraction_weight=Vec4f(1.0),
                        refraction_roughness=Vec4f(0.0), refraction_ior=Vec4f(1.5),
                        refraction_thin_surface=false, refraction_absorption_color=Vec4f(0.6, 0.8, 0.6, 0.0),
                        refraction_absorption_distance=Vec4f(150), refraction_caustics=true,
                        coating_weight=Vec4f(0), sheen_weight=Vec4f(0), emission_weight=Vec4f(0),
                        transparency=Vec4f(0), kw...)
end

function arithmetic_material(op, args...)
    idx = findfirst(x -> x isa Material, args)
    mat_arg = args[idx]
    mat = ArithmeticMaterial(mat_arg.matsys)
    for (i, arg) in enumerate(args)
        setproperty!(mat, Symbol("color$i"), arg)
    end
    mat.op = op
    return mat
end

function Base.:(*)(mat::Material, arg)
    return arithmetic_material(RPR.RPR_MATERIAL_NODE_OP_MUL, mat, arg)
end

function Texture(matsys::MaterialSystem, matrix::AbstractMatrix; uv=nothing)
    context = get_context(matsys)
    img = Image(context, matrix)
    tex = ImageTextureMaterial(matsys)
    tex.data = img
    if !isnothing(uv)
        tex.uv = uv
    end
    return tex
end

"""
    Matx(matsys::MaterialSystem, path::String)
Load a Material X xml definition from a `.mtlx` file.
Only works with Northstar plugin!
"""
function Matx(matsys::MaterialSystem, path::String)
    ctx = get_context(matsys)
    if ctx.plugin != Northstar
        error("Matx needs Northstar plugin")
    end
    mat = MatxMaterial(matsys)
    RPR.rprMaterialXSetFile(mat.node, path)
    return mat
end

function DielectricBrdfX(matsys::MaterialSystem)
    return Matx(matsys, assetpath("textures", "matx_dielectricBrdf.mtlx"))
end

function SurfaceGoldX(matsys::MaterialSystem)
    ctx = get_context(matsys)
    RPR.rprMaterialXAddDependencyMtlx(ctx, assetpath("textures", "matx_standard_surface.mtlx"))
    return Matx(matsys, assetpath("textures", "matx_standard_surface_gold.mtlx"))
end

function CarPaintX(matsys::MaterialSystem)
    ctx = get_context(matsys)
    RPR.rprMaterialXAddDependencyMtlx(ctx, assetpath("textures", "matx_standard_surface.mtlx"))
    return Matx(matsys, assetpath("textures", "matx_standard_surface_xi_carPaint.mtlx"))
end

function StandardMetalX(matsys::MaterialSystem)
    return Matx(matsys, assetpath("textures", "matx_standard_surface_metal_texture.mtlx"))
end

# begin
#     mesh!(ax, m, material=athmo)
#     athmo = RPR.UberMaterial(matsys)
#     athmo.color = Vec4f(1, 1, 1, 1)
#     athmo.diffuse_weight = Vec4f(0, 0, 0, 0)
#     athmo.diffuse_roughness = Vec4f(0)

#     athmo.reflection_ior = Vec4f(1.0)
#     athmo.refraction_color = Vec4f(0.5, 0.5, 0.7, 0)
#     athmo.refraction_weight = Vec4f(1)
#     athmo.refraction_roughness = Vec4f(0)
#     athmo.refraction_ior =Vec4f(1.0)
#     athmo.refraction_absorption_color = Vec4f(0.9, 0.3, 0.1, 1)
#     athmo.refraction_absorption_distance = Vec4f(1)
#     athmo.refraction_caustics = false

#     athmo.sss_scatter_color = Vec4f(1.4, 0.8, 0.3, 0)
#     athmo.sss_scatter_distance = Vec4f(0.3)
#     athmo.sss_scatter_direction = Vec4f(0)
#     athmo.sss_weight = Vec4f(1)
#     athmo.backscatter_weight = Vec4f(1)
#     athmo.backscatter_color = Vec4f(0.9, 0.1, 0.1, 1)

#     athmo.reflection_mode = UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR)
#     athmo.emission_mode = UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_DOUBLESIDED)
#     athmo.coating_mode = UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR)
#     athmo.sss_multiscatter = true
#     athmo.refraction_thin_surface = false
#     notify(refresh)
# end


text:
function sinh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    x < 0 && return -sinh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))

    y = exp(x) - exp(-x)
    z = DoubleFloat{T}(y.hi*0.5, y.lo*0.5)
    return z
end

function cosh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    x < 0 && return cosh(-x)
    iszero(x) && return one(x)
    !isfinite(x) && return nan(typeof(x))

    y = exp(x) + exp(-x)
    z = DoubleFloat{T}(y.hi*0.5, y.lo*0.5)
    return z
end

function tanh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    return sinh(x) / cosh(x)
end
#=
function tanh(x::DoubleFloat{T}) where {T<:IEEEFloat}
    x < 0 && return -tanh(-x)
    iszero(x) && return zero(x)
    !isfinite(x) && return nan(typeof(x))
    epos = exp(x)
    eneg = exp(-x)
    z = (epos - eneg) / (epos + eneg)
    return z
end
=#

csch(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    inv(sinh(x))

sech(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    inv(cosh(x))

coth(x::DoubleFloat{T}) where {T<:IEEEFloat} =
    inv(tanh(x))


text:
function jt05()
#
#  test polynomial suggested by Jenkins and Traub
#
    p = poly([0.1*ones(4);0.2*ones(3);0.3;0.3;0.4])
    z = [0.4  1; 0.3  2; 0.2  3; 0.1  4];
    p, PolyZeros(z)
end


text:
_left(x, _) = x
_right(_, y) = y

# TODO We should add the concept of alignment_base, i.e. an ancestor that provably has the
#   same alignment as a particular node. This can allow for extra pruning of the graph.

"""
    left(x, y[, alignment::Alignment; initial_values=nothing])

Construct a node that ticks according to `alignment` with the latest value of `x`.

It is "left", in the sense of picking the left-hand of the two arguments `x` and `y`.
"""
function left(x, y, alignment::Alignment=DEFAULT_ALIGNMENT; initial_values=nothing)
    return apply(_left, x, y, alignment; initial_values)
end

"""
    right(x, y[, alignment::Alignment; initial_values=nothing])

Construct a node that ticks according to `alignment` with the latest value of `y`.

It is "right", in the sense of picking the right-hand of the two arguments `x` and `y`.
"""
function right(x, y, alignment::Alignment=DEFAULT_ALIGNMENT; initial_values=nothing)
    return apply(_right, x, y, alignment; initial_values)
end

"""
    align(x, y)

Form a node that ticks with the values of `x` whenever `y` ticks.
"""
align(x, y) = right(y, x, LEFT)

# TODO support initial_values in coalign.
"""
    coalign(node_1, [node_2...; alignment::Alignment]) -> Node...

Given at least one node(s) `x`, or values that are convertible to nodes, align all of them.

We guarantee that all nodes that are returned will have the same alignment. The values of
each node will correspond to the values of the input nodes.

The choice of alignment is controlled by `alignment`, which defaults to [`UNION`](@ref).
"""
function coalign(x, x_rest...; alignment::Alignment=DEFAULT_ALIGNMENT)
    x = map(_ensure_node, [x, x_rest...])

    # Deal with simple case where we only have one input. There is no aligning to do.
    length(x) == 1 && return only(x)

    # Find a well-defined ordering of the inputs -- this is an optimisation designed to
    # avoid creating equivalent nodes if `coalign` is called multiple times.
    # As such we use objectid. Strictly this is a hash, and so there could be clashes. We
    # accept this, since if such a clash were to occur it would result only in sub-optimal
    # performance, and most likely in a very minor way.
    index = if isa(alignment, LeftAlignment)
        # For left alignment we must leave the first node in place.
        [1; 1 .+ sortperm(@view(x[2:end]); by=objectid)]
    else
        sortperm(x; by=objectid)
    end
    x, x_rest... = x[index]

    # Construct one node with the correct alignment. This will also have the correct values
    # for the first node to return.
    for node in x_rest
        x = left(x, node, alignment)
    end

    # For all of the remaining nodes, align them.
    new_nodes = (x, (align(node, x) for node in x_rest)...)

    # Convert nodes back to original order.
    return new_nodes[invperm(index)]
end

struct FirstKnot{T} <: NodeOp{T} end

mutable struct FirstKnotState <: NodeEvaluationState
    ticked::Bool
    FirstKnotState() = new(false)
end

create_evaluation_state(::Tuple{Node}, ::FirstKnot) = FirstKnotState()

function run_node!(
    ::FirstKnot{T},
    state::FirstKnotState,
    time_start::DateTime,
    time_end::DateTime,
    block::Block{T},
) where {T}
    # If we have already ticked, or the input is empty, we should not emit any knots.
    (state.ticked || isempty(block)) && return Block{T}()

    # We should tick, and record the fact that we have done so.
    state.ticked = true
    time = @inbounds first(block.times)
    value = @inbounds first(block.values)
    return Block(unchecked, [time], T[value])
end

"""
    first_knot(x::Node{T}) -> Node{T}

Get a node which ticks with only the first knot of `x`, and then never ticks again.
"""
function first_knot(node::Node{T}) where {T}
    # TODO There are various optimisations here
    #   * idempotent for constant nodes
    #   * ? eagerly filter block node
    return obtain_node((node,), FirstKnot{T}())
end

"""
    active_count(nodes...) -> Node{Int64}

Get a node of the number of the given `nodes` (at least one) which are active.
"""
function active_count(x, x_rest...)
    nodes = map(_ensure_node, [x, x_rest...])

    # Perform the same ordering optimisation that we use in coalign. This aims to give the
    # same node regardless of the order in which `nodes` were passed in.
    sort!(nodes; by=objectid)
    return reduce((x, y) -> +(x, y; initial_values=(0, 0)), align.(1, first_knot.(nodes)))
end


text:
import Base: +, -, ==, hash
using Plots

include("Hand.jl")

plot(hand, xlabel = "Card", ylabel = "Frequency", xlims = (1, 30))

hand = Hand([])
anim = @animate for i = 1:100
    hand = add_optimal(hand)
    plot(hand, xlabel = "Card", ylabel = "Frequency", xlims = (1, 30))
end

gif(anim, "optimal_hand.gif", fps = 120)

bad_hand = Hand(rand(1:10, 10))
good_hand = Hand([])
for _ = 1:sum(bad_hand.cards)
    good_hand = add_optimal(good_hand)
end

starting_distance = distance(good_hand, bad_hand)
dist = distance(good_hand, bad_hand)
dist_vect = Int[]
while dist != 0
    good_hand = add_optimal(good_hand)
    bad_hand = add_optimal(bad_hand)
    dist = distance(good_hand, bad_hand)
    push!(dist_vect, dist)
end

Plots.plot(dist_vect)


text:
@everywhere begin ## Define functions and variables on each processor
## Generate data for Monte Carlo based on West's (1996) IV simulation
function makedata!(y::Vector{Float64}, w::Matrix{Float64}, z::Matrix{Float64})
    # We're storing the initial errors v in the vector y to save
    # memory. The code matches the mathematics better using v as an
    # alias for y.
    v = y
    randn!(v); randn!(w); randn!(z)
    for i in 1:length(y)
        w[i,1] = z[i,1] + v[i]
        w[i,2] = z[i,2] + v[i]
        y[i] = w[i,1] + w[i,2] + v[i] # This line writes over v[i] as well.
    end
end

## oosstat! constructs the OOS test statistic corresponding to example
## 5.2 West's 1996 Econometrica paper. It also accomodates a data reordering
## a la the block bootstrap, which is what my paper proposes.
##
## βhat - An array that stores the recursive parameter estimates after the
##        function exectes. This array is overwritten.
## f -    The vector that will hold the out-of-sample statistics; this vector
##        is overwritten by the function.
## ZW_t - preallocated storage for the recursive window matrices Z[i,1:t]'*W[i,1:t];
##        this array is written over by the function (2 × 2 × k)
## ZY_t - preallocated storage for the recursive window vector Z[i,1:t]'*Y[1:t];
##        this matrix is written over by the function, (2 × k)
## l_t  - preallocated storage for the period t forecast loss, k-vector.
## y -    Data: the target varable, n-vector.
## w -    A matrix with the predictors; (n × k) each column corresponds to
##        a different forecasting model. The models will be estimated with IV
## z -    A matrix with the instruments for each model (n × k).
## bootindex - A vector of the bootstrap-generated index. Defaults to no bootstrap.

function oosstat!(βhat::Array{Float64,3}, f::Vector{Float64},
                  ZW_t::Array{Float64,3}, ZY_t::Matrix{Float64}, l_t::Vector{Float64},
                  y::Vector{Float64}, w::Matrix{Float64}, z::Matrix{Float64},
                  bootindex::Vector{Int})
    T = length(y)
    R = T - length(f)
    return oosstat_broken!(βhat, f, ZW_t, ZY_t, l_t, y, w, z,
                           bootindex[1:R], bootindex[(R+1):T])
end

function oosstat!(βhat::Array{Float64,3}, f::Vector{Float64},
                  ZW_t::Array{Float64,3}, ZY_t::Matrix{Float64}, l_t::Vector{Float64},
                  y::Vector{Float64}, w::Matrix{Float64}, z::Matrix{Float64})
    return oosstat!(βhat, f, ZW_t, ZY_t, l_t, y, w, z, [1:length(y)])
end

## oosstat_broken! constructs the OOS test statistic corresponding to
## example 5.2 West's 1996 Econometrica paper, but allowing the first
## R and the remaining P observations to be indexed independently. It
## also accomodates a data reordering a la the block bootstrap, which
## is what my paper proposes. Note that `oosstat!` calls
## `oosstat_broken!` for its functionality, but the paper recommends
## calling `oosstat_broken!` directly using independent bootstraps
## over 1:R and R+1:T.
##
## βhat - An array that stores the recursive parameter estimates after the
##        function exectes. This array is overwritten.
## f -    The vector that will hold the out-of-sample statistics; this vector
##        is overwritten by the function.
## ZW_t - preallocated storage for the recursive window matrices Z[i,1:t]'*W[i,1:t];
##        this array is written over by the function (2 × 2 × k)
## ZY_t - preallocated storage for the recursive window vector Z[i,1:t]'*Y[1:t];
##        this matrix is written over by the function, (2 × k)
## l_t  - preallocated storage for the period t forecast loss, k-vector.
## y -    Data: the target varable, n-vector.
## w -    A matrix with the predictors; (n × k) each column corresponds to
##        a different forecasting model. The models will be estimated with IV
## z -    A matrix with the instruments for each model (n × k).
## bootindexR - A vector of the bootstrap-generated index over the first R observations.
## bootindexP - A vector of the bootstrap-generated index over the last P observations.
function oosstat_broken!(βhat::Array{Float64,3}, f::Vector{Float64},
                  ZW_t::Array{Float64,3}, ZY_t::Matrix{Float64}, l_t::Vector{Float64},
                  y::Vector{Float64}, w::Matrix{Float64}, z::Matrix{Float64},
                  bootindexR::Vector{Int}, bootindexP::Vector{Int})
    R = length(bootindexR)
    _,_,k = size(ZW_t)
    # Initialize ZW_t and ZY_t
    for i in 1:k
        ZW_t[1, 1, i] = R
        ZW_t[1, 2, i] = sum(w[bootindexR,i])
        ZW_t[2, 1, i] = sum(z[bootindexR,i])
        ZW_t[2, 2, i] = sum(w[bootindexR, i] .* z[bootindexR, i])
        ZY_t[1, i] = sum(y[bootindexR])
        ZY_t[2, i] = sum(y[bootindexR] .* z[bootindexR, i])
    end
    # Update with next observation and produce forecasts
    for t = 1:(length(y) - R)
        boot_t = bootindexP[t]
        for i in 1:k
            # Construct the coefficient estimator from its previously
            # calculated components and forecast
            βhat[:,t,i] = ZW_t[:,:,i] \ ZY_t[:,i]
            l_t[i] = (y[boot_t] - (βhat[1,t,i] + βhat[2,t,i] * w[boot_t,i]))^2
            # Update estimates with the current observation.
            ZW_t[1,1,i] += 1.
            ZW_t[1,2,i] += w[boot_t,i]
            ZW_t[2,1,i] += z[boot_t,i]
            ZW_t[2,2,i] += w[boot_t,i] * z[boot_t,i]
            ZY_t[1,i] += y[boot_t]
            ZY_t[2,i] += z[boot_t,i] * y[boot_t]
        end
        f[t] = l_t[1] - l_t[2]
    end
    return mean(f)
end

## oosnaive constructs bootstrapped OOS test statistic that one
## gets by just bootstrapping the observed out-of-sample loss differences.
##
## f -  A vector holding the original out-of-sample statistics.
## bootindex - A vector of the bootstrap-generated index.

function oosnaive(f::Vector{Float64}, bootindex::Vector{Int})
    P = length(f)
    bootstat = 0.0
    for t in bootindex
        bootstat += f[t] / P
    end
    return bootstat
end

# newbootmean! calculates the centering term for our bootstrap using
# West's OOS statistic. The arguments are simliar to those in oosstat!
function newbootmean!(ZW::Array{Float64,3}, ZY::Matrix{Float64}, l::Vector{Float64},
                      y::Vector{Float64}, w::Matrix{Float64}, z::Matrix{Float64})
    _,k = size(z)
    for i in 1:k
        ZW[1,1,i] = 1.
        ZW[1,2,i] = mean(w[:,i])
        ZW[2,1,i] = mean(z[:,i])
        ZW[2,2,i] = mean(z[:,i] .* w[:,i])
        ZY[1,i] = mean(y)
        ZY[2,i] = mean(y .* z)
        coef = ZW[:,:,i] \ ZY[:,i]
        l[i] = mean((y - (coef[1] + coef[2] * w[:,i])).^2)
    end
    l[1] - l[2]
end

## runmc! runs the West (1996)-based Monte Carlo for a particular value of P, R, and α.
## - oosstat: preallocated vector that stores the oos test statistics for
##   each run. The length of this array is the number of simulations.
## - oostest: preallocated array that contains the result of each oos
##   bootstrap based test.
## - nboot: number of bootstrap replications (Integer).
## - P, R: number of oos and number of in-sample observations (both Integers).
## - α: nominal test size.

function runmc!(oosstat::Vector{Float64}, oostest::BitArray, nboot, P, R, α)
    n = P + R
    y = Array(Float64, n)
    w = Array(Float64, n, 2)
    z = Array(Float64, n, 2)
    ZW = Array(Float64, 2, 2, 2)
    ZY = Array(Float64, 2, 2)
    l = Array(Float64, 2)
    yboot = similar(y)
    wboot = similar(w)
    zboot = similar(z)
    f = Array(Float64,  P)
    βhat = Array(Float64, 2, P, 2)
    fullindex = Array(Int, n)
    initindex = Array(Int, R)
    naiveindex = Array(Int, P)
    oosboot = Array(Float64, 2, nboot)
    bootmean = Array(Float64, 2)
    for i in 1:length(oosstat)
        makedata!(y, w, z)
        oosstat[i] = oosstat!(βhat, f, ZW, ZY, l, y, w, z)
        ## Do "non-destructive" bootstraps
        for j in 1:nboot
            rand!(1:P, naiveindex)
            oosboot[1,j] = oosnaive(f, naiveindex)
        end
        bootmean[1] = mean(oosboot[1,:])
        ## Destructive bootstrap (overwrites f and βhat as well as ZW and ZY)
        for j in 1:nboot
            rand!(1:n, fullindex)
            oosboot[2,j] = oosstat!(βhat, f, ZW, ZY, l, y, w, z, fullindex)
        end
        bootmean[2] = mean(oosboot[2,:])
        for j in 1:length(bootmean)
            bootcrit = quantile(vec(oosboot[j,:]), [α/2, 1 - α/2]) - bootmean[j]
            oostest[j,i] = oosstat[i] < bootcrit[1] || oosstat[i] > bootcrit[2]
        end
    end
end

## allmcs! runs the West (1996)-based Monte Carlo for several
## different values of P and R; it basically wraps `runmc!`
## - nsim: number of simulations to run
## - nboot: number of bootstrap replications (Integer).
## - Ps, Rs: vectors with different values of P and R
## - α: nominal test size

function allmcs(nsim, nboot, RP, α)
    results = Array(Float64, size(RP)[1], 2)
    mcstat = Array(Float64, nsim)
    mctest = BitArray(2, nsim)
    for r in 1:size(RP)[1]
        runmc!(mcstat, mctest, nboot, RP[r,2], RP[r,1], α)
        results[r,:] = mean(mctest, 2)
    end
    return results
end
end

max_stars_repo_name: grayclhn-econ/oosbootstrap

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: grayclhn-econ/oosbootstrap

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: grayclhn-econ/oosbootstrap

max_forks_repo_licenses: MIT Unlicense


text:
# ---
# title: 211. Design Add and Search Words Data Structure
# id: problem211
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking, Depth-first Search, Design, Trie
# link: <https://leetcode.com/problems/design-add-and-search-words-data-structure/description/>
# hidden: true
# ---
#
# Design a data structure that supports adding new words and finding if a string
# matches any previously added string.
#
# Implement the `WordDictionary` class:
#
#   * `WordDictionary()` Initializes the object.
#   * `void addWord(word)` Adds `word` to the data structure, it can be matched later.
#   * `bool search(word)` Returns `true` if there is any string in the data structure that matches `word` or `false` otherwise. `word` may contain dots `'.'` where dots can be matched with any letter.
#
#
#
# **Example:**
#
#
#
#     **Input**
#     ["WordDictionary","addWord","addWord","addWord","search","search","search","search"]
#     [[],["bad"],["dad"],["mad"],["pad"],["bad"],[".ad"],["b.."]]
#     **Output**
#     [null,null,null,null,false,true,true,true]
#
#     **Explanation**
#     WordDictionary wordDictionary = new WordDictionary();
#     wordDictionary.addWord("bad");
#     wordDictionary.addWord("dad");
#     wordDictionary.addWord("mad");
#     wordDictionary.search("pad"); // return False
#     wordDictionary.search("bad"); // return True
#     wordDictionary.search(".ad"); // return True
#     wordDictionary.search("b.."); // return True
#
#
#
#
# **Constraints:**
#
#   * `1 <= word.length <= 500`
#   * `word` in `addWord` consists lower-case English letters.
#   * `word` in `search` consist of  `'.'` or lower-case English letters.
#   * At most `50000` calls will be made to `addWord` and `search`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TestExpressionsBroadcastForm

import Test: @testset, @test
import Volcanito: broadcast_form, ColumnName
import MacroTools: postwalk, rmlines

@testset "broadcast_form" begin
    input = broadcast_form(
        :(y + x + $(Expr(:$, :z))),
        (ColumnName(:x, false), ColumnName(:y, false)),
    )
    output = :((x, y) -> y + x + z)

    # @test postwalk(rmlines, input) == postwalk(rmlines, output)
end

end


text:
include("LowRankOperator.jl")
include("LowRankPertOperator.jl")