text:
include("/Users/ksb/Catlab.jl/src/atp/Chase.jl")

sch = Schema(Dict(:R=>2, :A=>1))
_, Isch = schema_to_csettype(sch)
I = Isch()
add_parts!(I, :V, 2, id_V=[101,202], sort=[:X, :X])
add_parts!(I, :R, 2, R_1=[1,2], R_2=[2,2], id_R=[1001,2002])
I2 = deepcopy(I)
add_parts!(I2, :V, 1, id_V=[303], sort=[:X])
add_parts!(I2, :A, 2, A_1=[2,3], id_A=[1000,2000])
I3 = deepcopy(I2)
add_parts!(I2, :R, 1, R_1=[2], R_2=[3], id_R=[3003]) # x2, y relation
add_parts!(I3, :R, 1, R_1=[1], R_2=[3], id_R=[3003]) # x1, y relation

t1 = @acset Isch begin
  V = 2
  R = 1
  R_1 = [1]
  R_2 = [2]
  sort = [:X,:X]
  id_V = [11, 22]
  id_R = [241]
end

t2 = @acset Isch begin
  V = 3
  R = 1
  A = 2
  R_1 = [1]
  R_2 = [3]
  A_1 = [3, 2]
  sort = [:X,:X,:X]
  id_V = [11,22,33]
  id_R = [799]
  id_A = [900, 901]
end

tgd3 = Dep(:d1, t1, t2)
egd3 = Dep(:d2, t1, Isch(), Set([Set([11,22])]))
te = Dep(:d3, t1, t2, Set([Set([11,22])]))

ch, _ = core_chase(I, Set([tgd3]))

redun = @acset Isch begin
  V = 2
  A = 2
  sort = [:X,:X]
  id_V = [11, 22]
  A_1 = [1,1]
  id_A = [1, 2]
end
rcore = findcore(redun)

xnyymx = trim(@program(C, (x::X,n::X,m::X),
              let nx = mul(n,x); (nx,[mul(nx,m),x]) end), 1,1)
s = [:e_intro, :sym, :mul_assoc, :cancel, :pos, :leftid]
h, newres, r, r2 = prove(T_pcs_mon, xnyymx, passx; strat=seq([[ss] for ss in s]), maxi=length(s))


text:
module BoundaryCellFieldsTests

using Test
using Gridap
using Gridap.CellArrayApply: IndexCellArrayFromKernel
import Gridap: ∇

tags = [23,24,25]
model = CartesianDiscreteModel(partition=(3,4,2))

trian = Triangulation(model)
btrian = BoundaryTriangulation(model,tags)

ufun(x) = 2*x[1] + x[2]
grad_ufun(x) = VectorValue(2.0,1.0)
∇(::typeof(ufun)) = grad_ufun

u = CellField(trian,ufun)

bu = restrict(u,btrian)
bu_grad = ∇(bu)

bquad = CellQuadrature(btrian,degree=0)
q = coordinates(bquad)

quad = CellQuadrature(trian,degree=2)
s = coordinates(quad)

bphi = CellGeomap(btrian)

x = evaluate(bphi,q)

buq = evaluate(bu,q)
@test isa(buq,IndexCellArrayFromKernel)
bu_grad_q = evaluate(bu_grad,q)

phi = CellGeomap(trian)
bphi2 = restrict(phi,btrian)

x2 = evaluate(bphi2,q)

@test x ≈ x2

#writevtk(btrian,"btrian",cellfields=["u"=>bu,"u_grad"=>bu_grad])

#writevtk(x,"x",pointdata=["u"=>buq])

cn = integrate(bu,btrian,bquad)
@test isa(cn,CellNumber)
_ = collect(cn)

v = CellBasis(trian)
bv = restrict(v,btrian)
cm = varinner(bv,bu)

ca = integrate(cm,btrian,bquad)
@test isa(ca,CellArray{Float64,1})
_ = collect(ca)

cm = varinner(bv,bv)
ca = integrate(cm,btrian,bquad)
@test isa(ca,CellArray{Float64,2})
_ = collect(ca)

end # module

ext: jl


text:
"""
AX2StageLeastSquares(y, D, instrument, control)
A simple two stage least squares implementation.
"""
struct AX2StageLeastSquares <: AXEstimator
    β::Array{Float64} # coefficient
    FS::Array{Float64} # first stage coefficients
    y::Matrix{Float64} # response
	Z::Matrix{Float64} # combined first stage variables
    X::Matrix{Float64} # combined second stage variables

	# Define constructor function
	function AX2StageLeastSquares(y::Matrix{Float64}, D::Matrix{Float64},
                                  instrument::Matrix{Float64}, control = nothing)
        # Add constant if no control is passed
        if isnothing(control) control = ones(length(y)) end

		# Define data matrices
		Z = hcat(control, instrument) # combined first stage variables
		X = hcat(D, control)          # combined second stage variables

		# Calculate TSLS coefficient
        PZ = Z * (Z' * Z)^(-1) * Z'
        β  = inv(X' * PZ * PZ * X) * (X' * PZ * y)
        FS = inv(Z' * Z) * (Z' * X)

		# Return output
		new(β, FS, y, Z, X)
	end #AX2StageLeastSquares
end #AX2StageLeastSquares

# Methods ======================================================================
"""
infer(fit::AX2StageLeastSquares; heteroskedastic, print_df)
A method to calculate standard errors of an AX2StageLeastSquares object.

"""
function infer(fit::AX2StageLeastSquares;
               heteroskedastic::Bool=false, cluster=nothing,
               print_df::Bool=true)
    # Retrieve necessary parameters
    N = length(fit.y)
    Kx = size(fit.X, 2)
    Kz = size(fit.Z, 2)

	# Calculate the variance-covariance matrix
    FShat = fit.Z * fit.FS
    FSinv = inv(FShat' * FShat)
    resid = fit.y - predict(fit)
	if heteroskedastic
        residFS = FShat .* resid
        PZuuZP  = (residFS' * residFS)
		vcv     = FSinv * PZuuZP * FSinv .* (N / (N - Kz))
	else
		vcv = FSinv .* (resid' * resid) ./ (N - Kz)
	end

	# Get standard errors, t-statistics and p-values
    se = sqrt.(vcv[diagind(vcv)])
    t_stat = fit.β ./ se
    p_val = 2 * cdf.(Normal(), -abs.(t_stat))

    # Print estimates
    if print_df
        out_df = DataFrame(hcat(fit.β, se, t_stat, p_val), :auto)
        rename!(out_df, ["coef", "se", "t-stat", "p-val"])
        display(out_df)
    end

    # Organize and return output
    output = (β = fit.β, se = se, t = t_stat, p = p_val)
    return output
end #infer.AX2StageLeastSquares


text:
###########################################################################################
#Define the marginal likelihood function
###########################################################################################
@everywhere function MargLik(x, y, g)
    N = size(y, 1)
    k = size(x, 2)
    P = eye(N) - x * inv(x'x) * x'
    yty = (y - mean(y))' * (y - mean(y))
    ymy = y'P * y
    g1 = g / (g + 1)
    g2 = 1 / (g + 1)
    ln_fy = .5 * k * log.(g1) - .5 * (N - 1) * log.(g2 * ymy + g1 * yty)
    return ln_fy[1]
end

###########################################################################################
#Use MCMCMC method for bayesian model averaging
###########################################################################################
@everywhere function bma_forc(Y, Y_lag, X_lag, vec_î, h, θ, pre_conv, post_conv)
    LHS = Y[1:vec_î - h]
    N = length(LHS)
    x = [Y_lag[1:N] X_lag[1:N, :]]
    x = x .- mean(x, 1)
    K = size(x, 2)
    τ_vec = zeros(Bool, K)
    k = sum(τ_vec)
    RHS = [ones(N) x[:, find(τ_vec)]]
    ln_mod_prior = log(θ ^ k * (1 - θ) ^ (K - k))
    # Compute marginal likelihood for initial model
    g = min(1 / N, 1 / K ^ 2)
    ln_fy = MargLik(RHS, LHS, g)
    y_forc = 0.
    # Begin MC3 simulations
    for s = 1:(pre_conv + post_conv)
        # Generate candidate Neighborhood Model
        neighbor_indx = rand(1:(1 + K))
        # Specify candidate model
        τ_vec_new = deepcopy(τ_vec)
        if neighbor_indx >= 2
            τ_vec_new[neighbor_indx - 1] = !τ_vec[neighbor_indx - 1]
        end
        k = sum(τ_vec_new)
        println("s = $(s), k = $(k): $(find(τ_vec_new))")
        RHS_new = [ones(N) x[:, find(τ_vec_new)]]
        ln_mod_prior_new = log(θ ^ k * (1 - θ) ^ (K - k))
        # Compute marginal likelihood for candidate model
        ln_fy_new = -Inf
        if k < K
            writedlm("RSM_new.csv", RHS_new, ',')
            writedlm("LHS.csv", LHS, ',')
            ln_fy_new = MargLik(RHS_new, LHS, g)
        end
        # MH Step
        prob_acc = min(exp(ln_fy_new + ln_mod_prior_new - (ln_fy + ln_mod_prior)), 1)
        if rand() <= prob_acc
            ln_fy = ln_fy_new
            τ_vec = deepcopy(τ_vec_new)
        end
        # If post-convergence, form forecasts
        if s > pre_conv
            x_selected = x[:, find(τ_vec)]
            post_β̄ = [mean(LHS); inv((1 + g) * (x_selected' * x_selected)) * (x_selected' * LHS)]
            x_forc = [Y_lag[1:vec_î] X_lag[1:vec_î, :]]
            x_forc = x_forc .- mean(x_forc, 1)
            x_forc = [ones(size(x_forc, 1)) x_forc[:, find(τ_vec)]]
            x_crnt = x_forc[end, :]
            y_forc = y_forc + x_crnt' * post_β̄
        end
    end
    return y_forc / (post_conv)
end

###########################################################################################
#The inner function
###########################################################################################
@everywhere function bma_inner_RMSE(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ, pre_conv, post_conv, oos_i)
    ## Vector Index for oos forecast goes from J_cv - (h_lag + h) to J - (1 + h_lag + h)
    #println("θ = $θ")
    ŷ_cv = zeros(Float64, size(J_cv - (h_lag + h):J + (oos_i - 1) - (1 + h_lag + h)))
    for (n, vec_î) in enumerate(J_cv - (h_lag + h):J + (oos_i - 1) - (1 + h_lag + h))
        println("n = $n")
        ŷ_cv[n] = bma_forc(Y, Y_lag, X_lag, vec_î, h, θ, pre_conv, post_conv)
    end
    return [sqrt(mean((Y[J_cv - (h_lag + h):J - (1 + h_lag + h)] - ŷ_cv) .^ 2)), θ]
end

###########################################################################################
#The outer function
###########################################################################################
function bma_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ_vec, pre_conv, post_conv, oos_i)
    # Step1 : pick the best θ
    RMSE_cv = zeros(length(θ_vec), 2)
    result = pmap(θ -> bma_inner_RMSE(J, h, h_lag, J_cv, Y, Y_lag, X_lag, θ, pre_conv, post_conv, oos_i),
                  θ_vec)::Array{Array{Float64,1},1}
    for i in 1:length(θ_vec)
        RMSE_cv[i, :] = result[i]
    end
    index = findmin(RMSE_cv[:, 1])[2]
    θ = θ_vec[index]
    # Step 2: oos forecasting
    ## Vector Index for oos forecast goes from J - (h_lag + h) to length(Y)
    vec_î = (J - (h_lag + h):length(Y))[oos_i]
    ŷ = bma_forc(Y, Y_lag, X_lag, vec_î, h, θ, pre_conv, post_conv)
    return θ, ŷ, RMSE_cv
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning


text:
__precompile__(true)

module BayesianNonparametricStatistics

export
# model.jl defines:
AbstractModel, SDEModel,
# samplepath.jl defines:
AbstractSamplePath, SamplePath, isincreasing,
# basisfunctions.jl defines:
fourier, faberschauder, faberschauderone,
# SDE.jl defines:
AbstractSDE, SDE,
# gaussianvectors.jl defines:
GaussianVector,
# gaussianprocess.jl defines:
AbstractGaussianProcess, GaussianProcess,
FaberSchauderExpansionWithGaussianCoefficients,
# defined in both gaussianvectors.jl and gaussianprocess.jl:
mean, var,
# calculateposterior.jl defines:
calculateposterior

if v"0.6" ≤ VERSION < v"0.7-"
	include(joinpath("version0dot6", "BNPversion0dot6.jl"))
elseif VERSION ≥ v"0.7-"
	include(joinpath("version0dot7and1dotX", "BNPversion0dot7and1dotX.jl"))
else
	throw(LoadError("BayesianNonparametricStatistics is only available for Julia version 0.6 and higher."))
end
end # End of module.


text:
using SymEngine

"""
Creates the necessary differentiated versions of base kernels required by GLOM.
You must pass it a SymEngine Basic object with the variables already declared with the @vars command.
The created function will look like this

    \$kernel_name(hyperparameters::Vector{<:Real}, δ::Real; dorder::Vector{<:Integer}=zeros(Int64, length(hyperparameters) + 2))

For example, you could define a kernel like so:

    "Radial basis function GP kernel (aka squared exonential, ~gaussian)"
    function se_kernel_base(λ::Number, δ::Number)
        return exp(-δ ^ 2 / (2 * λ ^ 2))
    end

And then calculate the necessary derivative versions like so:

    @vars δ λ
    kernel_coder(se_kernel_base(λ, δ), "se_kernel")

The function is saved in src/kernels/\$kernel_name.jl, so you can use it with a command akin to this:

    include("src/kernels/" * kernel_name * ".jl")

"""

function kernel_coder(
    symbolic_kernel_original::Basic,
    kernel_name::String;
    periodic_var::String="",
    cutoff_var::String="")
    # add_white_noise::Bool=false)

    δ = symbols("δ")
    periodic = periodic_var != ""
    if periodic
        P_sym_str = kernel_name * "_P"
        P_sym = symbols(P_sym_str)
        π_sym = symbols("π_sym")
        @assert occursin(periodic_var, SymEngine.toString(symbolic_kernel_original)) "can't find periodic variable"
        symbolic_kernel_original = subs(symbolic_kernel_original, abs(symbols(periodic_var))=>2*sin(π_sym*abs(δ)/P_sym))
        symbolic_kernel_original = subs(symbolic_kernel_original, symbols(periodic_var)=>2*sin(π_sym*δ/P_sym))
        # length(periodic_vars)==1 ? P_sym_str = kernel_name * "_P" : P_sym_str = kernel_name * "_P$i"
        # for i in 1:length(periodic_vars)
        #     length(periodic_vars)==1 ? P_sym_str = kernel_name * "_P" : P_sym_str = kernel_name * "_P$i"
        #     P_sym = symbols(P_sym_str)
        #     append!(symbols, [P_sym])
        #     append!(symbols_str, [P_sym_str])
        #     symbolic_kernel_original = subs(symbolic_kernel_original, periodic_vars[i]=>2*sin(π*δ/P_sym))
        # end
        # kernel_name *= "_periodic"
    end

    # if add_white_noise
    #     σ_white = symbols("σ_white")
    #     symbolic_kernel_original += σ_white * σ_white
    # end

    # get the symbols of the passed function
    symbs = free_symbols(symbolic_kernel_original)

    δ_inds = findall(x -> x==δ, symbs)
    @assert length(δ_inds) <= 1
    if length(δ_inds)==1
        deleteat!(symbs, δ_inds[1])
    end
    π_inds = findall(x -> x==symbols("π_sym"), symbs)
    if length(π_inds) > 0
        deleteat!(symbs, π_inds[1])
    end
    symbs_str = [string(symb) for symb in symbs]

    hyper_amount = length(symbs)

    # open the file we will write to
    kernel_name *= "_kernel"
    file_loc = "src/kernels/" * kernel_name * ".jl"
    io = open(file_loc, "w")

    # begin to write the function including assertions that the amount of hyperparameters are correct
    write(io, """import GPLinearODEMaker.powers_of_negative_one\n\n\"\"\"
$kernel_name function created by kernel_coder(). Requires $hyper_amount hyperparameters. Likely created using $(kernel_name)_base() as an input.
Use with include(\"src/kernels/$kernel_name.jl\").
hyperparameters == $symbs_str
\"\"\"
function $kernel_name(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==$hyper_amount+Int(include_shift) \"hyperparameters is the wrong length\"
    dorder_len = $(hyper_amount + 2) + Int(include_shift)
    @assert length(dorder)==dorder_len \"dorder is the wrong length\"
    @assert maximum(dorder) < 3 \"No more than two derivatives for each time or hyperparameter can be calculated\"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:$(hyper_amount+1) .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end\n\n""")

    # map the hyperparameters that will be passed to this function to the symbol names
    for i in 1:(hyper_amount)
        write(io, "    " * symbs_str[i] * " = hyperparameters_view[$i]\n")
    end
    # if add_white_noise
    #     write(io, "    σ_white *= Int(δ == 0)\n")  # accounting for white noise only being on the diagonal
    # end
    write(io, "\n")

    if cutoff_var!=""
        @assert cutoff_var in symbs_str
        write(io, "    if abs(δ) > $cutoff_var\n")
        write(io, "        dorder[2] = dorder2\n")
        write(io, "        return 0\n")
        write(io, "    end\n\n")
    end

    abs_vars = String[]
    for i in append!(["δ"], symbs_str)
        println("abs($i)")
        if occursin("abs($i)", SymEngine.toString(symbolic_kernel_original))
            append!(abs_vars, [i])
            println("found something")
            # write(io, "    dabsδ = sign(δ)  # store derivative of abs()\n")
            write(io, "    dabs_$i = powers_of_negative_one($i < 0)  # store derivative of abs()\n")
            write(io, "    abs_$i = abs($i)\n\n")
        end
    end

    max_δ_derivs = 7  # δf has 7 derivatives (0-(4+2), 4 from 2nd derivatives of GP and 2 from a possible time delay hyperparameter)
    max_hyper_derivs = 3  # other symbols have 3 (0-2)
    # calculate all of the necessary derivations we need for the GLOM model
    # for two symbols (δ and a hyperparameter), dorders is of the form:
    # [6 2; 4 2; 3 2; 2 2; 1 2; 0 2; ... ]
    # where the dorders[n, :]==[dorder of δ, dorder of symbol 2]
    # can be made for any number of symbols

    dorders = zeros(Int64, max_δ_derivs * (max_hyper_derivs ^ hyper_amount), hyper_amount + 1)
    amount_of_dorders = size(dorders,1)
    for i in 1:amount_of_dorders
          quant = amount_of_dorders - i
          dorders[i, 1] = rem(quant, max_δ_derivs)
          quant = div(quant, max_δ_derivs)
          for j in 1:(hyper_amount)
                dorders[i, j+1] = rem(quant, max_hyper_derivs)
                quant = div(quant, max_hyper_derivs)
          end
    end

    # for each differentiated version of the kernel
    for i in 1:amount_of_dorders

        # dorder = convert(Vector{Int64}, dorders[i, :])
        dorder = dorders[i, :]

        # only record another differentiated version of the function if we will actually use it
        # i.e. no instances where differentiations of multiple, non-time symbols are asked for
        if sum(dorder[2:end]) < max_hyper_derivs
            symbolic_kernel = copy(symbolic_kernel_original)

            # performing the differentiations
            symbolic_kernel = diff(symbolic_kernel, δ, dorder[1])
            for j in 1:hyper_amount
                # println(symbs[j], dorder[j+1])
                symbolic_kernel = diff(symbolic_kernel, symbs[j], dorder[j+1])
            end

            for j in abs_vars
                variable, abs_var, dabs_var = symbols("$j abs_$j dabs_$j")
                symbolic_kernel = subs(symbolic_kernel, diff(abs(variable), variable)=>dabs_var)
                for i in 2:6
                    symbolic_kernel = subs(symbolic_kernel, diff(abs(variable), variable, i)=>0)
                end
                for i in 1:2
                    symbolic_kernel = subs(symbolic_kernel, dabs_var^(i*2)=>1)
                end
                symbolic_kernel = subs(symbolic_kernel, abs(variable)=>abs_var)
            end

            symbolic_kernel_str = SymEngine.toString(symbolic_kernel)
            symbolic_kernel_str = replace(symbolic_kernel_str, "π_sym"=>"π")
            symbolic_kernel_str = " " * symbolic_kernel_str

            # some simplifications
            symbolic_kernel_str = replace(symbolic_kernel_str, "sqrt(δ^2)"=>"abs(δ)")
            symbolic_kernel_str = replace(symbolic_kernel_str, " 0.0 - "=>" -")
            symbolic_kernel_str = replace(symbolic_kernel_str, "(0.0 - "=>"(-")
            symbolic_kernel_str = replace(symbolic_kernel_str, "(0.0 + "=>"(")
            symbolic_kernel_str = replace(symbolic_kernel_str, " 0.0 + "=>" ")
            symbolic_kernel_str = replace(symbolic_kernel_str, " 1.0*"=>" ")
            symbolic_kernel_str = replace(symbolic_kernel_str, "-1.0*"=>"-")

            # println(symbolic_kernel_str)

            write(io, "    if dorder_view==" * string(dorder) * "\n")
            write(io, "        func =" * symbolic_kernel_str * "\n    end\n\n")
        end

    end

    write(io, "    dorder[2] = dorder2  # resetting dorder[2]\n")
    write(io, "    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives")
    write(io, "\n\nend\n\n\n")
    write(io, "return $kernel_name, $hyper_amount  # the function handle and the number of kernel hyperparameters\n")
    close(io)

    @warn "The order of hyperparameters in the function may be different from the order given when you made them in @vars. Check the created function file (or the print statement below) to see what order you need to actually use."
    println("$kernel_name() created at $file_loc")
    println("hyperparameters == ", symbs_str)
end


text:
using  Gtk

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

# Import binary definitions
const suffix = :Leaf
#=include("GLib/GLib.jl") #important =#
#= Import GLib - this is from staticfloat/GlibcBuilder imports  skippable
using .GLib
using .GLib.MutableTypes
import .GLib: set_gtk_property!, get_gtk_property, getproperty, FieldRef
import .GLib:
    signal_connect, signal_handler_disconnect,
    signal_handler_block, signal_handler_unblock, signal_handler_is_connected,
    signal_emit, unsafe_convert,
    AbstractStringLike, bytestring
=#
"""Guide

GTK3_jll: standard Gtk package Imports

#icon packages (2):
##1.hicolor_icon_theme_jll
##2. adwaita_icon_theme_jll

"""
#this is for MinGTK
using GTK3_jll, Glib_jll, Xorg_xkeyboard_config_jll, gdk_pixbuf_jll, adwaita_icon_theme_jll, hicolor_icon_theme_jll #done

using Librsvg_jll
using JLLWrappers
using Pkg.Artifacts #ok # ] add Artifacts
#const libgdk # = libgdk3
#const libgtk # = libgtk3
#const libgdk_pixbuf # = libgdkpixbuf

import Base: convert, show, run, size, resize!, length, getindex, setindex!,
             insert!, push!, append!, pushfirst!, pop!, splice!, delete!, deleteat!,
             parent, isempty, empty!, first, last, in, popfirst!,
             eltype, copy, isvalid, string, sigatomic_begin, sigatomic_end, (:), iterate #fix: adds c.jl fixes sigatomic_begin,  sigatomic_end ccalls #crucil

export showall, select!, start

using Reexport

#Import Graphics: width, height, getgc

@reexport using Graphics # lets you import the following from Graphics package:
import Graphics: width, height, getgc

using Cairo
import Cairo: destroy
using Serialization

const Index{I<:Integer} = Union{I, AbstractVector{I}}

export GAccessor

#= not found
include("basic_exports.jl")
include("long_exports.jl")
include("long_leaf_exports.jl")

=#

#=

=#
global const libgtk_version = VersionNumber(
      ccall((:gtk_get_major_version, libgtk), Cint, ()),
      ccall((:gtk_get_minor_version, libgtk), Cint, ()),
      ccall((:gtk_get_micro_version, libgtk), Cint, ()))

      """
        __init__

        ```input
        setting up Gtk required files (icons & schemas)
        sets a particular Environment variable for each OS(OS-specific):
        Windows:  ENV["XDG_DATA_DIRS"]

        ```
references to:
gdk-pixbuf-loaders-cache
gdk-pixbuf-loaders-dir

    """
function __init__()
    #for windows
    # Set XDG_DATA_DIRS so that Gtk can find its icons & schemas
    ENV["XDG_DATA_DIRS"] = join(filter(x -> x !== nothing, [
            dirname(adwaita_icons_dir),
            dirname(hicolor_icons_dir),
            joinpath(dirname(GTK3_jll.libgdk3_path::String), "..", "share"),
            get(ENV, "XDG_DATA_DIRS", nothing)::Union{String,Nothing}, # getting XDG_DATA_DIRS
        ]), Sys.iswindows() ? ";" : ":")

    # Next, ensure that gdk-pixbuf has its loaders.cache file; we generate a #review#1: what do you mean by this (whomever wrote this)
    # MutableArtifacts.toml file that maps in a loaders.cache we dynamically
    # generate by running `gdk-pixbuf-query-loaders:`
    mutable_artifacts_toml = joinpath(dirname(@__DIR__), "MutableArtifacts.toml")
    loaders_cache_name = "gdk-pixbuf-loaders-cache"
    loaders_cache_hash = artifact_hash(loaders_cache_name, mutable_artifacts_toml)
    loaders_dir_name = "gdk-pixbuf-loaders-dir"
    loaders_dir_hash = artifact_hash(loaders_dir_name, mutable_artifacts_toml)

    if loaders_cache_hash === nothing
        if Librsvg_jll.is_available()
            # Copy loaders into a directory
            loaders_dir_hash = create_artifact() do art_dir
                loaders_dir = mkdir(joinpath(art_dir, "loaders_dir"))
                pixbuf_loaders = joinpath.(gdk_pixbuf_loaders_dir, readdir(gdk_pixbuf_loaders_dir))
                push!(pixbuf_loaders, Librsvg_jll.libpixbufloader_svg)
                cp.(pixbuf_loaders, joinpath.(loaders_dir, basename.(pixbuf_loaders)))
            end

            loaders_dir = joinpath(artifact_path(loaders_dir_hash), "loaders_dir")
            # Pkg removes "execute" permissions on Windows
            Sys.iswindows() && chmod(artifact_path(loaders_dir_hash), 0o755; recursive = true)
            # Run gdk-pixbuf-query-loaders, capture output
            loader_cache_contents = gdk_pixbuf_query_loaders() do gpql
                withenv("GDK_PIXBUF_MODULEDIR" => loaders_dir, JLLWrappers.LIBPATH_env => Librsvg_jll.LIBPATH[]) do
                    return String(readchomp(`$gpql`))
                end
            end

            bind_artifact!(mutable_artifacts_toml,
                loaders_dir_name,
                loaders_dir_hash;
                force = true
            )
        else  # just use the gdk_pixbuf directory
            loader_cache_contents = gdk_pixbuf_query_loaders() do gpql
                withenv("GDK_PIXBUF_MODULEDIR" => gdk_pixbuf_loaders_dir) do
                    return String(read(`$gpql`))
                end
            end
        end
        # Write cache out to file in new artifact
        loaders_cache_hash = create_artifact() do art_dir
            open(joinpath(art_dir, "loaders.cache"), "w") do io
                write(io, loader_cache_contents)
            end
        end
        bind_artifact!(mutable_artifacts_toml,
            loaders_cache_name,
            loaders_cache_hash;
            force = true
        )
    end

    # Point gdk to our cached loaders
    ENV["GDK_PIXBUF_MODULE_FILE"] = joinpath(artifact_path(loaders_cache_hash), "loaders.cache")
    ENV["GDK_PIXBUF_MODULEDIR"] = Librsvg_jll.is_available() && loaders_dir_hash !== nothing ?
                                  joinpath(artifact_path(loaders_dir_hash), "loaders_dir") :
                                  gdk_pixbuf_loaders_dir

    if Sys.islinux() || Sys.isfreebsd()
        # Needed by xkbcommon:
        # https://xkbcommon.org/doc/current/group__include-path.html.  Related
        # to issue https://github.com/JuliaGraphics/Gtk.jl/issues/469
        ENV["XKB_CONFIG_ROOT"] = joinpath(Xorg_xkeyboard_config_jll.artifact_dir::String,
            "share", "X11", "xkb")
    end

    GError() do error_check # GError #not defined
        ccall((:gtk_init_with_args, libgtk), Bool,
            (Ptr{Nothing}, Ptr{Nothing}, Ptr{UInt8}, Ptr{Nothing}, Ptr{UInt8}, Ptr{Ptr{GError}}),
            C_NULL, C_NULL, "Julia Gtk Bindings", C_NULL, C_NULL, error_check)
    end

    # if g_main_depth > 0, a glib main-loop is already running.
    # unfortunately this call does not reliably reflect the state after the
    # loop has been stopped or restarted, so only use it once at the start
    gtk_main_running[] = ccall((:g_main_depth, GLib.libglib), Cint, ()) > 0

    # Given GLib provides `g_idle_add` to specify what happens during idle, this allows
    # that call to also start the eventloop
    GLib.gtk_eventloop_f[] = enable_eventloop

    auto_idle[] = get(ENV, "GTK_AUTO_IDLE", "true") == "true"

    # by default, defer starting the event loop until either `show`, `showall`, or `g_idle_add` is called
    enable_eventloop(!auto_idle[])
end

const auto_idle = Ref{Bool}(true) # control default via ENV["GTK_AUTO_IDLE"]
const gtk_main_running = Ref{Bool}(false)
const quit_task = Ref{Task}()
const enable_eventloop_lock = Base.ReentrantLock() #

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
##############################################################################
##
## Iteration: eachrow, eachcol
##
##############################################################################

# TODO: Reconsider/redesign eachrow -- ~100% overhead

# Iteration by rows
"""
    DataFrameRows{D<:AbstractDataFrame,S<:AbstractIndex} <: AbstractVector{DataFrameRow{D,S}}

Iterator over rows of an `AbstractDataFrame`,
with each row represented as a `DataFrameRow`.

A value of this type is returned by the [`eachrow`](@ref) function.
"""
struct DataFrameRows{D<:AbstractDataFrame,S<:AbstractIndex} <: AbstractVector{DataFrameRow{D,S}}
    df::D
    index::S
end

Base.summary(dfrs::DataFrameRows) = "$(length(dfrs))-element DataFrameRows"
Base.summary(io::IO, dfrs::DataFrameRows) = print(io, summary(dfrs))

"""
    eachrow(df::AbstractDataFrame)

Return a `DataFrameRows` that iterates a data frame row by row,
with each row represented as a `DataFrameRow`.

**Examples**

```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> eachrow(df)
4-element DataFrameRows:
 DataFrameRow (row 1)
x  1
y  11
 DataFrameRow (row 2)
x  2
y  12
 DataFrameRow (row 3)
x  3
y  13
 DataFrameRow (row 4)
x  4
y  14

julia> copy.(eachrow(df))
4-element Array{NamedTuple{(:x, :y),Tuple{Int64,Int64}},1}:
 (x = 1, y = 11)
 (x = 2, y = 12)
 (x = 3, y = 13)
 (x = 4, y = 14)

julia> eachrow(view(df, [4,3], [2,1]))
2-element DataFrameRows:
 DataFrameRow (row 4)
y  14
x  4
 DataFrameRow (row 3)
y  13
x  3
```
"""
eachrow(df::AbstractDataFrame) = DataFrameRows(df, index(df))

Base.IndexStyle(::Type{<:DataFrameRows}) = Base.IndexLinear()
Base.size(itr::DataFrameRows) = (size(itr.df, 1), )

Base.@propagate_inbounds Base.getindex(itr::DataFrameRows, i::Int) =
    DataFrameRow(itr.df, itr.index, i)
Base.@propagate_inbounds Base.getindex(itr::DataFrameRows{<:SubDataFrame}, i::Int) =
    DataFrameRow(parent(itr.df), itr.index, rows(itr.df)[i])

# Iteration by columns
"""
    DataFrameColumns{<:AbstractDataFrame, V} <: AbstractVector{V}

Iterator over columns of an `AbstractDataFrame` constructed using
[`eachcol(df, true)`](@ref) if `V` is a `Pair{Symbol,AbstractVector}`. Then each
returned value is a pair consisting of column name and column vector.
If `V` is an `AbstractVector` (a value returned by [`eachcol(df, false)`](@ref))
then each returned value is a column vector.
"""
struct DataFrameColumns{T<:AbstractDataFrame, V} <: AbstractVector{V}
    df::T
end

"""
    eachcol(df::AbstractDataFrame, names::Bool=true)

Return a `DataFrameColumns` that iterates an `AbstractDataFrame` column by column.
If `names` is equal to `true` (currently the default, in the future the default
will be set to `false`) iteration returns a pair consisting of column name
and column vector.
If `names` is equal to `false` then column vectors are yielded.

**Examples**

```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> collect(eachcol(df, true))
2-element Array{Pair{Symbol,AbstractArray{T,1} where T},1}:
 :x => [1, 2, 3, 4]
 :y => [11, 12, 13, 14]

julia> collect(eachcol(df, false))
2-element Array{AbstractArray{T,1} where T,1}:
 [1, 2, 3, 4]
 [11, 12, 13, 14]

julia> sum.(eachcol(df, false))
2-element Array{Int64,1}:
 10
 50

julia> map(eachcol(df, false)) do col
           maximum(col) - minimum(col)
       end
2-element Array{Int64,1}:
 3
 3
```
"""
@inline function eachcol(df::T, names::Bool) where T<: AbstractDataFrame
    if names
        DataFrameColumns{T, Pair{Symbol, AbstractVector}}(df)
    else
        DataFrameColumns{T, AbstractVector}(df)
    end
end

# TODO: remove this method after deprecation
# and add default argument value above
function eachcol(df::AbstractDataFrame)
    Base.depwarn("In the future eachcol will have names argument set to false by default", :eachcol)
    eachcol(df, true)
end

# TODO: remove this method after deprecation
# this is left to make sure we do not forget to properly fix columns calls
columns(df::AbstractDataFrame) = eachcol(df, false)

Base.size(itr::DataFrameColumns) = (size(itr.df, 2),)
Base.IndexStyle(::Type{<:DataFrameColumns}) = Base.IndexLinear()

@inline function Base.getindex(itr::DataFrameColumns{<:AbstractDataFrame,
                                                     Pair{Symbol, AbstractVector}}, j::Int)
    @boundscheck checkbounds(itr, j)
    @inbounds _names(itr.df)[j] => itr.df[j]
end

@inline function Base.getindex(itr::DataFrameColumns{<:AbstractDataFrame, AbstractVector}, j::Int)
    @boundscheck checkbounds(itr, j)
    @inbounds itr.df[j]
end

"""
    mapcols(f::Union{Function,Type}, df::AbstractDataFrame)

Return a `DataFrame` where each column of `df` is transformed using function `f`.
`f` must return `AbstractVector` objects all with the same length or scalars.

**Examples**

```jldoctest
julia> df = DataFrame(x=1:4, y=11:14)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 11    │
│ 2   │ 2     │ 12    │
│ 3   │ 3     │ 13    │
│ 4   │ 4     │ 14    │

julia> mapcols(x -> x.^2, df)
4×2 DataFrame
│ Row │ x     │ y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 121   │
│ 2   │ 4     │ 144   │
│ 3   │ 9     │ 169   │
│ 4   │ 16    │ 196   │
```
"""
function mapcols(f::Union{Function,Type}, df::AbstractDataFrame)
    # note: `f` must return a consistent length
    res = DataFrame()
    for (n, v) in eachcol(df, true)
        res[n] = f(v)
    end
    res
end


text:
"""
    @maybe_inbounds P blk

yields code that may compile block of code `blk` without bound checking if
allowed by `P`, that is if `P` inherits from `InBounds`.  This is equivalent to
(without escaping for clarity):

```julia
if \$P <: InBounds
     @inbounds \$blk
else
     \$blk
end
```

In words, bound checking is turned on/off at compilation time depending on the
type `P` without the needs to explicitly duplicate the code.

A typical usage is to write a method like:

```julia
function vsum(::Type{P}, x::AbstractArray{T}) where {T<:Real, P<:OptimLevel}
    s = zero(T)
    @maybe_inbounds P for i in eachindex(x)
        s += x[i]
    end
    return s
end
```

then `vsum(Debug,x)` will execute a version of the code that performs bound
checking; while `vsum(InBounds,x)` or `vsum(Vectorize,x)` will execute a
version of the code that avoids bound checking.

See [`@maybe_vectorized`], [`OptimLevel`].

"""
macro maybe_inbounds(P, blk)
    # See comments in `@maybe_vectorized`.
    opt = esc(:($P))
    blk0 = esc(:($blk))
    blk1 = esc(:(@inbounds $blk))
    quote
        if $opt <: InBounds
            $blk1
        elseif $opt <: Debug
            $blk0
        else
            error("argument `P` of macro `@maybe_inbounds` is not a valid optimization level")
        end
    end
end

"""
    @maybe_vectorized P blk

yields code that may compile block of code `blk` with vectorization and/or
without bound checking as allowed by the optimization level `P`.  This is
equivalent to (without escaping for clarity):

```julia
if \$P <:Vectorize
     @inbounds @simd \$blk
elseif \$P <: InBounds
     @inbounds \$blk
else
     \$blk
end
```

In words, vectorization and/or bound checking are turned on/off at compilation
time depending on the type `P` without the needs to explicitly duplicate the
code.

A typical usage:

```julia
function vsum(::Type{P}, x::AbstractArray{T}) where {T<:Real, P<:OptimLevel}
    s = zero(T)
    @maybe_vectorized P for i in eachindex(x)
        s += x[i]
    end
    return s
end
```

then `vsum(Debug,x)` will execute a version of the code that performs bound
checking and no vectorization; `vsum(InBounds,x)` will execute a version of the
code that avoids bound checking (and may vectorize code); finally
`vsum(Vectorize,x)` will execute a version of the code that avoids bound
checking and requires vectorization.

See [`@maybe_inbounds`], [`OptimLevel`].

"""
macro maybe_vectorized(P, blk)
    # Type to compare `P` with is to be interpreted in this module context
    # and thus must not be escaped.  But `opt` and `blk` must be escaped.
    # I empirically found that `esc(:(@macro1 @macro2 ... $expr))` was the
    # correct way to escape expression `expr`, possibly preceded by some
    # other macro calls; the alternative `:(esc($expr))` works for all
    # expressions below but not for the one with `@simd` because this macro
    # expects a `for` loop.
    opt = esc(:($P))
    blk0 = esc(:($blk))
    blk1 = esc(:(@inbounds $blk))
    blk2 = esc(:(@inbounds @simd $blk))
    quote
        if $opt <: Vectorize
            $blk2
        elseif $opt <: InBounds
            $blk1
        elseif $opt <: Debug
            $blk0
        else
            error("argument `P` of macro `@maybe_vectorized` is not a valid optimization level")
        end
    end
end

max_stars_repo_name: emmt/ConditionallyOptimize


text:
import NbodyGradient: set_state!, zero_out!, amatrix, initialize!

@testset "Transit Timing" begin
    N = 3
    t0 = 7257.93115525 - 7300.0 - 0.5 # Initialize IC before first transit
    h = 0.04
    itime = 10.0 # Amount of time the integrator will run
    tmax = itime

    # Setup initial conditions:
    elements = readdlm("elements.txt", ',')[1:N,:]
    elements[2:end,3] .-= 7300.0
    elements[:,end] .= 0.0 # Set Ωs to 0
    elements[2,1] *= 10.0
    elements[3,1] *= 10.0
    ic = ElementsIC(t0, N, elements)

    function calc_times(h, ti)
        intr = Integrator(ahl21!, h, 0.0, tmax)
        s = State(ic)
        tt = TransitTiming(itime, ic, ti)
        intr(s, tt)
        return tt
    end

    function calc_finite_diff(h, t0, itime, tmax, elements, ti)
        # Now do finite difference with big float
        dq0 = big(1e-10)
        ic_big = ElementsIC(big(t0), N, big.(elements))
        elements_big = copy(ic_big.elements)
        sp, dp = dState(ic_big)
        sm, dm = dState(ic_big)
        ttp = TransitTiming(big(itime), ic_big, ti);
        ttm = TransitTiming(big(itime), ic_big, ti);
        dtde_num = zeros(BigFloat, size(ttp.dtdelements));
        intr_big = Integrator(big(h), zero(BigFloat), big(tmax))

        for jq in 1:N
            for iq in 1:7
                zero_out!(ttp); zero_out!(ttm); zero_out!(dp); zero_out!(dm);
                ic_big.elements .= elements_big
                ic_big.m .= elements_big[:,1]
                if iq == 7; ivary = 1; else; ivary = iq + 1; end
                ic_big.elements[jq,ivary] += dq0
                if ivary == 1; ic_big.m[jq] += dq0; amatrix(ic_big); end # Masses don't update with elements array
                initialize!(sp, ic_big)
                intr_big(sp, ttp, dp; grad=false)
                ic_big.elements[jq,ivary] -= 2dq0
                if ivary == 1; ic_big.m[jq] -= 2dq0; amatrix(ic_big); end
                initialize!(sm, ic_big)
                intr_big(sm, ttm, dm; grad=false)
                for i in ttm.occs
                    for k in 1:ttm.count[i]
                    # Compute double-sided derivative for more accuracy:
                        dtde_num[i,k,iq,jq] = (ttp.tt[i,k] - ttm.tt[i,k]) / (2dq0)
                    end
                end
            end
        end
        return dtde_num
    end

    for ti in 1:N
        # Calculate transit times
        hs = h ./ [1.0,2.0,4.0,8.0]
        tts = TransitTiming[]
        for h in hs
            push!(tts, calc_times(h, ti))
        end

        mask = zeros(Bool, size(tts[2].dtdq0));
        for jq = 1:N
            for iq = 1:7
                if iq == 7; ivary = 1; else; ivary = iq + 1; end  # Shift mass variation to end
                for i = tts[1].occs
                    for k = 1:tts[1].count[i]
                    # Ignore inclination & longitude of nodes variations:
                        if iq != 5 && iq != 6 && ~(jq == 1 && iq < 7) && ~(jq == i && iq == 7)
                            mask[i,k,iq,jq] = true
                        end
                    end
                end
            end
        end

        dtde_num = calc_finite_diff(h, t0, itime, tmax, elements, ti)
        @test isapprox(asinh.(tts[1].dtdelements[mask]), asinh.(dtde_num[mask]);norm=maxabs)
        @test isapprox(asinh.(tts[1].dtdelements), asinh.(dtde_num);norm=maxabs)
    end

    # Test that the transit times for the derivatives and non-derivatives match
    tmax = 2000.0
    s_nograd = State(ic); tt_nograd = TransitTiming(tmax, ic);
    s_grad = State(ic); tt_grad = TransitTiming(tmax, ic);
    Integrator(h, tmax)(s_nograd, tt_nograd, grad=false)
    Integrator(h, tmax)(s_grad, tt_grad, grad=true)
    @test tt_nograd.tt == tt_grad.tt
end


text:
include("MatrixProductState.jl")
include("SparseTensor.jl")

#identity activation
self(x) = x
self_deriv(x) = 1
#sigmoid activation
sigmoid(x) = 1 / (1+exp(-x))
sigmoid_deriv(x) = sigmoid(x)*(1-sigmoid(x))
#binaray crossentropy loss
bce(y, f) = -y*log(f+eps()) - (1-y)*log(1-f+eps())
bce_deriv(y,f) = -y/(f+eps()) + (1-y)/(1-f+eps())
#squared error loss
se(y,f) = (f-y)^2
se_deriv(y,f) = 2*(f-y)
# l2 regularization
l2(A,K) = sum(A .* (A*K))
l2_deriv(A,K) = 2*(A*K)

# Loss: generic for all combinations in losses and regs: F and G are value and derivative w.r.t. A
# l = given as argument = l(y,f)  with f = a(A*K)
# r = given as argumnet = r(A,K)
# note that Y is a dense label tensor. inds denote what labels are really observed to evaluate the lossfunction.
function Loss_reg!(F,G,A,K,act,actderiv,λ,Y,inds, loss,lossderiv, reg, regderiv)
    z = A*K
    f = act.(z)
    if G!=nothing
        dldf_dfdz = zeros(size(f))
        for I in inds
            dldf_dfdz[I] = lossderiv(Y[I], f[I]) * actderiv(z[I])
        end
        dldf_dfdz_dzdA = dldf_dfdz * K
        dr_dA = regderiv(A,K)
        G[:] = dldf_dfdz_dzdA + λ * dr_dA
    end
    if F != nothing
        l = 0.
        for I in inds
            l+= loss(Y[I],f[I])
        end
        return l + λ * reg(A,K)
    end
end

# Loss: lossfunction for n-step kernel ridge regression: F and G are value and derivative w.r.t. A
# l = SE
# r = f * (H^(-1) - 1) * f
function NSKRR_loss!(F,G,A,K,inverseH, Y, inds)
    f = A*K
    if G!=nothing
        dldf = zeros(size(Y))
        for I in inds
            dldf[I] = 2*(f[I] - Y[I])
        end
        dldf_dfdA = dldf*K
        drdA = 2*(inverseH*f -f)*K
        G[:] = dldf_dfdA + drdA
    end
    if F!=nothing
        SE = 0.
        for I in inds
            SE+= (f[I]-Y[I])^2
        end
        return SE + sum(f.*(inverseH*f -f))
    end
end

## Below losses specific implemented for rigde regression and logistic ridge regression
## However not really needed anymore due to more generic formulation above, which as nearly as efficient
"""
# Loss: squared error + L2 regularization: F and G are value and derivative w.r.t. A
# l = SE = (f -Y)^2  with f = A*K
# r = L2 = A*K*A
function SE_reg!(F,G,A,K,λ,Y,inds)
    f = A*K
    if G != nothing
        dldf = zeros(size(f))
        for I in inds
            dldf[I] = 2*(f[I]-Y[I])
        end
        dldf_dfdA = dldf*K
        dr_dA = 2*A*K
        G[:] = dldf_dfdA + λ * dr_dA
    end
    if F != nothing
        SE = 0.
        for I in inds
            SE+= (f[I]-Y[I])^2
        end
        return SE + λ * sum(A.*(A*K))
    end
end

bce(1,1)
bce(1,0)
bce_deriv(1,0)

# Loss: binary crossentropy + L2 regularization: F and G are value and derivative w.r.t. A
# l = BCE = -[ y*log(f) + (1-y)* log(1-f) ] with f = a(A*K)
# r = L2 = A*K*A
function BCE_reg!(F,G,A,K,a,aderiv,λ,Y,inds)
    z = A*K
    f = a.(z)
    if G!=nothing
        dldf_dfdz = zeros(size(f))
        for I in inds
            dldf_dfdz[I] = bce_deriv(Y[I], f[I]) * aderiv(z[I])
        end
        dldf_dfdz_dzdA = dldf_dfdz * K
        dr_dA = 2*A*K
        G[:] = dldf_dfdz_dzdA + λ * dr_dA
    end
    if F != nothing
        BCE = 0.
        for I in inds
            BCE+= bce(Y[I],f[I])
        end
        return BCE + λ * sum(A.*(A*K))
    end
end
"""

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
include("../domain/domain.jl")
include("./doc.jl")

# https://www.investopedia.com/terms/j/journal.asp

create_journal_entry(
    id::String,
    customer_id::String,
    invoice_nbr::String,
    from::Int64,
    to::Int64,
    debit::Float64,
    credit::Float64,
    vat::Float64,
    descr::String) = JournalEntry(
        id,
        Dates.month(now()),
        now(),
        customer_id,
        invoice_nbr,
        from,
        to,
        debit,
        credit,
        vat,
        descr)

create_account(id::String, name::String, group::String) = Account(id, name)

create_record_from(entry::JournalEntry) = Record(
    entry.id,
    entry.from,
    entry.date,
    entry.customer_id,
    entry.invoice_nbr,
    entry.debit + entry.vat,
    0,
    entry.descr)

create_record_to(entry::JournalEntry) = Record(
    entry.id,
    entry.to,
    entry.date,
    entry.customer_id,
    entry.invoice_nbr,
    0,
    entry.debit,
    entry.descr)

create_record_vat(entry::JournalEntry, vat_account_nbr::Int) = Record(
    entry.id,
    vat_account_nbr,
    entry.date,
    entry.customer_id,
    entry.invoice_nbr,
    0,
    entry.vat,
    entry.descr)


text:
# Checkers #
# ========================== #
check_generators(stochasticprogram::StochasticProgram{N}) where N = _check_generators(stochasticprogram, Val(N))
function _check_generators(stochasticprogram::StochasticProgram, ::Val{N}) where N
    _check_stage_generator(stochasticprogram, N)
    _check_generators(stochasticprogram, Val(N-1))
    return nothing
end
function _check_generators(stochasticprogram::StochasticProgram, ::Val{1})
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    return nothing
end
function _check_stage_generator(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    stage_key = Symbol(:stage_, s)
    has_generator(stochasticprogram, stage_key) || error("Stage problem $stage not defined in stochastic program. Consider @stage $stage.")
    return nothing
end
# Model generation #
# ========================== #
"""
    stage_one_model(stochasticprogram::StochasticProgram; optimizer = nothing)

Return a generated copy of the first stage model in `stochasticprogram`. Optionally, supply a capable `optimizer` to the stage model.
"""
function stage_one_model(stochasticprogram::StochasticProgram; optimizer = nothing)
    # Return possibly cached model
    cache = problemcache(stochasticprogram)
    if haskey(cache, :stage_1)
        m = cache[:stage_1]
        optimizer != nothing && set_optimizer(m, optimizer)
        return m
    end
    # Check that the required generators have been defined
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    model = optimizer == nothing ? Model() : Model(optimizer)
    # Prepare decisions
    model.ext[:decisions] = (Decisions(),)
    add_decision_bridges!(model)
    # Generate and cache first-stage model
    generator(stochasticprogram, :stage_1)(model, stage_parameters(stochasticprogram, 1))
    cache[:stage_1] = model
    return model
end
"""
    stage_model(stochasticprogram::StochasticProgram, stage::Integer, scenario::AbstractScenario; optimizer = nothing)

Return a generated stage model corresponding to `scenario`, in `stochasticprogram`. Optionally, supply a capable `optimizer` to the stage model.
"""
function stage_model(stochasticprogram::StochasticProgram{N},
                     stage::Integer,
                     scenario::AbstractScenario;
                     optimizer = nothing) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N")
    stage == 1 && return stage_one_model(stochasticprogram, optimizer)
    # Check generators
    stage_key = Symbol(:stage_, stage)
    decision_key = Symbol(:stage_, stage - 1, :_decisions)
    has_generator(stochasticprogram, stage_key) || error("Stage problem $stage not defined in stochastic program. Consider @stage $stage")
    has_generator(stochasticprogram, decision_key) || error("No decision variables defined in stage problem $(stage-1).")
    # Create stage model
    stage_model = optimizer == nothing ? Model() : Model(optimizer)
    # Prepare decisions
    stage_model.ext[:decisions] = Decisions()
    add_decision_bridges!(stage_model)
    # Generate and return the stage model
    generator(stochasticprogram, decision_key)(stage_model, decision_params)
    generator(stochasticprogram, stage_key)(stage_model, stage_params, scenario)
    return stage_model
end
function generate_stage_one!(stochasticprogram::StochasticProgram)
    haskey(stochasticprogram.problemcache, :stage_1) && return nothing
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    stochasticprogram.problemcache[:stage_1] = JuMP.Model()
    generator(stochasticprogram, :stage_1)(stochasticprogram.problemcache[:stage_1], stage_parameters(stochasticprogram, 1))
    return nothing
end

function generate_proxy!(stochasticprogram::StochasticProgram{N}) where N
    # First-stage decisions are unique (reuse)
    stochasticprogram.proxy[1].ext[:stage_map] = Dict{MOI.VariableIndex, Int}()
    stochasticprogram.proxy[1].ext[:decisions] = (stochasticprogram.decisions[1],)
    # Generate first stage
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    generator(stochasticprogram, :stage_1)(stochasticprogram.proxy[1], stage_parameters(stochasticprogram, 1))
    # Generate remaining stages
    for s in 2:N
        # Initialize decisions
        stochasticprogram.proxy[s].ext[:stage_map] = Dict{MOI.VariableIndex, Int}()
        stochasticprogram.proxy[s].ext[:decisions] = ntuple(Val{s}()) do i
            if i == s - 1
                # Known decisions from the previous stages are
                # the same everywhere.
                return stochasticprogram.decisions[s]
            end
            return Decisions()
        end
        # Check generator
        stage_key = Symbol(:stage_, s)
        decision_key = Symbol(:stage_, s - 1, :_decisions)
        has_generator(stochasticprogram, stage_key) || error("Stage problem $stage not defined in stochastic program. Consider @stage $stage.")
        has_generator(stochasticprogram, decision_key) || error("No decision variables defined in stage problem $(stage-1).")
        # Generate
        generator(stochasticprogram, decision_key)(stochasticprogram.proxy[s], stage_parameters(stochasticprogram, s-1))
        generator(stochasticprogram, stage_key)(stochasticprogram.proxy[s], stage_parameters(stochasticprogram, s), scenario(stochasticprogram, s, 1))
    end
    return nothing
end
"""
    clear!(stochasticprogram::StochasticProgram)

Clear the `stochasticprogram`, removing all model definitions.
"""
function clear!(stochasticprogram::StochasticProgram{N}) where N
    # Clearing invalidates the cache
    invalidate_cache!(stochasticprogram)
    # Dispatch clearup on stochastic structure
    clear!(structure(stochasticprogram))
    # Clear proxy
    for s in 1:N
        proxy_ = proxy(stochasticprogram, s)
        # Clear decisions
        if haskey(proxy_.ext, :decisions)
            map(clear!, proxy_.ext[:decisions])
        end
        # Clear model
        empty!(proxy_)
    end
    return nothing
end
"""
    generate!(stochasticprogram::StochasticProgram)

Generate the `stochasticprogram` using the model definitions from @stage and available data.
"""
function generate!(stochasticprogram::StochasticProgram{N}) where N
    # Clear stochasticprogram before re-generation
    clear!(stochasticprogram)
    # Abort early if there are no scenarios
    num_scenarios(stochasticprogram) == 0 && return nothing
    # Check generators
    check_generators(stochasticprogram::StochasticProgram{N})
    # Generate proxy
    generate_proxy!(stochasticprogram)
    # Dispatch generation on stochastic structure
    generate!(stochasticprogram, structure(stochasticprogram))
    return nothing
end
# Outcome model generation #
# ========================== #
function _outcome_model!(outcome_model::JuMP.Model,
                         decision_generator::Function,
                         generator::Function,
                         decision_params::Any,
                         stage_params::Any,
                         stage_one_model::JuMP.Model,
                         decisions::AbstractVector,
                         scenario::AbstractScenario)
    # Prepare decisions
    outcome_model.ext[:decisions] = (Decisions(), Decisions())
    add_decision_bridges!(outcome_model)
    # Generate the outcome model
    decision_generator(outcome_model, decision_params)
    generator(outcome_model, stage_params, scenario)
    # Copy first-stage objective
    copy_decision_objective!(stage_one_model,
                             outcome_model,
                             all_known_decision_variables(outcome_model)[1])
    # Update the known decision values
    update_known_decisions!(outcome_model, decisions)
    return nothing
end
"""
    outcome_model(stochasticprogram::TwoStageStochasticProgram,
                  decision::AbstractVector,
                  scenario::AbstractScenario;
                  optimizer = nothing)

Return the resulting second stage model if `decision` is the first-stage decision in the provided `scenario`, in `stochasticprogram`. The supplied `decision` must match the defined decision variables in `stochasticprogram`. Optionally, supply a capable `optimizer` to the outcome model.
"""
function outcome_model(stochasticprogram::TwoStageStochasticProgram,
                       decision::AbstractVector,
                       scenario::AbstractScenario;
                       optimizer = nothing)
    # Check generators
    has_generator(stochasticprogram,:stage_1_decisions) || error("First-stage not defined in stochastic program. Consider @first_stage or @stage 1.")
    has_generator(stochasticprogram,:stage_2) || error("Second-stage problem not defined in stochastic program. Consider @second_stage.")
    # Sanity checks on given decision vector
    length(decision) == num_decisions(stochasticprogram) || error("Incorrect length of given decision vector, has ", length(decision), " should be ", num_decisions(stochasticprogram))
    all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
    # Create outcome model
    outcome_model = optimizer == nothing ? Model() : Model(optimizer)
    _outcome_model!(outcome_model,
                    generator(stochasticprogram,:stage_1_decisions),
                    generator(stochasticprogram,:stage_2),
                    stage_parameters(stochasticprogram, 1),
                    stage_parameters(stochasticprogram, 2),
                    stage_one_model(stochasticprogram),
                    decision,
                    scenario)
    return outcome_model
end
"""
    outcome_model(stochasticprogram::StochasticProgram{N},
                  decisions::NTuple{N-1,AbstractVector}
                  scenario_path::NTuple{N-1,AbstractScenario},
                  optimizer = nothing)

Return the resulting `N`:th stage model if `decisions` are the decisions taken in the previous stages and `scenario_path` are the realized scenarios up to stage `N` in `stochasticprogram`. Optionally, supply a capable `solver` to the outcome model.
"""
function outcome_model(stochasticprogram::StochasticProgram{N},
                       decisions::NTuple{M,AbstractVector},
                       scenario_path::NTuple{M,AbstractScenario};
                       optimizer = nothing) where {N,M}
    N == M - 1 || error("Inconsistent number of stages $N and number of decisions and scenarios $M")
    # TODO
end
# ========================== #


text:
#==============================================================================#
#       ___ _   _ ___  ____ ____      /    ____ ____ _  _ ____ ___ ____        #
#        |   \_/  |__] |___ [__      /     |    |  | |\ | [__   |  |__/        #
#        |    |   |    |___ ___]    /      |___ |__| | \| ___]  |  |  \ .      #
#==============================================================================#

PFB{T} = Matrix{T}          # polyphase filter bank
PNFB{T} = Vector{Poly{T}}   # polynomial filter bank (used for farrow filter)

abstract type Filter end
abstract type FIRKernel end
# TODO: all kernels: add field original taps

# Single rate FIR kernel
type FIRStandard{T} <: FIRKernel
    h::Vector{T}
    hLen::Int
end

function FIRStandard( h::Vector )
    h    = flipdim( h, 1 )
    hLen = length( h )
    FIRStandard( h, hLen )
end

# Interpolator FIR kernel
type FIRInterpolator{T} <: FIRKernel
    pfb::PFB{T}
    interpolation::Int
    N𝜙::Int
    tapsPer𝜙::Int
end

function FIRInterpolator( h::Vector, interpolation::Integer )
    pfb           = taps2pfb( h, interpolation )
    tapsPer𝜙      = size( pfb )[1]
    N𝜙            = size( pfb )[2]
    interpolation = interpolation
    FIRInterpolator( pfb, interpolation, N𝜙, tapsPer𝜙 )
end

# Decimator FIR kernel
type FIRDecimator{T} <: FIRKernel
    h::Vector{T}
    hLen::Int
    decimation::Int
    inputDeficit::Int
end

function FIRDecimator( h::Vector, decimation::Integer )
    h            = flipdim( h, 1 )
    hLen         = length( h )
    decimation   = decimation
    inputDeficit = 1
    FIRDecimator( h, hLen, decimation, inputDeficit )
end

# Rational resampler FIR kernel
type FIRRational{T}  <: FIRKernel
    pfb::PFB{T}
    ratio::Rational{Int}
    N𝜙::Int
    tapsPer𝜙::Int
    criticalYidx::Int
    𝜙Idx::Int
    inputDeficit::Int
end

function FIRRational( h::Vector, ratio::Rational )
    pfb          = taps2pfb( h, numerator(ratio) )
    N𝜙           = size( pfb )[2]
    tapsPer𝜙     = size( pfb )[1]
    criticalYidx = floor( Int, tapsPer𝜙 * ratio )
    𝜙Idx         = 1
    inputDeficit = 1
    FIRRational( pfb, ratio, N𝜙, tapsPer𝜙, criticalYidx, 𝜙Idx, inputDeficit )
end

# Arbitrary resampler FIR kernel
# This kernel is different from the others in that it has two polyphase filtlter banks.
# The the second filter bank, dpfb, is the derivative of pfb. The purpose of this is to
# allow us to compute two y values, yLower & yUpper, whitout having to advance the input
# index by 1. It makes the kernel simple by not having to store extra state in the case
# when where's at the last polphase branch and the last available input sample. By using
# a derivitive filter, we can always compute the output in that scenario.
# See section 7.6.1 in [1] for a better explanation.
type FIRArbitrary{T} <: FIRKernel # TODO: since farrow is also arbitrary, find a new name
    rate::Float64
    pfb::PFB{T}
    dpfb::PFB{T}
    N𝜙::Int
    tapsPer𝜙::Int
    𝜙Accumulator::Float64
    𝜙Idx::Int
    α::Float64
    Δ::Float64
    inputDeficit::Int
    xIdx::Int
end

function FIRArbitrary{T <: Number}( h::Vector{T}, rate::Real, N𝜙::Integer )
    dh           = [ diff( h ); zero(T) ]
    pfb          = taps2pfb( h,  N𝜙 )
    dpfb         = taps2pfb( dh, N𝜙 )
    tapsPer𝜙     = size( pfb )[1]
    𝜙Accumulator = 1.0
    𝜙Idx         = 1
    α            = 0.0
    Δ            = N𝜙/rate
    inputDeficit = 1
    xIdx         = 1
    FIRArbitrary( rate, pfb, dpfb, N𝜙, tapsPer𝜙, 𝜙Accumulator, 𝜙Idx, α, Δ, inputDeficit, xIdx )
end

# Farrow filter kernel.
# Takes a polyphase filterbank and converts each row of taps into a polynomial.
# That we can calculate filter tap values for any arbitrary 𝜙Idx, not just integers between 1 and N𝜙
type FIRFarrow{T} <: FIRKernel
    rate::Float64
    pfb::PFB{T}
    pnfb::PNFB{T}
    polyorder::Int
    currentTaps::Vector{T}
    N𝜙::Int
    tapsPer𝜙::Int
    𝜙Idx::Float64
    Δ::Float64
    inputDeficit::Int
    xIdx::Int
end

function FIRFarrow{T}( h::Vector{T}, rate::Real, N𝜙::Integer, polyorder::Integer )
    pfb          = taps2pfb( h,  N𝜙 )
    pnfb         = pfb2pnfb( pfb, polyorder )
    tapsPer𝜙     = size( pfb )[1]
    𝜙Idx         = 1.0
    Δ            = N𝜙/rate
    inputDeficit = 1
    xIdx         = 1
    currentTaps  = T[ polyval( pnfb[tapIdx], 𝜙Idx ) for tapIdx in 1:tapsPer𝜙 ]
    FIRFarrow( rate, pfb, pnfb, polyorder, currentTaps, N𝜙, tapsPer𝜙, 𝜙Idx, Δ, inputDeficit, xIdx )
end

# FIRFilter - the kernel does the heavy lifting
type FIRFilter{Tk<:FIRKernel} <: Filter
    kernel::Tk
    history::Vector
    historyLen::Int
end

# Constructor for single-rate, decimating, interpolating, and rational resampling filters
function FIRFilter( h::Vector, resampleRatio::Rational = 1//1 )
    interpolation = numerator( resampleRatio )
    decimation    = denominator( resampleRatio )
    historyLen    = 0

    if resampleRatio == 1                                     # single-rate
        kernel     = FIRStandard( h )
        historyLen = kernel.hLen - 1
    elseif interpolation == 1                                 # decimate
        kernel     = FIRDecimator( h, decimation )
        historyLen = kernel.hLen - 1
    elseif decimation == 1                                    # interpolate
        kernel     = FIRInterpolator( h, interpolation )
        historyLen = kernel.tapsPer𝜙 - 1
    else                                                      # rational
        kernel     = FIRRational( h, resampleRatio )
        historyLen = kernel.tapsPer𝜙 - 1
    end

    history = zeros( historyLen )

    FIRFilter( kernel, history, historyLen )
end

# Constructor for arbitrary resampling filter (polyphase interpolator w/ intra-phase linear interpolation )
function FIRFilter( h::Vector, rate::AbstractFloat, N𝜙::Integer = 32 )
    rate > 0.0 || error( "rate must be greater than 0" )
    kernel     = FIRArbitrary( h, rate, N𝜙 )
    historyLen = kernel.tapsPer𝜙 - 1
    history    = zeros( historyLen )
    FIRFilter( kernel, history, historyLen )
end

# Constructor for farrow filter (polyphase interpolator w/ polynomial genrated intra-phase taps )
function FIRFilter( h::Vector, rate::AbstractFloat, N𝜙::Integer, polyorder::Integer )
    rate > 0.0 || error( "rate must be greater than 0" )
    kernel     = FIRFarrow( h, rate, N𝜙, polyorder )
    historyLen = kernel.tapsPer𝜙 - 1
    history    = zeros( historyLen )
    FIRFilter( kernel, history, historyLen )
end

#==============================================================================#
#                    ____ ____ ___ ___  _  _ ____ ____ ____                    #
#                    [__  |___  |  |__] |__| |__| [__  |___                    #
#                    ___] |___  |  |    |  | |  | ___] |___                    #
#==============================================================================#
# Sets the kernel's phase (𝜙Idx+α).
#   Valid input is [0, 1]

function setphase( kernel::Union{FIRInterpolator, FIRRational}, 𝜙::Number )
    @assert zero(𝜙) <= 𝜙 <= one(𝜙)
    kernel.𝜙Idx = Int(𝜙Idx)
    return kernel.𝜙Idx
end

function setphase( kernel::FIRArbitrary, 𝜙::Number )
    @assert zero(𝜙) <= 𝜙 <= one(𝜙)
    (α, 𝜙Idx)   = modf( 𝜙 * kernel.N𝜙 )
    kernel.𝜙Idx = Int(𝜙Idx)
    kernel.α    = α
    return 𝜙Idx, α
end

function setphase( kernel::FIRFarrow, 𝜙::Number )
    @assert zero(𝜙) <= 𝜙 <= one(𝜙)
    kernel.𝜙Idx = 𝜙*(kernel.N𝜙-1)+1
    tapsforphase!( kernel.currentTaps, kernel, kernel.𝜙Idx  )
    return kernel.𝜙Idx
end

setphase( self::FIRFilter, 𝜙::Number ) = setphase( self.kernel, 𝜙 )

#==============================================================================#
#                            ____ ____ ____ ____ ___                           #
#                            |__/ |___ [__  |___  |                            #
#                            |  \ |___ ___] |___  |                            #
#==============================================================================#

# Resets filter and its kernel to an initial state

# Does nothing for non-rational kernels
reset( self::FIRKernel ) = self

# For rational kernel, set 𝜙Idx back to 1
reset( self::FIRRational ) = self.𝜙Idx = 1

# For rational kernel, set 𝜙Idx back to 1
function reset( self::FIRArbitrary )
    self.yCount = 0
    update( self )
end

# For FIRFilter, set history vector to zeros of same type and required length
function reset( self::FIRFilter )
    self.history = zeros( eltype( self.history ), self.historyLen )
    reset( self.kernel )
    return self
end

#==============================================================================#
#                      _  _    ___ ____    ___  ____ ___                       #
#                      |__|     |  |  |    |__] |___ |__]                      #
#                      |  |     |  |__|    |    |    |__]                      #
#==============================================================================#

# Converts a vector of coefficients to a matrix. Each column is a filter.
# NOTE: also flips the matrix up/down so computing the dot product of a
#       column with a signal vector is more efficient (since filter is convolution)
# Appends zeros if necessary.
# Example:
#   julia> taps2pfb( [1:9], 4 )
#   3x4 Array{Int64,2}:
#    9  0  0  0
#    5  6  7  8
#    1  2  3  4
#
#  In this example, the first phase, or 𝜙, is [9, 5, 1].

function taps2pfb{T}( h::Vector{T}, N𝜙::Integer )
    hLen     = length( h )
    tapsPer𝜙 = ceil( Int, hLen/N𝜙 )
    pfbSize  = tapsPer𝜙 * N𝜙
    pfb      = Array{T}( tapsPer𝜙, N𝜙 )
    hIdx     = 1

    for rowIdx in tapsPer𝜙:-1:1, colIdx in 1:N𝜙
        tap = hIdx > hLen ? zero(T) : h[hIdx]
        @inbounds pfb[rowIdx,colIdx] = tap
        hIdx += 1
    end

    return pfb
end

#==============================================================================#
#        ___  ____ _    _   _ _  _ ____ _  _ _ ____ _       ____ ___           #
#        |__] |  | |     \_/  |\ | |  | |\/| | |__| |       |___ |__]          #
#        |    |__| |___   |   | \| |__| |  | | |  | |___    |    |__]          #
#==============================================================================#

# Convert a polyphase filterbank into a polynomial filterbank

function pfb2pnfb{T}( pfb::PFB{T}, polyorder::Integer )
    (tapsPer𝜙, N𝜙) = size( pfb )
    result         = Array{Poly{T}}( tapsPer𝜙 )

    for i in 1:tapsPer𝜙
        row = vec( pfb[i,:] )
        result[i] = polyfit( row, polyorder )
    end

    return result
end

function taps2pnfb{T}( h::Vector{T}, N𝜙::Integer, polyorder::Integer )
    hLen     = length( h )
    tapsPer𝜙 = ceil( Int, hLen/N𝜙 )
    pnfb     = Array{Poly{T}}( tapsPer𝜙 )
    pfbSize  = N𝜙 * tapsPer𝜙
    h        = hLen < pfbSize + 1 ? [ h, zeros( T, pfbSize+1-hLen ) ] : h

    pnfbIdx = tapsPer𝜙
    for startIdx in 0:N𝜙:hLen-N𝜙
        row           = h[startIdx+1:startIdx+N𝜙+1]
        pnfb[pnfbIdx] = polyfit( row, polyorder )
        pnfbIdx      -= 1
    end

    return pnfb
end

#==============================================================================#
#               ____ _  _ ___ ___  _  _ ___    _    ____ _  _                  #
#               |  | |  |  |  |__] |  |  |     |    |___ |\ |                  #
#               |__| |__|  |  |    |__|  |     |___ |___ | \|                  #
#==============================================================================#

# Calculates the resulting length of a multirate filtering operation, given a
#   FIRFilter{FIRRational} object and an input vector
#
# ( It's hard to explain how this works without a diagram )

function outputlength( inputlength::Integer, ratio::Rational, initial𝜙::Integer )
    interpolation = numerator( ratio )
    decimation    = denominator( ratio )
    outLen        = (( inputlength * interpolation ) - initial𝜙 + 1 ) / decimation
    ceil( Int, outLen )
end

function outputlength( kernel::FIRStandard, inputlength::Integer )
    inputlength
end

function outputlength( kernel::FIRInterpolator, inputlength::Integer )
    kernel.interpolation * inputlength
end

function outputlength( kernel::FIRDecimator, inputlength::Integer )
    outputlength( inputlength-kernel.inputDeficit+1, 1//kernel.decimation, 1 )
end

function outputlength( kernel::FIRRational, inputlength::Integer )
    outputlength( inputlength-kernel.inputDeficit+1, kernel.ratio, kernel.𝜙Idx )
end

function outputlength( kernel::FIRArbitrary, inputlength::Integer )
    ceil( Int, (inputlength-kernel.inputDeficit+1) * kernel.rate )
end

function outputlength( kernel::FIRFarrow, inputlength::Integer )
    ceil( Int, (inputlength-kernel.inputDeficit+1) * kernel.rate )
end

function outputlength( self::FIRFilter, inputlength::Integer )
    outputlength( self.kernel, inputlength )
end

#==============================================================================#
#                 _ _  _ ___  _  _ ___    _    ____ _  _                       #
#                 | |\ | |__] |  |  |     |    |___ |\ |                       #
#                 | | \| |    |__|  |     |___ |___ | \|                       #
#==============================================================================#

function inputlength( outputlength::Int, ratio::Rational, initial𝜙::Integer )
    interpolation = numerator( ratio )
    decimation    = denominator( ratio )
    inLen         = ( outputlength * decimation + initial𝜙 - 1 ) / interpolation
    ceil( Int, inLen )
end

function inputlength( self::FIRFilter{FIRStandard}, outputlength::Integer )
    outputlength
end

function inputlength( self::FIRFilter{FIRInterpolator}, outputlength::Integer )
    kernel = self.kernel
    inputlength( outputlength, kernel.interpolation//1, 1 )
end

function inputlength( self::FIRFilter{FIRDecimator}, outputlength::Integer )
    kernel = self.kernel
    inLen  = inputlength( outputlength, 1//kernel.decimation, 1 )
    inLen  = inLen + kernel.inputlength - 1
end

function inputlength( self::FIRFilter{FIRRational}, outputlength::Integer )
    kernel = self.kernel
    inLen = inputlength( outputlength, kernel.ratio, kernel.𝜙Idx )
    inLen = inLen + kernel.inputDeficit - 1
end

#==============================================================================#
#              _  _ ____ _  _ ___    ___  _  _ ____ ____ ____                  #
#              |\ | |___  \/   |     |__] |__| |__| [__  |___                  #
#              | \| |___ _/\_  |     |    |  | |  | ___] |___                  #
#==============================================================================#

function nextphase( currentphase::Integer, ratio::Rational )
    interpolation = numerator( ratio )
    decimation    = denominator( ratio )
    𝜙Step         = mod( decimation, interpolation )
    𝜙Next         = currentphase + 𝜙Step
    𝜙Next         = 𝜙Next > interpolation ? 𝜙Next - interpolation : 𝜙Next
end

#==============================================================================#
#               ____ _ _  _ ____ _    ____    ____ ____ ___ ____               #
#               [__  | |\ | | __ |    |___    |__/ |__|  |  |___               #
#               ___] | | \| |__] |___ |___    |  \ |  |  |  |___               #
#==============================================================================#

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRStandard{Th}}, x::Vector{Tx} )
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    hLen                = kernel.hLen
    historyLen          = self.historyLen
    bufLen              = length( buffer )
    xLen                = length( x )
    outLen              = xLen
    criticalYidx        = min( hLen, outLen )

    bufLen >= xLen || error( "buffer length must be >= x length" )

    for yIdx in 1:criticalYidx        # this first loop takes care of filter ramp up and previous history
        @inbounds buffer[yIdx] = unsafedot( kernel.h, history, x, yIdx )
    end

    for yIdx in criticalYidx+1:xLen
        @inbounds buffer[yIdx] = unsafedot( kernel.h, x, yIdx )
    end

    self.history = shiftin!( history, x )

    return buffer
end

function filt{Th,Tx}( self::FIRFilter{FIRStandard{Th}}, x::Vector{Tx} )
    buffer = Array{promote_type(Th, Tx)}( length(x) )
    filt!( buffer, self, x )
end

#==============================================================================#
#               _ _  _ ___ ____ ____ ___  _    ____ ____ ___ ____              #
#               | |\ |  |  |___ |__/ |__] |    |  | |__|  |  |___              #
#               | | \|  |  |___ |  \ |    |___ |__| |  |  |  |___              #
#==============================================================================#

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRInterpolator{Th}}, x::Vector{Tx} )
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    interpolation       = kernel.interpolation
    N𝜙                  = kernel.N𝜙
    tapsPer𝜙            = kernel.tapsPer𝜙
    xLen                = length( x )
    bufLen              = length( buffer )
    historyLen          = self.historyLen
    outLen              = outputlength( self, xLen )
    criticalYidx        = min( historyLen*interpolation, outLen )
    inputIdx            = 1
    𝜙                   = 1

    bufLen >= outLen || error( "length( buffer ) must be >= interpolation * length(x)")

    for yIdx in 1:criticalYidx
        @inbounds buffer[yIdx] = unsafedot( kernel.pfb, 𝜙, history, x, inputIdx )
        (𝜙, inputIdx) = 𝜙 == N𝜙 ? ( 1, inputIdx+1 ) : ( 𝜙+1, inputIdx )
    end
    for yIdx in criticalYidx+1:outLen
        @inbounds buffer[yIdx] = unsafedot( kernel.pfb, 𝜙, x, inputIdx )
        (𝜙, inputIdx) = 𝜙 == N𝜙 ? ( 1, inputIdx+1 ) : ( 𝜙+1, inputIdx )
    end

    self.history = shiftin!( history, x )

    return buffer
end

function filt{Th,Tx}( self::FIRFilter{FIRInterpolator{Th}}, x::Vector{Tx} )
    xLen   = length( x )
    outlen = outputlength( self, xLen )
    buffer = Array{promote_type(Th,Tx)}( outlen )
    filt!( buffer, self, x )
    return buffer
end

#==============================================================================#
#           ____ ____ ___     ____ ____ ____ ____ _  _ ___  _    ____          #
#           |__/ |__|  |      |__/ |___ [__  |__| |\/| |__] |    |___          #
#           |  \ |  |  |  .   |  \ |___ ___] |  | |  | |    |___ |___          #
#==============================================================================#

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRRational{Th}}, x::Vector{Tx} )
    kernel              = self.kernel
    history::Vector{Tx} = self.history
    xLen                = length( x )
    bufLen              = length( buffer )
    bufIdx              = 0

    if xLen < kernel.inputDeficit
        self.history = shiftin!( history, x )
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    outLen = outputlength( xLen-kernel.inputDeficit+1, kernel.ratio, kernel.𝜙Idx )
    bufLen >= outLen || error( "buffer is too small" )

    interpolation       = numerator( kernel.ratio )
    decimation          = denominator( kernel.ratio )
    𝜙IdxStepSize        = mod( decimation, interpolation )
    critical𝜙Idx        = kernel.N𝜙 - 𝜙IdxStepSize
    inputIdx            = kernel.inputDeficit

    while inputIdx <= xLen
        bufIdx += 1
        if inputIdx < kernel.tapsPer𝜙
            accumulator = unsafedot( kernel.pfb, kernel.𝜙Idx, history, x, inputIdx )
        else
            accumulator = unsafedot( kernel.pfb, kernel.𝜙Idx, x, inputIdx )
        end

        buffer[ bufIdx ] = accumulator
        inputIdx      += floor( Int, ( kernel.𝜙Idx + decimation - 1 ) / interpolation )
        kernel.𝜙Idx    = nextphase( kernel.𝜙Idx, kernel.ratio )
    end

    kernel.inputDeficit = inputIdx - xLen
    self.history        = shiftin!( history, x )

    return bufIdx
end

function filt{Th,Tx}( self::FIRFilter{FIRRational{Th}}, x::Vector{Tx} )
    kernel         = self.kernel
    xLen           = length( x )
    bufLen         = outputlength( self, xLen )
    buffer         = Array{promote_type(Th,Tx)}( bufLen )
    samplesWritten = filt!( buffer, self, x )

    samplesWritten == bufLen || resize!( buffer, samplesWritten)

    return buffer
end

#==============================================================================#
#                      ___  ____ ____ _ _  _ ____ ___ ____                     #
#                      |  \ |___ |    | |\/| |__|  |  |___                     #
#                      |__/ |___ |___ | |  | |  |  |  |___                     #
#==============================================================================#

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRDecimator{Th}}, x::Vector{Tx} )
    kernel = self.kernel
    xLen   = length( x )

    if xLen < kernel.inputDeficit
        self.history = shiftin!( history, x )
        kernel.inputDeficit -= xLen
        return Tx[]
    end

    outLen              = outputlength( self, xLen )
    history::Vector{Tx} = self.history
    inputIdx            = kernel.inputDeficit
    yIdx                = 0

    while inputIdx <= xLen
        accumulator = zero( Tb )
        yIdx       += 1

        if inputIdx < kernel.hLen
            accumulator = unsafedot( kernel.h, history, x, inputIdx )
        else
            accumulator = unsafedot( kernel.h, x, inputIdx )
        end

        buffer[ yIdx ] = accumulator
        inputIdx      += kernel.decimation
    end

    kernel.inputDeficit = inputIdx - xLen
    self.history        = shiftin!( history, x )

    return yIdx
end

function filt{Th,Tx}( self::FIRFilter{FIRDecimator{Th}}, x::Vector{Tx} )
    kernel = self.kernel
    xLen   = length( x )
    Tb     = promote_type( Th, Tx)

    if xLen < kernel.inputDeficit
        history::Vector{Tx} = self.history
        self.history       = shiftin!( history, x )
        kernel.inputDeficit -= xLen
        return Tb[]
    end

    outLen = outputlength( self, xLen )
    buffer = Array{Tb}( outLen )
    filt!( buffer, self, x )

    return buffer
end

#==============================================================================#
#        ____ ____ ___      ____ ____ ____ ____ _  _ ___  _    ____ ____       #
#        |__| |__/ |__]     |__/ |___ [__  |__| |\/| |__] |    |___ |__/       #
#        |  | |  \ |__] .   |  \ |___ ___] |  | |  | |    |___ |___ |  \       #
#==============================================================================#

# Updates FIRArbitrary state. See Section 7.5.1 in [1].
#   [1] uses a phase accumilator that increments by Δ (N𝜙/rate)
function update( kernel::FIRArbitrary )
    kernel.𝜙Accumulator += kernel.Δ

    if kernel.𝜙Accumulator > kernel.N𝜙
        kernel.xIdx        += floor( Int, (kernel.𝜙Accumulator-1) / kernel.N𝜙 )
        kernel.𝜙Accumulator = mod( (kernel.𝜙Accumulator-1), kernel.N𝜙 ) + 1
    end

    kernel.𝜙Idx = floor( Int, kernel.𝜙Accumulator )
    kernel.α    = kernel.𝜙Accumulator - kernel.𝜙Idx
end

# Generates a vector of filter taps for an arbitrary phase index.
function tapsforphase!{T}( buffer::Vector{T}, kernel::FIRArbitrary{T}, phase::Real )
    0 <= phase <= kernel.N𝜙 + 1         || error( "phase must be >= 0 and <= N𝜙+1" )
    length( buffer ) >= kernel.tapsPer𝜙 || error( "buffer is too small" )

    (α, 𝜙Idx) = modf( phase )
    𝜙Idx      = Int( 𝜙Idx )

    for tapIdx in 1:kernel.tapsPer𝜙
        buffer[tapIdx] = kernel.pfb[tapIdx,𝜙Idx] + α*kernel.dpfb[tapIdx,𝜙Idx]
    end
    buffer
end

tapsforphase{T}( kernel::FIRArbitrary{T}, phase::Real ) = tapsforphase!( Array{T}(kernel.tapsPer𝜙), kernel, phase )

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRArbitrary{Th}}, x::Vector{Tx} )
    kernel              = self.kernel
    pfb                 = kernel.pfb
    dpfb                = kernel.dpfb
    xLen                = length( x )
    bufIdx              = 0
    history::Vector{Tx} = self.history
    # TODO: Remove when arb and farrow filters are rock-solid.
    # db_vec_phi          = Array{Float64}(bufLen)
    # db_vec_xidx         = Array{Int}(bufLen)

    # Do we have enough input samples to produce one or more output samples?
    if xLen < kernel.inputDeficit
        self.history = shiftin!( history, x )
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    # Skip over input samples that are not needed to produce output results.
    # We do this by seting inputIdx to inputDeficit which was calculated in the previous run.
    # InputDeficit is set to 1 when instantiation the FIRArbitrary kernel, that way the first
    #   input always produces an output.
    kernel.xIdx = kernel.inputDeficit

    while kernel.xIdx <= xLen
        # TODO: Remove when arb and farrow filters are rock-solid.
        # db_vec_xidx[bufIdx] = kernel.xIdx
        # db_vec_phi[bufIdx]  = kernel.𝜙Idx + kernel.α
        bufIdx += 1

        if kernel.xIdx < kernel.tapsPer𝜙
            yLower = unsafedot( pfb,  kernel.𝜙Idx, history, x, kernel.xIdx )
            yUpper = unsafedot( dpfb, kernel.𝜙Idx, history, x, kernel.xIdx )
        else
            yLower = unsafedot( pfb,  kernel.𝜙Idx, x, kernel.xIdx )
            yUpper = unsafedot( dpfb, kernel.𝜙Idx, x, kernel.xIdx )
        end
        buffer[bufIdx] = yLower + yUpper * kernel.α
        update( kernel )
    end

    kernel.inputDeficit = kernel.xIdx - xLen
    self.history        = shiftin!( history, x )

    # TODO: Remove when arb and farrow filters are rock-solid.
    # resize!( db_vec_phi, length(buffer) )
    # resize!( db_vec_xidx, length(buffer) )
    # return buffer, db_vec_xidx, db_vec_phi
    return bufIdx
end

function filt{Th,Tx}( self::FIRFilter{FIRArbitrary{Th}}, x::Vector{Tx} )
    bufLen         = outputlength( self, length(x) )
    buffer         = Array{promote_type(Th,Tx)}( bufLen )
    samplesWritten = filt!( buffer, self, x )

    samplesWritten == bufLen || resize!( buffer, samplesWritten)

    return buffer
end

#==============================================================================#
#              ____ ____ ____ ____ ____ _ _ _    ____ _ _    ___               #
#              |___ |__| |__/ |__/ |  | | | |    |___ | |     |                #
#              |    |  | |  \ |  \ |__| |_|_|    |    | |___  |                #
#==============================================================================#

# Generates a vector of filter taps for an arbitray (non-integer) phase index using polynomials
function tapsforphase!{T}( buffer::Vector{T}, kernel::FIRFarrow{T}, phase::Real )
    0 <= phase <= kernel.N𝜙 + 1         || error( "phase must be >= 0 and <= N𝜙+1" )
    length( buffer ) >= kernel.tapsPer𝜙 || error( "buffer is too small" )

    for tapIdx in 1:kernel.tapsPer𝜙
        buffer[tapIdx] = polyval( kernel.pnfb[tapIdx], phase  )
    end

    return buffer
end

tapsforphase{T}( kernel::FIRFarrow{T}, phase::Real ) = tapsforphase!( Array{T}(kernel.tapsPer𝜙), kernel, phase )

# Updates farrow filter state.
# Generates new taps.
function update( kernel::FIRFarrow )
    kernel.𝜙Idx += kernel.Δ

    if kernel.𝜙Idx > kernel.N𝜙
        kernel.xIdx += floor( Int, (kernel.𝜙Idx-1) / kernel.N𝜙 )
        kernel.𝜙Idx  = mod( (kernel.𝜙Idx-1), kernel.N𝜙 ) + 1
    end

    # tapsforphase!( kernel.currentTaps, kernel, kernel.𝜙Idx ) # TODO: why does this produce worse results than below?
    for tapIdx in 1:kernel.tapsPer𝜙
        @inbounds kernel.currentTaps[tapIdx] = polyval( kernel.pnfb[tapIdx], kernel.𝜙Idx )
    end
end

function filt!{Tb,Th,Tx}( buffer::Vector{Tb}, self::FIRFilter{FIRFarrow{Th}}, x::Vector{Tx} )
    kernel              = self.kernel
    xLen                = length( x )
    bufIdx              = 0
    history::Vector{Tx} = self.history
    # TODO: Remove when arb and farrow filters are rock-solid.
    # db_vec_phi          = Array{Float64}(bufLen)
    # db_vec_xidx         = Array{Int}(bufLen)

    # Do we have enough input samples to produce one or more output samples?
    if xLen < kernel.inputDeficit
        self.history = shiftin!( history, x )
        kernel.inputDeficit -= xLen
        return bufIdx
    end

    # Skip over input samples that are not needed to produce output results.
    kernel.xIdx = kernel.inputDeficit

    while kernel.xIdx <= xLen
        bufIdx        += 1
        # TODO: Remove when arb and farrow filters are rock-solid.
        # db_vec_xidx[bufIdx] = kernel.xIdx
        # db_vec_phi[bufIdx]  = kernel.𝜙Idx
        if kernel.xIdx < kernel.tapsPer𝜙
            y = unsafedot( kernel.currentTaps, history, x, kernel.xIdx )
        else
            y = unsafedot( kernel.currentTaps, x, kernel.xIdx )
        end
        buffer[bufIdx] = y
        update( kernel )
    end

    kernel.inputDeficit = kernel.xIdx - xLen
    self.history        = shiftin!( history, x )

    # TODO: Remove when arb and farrow filters are rock-solid.
    # resize!( db_vec_phi, length(buffer) )
    # resize!( db_vec_xidx, length(buffer) )
    # return buffer, db_vec_xidx, db_vec_phi
    return bufIdx
end

function filt{Th,Tx}( self::FIRFilter{FIRFarrow{Th}}, x::Vector{Tx} )
    bufLen         = outputlength( self, length(x) )
    buffer         = Array{promote_type(Th,Tx)}( bufLen )
    samplesWritten = filt!( buffer, self, x )

    samplesWritten == bufLen || resize!( buffer, samplesWritten)

    return buffer
end

#==============================================================================#
#              ____ _  _ ____ _  _ _  _ ____ _    _ ___  ____ ____             #
#              |    |__| |__| |\ | |\ | |___ |    |   /  |___ |__/             #
#              |___ |  | |  | | \| | \| |___ |___ |  /__ |___ |  \             #
#==============================================================================#
include( "Channelizer.jl" )

#==============================================================================#
#       ____ ___ ____ ___ ____ _    ____ ____ ____    ____ _ _    ___          #
#       [__   |  |__|  |  |___ |    |___ [__  [__     |___ | |     |           #
#       ___]  |  |  |  |  |___ |___ |___ ___] ___]    |    | |___  |           #
#==============================================================================#

# Single-rate, decimation, interpolation, and rational resampling.
function filt( h::Vector, x::Vector, ratio::Rational = 1//1 )
    self = FIRFilter( h, ratio )
    filt( self, x )
end

# Arbitrary resampling with polyphase interpolation and two neighbor lnear interpolation.
function filt( h::Vector, x::Vector, rate::AbstractFloat, N𝜙::Integer = 32 )
    self = FIRFilter( h, rate, N𝜙 )
    filt( self, x )
end

# Arbitrary resampling with polyphase interpolation and polynomial generated intra-phase taps.
function filt( h::Vector, x::Vector, rate::AbstractFloat, N𝜙::Integer, polyorder::Integer )
    self = FIRFilter( h, rate, N𝜙, polyorder )
    filt( self, x )
end

#==============================================================================#
#               ____ ____ ____ ____ ____ ____ _  _ ____ ____ ____              #
#               |__/ |___ |___ |___ |__/ |___ |\ | |    |___ [__               #
#               |  \ |___ |    |___ |  \ |___ | \| |___ |___ ___]              #
#==============================================================================#

# [1] F.J. Harris, *Multirate Signal Processing for Communication Systems*. Prentice Hall, 2004
# [2] Dick, C.; Harris, F., "Options for arbitrary resamplers in FPGA-based modulators," Signals, Systems and Computers, 2004. Conference Record of the Thirty-Eighth Asilomar Conference on , vol.1, no., pp.777,781 Vol.1, 7-10 Nov. 2004
# [3] Kim, S.C.; Plishker, W.L.; Bhattacharyya, S.S., "An efficient GPU implementation of an arbitrary resampling polyphase channelizer," Design and Architectures for Signal and Image Processing (DASIP), 2013 Conference on, vol., no., pp.231,238, 8-10 Oct. 2013
# [4] Horridge, J.P.; Frazer, Gordon J., "Accurate arbitrary resampling with exact delay for radar applications," Radar, 2008 International Conference on , vol., no., pp.123,127, 2-5 Sept. 2008
# [5] Blok, M., "Fractional delay filter design for sample rate conversion," Computer Science and Information Systems (FedCSIS), 2012 Federated Conference on , vol., no., pp.701,706, 9-12 Sept. 2012


text:
# TODO: DESIGN PROPOSALS

"""
torus(x, a, b)
# Description
- Takes x ∈ ℝ and outputs torus(x) ∈ [a, b] in a periodic way.
- If a particle is moving to the right then it will pop from b to the point a
# Arguments: x, a, b
- `x`: (scalar). Current location of particle
- `a`: (scalar). left endpoint of interval
- `b`: (scalar). right endpoint of interval
# Output
-  `y`: (scalar). a value in the interval [a,b]
"""
torus(x::Number, a::Number, b::Number) = (((x-a)/(b-a))%1 - 0.5 * (sign((x-a)/(b-a)) - 1) )*(b-a) + a

"""
torus(x, a, b)
# Description
- Takes x ∈ ℝⁿ and outputs torus(x) ∈ ∏[aⁿ, bⁿ] in a periodic way.
- If a particle is moving to the right then it will pop from one part of the box to the oher
# Arguments: x, a, b
- `x`: (array). Current location of particle
- `a`: (array). left endpoint of tensor product interval
- `b`: (array). right endpoint of tensor product interval
# Output
-  `y`: (array). a value in the interval ∏[aⁿ, bⁿ]
"""
function torus(x::AbstractArray, a::AbstractArray, b::AbstractArray)
    N = length(x)
    y = zeros(N)
    for i in 1:N
        y[i] = torus(x[i], a[i], b[i])
    end
    return y
end

"""
closure_proprosal(covariance = Σ; left_bounds = [], right_bounds = []))
# Description
- Constructs a proposal for the Monte Carlo method.
# Arguments
- `covariance`: (vector) proposal parameter
# Keyword Arguments
- `left_bounds`: (array), left bounds for parameters
- `right_bounds`: (array), right bounds for parameters
# Output:
- `proposal`: (function), a function that outputs the proposal parameter
"""
function closure_proposal(Σ; left_bounds = [], right_bounds = [])
    perturbation = MvNormal(Σ)
    function proposal(C)
        proposal_C = C + rand(perturbation)
        # limit ranges for the parameters
        if isempty(left_bounds)
            return proposal_C
        else
            return torus(proposal_𝑪, left_bounds, right_bounds)
        end
        return proposal_C
    end
    return proposal
end

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
module DummyPool

# dummy allocator that passes through any requests, calling into the GC if that fails.

import ..@pool_timeit, ..actual_alloc, ..actual_free

using CUDAdrv

init() = return

deinit() = @assert isempty(allocated) "Cannot deinitialize memory pool with outstanding allocations"

const allocated = Set{Mem.Buffer}()

function alloc(sz)
    buf = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc(false)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc(true)" GC.gc(true)
        end

        @pool_timeit "$phase.1 alloc" begin
            buf = actual_alloc(sz)
        end
        buf === nothing || break
    end

    if buf !== nothing
        push!(allocated, buf)
        return buf
    end
end

function free(buf)
    delete!(allocated, buf)
    actual_free(buf)
    return
end

used_memory() = isempty(allocated) ? 0 : sum(sizeof, allocated)

cached_memory() = 0

dump() = return

end