text:
# '''
# Athesis.jl
# A Thramework for Efficient Simulations in Software
# A flexible and portable toolbox for scientific computing on CPUs and GPUs
# '''

using Test
using CuArrays
using CUDAdrv, CUDAnative

include("reduce.jl")

"""
Get model input:
grid, model type, physcial and numerical parameters
"""
function model_input(GPU)
    model_data = []
    grid_data  = []
    time_data  = []
    parameters = []
    err        = ' '

    println("Reading model input...")

    # Grid boundaries or south-west (lower-left) corner
    x0 = 0.0
    y0 = 0.0

    # grid dimensions
    nx = 101
    ny = 101

    # grid step size
    dx = 10.0
    dy = dx

    # Model components
    advection = false
    diffusion = true
    friction  = false

    # permeability?
    # 1.0e-06 m/s for some kind of sand
    # isotropic
    K = 1.0e-06

    # Courant number
    C = 0.98

    # time step size
    maxdx = max(dx, dy)
    dx2   = maxdx * maxdx
    dt    = C * 0.5 * dx2 / K

    # Maximum number of time steps
    maxsteps = 10000

    # Initialize hydraulic head
    h0 = 15.0 * ones(Float64, nx, ny)

    # boundary conditions
    bcleft = 10.0 * ones(Float64, nx)
    bcright = 20.0 * ones(Float64, nx)

    # solver parameters
    conv_limit = 1.0e-06

    # time integtation
    time_integration = "forward_euler"

    # Hardware/computing information
    # GPU = false

    grid_data  = Dict("x0" => x0, "y0" => y0, "nx" => nx, "ny" => ny, "dx" => dx, "dy" => dy)
    model_data = Dict("advection" => advection, "diffusion" => diffusion, "friction" => friction, "h0" => h0, "bcleft" => bcleft, "bcright" => bcright)
    time_data  = Dict("dt" => dt, "maxsteps" => maxsteps, "Courant" => C, "time_integration" => time_integration)
    parameters = Dict("K" => K, "conv_limit" => conv_limit, "GPU" => GPU)

    return model_data, grid_data, time_data, parameters, err
end

"""Setup the grid: structured / unstructured"""
function grid_setup(grid_data)

    println("Setting up the grid ...")
    grid = []
    err  = ' '

    x0 = grid_data["x0"]
    y0 = grid_data["y0"]
    nx = grid_data["nx"]
    ny = grid_data["ny"]
    dx = grid_data["dx"]
    dy = grid_data["dy"]

    xc = zeros(Float64, nx, ny)
    yc = zeros(Float64, nx, ny)
    for j = 1:ny
        for i = 1:nx
            xc[i,j] = x0 + (i - 0.5) * dx
            yc[i,j] = y0 + (j - 0.5) * dy
        end
    end

    grid = Dict("nx" => nx, "ny" => ny, "dx" => dx, "dy" => dy, "xc" => xc, "yc" => yc)

    return grid, err
end

"""
Setup the model:
discretization in space and time, initializations
"""
function model_setup(model_data, time_data, grid, parameters)

    println("Setting up the model ...")

    model = []
    pars  = []
    err   = ' '

    advection = model_data["advection"]
    diffusion = model_data["diffusion"]
    friction  = model_data["friction"]
    h0        = model_data["h0"]
    bcleft    = model_data["bcleft"]
    bcright   = model_data["bcright"]

    bcs       = [bcleft, bcright]

    dt               = time_data["dt"]
    maxsteps         = time_data["maxsteps"]
    C                = time_data["Courant"]
    # time_integration = time_data["time_integration"]

    conv_limit = parameters["conv_limit"]
    GPU        = parameters["GPU"]

    nx = grid["nx"]
    ny = grid["ny"]
    # dx = grid["dx"]
    # dy = grid["dy"]

    # Set initial h
    h = h0

    # Special GPU initializations and copying
    if GPU
        # copy to GPU
        h = CuArray(h)
        hnew = copy(h)
        hchange = 2.0 * conv_limit * CuArrays.ones(eltype(h), nx * ny)
        state = [h, hnew, hchange]

        # GPU stuff
        ths = (16, 16)
        bls = (Int(ceil(nx / ths[1])), Int(ceil(ny / ths[2])))
        pars  = [ths, bls]
    else
        hnew = copy(h)
        hchange = 2.0 * conv_limit * ones(eltype(h), nx * ny)
    end

    state = [h, hnew, hchange]
    model = Dict("state" => state, "bcs" => bcs)

    # For now the simple diffusion model, with forward euler (explicit) time integration
    # TO DO

    return model, pars, err
end

"""Kernel for change in 2D, result in u"""
function calculateChange_gpu!(u::AbstractArray, v::AbstractArray, w::AbstractArray)
    ix = blockDim().x * (blockIdx().x - 1) + threadIdx().x;
    iy = blockDim().y * (blockIdx().y - 1) + threadIdx().y;
    idx = size(u, 1) * (iy - 1) + ix
    if ix < size(u, 1) + 1 && iy < size(u, 2) + 1
        @inbounds w[idx] = abs(v[ix,iy] - u[ix,iy])
    end

    return nothing
end

"""Kernel for diffusion in 2D"""
function calculateDiffusion_gpu!(unew::AbstractArray, u::AbstractArray, dx, dy, dt, K)
    ix = blockDim().x * (blockIdx().x - 1) + threadIdx().x;
    iy = blockDim().y * (blockIdx().y - 1) + threadIdx().y;

    if ix > 1 && iy > 1 && ix < size(unew, 1) && iy < size(unew, 2)
        @inbounds F = K * ( (u[ix + 1,iy] + u[ix - 1,iy] - 2.0 * u[ix,iy]) / (dx * dx) + (u[ix,iy + 1] + u[ix,iy - 1] - 2.0 * u[ix,iy]) / (dy * dy) )
        @inbounds unew[ix,iy] = u[ix,iy] + dt * F
    end

    return nothing
end

function update_boundary_conditions!(model)
    bcs = []
    err = ' '

    # Boundary conditions
    bcleft    = model["bcs"][1]
    bcright   = model["bcs"][2]

    # Unpack present state
    h         = model["state"][1]
    hnew      = model["state"][2]
    hchange   = model["state"][3]

    # set boundaries
    h[:,1] = bcleft
    h[:,end] = bcright

    state = [h,hnew, hchange]
    model["state"] = state

    return model, err
end

"""Time loop: integrate the solution in time"""
function timeloop!(model, pars, grid, time_data, parameters)

    println("Starting time loop ...")
    # solution = []
    err      = ' '

    dt       = time_data["dt"]
    maxsteps = time_data["maxsteps"]
    # state    = model["state"]
    GPU      = parameters["GPU"]

    # niter = 0
    conv = false

    for n = 1:maxsteps
        if GPU
            model, conv, err = single_timestep_GPU!(n, dt, model, pars, grid, parameters, conv)
        else
            model, conv, err = single_timestep!(n, dt, model, grid, parameters, conv)
        end

        if conv
            break
        end
    end

    # println("convergence in ", niter, " iterations (", maxchange, ")")
    # display(hnew)
    display(model["state"][2])

    return model, err
end

"""Perform a single time step on the GPU"""
function single_timestep_GPU!(niter, dt, model, pars, grid, parameters, conv)
    err = ' '
    maxchange = 0.0

    # set boundaries for present time steps
    model, err = update_boundary_conditions!(model)

    # Unpack present state
    h       = model["state"][1]
    hnew    = model["state"][2]
    hchange = model["state"][3]

    conv_limit = parameters["conv_limit"]
    K          = parameters["K"]

    # GPU parameters
    ths = pars[1]
    bls = pars[2]

    # Unpack grid
    dx = grid["dx"]
    dy = grid["dy"]
    # nx = grid["nx"]
    # ny = grid["ny"]

    # niter = niter + 1

    # run CUDA kernels
    @cuda threads = ths blocks = bls calculateDiffusion_gpu!(hnew, h, dx, dy, dt, K)
    @cuda threads = ths blocks = bls calculateChange_gpu!(h, hnew, hchange)

    # reduce to find maximum change in h
    gpu_reduce(max, hchange, hchange)
    maxchange = hchange[1]

    # set no-flow boundaries
    hnew[1,:] = hnew[2,:]
    hnew[end,:] = hnew[end - 1,:]

    # swap
    htemp = hnew
    hnew  = h
    h     = htemp

    state = [h, hnew, hchange]

    if maxchange < conv_limit
        conv = true
    end
    # maxchange = 0.0

    if (mod(niter, 10) == 0)
        println("iteration: ", niter, ". Maximum change: ", maxchange)
    end
    model["state"] = state

    return model, conv, err
end

"""Perform a single time step"""
function single_timestep!(niter, dt, model, grid, parameters, conv)
    err = ' '
    maxchange = 0.0

    # set boundaries for present time steps
    model, err = update_boundary_conditions!(model)

    # Unpack present state
    h       = model["state"][1]
    hnew    = model["state"][2]
    hchange = model["state"][3]

    conv_limit = parameters["conv_limit"]
    K = parameters["K"]

    dx = grid["dx"]
    dy = grid["dy"]
    nx = grid["nx"]
    ny = grid["ny"]
    # niter = niter + 1

    # Compute the spatial operator
    for j = 2:ny - 1
        for i = 2:nx - 1
            @inbounds F = K * ( (h[i + 1,j] + h[i - 1,j] - 2.0 * h[i,j]) / (dx * dx) + (h[i,j + 1] + h[i,j - 1] - 2.0 * h[i,j]) / (dy * dy) )
            @inbounds hnew[i,j] = h[i,j] + dt * F

            change = abs(h[i,j] - hnew[i,j])
            # replace by max(..)?
            if change > maxchange
                maxchange = change
            end
        end
    end

    # set no-flow boundaries
    hnew[1,:] = hnew[2,:]
    hnew[end,:] = hnew[end - 1,:]

    # swap
    htemp = hnew
    hnew  = h
    h     = htemp

    model["state"] = [h, hnew, hchange]

    if maxchange < conv_limit
        conv = true
    end
    # maxchange = 0.0
    if (mod(niter, 10) == 0)
        println("iteration: ", niter, ". Maximum change: ", maxchange)
    end

    # solution = state

    return model, conv, err
end

"""Finalize the program"""
function finalization!(model, grid)

    println("Finalizing model ...")
    err = ' '

    model = nothing
    grid  = nothing
    return err
end

"""Main program"""
function athesis(; GPU::Bool=false)

    println("Running Athesis ...")

    # Model input
    model_data, grid_data, time_data, parameters, err = model_input(GPU)

    # Grid setup
    grid, err = grid_setup(grid_data)

    # Model setup
    model, pars, err = model_setup(model_data, time_data, grid, parameters)

    # Start of the time loop
    solution, err = timeloop!(model, pars, grid, time_data, parameters)

    # Finalization
    err = finalization!(model, grid)

    return nothing

end

function time_athesis(; GPU::Bool=false)

    @time athesis(GPU=GPU)

end

time_athesis(GPU=true)


text:
using juliaPackage, Test

@testset "powell" begin

    x, f = 10, x->x^4+x^3+x^2+sin(x)
    @test powell(f, x) == (-0.5480747372690095, -0.29506124329088085)

    x, f = -24.5, x->x^2
    @test powell(f, x) == (0.0, 0.0)

end

ext: jl

max_stars_repo_name: Haloxxx/juliaPackage

max_issues_repo_name: Haloxxx/juliaPackage

max_forks_repo_name: Haloxxx/juliaPackage


text:
"""
get_executable(rulenode::RuleNode, grammar::Grammar)

Returns the executable julia expression represented in the expression tree with root rulenode.  The returned expression can be evaluated using eval().
"""
function get_executable(rulenode::RuleNode, grammar::Grammar)
	root = (rulenode._val != nothing) ?
		rulenode._val : deepcopy(grammar.rules[rulenode.ind])
	if !grammar.isterminal[rulenode.ind] # not terminal
		root,j = _get_executable(root, rulenode, grammar)
	end
	return root
end

function _get_executable(expr::Expr, rulenode::RuleNode, grammar::Grammar, j=0)
	for (k,arg) in enumerate(expr.args)
		if isa(arg, Expr)
			expr.args[k],j = _get_executable(arg, rulenode, grammar, j)
		elseif haskey(grammar.bytype, arg)
			child = rulenode.children[j+=1]
			expr.args[k] = (child._val != nothing) ?
			child._val : deepcopy(grammar.rules[child.ind])
			if !isterminal(grammar, child)
				expr.args[k],_ = _get_executable(expr.args[k], child, grammar, 0)
			end
		end
	end
	return expr, j
end

function _get_executable(typ::Symbol, rulenode::RuleNode, grammar::Grammar, j=0)
	retval = typ
		if haskey(grammar.bytype, typ)
			child = rulenode.children[1]
			retval = (child._val !== nothing) ?
				child._val : deepcopy(grammar.rules[child.ind])
			if !grammar.isterminal[child.ind]
				retval,_ = _get_executable(retval, child, grammar, 0)
			end
		end
	retval, j
end

"""
Core.eval(rulenode::RuleNode, grammar::Grammar)

Evaluate the expression tree with root rulenode.
"""
Core.eval(rulenode::RuleNode, grammar::Grammar) = Core.eval(Main, get_executable(rulenode, grammar))
Core.eval(grammar::Grammar, index::Int) = Core.eval(Main, grammar.rules[index].args[2])

"""
Core.eval(tab::SymbolTable, ex::Expr)

Evaluate the expression ex using symbol table tab
"""
Core.eval(tab::SymbolTable, ex) = interpret(tab, ex)

function Base.display(rulenode::RuleNode, grammar::Grammar)
	root = get_executable(rulenode, grammar)
	if isa(root, Expr)
	    walk_tree(root)
	else
	    root
	end
    end


text:
using Roots
using Statistics
using FileIO
using Parameters
include("random_matrices.jl")
include("network_sims_gpu.jl")

N = 1000
g = 10.0f0
eta = 0.0f0
s = 0.0f0

J_seeds = collect(5000:5019)
h0_seed = 6000

h0 = generate_nv(N, Float32, seed=h0_seed)

wait = 200.0f0
t_bound = 2000.0f0
dt = 0.01f0

startstep = round(Int, wait / dt) + 1

# Global array of the average temporal variances of network units
vanilla_network_avg_vars = Vector{Float32}(undef, length(J_seeds))

starttime = time()
for (i, J_seed) in enumerate(J_seeds)
    J = generate_J(N, eta, seed=J_seed)
    A = sim_activity_gpu(J, g, h0; s=s, dt=dt, t_bound=t_bound,
            minimum_norm=1e-90)
    if size(A, 2) < startstep
        avg_var = 0
    else
        avg_var = mean(A[:,startstep:end].^2)
    end
    vanilla_network_avg_vars[i] = avg_var
end
println("Ensemble integration done in ", round(time() - starttime, digits=3),
        " seconds")

target_avg_var = mean(vanilla_network_avg_vars)

out = Dict{String,Any}()
@pack! out = N, g, eta, s, J_seeds, h0_seed, wait, t_bound, dt, target_avg_var
save("vanilla_target_avg_var.jld2", out)


text:
using Test
import Pluto
import Pluto: update_save_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell

import UUIDs: UUID, uuid1

function get_unique_short_id()
    string(uuid1())[1:8]
end

function stringify_keys(d::Dict)
    Dict(string(k) => stringify_keys(v) for (k, v) in d)
end
stringify_keys(x::Any) = x

import Pluto.Firebasey

function await_with_timeout(check::Function, timeout::Real=60.0, interval::Real=.05)
    starttime = time()
    while !check()
        sleep(interval)
        if time() - starttime >= timeout
            error("Timeout after $(timeout) seconds")
        end
    end
end

@testset "Communication protocol" begin

    @testset "Functionality sweep" begin
        u = [uuid1() for _ in 1:100]

        buffer = IOBuffer()
        client = ClientSession(:buffery, buffer)

        🍭 = ServerSession()
        # 🍭.connected_clients[client.id] = client

        notebook = Notebook([
            Cell(
                u[1],
                ""
            ),
        ])
        🍭.notebooks[notebook.notebook_id] = notebook
        update_save_run!(🍭, notebook, notebook.cells)
        client.connected_notebook = notebook
        read(buffer)

        local_state = Pluto.notebook_to_js(notebook)

        function send(type, body, metadata = Dict(:notebook_id => string(notebook.notebook_id)))
            request_id = get_unique_short_id()
            Pluto.process_ws_message(🍭, Dict(
                "type" => string(type),
                "client_id" => string(client.id),
                "request_id" => request_id,
                "body" => body,
                metadata...
            ) |> stringify_keys, client.stream)
        end

        function send_new_state(new_state)
            patches::Array{Dict} = Firebasey.diff(new_state, local_state)
            # @info "patches" patches
            send(:update_notebook, Dict("updates" => patches))

            new_state
        end

        # function update_local_notebook(mutate_fn::Function)
        #     mutable_notebook = deepcopy(notebook)

        #     mutate_fn(mutable_notebook)

        #     new_state = Pluto.notebook_to_js(mutable_notebook)
        #     send_new_state(new_state)
        # end

        last_position = position(buffer)
        function wait_for_updates(process=true)
            # @info "Waiting for updates"
            await_with_timeout() do
                position(buffer) != last_position
            end
            seek(buffer, last_position)
            response = Pluto.unpack(buffer)
            last_position = position(buffer)

            if process
                if response["type"] == "notebook_diff"
                    message = response["message"]
                    patches = [Base.convert(Firebasey.JSONPatch, update) for update in message["patches"]]
                    # @show patches
                    for patch in patches
                        Firebasey.applypatch!(local_state, patch)
                    end
                end
            end

            # @info "Response received" response

            local_state
        end

        # function patch_and_check(mutate_fn::Function)
        #     desired = update_local_notebook(mutate_fn)
        #     result = wait_for_updates()
        #     desired == result
        # end

        @test_nowarn send(:connect, Dict())
        wait_for_updates()

        @test_nowarn send_new_state(local_state)
        wait_for_updates(false)

        read(buffer)

        #=
        We would also like to test:
        - add cell
        - set code and run
        - fold cell
        - move cell
        - delete cell

        - run multiple cells
        - move cells
        - set bond

        - move notebook file
        - search for docs
        - show more items of an array =#

        send(:shutdown_notebook, Dict("keep_in_session" => false))

        @test_nowarn await_with_timeout() do
            !haskey(🍭.notebooks, notebook.notebook_id)
        end
    end

    @testset "Docs" begin
        @test occursin("square root", Pluto.PlutoRunner.doc_fetcher("sqrt")[1])
        @test occursin("square root", Pluto.PlutoRunner.doc_fetcher("Base.sqrt")[1])
        @test occursin("No documentation found", Pluto.PlutoRunner.doc_fetcher("Base.findmeta")[1])
    end
end


text:
using Localize

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

replace!(Base.Math, Docs.@var(sin), Docs.docstr("正弦函数"), Tuple{Number})


text:
using CSV
using DataFrames
using Plots

### Import data

# Generate a temporary file path
tmp = tempname()

# Download data in the temporary file path from the UCI website
download("https://archive.ics.uci.edu/ml/machine-learning-databases/00236/seeds_dataset.txt",
        tmp)

# Read the seeds dataset
# Dropmissing used to drop pseudo-missing values in the raw dataset
seeds = dropmissing(CSV.read(tmp; header=0, delim='\t'))

# Name the variables (measures of wheat kernels = grains)
# 1 = area (A)
# 2 = perimeter (P)
# 3 = compactness (C = 4*pi*A/P^2)
# 4 = length of kernel (Fr: graine)
# 5 = width of kernel
# 6 = asymmetry coefficient
# 7 = length of kernel groove (Fr: rainure)
# 8 = cultivar (1, 2 or 3) : variety of wheat
rename!(seeds,
  [:Column1 => :area, :Column2 => :perimeter,
  :Column3 => :compactness, :Column4 => :kernel_length,
  :Column5 => :kernel_width, :Column6 => :asymmetry,
  :Column7 => :kernel_groove, :Column8 => :cultivar]
  )

### Scatter plots

## Length of kernel groove as a function of kernel length

# Basic plot
scatter(seeds.kernel_length, seeds.kernel_groove)

# Color palette
pal = (
    C1=RGB(230/255,159/255,0/255),
    C2=RGB(86/255,190/255,233/255),
    C3=RGB(0/255,158/255,115/255),
    lr=RGB(204/255,121/255,167/255)
)

# Make it pretty
plotA = scatter(seeds.kernel_length, seeds.kernel_groove,
                color=:grey, alpha=0.8,
                label="",
                frame=:box, dpi=120,
                smooth=true, linecolor=pal.lr, linewidth=3)
xaxis!(xlabel="Kernel length (cm)", xlims=(4.8,6.8))
yaxis!(ylabel="Groove length (cm)", ylims=(4.4,6.7))

# Add 1:1 line
x=minimum(seeds.kernel_length):0.1:maximum(seeds.kernel_length)
plot!(x, x,
     label="1:1",
     linestyle=:dot, linecolor=:black, linewidth=1,
     legend=:right, foreground_color_legend=nothing)

### Density plots
using StatsPlots

# Density of Kernel length
plotB = density(seeds.kernel_length, groups=seeds.cultivar,
                linewidth=2,
                xlabel="Kernel or groove length (cm)", ylabel="Density",
                xlims=(4, 7.5),
                label=["C1" "C2" "C3"],
                linecolor=[pal.C1 pal.C2 pal.C3],
                legend=:topright, foreground_color_legend=nothing,
                frame=:box, dpi=120)

# Density of Groove length
density!(seeds.kernel_groove, groups=seeds.cultivar,
         linewidth=2, linestyle=:dot,
         linecolor=[pal.C1 pal.C2 pal.C3], lab="")

### Violin plots
cultivar_int = string.(convert.(Int64, seeds.cultivar))

plotC = violin(cultivar_int, seeds.kernel_length,
       color=[pal.C1, pal.C2, pal.C3],
       label="",
       xlabel="Cultivar",
       ylabel="Kernel length (cm)",
       frame=:box, dpi=120)
boxplot!(cultivar_int, seeds.kernel_length,
        color=:grey, alpha=0.5,
        label="")
dotplot!(cultivar_int, seeds.kernel_length,
        color=:grey, alpha=0.5, markersize=1,
        label="")

### Plots together
l = @layout [a{0.6h}; b c]
plot(plotA, plotB, plotC, layout=l,
    legend = false,
    title = ["(1)" "(2)" "(3)"], titleloc = :right, titlefont = font(8))
savefig(joinpath("tutorials", "Julia_tutorial", "fig.png"))

### Resources
# http://docs.juliaplots.org/latest/

### To go further
# Gadfly.jl: https://github.com/GiovineItalia/Gadfly.jl
# PyPlot.jl: https://github.com/JuliaPy/PyPlot.jl
# StatsPlots.jl: https://github.com/JuliaPlots/StatsPlots.jl


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

using DataFrames
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame);

mean_ma = mean(df[!, :MedianAgeMarriage])
df[!, :MedianAgeMarriage_s] =
    convert(Vector{Float64},  (df[!, :MedianAgeMarriage]) .-
        mean_ma)/std(df[!, :MedianAgeMarriage]);
df.Divorce_SE = df[:, "Divorce SE"]

stan14_1 = "
  data {
    int N;
    vector[N] A;
    vector[N] R;
    vector[N] Dobs;
    vector[N] Dsd;
  }
  parameters {
    real a;
    real ba;
    real br;
    real<lower=0> sigma;
    vector[N] Dest;
  }
  model {
    vector[N] mu;
    // priors
    target += normal_lpdf(a | 0, 10);
    target += normal_lpdf(ba | 0, 10);
    target += normal_lpdf(br | 0, 10);
    target += cauchy_lpdf(sigma | 0, 2.5);

    // linear model
    mu = a + ba * A + br * R;

    // likelihood
    target += normal_lpdf(Dest | mu, sigma);

    // prior for estimates
    target += normal_lpdf(Dobs | Dest, Dsd);
  }
  generated quantities {
    vector[N] log_lik;
    {
      vector[N] mu;
      mu = a + ba * A + br * R;
      for(i in 1:N) log_lik[i] = normal_lpdf(Dest[i] | mu[i], sigma);
    }
  }
";

m14_1s = SampleModel("m14.1s", stan14_1)

m14_1_data = Dict(
    "N" => size(df, 1),
    "A" => df[!, :MedianAgeMarriage],
    "R" => df[!, :Marriage],
    "Dobs" => df[!, :Divorce],
    "Dsd" => df[!, :Divorce_SE]
)

rc14_1s = stan_sample(m14_1s, data=m14_1_data)

if success(rc14_1s)
    chns14_1s = read_samples(m14_1s, :dataframe)
    chns14_1s |> display
end


text:
# This file was generated, do not modify it. # hide
rand(rng, Normal())

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
using BaseBenchmarks
using BenchmarkTools
if VERSION >= v"0.7.0-DEV.2954"
    using Distributed
end

addprocs(1)

# re-tune the entire suite

BaseBenchmarks.loadall!(tune = false)
warmup(BaseBenchmarks.SUITE)
tune!(BaseBenchmarks.SUITE; verbose = true)
BenchmarkTools.save(BaseBenchmarks.PARAMS_PATH, params(BaseBenchmarks.SUITE))


text:
using SHA

bytes(n::Integer, l::Int) = collect(UInt8, (n >> 8i) & 0xFF for i in l-1:-1:0)

function decodebase58(bc::String, l::Int)
    digits = "123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz"
    num = big(0)
    for c in bc
        num = num * 58 + search(digits, c) - 1
    end
    return bytes(num, l)
end

function checkbcaddress(addr::String)
    if !(25 ≤ length(addr) ≤ 35) return false end
    bcbytes = decodebase58(addr, 25)
    sha = sha256(sha256(bcbytes[1:end-4]))
    return bcbytes[end-3:end] == sha[1:4]
end

const addresses = Dict(
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62i" => true,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62j" => false,
    "1Q1pE5vPGEEMqRcVRMbtBK842Y6Pzo6nK9" => true,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62X" => true,
    "1ANNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62i" => false,
    "1A Na15ZQXAZUgFiqJ2i7Z2DPU2J6hW62i" => false,
    "BZbvjr" => false,
    "i55j" => false,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62!" => false,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62iz" => false,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62izz" => false,
    "1Q1pE5vPGEEMqRcVRMbtBK842Y6Pzo6nJ9" => false,
    "1AGNa15ZQXAZUgFiqJ2i7Z2DPU2J6hW62I" => false)

for (addr, corr) in addresses
    println("Address: $addr\nExpected: $corr\tChecked: ", checkbcaddress(addr))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
###########
# Helpers #
###########

unref(value::T) where T = [value]
unref(node::Node, ::Type{T}) where T<:Rule = [node]
unref(node::Node, ::Type{ReferencedRule}) = node.children
unref(node::Node) = unref(node, node.ruleType)

function make_node(rule, value, first, last, children::Array)
  #println("make_node: $(rule.action)($rule, $value, first=$first, last=$last, children=$children)")
  return rule.action(rule, value, first, last, children)
end

###############
# ParserCache #
###############

abstract type ParserCache end

struct StandardCache <: ParserCache
  values::Dict{AbstractString, Node}

  function StandardCache()
    return new(Dict{AbstractString, Any}())
  end
end

###################
# parse (generic) #
###################

"""
    parse(grammar, text; cache::ParserCache=nothing, start=:start)
parses `text` according to `grammar` to yield a tuple consisting of Abstract Syntax Tree, final matched position and error (ast, pos, error) depending on outcome.
`start` specifies the symbol associated to the rule at the top of the AST. Specifying a `cache` different than `nothing` allows to reuse previous work, whenever the same rule is evaluated at the same position again.
"""
function parse(grammar::Grammar, text::AbstractString; cache=nothing, start=:start)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)

  rule = grammar.rules[start]
  (ast, pos, error) = parse(grammar, rule, text, 1, cache)

  if pos < length(text) + 1
    sequence = text[utf8ind[pos]:end]
    if length(sequence) > 15
      sequence = text[utf8ind[pos]:utf8ind[pos+15]]*"..."
    end
    error = Meta.ParseError("Entire string did not match at pos: $pos ($sequence)")
  end

  return (ast, pos, error)
end

"""
    parse(grammar, rule, text, pos, cache)
parses `text` according to `rule` within `grammar` starting with position `pos`. Specifying a `cache` different than `nothing` allows to reuse previous work, whenever the same rule has been matched at the same position before. If no `cache` is specified or no match in `cache` is found `parse` resorts to `parse_newcachekey`, because in an non-existent cache every cachekey is new.
"""
function parse(grammar::Grammar, rule::Rule, text::AbstractString, pos::Int, cache::Nothing)
  return parse_newcachekey(grammar, rule, text, pos, cache)
end

function parse(grammar::Grammar, rule::Rule, text::AbstractString, pos::Int, cache::StandardCache)
  cachekey::AbstractString = "$(object_id(rule))$pos"
  if haskey(cache.values, cachekey)
    # lookup cachekey
    cachedresult = cache.values[cachekey]
    (node, pos, error) = (cachedresult, cachedresult.last, nothing)
  else
    # parse new cachekey
    (node, pos, error) = parse_newcachekey(grammar, rule, text, pos, cache)

    # store new cachekey
    if node !== nothing
      cache.values[cachekey] = node
    end
  end

  return (node, pos, error)
end

function parse(grammar::Grammar, symbol::Symbol, text::AbstractString, pos::Int, cache::Union{StandardCache,Nothing})
  parse(grammar, grammar.rules[symbol], text, pos, cache)
end

##############################################
# parse_newcachekey (specific for each rule) #
##############################################

"""
    parse_newcachekey(grammar, rule, text, pos, cacheforsubnodes)
parses `text` according to `rule` within `grammar` starting with position `pos` without trying to lookup the complete match in the specified `cache`. Matches of children of the current `rule` along the `text` to parse will however be looked-up in `cache`.
"""
function parse_newcachekey(grammar::Grammar, rule::ReferencedRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  refrule = grammar.rules[rule.symbol]

  firstPos = pos
  (childNode, pos, error) = parse(grammar, refrule, text, pos, cache)

  if childNode !== nothing
    node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, [childNode])
    return (node, pos, error)
  else
    return (nothing, pos, error)
  end
end

function parse_newcachekey(grammar::Grammar, rule::OrRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  # Try branches in order (left to right). The first branch to match will be marked
  # as a success. If no branches match, then return an error.
  firstPos = pos
  for branch in rule.values
    (child, pos, error) = parse(grammar, branch, text, pos, cache)

    if error == nothing
      node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, unref(child))
      return (node, pos, error)
    end
  end

  # give error
  return (nothing, pos, Meta.ParseError("No match (OrRule) at pos: $pos"))
end

function parse_newcachekey(grammar::Grammar, rule::AndRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  firstPos = pos;

  # All items in sequence must match, otherwise give an error
  value = Any[]
  for item in rule.values
    (child, pos, error) = parse(grammar, item, text, pos, cache)

    # check for error
    if error !== nothing
      return (nothing, firstPos, error)
    end

    if child !== nothing
      append!(value, unref(child))
    end
  end

  node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, value)
  return (node, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::Terminal, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  local size::Int = length(rule.value)

  if startswith(text[utf8ind[pos]:end],rule.value)
    node = make_node(rule, text[utf8ind[pos]:utf8ind[pos+size-1]], pos, pos+size, [])
    return (node, pos+size, nothing)
  end

  len = min(pos+size-1, length(text))
  return (nothing, pos, Meta.ParseError("'$(text[utf8ind[pos]:utf8ind[len]])' does not match '$(rule.value)'. At pos: $pos"))
end

function parse_newcachekey(grammar::Grammar, rule::OneOrMoreRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  firstPos = pos
  (child, pos, error) = parse(grammar, rule.value, text, pos, cache)

  # make sure there is at least one
  if child === nothing
    return (nothing, pos, Meta.ParseError("No match (OneOrMoreRule) at pos: $pos"))
  end

  # and continue making matches for as long as we can
  children = unref(child)
  while error == nothing
    (child, pos, error) = parse(grammar, rule.value, text, pos, cache)

    if error === nothing && child !== nothing
      children = [children;collect(unref(child))]
    end
  end

  node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, children)
  return (node, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::ZeroOrMoreRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  firstPos::Int = pos
  children::Array = Any[]

  error = nothing
  prepos = -1
  while error == nothing && pos != prepos
    prepos = pos

    (child, pos, error) = parse(grammar, rule.value, text, pos, cache)

    if error === nothing && child !== nothing
      append!(children, unref(child))
    end
  end

  if length(children) > 0
    node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, children)
  else
    node = nothing
  end

  return (node, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::RegexRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  firstPos = pos

  # use regex match
  pattern = Regex("^$(rule.value.pattern)")
  if occursin(pattern, text[utf8ind[firstPos]:end])
    value = match(pattern, text[utf8ind[firstPos]:end])

    if length(value.match) == 0
      # this means that we didn't match, but the regex was optional, so we don't want to give an
      # error
      return (nothing, firstPos, nothing)
    else
      pos += length(value.match)
      node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, [])

      return (node, pos, nothing)
    end
  else
    return (nothing, firstPos, Meta.ParseError("Could not match RegEx at pos: $pos"))
  end
end

function parse_newcachekey(grammar::Grammar, rule::OptionalRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  (child, pos, error) = parse(grammar, rule.value, text, pos, cache)
  firstPos = pos

  if child !== nothing
    node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, unref(child))
    return (node, pos, error)
  end

  # no error, but we also don't move the position or return a valid node
  return (nothing, firstPos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::ListRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  firstPos = pos

  # number of occurances
  count = 0

  error = nothing
  children = Any[]

  # continue making matches for as long as we can
  while error === nothing
    (child, pos, error) = parse(grammar, rule.entry, text, pos, cache)

    if child !== nothing
      append!(children, unref(child))
      (dchild, pos, error) = parse(grammar, rule.delim, text, pos, cache)
    else
      break
    end

    count += 1
  end

  if count < rule.min
    return (nothing, pos, Meta.ParseError("No match (ListRule) at pos: $pos"))
  end

  node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, children)
  return (node, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::SuppressRule, text::AbstractString, pos::Int, cache)
  # use rule contained in the SuppressRule to parse, but don't return anything
  (_, pos, error) = parse_newcachekey(grammar, rule.value, text, pos, cache)
  return (nothing, pos, error)
end

function parse_newcachekey(grammar::Grammar, rule::LookAheadRule, text::AbstractString, pos::Int, cache)
    (_, newPos, error) = parse_newcachekey(grammar, rule.value, text, pos, cache)
    if error !== nothing
        return (nothing, newPos, error)
    else
        return (nothing, pos, nothing)
    end
end

function parse_newcachekey(grammar::Grammar, rule::NotRule, text::AbstractString, pos::Int, cache)
  # try to parse rule
  (child, newpos, error) = parse(grammar, rule.entry, text, pos, cache)

  # if we match, it's an error
  if error == nothing
    error = Meta.ParseError("No match (NotRule) at pos: $pos")
  else
    # otherwise, return a success
    error = nothing
  end

  return (nothing, pos, error)
end

function parse_newcachekey(grammar::Grammar, rule::EmptyRule, text::AbstractString, pos::Int, cache)
  # need to explicitely call rule's action because nothing is consumed
  if rule.action != nothing
    rule.action(rule, "", pos, pos, [])
  end

  return (nothing, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::EndOfFileRule, text::AbstractString, pos::Int, cache)
  # need to explicitely call rule's action because nothing is consumed
  if pos == length(text)
    #rule.action(rule, value, first, last, children)
    rule.action(rule, "", length(text), length(text), [])
  end

  return (nothing, pos, nothing)
end

function parse_newcachekey(grammar::Grammar, rule::IntegerRule, text::AbstractString, pos::Int, cache)
  utf8ind = collect(eachindex(text)); push!(utf8ind,ncodeunits(text)+1)
  #rexpr = r"^[-+]?[0-9]+([eE][-+]?[0-9]+)?"
  # Julia treats anything with 'e' to be a float, so for now follow suit
  rexpr = r"^[-+]?[0-9]+"
  firstPos = pos

  # use regex match
  if occursin(rexpr, text[utf8ind[firstPos]:end])
    value = match(rexpr, text[utf8ind[firstPos]:end])

    if length(value.match) != 0
      pos += length(value.match)
      node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, [])

      return (node, pos, nothing)
    end
  else
    return (nothing, firstPos, Meta.ParseError("Could not match IntegerRule at pos: $pos"))
  end
end

function parse_newcachekey(grammar::Grammar, rule::FloatRule, text::AbstractString, pos::Int, cache)
  rexpr = r"^[-+]?[0-9]*\.[0-9]+([eE][-+]?[0-9]+)?"
  firstPos = pos

  # use regex match
  if occursin(rexpr, text[utf8ind[firstPos]:end])
    value = match(rexpr, text[utf8ind[firstPos]:end])

    if length(value.match) != 0
      pos += length(value.match)
      node = make_node(rule, text[utf8ind[firstPos]:utf8ind[pos-1]], firstPos, pos, [])

      return (node, pos, nothing)
    end
  else
    return (nothing, firstPos, Meta.ParseError("Could not match FloatRule at pos: $pos"))
  end
end


text:
layers(::Type{WorldClim{Climate}}) = (:tmin, :tmax, :tavg, :prec, :srad, :wind, :vapr)

"""
    getraster(T::Type{WorldClim{Climate}}, [layer::Union{Symbol,Tuple}]; month=1:12, res::String="10m") => Vector{String}
    getraster(T::Type{WorldClim{Climate}}, layer::Symbol, month::Integer, res::String)

Download WorldClim weather data, choosing `layer` from $(layers(WorldClim{Climate})),
and `res` from $(resolutions(WorldClim{Climate})), and months from `1:12`.

Without a layer argument, all layers will be downloaded, and a tuple of paths is returned.
By default all months are downloaded , but can also be downloaded individually.
If the data is already downloaded the path will be returned.
"""
function getraster(T::Type{WorldClim{Climate}}, layer; month=1:12, res::String=defres(T))
    getraster(T, layer, month, res)
end
function getraster(T::Type{WorldClim{Climate}}, layers::Tuple, month, res::String)
    map(l -> getraster(T, l, month, res), layers)
end
function getraster(T::Type{WorldClim{Climate}}, layer::Symbol, month::AbstractArray, res::String)
    getraster.(T, layer, month, Ref(res))
end
function getraster(T::Type{WorldClim{Climate}}, layer::Symbol, month::Integer, res::String)
    _check_layer(T, layer)
    _check_res(T, res)
    raster_path = rasterpath(T, layer; res, month)
    if !isfile(raster_path)
        zip_path = zippath(T, layer; res, month)
        _maybe_download(zipurl(T, layer; res, month), zip_path)
        zf = ZipFile.Reader(zip_path)
        mkpath(dirname(raster_path))
        raster_name = rastername(T, layer; res, month)
        write(raster_path, read(_zipfile_to_read(raster_name, zf)))
        close(zf)
    end
    return raster_path
end

# Climate layers don't get their own folder
rasterpath(T::Type{<:WorldClim{Climate}}, layer; res, month) =
    joinpath(_rasterpath(T, layer), rastername(T, layer; res, month))
_rasterpath(T::Type{<:WorldClim{Climate}}, layer) = joinpath(rasterpath(T), string(layer))
rastername(T::Type{<:WorldClim{Climate}}, layer; res, month) =
    "wc2.1_$(res)_$(layer)_$(_pad2(month)).tif"
zipname(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    "wc2.1_$(res)_$(layer).zip"
zipurl(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    joinpath(WORLDCLIM_URI, "base", zipname(T, layer; res, month))
zippath(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    joinpath(rasterpath(T), "zips", zipname(T, layer; res, month))

_pad2(month) = lpad(month, 2, '0')


text:
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    kalmanFilter(data, ssmodel)

Description:
Apply Kalman filter to observed data.
Measurement Equation:
    y_{t} = H_{t} β_{t} + A z_{t} + e_{t} .
Transition Equation:
    β_{t} = μ + F β_{t-1} + v_{t};
    e_{t} ~ i.i.d.N(0,R);
    v_{t} ~ i.i.d.N(0,Q);
    z_{t} ~ i.i.d.N(0,Z);
    E(e_t v_s') = 0.

Inputs:
- data      = observed data
- H         = measurement eq. state coef. matrix
- A         = measurement eq. exogenous coef. matrix
- F         = state eq. companion matrix
- μ         = state eq. intercept term
- R         = covariance matrix on measurement disturbance
- Q         = covariance matrix on state disturbance
- Z         = covariance matrix on predetermined var vector
"""
function kalmanFilter(data_y, ssmodel)

    @unpack H, A, F, μ, R, Q, Z = ssmodel

    # Save number of observations
    num_obs = size(data_y)[1]

    # Empty filtered data matrices
    data_filtered_y = similar(data_y)
    data_filtered_β = zeros(num_obs, size(Q)[1])

    # Create empty lists for P_{t}, P_{t|t-1}
    Ptt = Any[]
    Pttlag = Any[]

    # Initialize β_pred and P_pred
    β_pred_laglag = inv(I - F) * μ
    P_pred_laglag = ones(size(Q)[1], size(Q)[1])

    for t = 1:num_obs
        # Save current obs of y
        y = data_y[t, :]

        # Prediction
        β_pred_lag = μ + F * β_pred_laglag
        P_pred_lag = F * P_pred_laglag * transpose(F) + Q
        y_pred_lag = H * β_pred_lag
        η_pred_lag = y - y_pred_lag
        f_pred_lag = H * P_pred_lag * transpose(H) + R

        # Save P_{t|t-1}
        push!(Pttlag, P_pred_lag)

        # Updating
        K = P_pred_lag * transpose(H) * pinv(f_pred_lag)
        β_pred = β_pred_lag + K * η_pred_lag
        P_pred = P_pred_lag - K * H * P_pred_lag

        # Save P_{t|t}
        push!(Ptt, P_pred)

        # Save data
        data_filtered_y[t, :] = y_pred_lag
        data_filtered_β[t, :] = β_pred

        # Lag the predictions
        β_pred_laglag = β_pred
        P_pred_laglag = P_pred
    end

    # Returned filtered series
    # for obs variable and state
    return data_filtered_y, data_filtered_β, Pttlag, Ptt
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    kalmanSmoother(data, ssmodel)

Description:
Apply Kalman smoother to observed data.
Measurement Equation:
    y_{t} = H_{t} β_{t} + A z_{t} + e_{t}.
Transition Equation:
    β_{t} = μ + F β_{t-1} + v_{t};
    e_{t} ~ i.i.d.N(0,R);
    v_{t} ~ i.i.d.N(0,Q);
    z_{t} ~ i.i.d.N(0,Z);
    E(e_t v_s') = 0.

Inputs:
- data      = observed data
- H         = measurement eq. state coef. matrix
- A         = measurement eq. exogenous coef. matrix
- F         = state eq. companion matrix
- μ         = state eq. intercept term
- R         = covariance matrix on measurement disturbance
- Q         = covariance matrix on state disturbance
- Z         = covariance matrix on predetermined var vector
"""
function kalmanSmoother(data_y, ssmodel)

    @unpack H, A, F, μ, R, Q, Z = ssmodel

    # Save number of observations
    num_obs = size(data_y)[1]

    # Empty filtered data matrices
    data_smoothed_y = similar(data_y)
    data_smoothed_β = zeros(num_obs, size(Q)[1])

    # Create empty list for P_{t|T}
    PtT = Any[]

    # Run Kalman filter
    data_filtered_y, data_filtered_β, Pttlag, Ptt = kalmanFilter(data_y, ssmodel)

    # Initialize β_{t+1|T} (β_{T|T})
    βtflagT = data_filtered_β[end]
    data_smoothed_β[end] = βtflagT

    # Initialize P_{t+1|T} (P_{T|T})
    Ptflag_T = Ptt[end]
    push!(PtT, Ptflag_T)

    # Initialize y_{t|T} (y_{T|T})
    data_smoothed_y[end] = data_filtered_y[end]

    # Run Kalman smoother
    for i = 1:(num_obs-1)

        # Retrieve β_{t|t}
        βtt = data_filtered_β[end-i]

        # Compute β_{t|T} using β_{t+1|T}, β_{t|t}, P_{t|t}, and P_{t+1|t}
        βtT = βtt +
              Ptt[end-i] * transpose(F) * inv(Pttlag[end-i+1]) *
              (βtflagT - F * βtt - μ)

        # Store β_{t+1|T} in smoothed data
        data_smoothed_β[end-i] = βtT

        # Set β_{t|T} as new β_{t+1|T} for next iteration
        βtflagT = βtT

        # Compute P_{t|T} using P_{t|t}, P_{t+1|t}, and P_{t+1|T}
        Pt_T = Ptt[end-i] +
               Ptt[end-i] * transpose(F) * inv(Pttlag[end-i+1]) *
               (Ptflag_T - Pttlag[end-i+1]) *
               transpose(Ptt[end-i] * transpose(F) * inv(Pttlag[end-i+1]))

        # Store P_{t|T}
        push!(PtT, Pt_T)

        # Set P_{t|T} as new P_{t+1|T} for next iteration
        Ptflag_T = Pt_T

        # Generate y_{t|T} (smoothed obs.)
        ytT = H * βtT

        # Store smoothed obs.
        data_smoothed_y[end-i] = ytT
    end

    # Flip P_{t|T} list
    PtT = reverse(PtT)

    # Returned filtered series
    # for obs variable and state
    return data_smoothed_y, data_smoothed_β, PtT
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    KNFactorSampler(data_y, ssmodel)

Description:
Draw a sample series of dynamic factor from conditional distribution in Ch 8, Kim & Nelson (1999).
Measurement Equation:
    y_{t} = H_{t} β_{t} + A z_{t} + e_{t}.
Transition Equation:
    β_{t} = μ + F β_{t-1} + v_{t};
    e_{t} ~ i.i.d.N(0,R);
    v_{t} ~ i.i.d.N(0,Q);
    z_{t} ~ i.i.d.N(0,Z);
    E(e_t v_s') = 0.

Inputs:
- data      = observed data
- H         = measurement eq. state coef. matrix
- A         = measurement eq. exogenous coef. matrix
- F         = state eq. companion matrix
- μ         = state eq. intercept term
- R         = covariance matrix on measurement disturbance
- Q         = covariance matrix on state disturbance
- Z         = covariance matrix on predetermined var vector
"""
function KNFactorSampler(data_y, ssmodel)

    @unpack H, A, F, μ, R, Q, Z = ssmodel

    # Run Kalman filter
    data_filtered_y, data_filtered_β, Pttlag, Ptt = kalmanFilter(data_y, ssmodel)

    # Create placeholders for factor distr.
    # mean vector and covariance matrix for all t
    β_t_mean = Any[]
    β_t_var = Any[]

    # Record number of time periods
    T = size(data_y)[1]

    # Create empty vector for factor realizations
    β_realized = similar(data_filtered_β)

    # Initialize β_realized
    push!(β_t_mean, data_filtered_β[T, :])
    push!(β_t_var, Ptt[T])
    β_realized[T, :] = sim_MvNormal_alt(β_t_mean[1], β_t_var[1])

    # Generate `β_t_mean` and `β_t_var`
    # for all time periods
    if isposdef(Q) == false
        ## IF Q IS SINGULAR

        # Determine number of rows
        # of state factor used
        num_use_rows = 0
        for i = 1:size(Q)[1]
            if Q[i, i] != 0
                num_use_rows += 1
            end
        end

        # Create modified F and Q matrices
        F_star = F[1:num_use_rows, :]
        Q_star = Q[1:num_use_rows, 1:num_use_rows]

        # Iteratively generate conditional draws
        # of state vector
        for j = 1:(T-1)

            # β_{t|t,β*_{t+1}}
            β_t_mean_temp = data_filtered_β[T-j, :]
            +Ptt[T-j] * transpose(F_star) * inv(F_star * Ptt[T-j] * transpose(F_star) + Q_star) * (β_realized[T+1-j, :][1:num_use_rows] - μ[1:num_use_rows] - F_star * data_filtered_β[T-j, :])
            push!(β_t_mean, β_t_mean_temp)

            # P_{t|t,β*_{t+1}}
            β_t_var_temp = Ptt[T-j] - Ptt[T-j] * transpose(F_star) * inv(F_star * Ptt[T-j] * transpose(F_star) + Q_star) * F_star * Ptt[T-j]
            push!(β_t_var, β_t_var_temp)

            # Draw new β_t
            β_realized[T-j, :] = sim_MvNormal_alt(β_t_mean[1+j], β_t_var[1+j])
        end
    else
        ## IF Q IS NOT SINGULAR (redundant, but potentially faster)
        for j = 1:(T-1)

            # β_{t|t,β_{t+1}}
            β_t_mean_temp = data_filtered_β[T-j, :]
            +Ptt[T-j] * transpose(F) * inv(F * Ptt[T-j] * transpose(F) + Q) * (β_realized[T+1-j, :] - μ - F * data_filtered_β[T-j, :])
            push!(β_t_mean, β_t_mean_temp)

            # P_{t|t,β_{t+1}}
            β_t_var_temp = Ptt[T-j] - Ptt[T-j] * transpose(F) * inv(F * Ptt[T-j] * transpose(F) + Q) * F * Ptt[T-j]
            push!(β_t_var, β_t_var_temp)

            # Draw new β_t
            β_realized[T-j, :] = sim_MvNormal_alt(β_t_mean[1+j], β_t_var[1+j])
        end
    end

    # Return sampled factor series
    # fot t = 1,...,T
    return β_realized
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    linearRegressionSampler(Y,X)

Description:
Estimate β and σ^2 in Y = Xβ + e, e ~ N(0,σ^2 I_T).
Generate samples of β and σ^2.

Inputs:
- Y     = Dependent data matrix
- X     = Independent data matrix
"""
function linearRegressionSampler(Y, X)

    # Save number of obs
    T = size(X)[1]

    # Initialize σ2
    σ2 = 1.0

    ##################################
    ##################################
    # Generate new β

    ## Prior parameters in N(β0,Σ0)
    β0 = zeros(size(X)[2])
    Σ0 = Matrix(I, size(β0)[1], size(β0)[1]) .* 1000.0

    ## Posterior parameters in N(β1,Σ1)
    β1 = inv(inv(Σ0) + inv(σ2) * transpose(X) * X) * (inv(Σ0) * β0 + inv(σ2) * transpose(X) * Y)
    β1 = vec(β1)
    Σ1 = inv(inv(Σ0) + inv(σ2) * transpose(X) * X)
    Σ1 = Hermitian(Σ1)

    ## Generate new β
    β = sim_MvNormal(β1, Σ1)

    ##################################
    ##################################
    # Update σ2

    ## Prior parameters in IG(ν0/2, δ0/2)
    ν0 = 0.002
    δ0 = 0.002

    ## Posterior parameters in IG(ν1/2, δ1/2)
    ν1 = ν0 + T
    δ1 = δ0 + norm(Y - X * β)^2

    ## Generate new σ2
    σ2 = rand(InverseGamma(ν1 / 2, δ1 / 2))

    ##################################
    ##################################
    # Return parameters
    return β, σ2
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    linearRegressionSamplerRestrictedVariance(Y, X, σ2)

Description:
Estimate β and σ^2 in Y = Xβ + e, e ~ N(0,σ^2 I_T),
where σ^2 is restricted to some chosen value.
Generate samples of β and σ^2.

Inputs:
- Y     = Dependent data matrix
- X     = Independent data matrix
- σ2    = Restricted error variance
"""
function linearRegressionSamplerRestrictedVariance(Y, X, σ2)

    # Save number of obs
    T = size(X)[1]

    ##################################
    ##################################
    # Generate new β

    ## Prior parameters in N(β0,Σ0)
    β0 = zeros(size(X)[2])
    Σ0 = Matrix(I, size(β0)[1], size(β0)[1]) .* 1000.0

    ## Posterior parameters in N(β1,Σ1)
    β1 = inv(inv(Σ0) + inv(σ2) * transpose(X) * X) * (inv(Σ0) * β0 + inv(σ2) * transpose(X) * Y)
    β1 = vec(β1)
    Σ1 = inv(inv(Σ0) + inv(σ2) * transpose(X) * X)
    Σ1 = Hermitian(Σ1)

    ## Generate new β
    β = sim_MvNormal(β1, Σ1)

    ##################################
    ##################################
    # Return parameters
    return β
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """

    autocorrErrorLinearRegressionSampler(Y, X, error_lag_num)

Description:
Estimate β, σ^2, and ϕ in Y = Xβ + e, e = Eϕ + ν, ν_t ~ i.i.d.N(0,σ^2).
Generate samples of β, σ^2, and ϕ.
Procedure descried in Section 7.4.2 in Kim & Nelson.

Inputs:
- Y             = Dependent data matrix
- X             = Independent data matrix
- error_lag_num = Number of lags in the disturbance DGP
"""

function autocorrErrorLinearRegressionSampler(Y, X, ϕold, σ2old, error_lag_num)

    # Save number of obs
    T = size(X)[1]

    # Initialize σ2
    ϕ = ϕold
    σ2 = σ2old

    ##################################
    ##################################
    # Generate β

    ## Prior parameters in N(b0,A0)
    b0 = zeros(size(X)[2])
    A0 = Matrix(I, size(b0)[1], size(b0)[1]) .* 1000.0

    ## Generate X^⋆
    X_star = similar(X[(1+length(ϕ)):end, :])
    for i = 1:size(X)[2] # iterate over variables in X
        x_temp = X[:, i]
        for p = 1:length(ϕ) # iterate over lag params in ϕ
            x_temp[(1+length(ϕ)):end, :] = x_temp[(1+length(ϕ)):end, :] - ϕ[p] .* lag(x_temp, p)[(1+length(ϕ)):end, :]
        end
        x_temp = x_temp[(1+length(ϕ)):end, :]
        X_star[:, i] = x_temp
    end

    ## Generate Y^⋆
    Y_star = similar(Y[(1+length(ϕ)):end, :])
    y_temp = Y
    for p = 1:length(ϕ)
        y_temp = Y
        y_temp = y_temp[(1+length(ϕ)):end, :] - ϕ[p] .* lag(y_temp, p)[(1+length(ϕ)):end, :]
    end
    Y_star = y_temp

    ## Posterior parameters in N(b1,A1)
    b1 = inv(inv(A0) + inv(σ2) * transpose(X_star) * X_star) * (inv(A0) * b0 + inv(σ2) * transpose(X_star) * Y_star)
    b1 = vec(b1)
    A1 = inv(inv(A0) + inv(σ2) * transpose(X_star) * X_star)
    A1 = Hermitian(A1)

    ## Draw new β
    β = sim_MvNormal(b1, A1)

    ##################################
    ##################################
    # Generate ϕ

    ## Prior parameters in N(c0,B0)
    c0 = zeros(size(ϕ)[1])
    B0 = Matrix(I, size(c0)[1], size(c0)[1]) .* 1000.0

    ## Generate e^⋆
    e_star = Y - X * β

    ## Generate E^⋆
    E_star = zeros(T, length(ϕ))
    E_star = E_star[(1+length(ϕ)):end, :]
    for i = 1:(size(E_star)[2]) # iterate over variables in X
        e_temp = lag(e_star, i)
        E_star[:, i] = e_temp[(1+length(ϕ)):end]
    end
    e_star = e_star[(1+length(ϕ)):end]

    ## Posterior parameters in N(c1,B1)
    c1 = inv(inv(B0) + inv(σ2) * transpose(E_star) * E_star) * (inv(B0) * c0 + inv(σ2) * transpose(E_star) * e_star)
    c1 = vec(c1)
    B1 = inv(inv(B0) + inv(σ2) * transpose(E_star) * E_star)
    B1 = Hermitian(B1)

    ## Generate new ϕ
    ###=
    ind = 0
    accept = 0
    ϕ = similar(ϕold)
    while accept == 0

        ind += 1

        ## Draw ψ
        ϕ = sim_MvNormal(c1, B1)

        ## Check for stationarity
        coef = [-reverse(vec(ϕ), dims = 1); 1]                      # check stationarity
        root = roots(Polynomial(reverse(coef)))
        rootmod = abs.(root)
        accept = min(rootmod...) >= 1.01

        ## If while loop goes on for too long
        if ind > 100
            ϕ = ϕold
            coef = [-reverse(vec(ϕ), dims = 1); 1]                      # check stationarity
            root = roots(Polynomial(reverse(coef)))
            rootmod = abs.(root)
            accept = min(rootmod...) >= 1.01
        end
    end
    ##=#

    #ϕ = sim_MvNormal(c1, B1)

    ##################################
    ##################################
    # Generate σ2

    ## Prior parameters in IG(ν0/2, δ0/2)
    ν0 = 0.002
    δ0 = 0.002

    ## Posterior parameters in IG(ν1/2, δ1/2)
    ν1 = ν0 + T
    δ1 = δ0 .+ norm(Y_star - X_star * β)^2
    δ1 = δ1[1, 1]

    ## Generate new σ2
    σ2 = rand(InverseGamma(ν1 / 2, δ1 / 2))

    ##################################
    ##################################
    # Return parameters
    return β, σ2, ϕ
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################


text:
abstract type AbstractCoefficientCache{T} end

function Base.checkbounds(::Type{Bool}, u::AbstractVector, cache::AbstractCoefficientCache)
    length(u) > length(cache.inv_left_weights)+length(cache.inv_right_weights)
end

"""
    VarCoefDerivativeCoefficients

The coefficients of a variable coefficient derivative operator on a nonperiodic grid.
"""
struct VarCoefDerivativeCoefficients{T,CoefficientCache<:AbstractCoefficientCache{T},
                                     LeftWidth,RightWidth,Parallel,SourceOfCoefficients} <: AbstractDerivativeCoefficients{T}
    # coefficients defining the operator and its action
    coefficient_cache::CoefficientCache
    left_weights::SVector{LeftWidth, T}
    right_weights::SVector{RightWidth, T}
    parallel::Parallel
    # corresponding orders etc.
    derivative_order::Int
    accuracy_order  ::Int
    source_of_coeffcients::SourceOfCoefficients

    function VarCoefDerivativeCoefficients(coefficient_cache::CoefficientCache,
                                           left_weights::SVector{LeftWidth, T},
                                           right_weights::SVector{RightWidth, T},
                                           parallel::Parallel,
                                           derivative_order::Int, accuracy_order::Int,
                                           source_of_coeffcients::SourceOfCoefficients) where {T,CoefficientCache<:AbstractCoefficientCache{T},
                                                                                LeftWidth,RightWidth,Parallel,SourceOfCoefficients}
        new{T,CoefficientCache,LeftWidth,RightWidth,Parallel,SourceOfCoefficients}(
            coefficient_cache, left_weights, right_weights, parallel, derivative_order, accuracy_order, source_of_coeffcients)
    end
end

function Base.show(io::IO, coefficients::VarCoefDerivativeCoefficients)
    print(io, "Coefficients of the variable coefficient ")
    if  derivative_order(coefficients) == 1
        print(io, "1st")
    elseif  derivative_order(coefficients) == 2
        print(io, "2nd")
    elseif  derivative_order(coefficients) == 3
        print(io, "3rd")
    else
        print(io, derivative_order(coefficients), "th")
    end
    print(io, "derivative operator with order of accuracy ",
            accuracy_order(coefficients), " given in \n")
    print(io, source_of_coeffcients(coefficients))
end

"""
    mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients, u::AbstractVector, b::AbstractVector, α, β)

Compute `α*D*u + β*dest` using the coefficients `b` and store the result in `dest`.
"""
function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                u::AbstractVector, b::AbstractVector, α, β)
    @unpack coefficient_cache, parallel = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(dest) == length(b) DimensionMismatch
        @argcheck checkbounds(Bool, dest, coefficient_cache) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, coefficient_cache, u, b, α, β)
    convolve_interior_coefficients!(dest, coefficient_cache, u, b, α, β, parallel)
end

"""
    mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients, u::AbstractVector, b::AbstractVector, α)

Compute `α*D*u` using the coefficients `b` and store the result in `dest`.
"""
function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                u::AbstractVector, b::AbstractVector, α)
    @unpack coefficient_cache, parallel = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(dest) == length(b) DimensionMismatch
        @argcheck checkbounds(Bool, dest, coefficient_cache) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, coefficient_cache, u, b, α)
    convolve_interior_coefficients!(dest, coefficient_cache, u, b, α, parallel)
end

@inline function left_length(cache::AbstractCoefficientCache)
    length(cache.inv_left_weights)
end

@inline function right_length(cache::AbstractCoefficientCache)
    length(cache.inv_right_weights)
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, parallel)
    for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, ::Val{:threads})
    Threads.@threads for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, β, parallel)
    for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval + β*dest[i]
    end end
end

function convolve_interior_coefficients!(dest::AbstractVector, cache::AbstractCoefficientCache,
                                         u::AbstractVector, b::AbstractVector, α, β, ::Val{:threads})
    Threads.@threads for i in (left_length(cache)+1):(length(dest)-right_length(cache)) @inbounds begin
        retval = convolve_interior_coefficients_loopbody(i, cache, u, b)
        dest[i] = α*retval + β*dest[i]
    end end
end

abstract type AbstractVariableCoefficientDerivativeOperator{T} <: AbstractDerivativeOperator{T} end

@inline source_of_coeffcients(D::AbstractVariableCoefficientDerivativeOperator) = source_of_coeffcients(D.coefficients)

@inline function lower_bandwidth(D::AbstractVariableCoefficientDerivativeOperator)
    lower_bandwidth(D.coefficients.coefficient_cache)
end

@inline function upper_bandwidth(D::AbstractVariableCoefficientDerivativeOperator)
    upper_bandwidth((D.coefficients.coefficient_cache))
end

"""
    DissipationOperator

A dissipation operator on a nonperiodic finite difference grid.
"""
struct DissipationOperator{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid} <: AbstractVariableCoefficientDerivativeOperator{T}
    factor::T
    coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients}
    grid::Grid
    b::Vector{T}

    function DissipationOperator(factor::T,
                                 coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients},
                                 grid::Grid, b::Vector{T}) where {T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}
        @argcheck checkbounds(Bool, grid, coefficients.coefficient_cache) DimensionMismatch
        @argcheck length(grid) == length(b)
        factor < 0 && @warn("Negative dissipation strength shouldn't be used.")

        new{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}(
            factor,coefficients, grid, b)
    end
end

function Base.show(io::IO, D::DissipationOperator{T}) where {T}
    if  derivative_order(D) == 2
        print(io, "SBP 2nd derivative dissipation operator of order ")
    else
        print(io, "SBP ", derivative_order(D), "th derivative dissipation operator of order ")
    end
    print(io, accuracy_order(D), " {T=", T, ", Parallel=", typeof(D.coefficients.parallel), "} \n")
    print(io, "on a grid in [", first(grid(D)), ", ", last(grid(D)),
                "] using ", length(grid(D)), " nodes \n")
    print(io, "and coefficients given in \n")
    print(io, source_of_coeffcients(D))
end

"""
    mul!(dest::AbstractVector, D::DissipationOperator, u::AbstractVector, α, β)

Compute `α*D*u + β*dest` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                       u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α, β)
end

"""
    mul!(dest::AbstractVector, D::DissipationOperator, u::AbstractVector, α)

Compute `α*D*u` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                       u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α)
end

"""
    dissipation_operator(source_of_coefficients, order, xmin, xmax, N,
                         left_weights, right_weights, parallel=Val{:serial}())

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a weighted `order`-th derivative on a grid between `xmin` and
`xmax` with `N` grid points up to order of accuracy 2 with coefficients given
by `source_of_coefficients`.
The norm matrix is given by `left_weights` and `right_weights`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function dissipation_operator(source_of_coefficients, order, xmin, xmax, N,
                              left_weights, right_weights,
                              strength=one(xmin+xmax),
                              parallel=Val{:serial}())
    grid = construct_grid(source_of_coefficients, order, xmin, xmax, N)
    coefficients, b = dissipation_coefficients(source_of_coefficients, order, grid, left_weights, right_weights, parallel)
    DissipationOperator(strength, coefficients, grid, b)
end

"""
    dissipation_operator(source_of_coefficients, D::DerivativeOperator{T};
                         strength=one(T),
                         order::Int=accuracy_order(D),
                         parallel=D.coefficients.parallel)

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a weighted `order`-th derivative adapted to the derivative
operator `D` with coefficients given in `source_of_coefficients`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function dissipation_operator(source_of_coefficients, D::DerivativeOperator{T};
                              strength=one(T),
                              order::Int=accuracy_order(D),
                              parallel=D.coefficients.parallel) where {T}
    x = grid(D)
    dissipation_operator(source_of_coefficients, order, first(x), last(x), length(x), D.coefficients.left_weights, D.coefficients.right_weights, strength, parallel)
end

"""
    dissipation_operator(D::DerivativeOperator; kwargs...)

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a weighted `order`-th derivative adapted to the derivative
operator `D`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function dissipation_operator(D::DerivativeOperator; kwargs...)
    source = MattssonSvärdNordström2004()
    dissipation_operator(source, D; kwargs...)
end


text:
"""
    struct MethodMatchInfo

Captures the result of a `method_matches` lookup for the given call. This
info may then be used by the optimizer to inline the matches, without having
to re-consult the method table. This info is illegal on any statement that is
not a call to a generic function.
"""
struct MethodMatchInfo
    results::Union{Missing, MethodLookupResult}
end

"""
    struct UnionSplitInfo

If inference decides to partition the method search space by splitting unions,
it will issue a method lookup query for each such partition. This info indicates
that such partitioning happened and wraps the corresponding MethodMatchInfo for
each partition. This info is illegal on any statement that is not a call to a
generic function.
"""
struct UnionSplitInfo
    matches::Vector{MethodMatchInfo}
end

"""
    struct CallMeta

A simple struct that captures both the return type any any additional `info`
for a given generic call.
"""
struct CallMeta
    rt::Any
    info::Any
end

"""
    struct AbstractIterationInfo

Captures all the information for abstract iteration analysis of a single value.
Each (abstract) call to `iterate`, corresponds to one entry in `each`.
"""
struct AbstractIterationInfo
    each::Vector{CallMeta}
end

"""
    struct ApplyCallInfo

This info applies to any call of _apply_iterate(...) and captures both the
info of the actual call being applied and the info for any implicit call
to the `iterate` function. Note that it is possible for the call itself
to be yet another `_apply_iterate`, in which case the `.call` field will
be another ApplyCallInfo. This info is illegal on any statement that is
not an _apply_iterate call.
"""
struct ApplyCallInfo
    # The info for the call itself
    call::Any
    # AbstractIterationInfo for each argument, if applicable
    arginfo::Vector{Union{Nothing, AbstractIterationInfo}}
end

"""
    struct UnionSplitApplyCallInfo

Like `UnionSplitInfo`, but for `ApplyCallInfo` rather than MethodMatchInfo.
This info is illegal on any statement that is not an _apply_iterate call.
"""
struct UnionSplitApplyCallInfo
    infos::Vector{ApplyCallInfo}
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 0b19c578-762b-11eb-34b6-01e80cef1406
using Pkg, DrWatson

# ╔═╡ 3f814e9e-762b-11eb-1340-91617ca7b58a
begin
    #@quickactivate "StatisticalRethinkingStan"
    using StanSample
    using StatisticalRethinking
end

# ╔═╡ af1efcf4-779a-11eb-0d62-e3c63876a09b
md" ## Clip-09-11-21s.jl"

# ╔═╡ cb88768a-779b-11eb-1a4e-99f92ded29ca
begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	dropmissing!(df, :rgdppc_2000)
	dropmissing!(df, :rugged)
	df.log_gdp = log.(df[:, :rgdppc_2000])
	df.log_gdp_s = df.log_gdp / mean(df.log_gdp)
	df.rugged_s = df.rugged / maximum(df.rugged)
	df.cid = [df.cont_africa[i] == 1 ? 1 : 2 for i in 1:size(df, 1)]

	data = (N = size(df, 1), K = length(unique(df.cid)),
		G = df.log_gdp_s, R = df.rugged_s, cid=df.cid)

	PRECIS(df[:, [:rgdppc_2000, :log_gdp, :log_gdp_s, :rugged, :rugged_s, :cid]])
end

# ╔═╡ e37331ea-779b-11eb-1b02-1fdb61aae124
stan9_1 = "
data {
	int N;
	int K;
	vector[N] G;
	vector[N] R;
	int cid[N];
}
parameters {
	vector[K] a;
	vector[K] b;
	real<lower=0> sigma;
}
transformed parameters {
	vector[N] mu;
	for (i in 1:N)
		mu[i] = a[cid[i]] + b[cid[i]] * (R[i] - 0.215);
}
model {
	a ~ normal(1, 0.1);
	b ~ normal(0, 0.3);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
generated quantities {
	vector[N] log_lik;
	for (i in 1:N)
		log_lik[i] = normal_lpdf(G[i] | mu[i], sigma);
}
";

# ╔═╡ 3fafc7fe-762b-11eb-37da-254b6d75f5ca
begin
	m9_1s = SampleModel("m9.1s", stan9_1;
		method=StanSample.Sample(num_samples=500))

	rc9_1s = stan_sample(m9_1s; data)

	if success(rc9_1s)
		m9_1s_df = read_samples(m9_1s, :dataframe)
		params = [Symbol("a.1"), Symbol("a.2"), Symbol("b.1"), Symbol("b.2")]
		PRECIS(m9_1s_df[:, append!(params, [:sigma])])
	end
end

# ╔═╡ 15d035d6-779d-11eb-0393-e54390b408dd
begin
	dfs = read_summary(m9_1s)
	dfs[8:12, :]
end

# ╔═╡ b4d020b8-762d-11eb-167a-9fd79d948647
begin
	figs_a = trankplot(m9_1s, :a)
	plot(figs_a..., layout=(1, 2))
end

# ╔═╡ 28a050e8-7857-11eb-18c8-bd4ef6695f24
begin
	figs_b = trankplot(m9_1s, :b)
	plot(figs_b..., layout=(1, 2))
end

# ╔═╡ 369b3044-7935-11eb-3ae7-abe15e53c060
begin
	fig = trankplot(m9_1s, :sigma)
	plot(fig...)
end

# ╔═╡ ca61c622-779a-11eb-1292-3f272953867f
md" ## End of clip-09-11-21s.jl"

# ╔═╡ Cell order:
# ╠═af1efcf4-779a-11eb-0d62-e3c63876a09b
# ╠═0b19c578-762b-11eb-34b6-01e80cef1406
# ╠═3f814e9e-762b-11eb-1340-91617ca7b58a
# ╠═cb88768a-779b-11eb-1a4e-99f92ded29ca
# ╠═e37331ea-779b-11eb-1b02-1fdb61aae124
# ╠═3fafc7fe-762b-11eb-37da-254b6d75f5ca
# ╠═15d035d6-779d-11eb-0393-e54390b408dd
# ╠═b4d020b8-762d-11eb-167a-9fd79d948647
# ╠═28a050e8-7857-11eb-18c8-bd4ef6695f24
# ╠═369b3044-7935-11eb-3ae7-abe15e53c060
# ╟─ca61c622-779a-11eb-1292-3f272953867f


text:
include("resTime.jl")

resTime(spawnRateSet=30,
        CASType=:correctedSector,
        policyFile="LOWFI_mid3_maxdist_2km_nmac_150m_Discount_0.95_ActWeight_0.5_IfSample_false_PenAction_3.0_PenConflict_50.0_PenCloseness_500.0_PenNmac_7000.0_Sigma_v_5.0_alert_5.0_theta_3.0_phi_7.0_coc_3.0_2018-11-13T18\:55\:43.491",
        numEvalEposide=1,
        maxEpisodeLen=1000,
        render=false)

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using DiffEqBase, RecursiveArrayTools, Test

mutable struct VectorType{T} <: DEDataVector{T}
    x::Vector{T}
    f::T
end

mutable struct MatrixType{T,S} <: DEDataMatrix{T}
    f::S
    x::Matrix{T}
end

A = [0.0; 1.0]
B = [1.0  2; 4  3]

a = VectorType{Float64}(copy(A), 1.0)
b = MatrixType{Float64,Float64}(2.0, copy(B))

# basic methods of AbstractArray interface
@test eltype(a) == Float64 && eltype(b) == Float64

# size
@test length(a) == 2 && length(b) == 4
@test size(a) == (2,) && size(b) == (2,2)

# iteration
@test first(a) == 0.0 && first(b) == 1
@test last(a) == 1.0 && last(b) == 3
for (i, (ai, Ai)) in enumerate(zip(a, A))
    @test ai == Ai
end
for (i, (bi, Bi)) in enumerate(zip(b, B))
    @test bi == Bi
end

# indexing
@test eachindex(a) == Base.LinearIndices(A) && eachindex(b) == vec(Base.LinearIndices(B))
@test a[2] == a[end] == 1.0
@test b[3] == b[1,2] == 2 && b[:, 1] == [1; 4]

a[1] = 3; b[2,1] = 1
@test a[1] == 3.0 && b[2] == 1

a[:] = A; b[:, 1] = B[1:2]
@test a.x == A && b.x == B

# simple broadcasts
@test [1//2] .* a == [0.0; 0.5]
@test b ./ 2 ==  [0.5  1.0; 2.0  1.5]
@test b .+ a  == [1.0  2.0; 5.0  4.0]
@test_broken a .+ b  == [1.0  2.0; 5.0  4.0] # Doesn't find the largest

# similar data arrays
a2 = similar(a); b2 = similar(b, (1,4)); b3 = similar(b, Float64, (1, 4))
@test typeof(a2.x) == Vector{Float64} && a2.f == 1.0
@test typeof(b2.x) == Matrix{Float64} && b2.f == 2.0
@test typeof(b3.x) == Matrix{Float64} && b3.f == 2.0

# copy all fields of data arrays
recursivecopy!(a2, a)
recursivecopy!(b2, b)
@test a2.x == A && vec(b2.x) == vec(B)
recursivecopy!(b3, b)

# copy all fields except of field x
a2.f = 3.0; b2.f = -1; b3.f == 0
DiffEqBase.copy_fields!(a, a2)
DiffEqBase.copy_fields!(b, b2)
@test a.f == 3.0 && b.f == -1.0
DiffEqBase.copy_fields!(b, b3)

# create data array with field x replaced by new array
a3 = DiffEqBase.copy_fields([1.0; 0.0], a)
@test a3 == VectorType{Float64}([1; 0], 3.0)

# broadcast assignments
a.f = 0.0
a .= a2 .+ a3
@test a == VectorType{Float64}([1.0; 1.0], 0.0)
@test a == a2 .+ a3
@test (a2 .+ a3) isa VectorType

old_b = copy(b)
b .= b .^ 2 .+ a3
@test b == MatrixType{Int,Float64}(-1.0, [2  5; 16  9])
@test b == old_b .^ 2 .+ a3
@test (b .^ 2 .+ a3) isa MatrixType

using StaticArrays

# Test ability to use MVectors in the solvers
mutable struct SimWorkspace{T} <: DEDataVector{T}
  x::MVector{2,T}
  a::T
end
s0   = SimWorkspace{Float64}(MVector{2,Float64}(1.0,4.0),1.)
similar(s0,Float64,size(s0))
s1   = SimWorkspace{Float64}(MVector{2,Float64}(2.0,1.0),1.)
s0 .+ s1 == SimWorkspace{Float64}(MVector{2,Float64}(3.0,5.0),1.)

mutable struct SimWorkspace2{T} <: DEDataVector{T}
  x::SVector{2,T}
  a::T
end
s0   = SimWorkspace2{Float64}(SVector{2,Float64}(1.0,4.0),1.)
s1   = SimWorkspace2{Float64}(SVector{2,Float64}(2.0,1.0),1.)
s0 .+ s1 == SimWorkspace2{Float64}(SVector{2,Float64}(3.0,5.0),1.)


text:
# Required for REPL.it
include("./The_Austroporean_Cards/main.jl")

max_stars_repo_name: NarutoXY/austroporean-cards

max_issues_repo_name: NarutoXY/austroporean-cards

max_forks_repo_name: NarutoXY/austroporean-cards


text:
module FixedPointSinCosApproximations

    using LoopVectorization, DocStringExtensions
    using LoopVectorization: vifelse, SVec

    export fpsin,
        fpcos,
        fpsincos

    const VInt16 = Union{Int16,SVec{<:Any,Int16}}
    const VInt32 = Union{Int32,SVec{<:Any,Int32}}

    include("second_order.jl")
    include("third_order.jl")
    include("fourth_order.jl")
    include("fifth_order.jl")
    include("sixth_order.jl")
    include("approximations.jl")

end # module


text:
using ShaderAbstractions: Buffer, Sampler, VertexArray

# Mesh
mesh(Sphere(Point3f0(0), 1f0)) |> display
mesh(Sphere(Point3f0(0), 1f0), color=:red, ambient=Vec3f0(0.9))

tocolor(x) = RGBf0(x...)
positions = Observable(decompose(Point3f0, Sphere(Point3f0(0), 1f0)))
triangles = Observable(decompose(GLTriangleFace, Sphere(Point3f0(0), 1f0)))
uv = Observable(GeometryBasics.decompose_uv(Sphere(Point3f0(0), 1f0)))
xyz_vertex_color = Observable(tocolor.(positions[]))
texture = Observable(rand(RGBAf0, 10, 10))

pos_buff = Buffer(positions)
triangles_buff = Buffer(triangles)
vert_color_buff = Buffer(xyz_vertex_color)
uv_buff = Buffer(uv)
texture_buff = Sampler(texture)
texsampler = Makie.sampler(:viridis, rand(length(positions)))

coords = VertexArray(pos_buff, triangles_buff, color=vert_color_buff)
mesh = GeometryBasics.Mesh(coords)
GeometryBasics.coordinates(mesh);

using ShaderAbstractions: data

posmeta = Buffer(meta(data(pos_buff); color = data(vert_color_buff)))

program = disp.renderlist[1][3].vertexarray.program

struct OGLContext <: ShaderAbstractions.AbstractContext end

instance = ShaderAbstractions.VertexArray(posmeta, triangles_buff)

ShaderAbstractions.type_string(OGLContext(), Vec2f0)

p = ShaderAbstractions.Program(
    OGLContext(),
    read(loadshader("mesh.vert"), String),
    read(loadshader("mesh.frag"), String),
    instance;
)

println(p.vertex_source)

uniforms = Dict{Symbol, Any}(
    :texturecoordinates => Vec2f0(0),
    :image => nothing
)
rshader = GLMakie.GLAbstraction.gl_convert(shader, uniforms)

vbo = GLMakie.GLAbstraction.GLVertexArray(program, posmeta, triangles_buff)

m = GeometryBasics.Mesh(posmeta, triangles_buff)
disp = display(Makie.mesh(m, show_axis=false));

mesh_normals = GeometryBasics.normals(positions, triangles)
coords = meta(positions, color=xyz_vertex_color, normals=mesh_normals)
vertexcolor_mesh = GeometryBasics.Mesh(coords, triangles)
scren = mesh(vertexcolor_mesh, show_axis=false) |> display

function getter_function(io::IO, ::Fragment, sampler::Sampler, name::Symbol)
    index_t = type_string(context, sampler.values)
    sampler_t = type_string(context, sampler.colors)

    println(io, """
    in $(value_t) fragment_$(name)_index;
    uniform $(sampler_t) $(name)_texture;

    vec4 get_$(name)(){
        return texture($(name)_texture, fragment_$(name)_index);
    }
    """)
end

function getter_function(io::IO, ::Vertex, sampler::Sampler, name::Symbol)
    index_t = type_string(context, sampler.values)
    println(io, """
    in $(index_t) $(name)_index;
    out $(index_t) fragment_$(name)_index;

    vec4 get_$(name)(){
        fragment_uv = uv;
        // color gets calculated in fragment!
        return vec4(0);
    }
    """)
end

function getter_function(io::IO, ::Fragment, ::AbstractVector{T}, name) where T
    t_str = type_string(context, T)
    println(io, """
    in $(t_str) fragment_$(name);
    $(t_str) get_$(name)(){
        return fragment_$(name);
    }
    """)
end

function getter_function(io::IO, ::Vertex, ::AbstractVector{T}, name) where T
    t_str = type_string(context, T)
    println(io, """
    in $(t_str) $(name);
    out $(t_str) fragment_$(name);

    $(t_str) get_$(name)(){
        fragment_$(name) = $(name);
        return $(name);
    }
    """)
end

texsampler = Makie.sampler(rand(RGBf0, 4, 4), uv)
coords = meta(positions, color=texsampler, normals=mesh_normals)
texture_mesh = GeometryBasics.Mesh(coords, triangles)

scren = mesh(texture_mesh, show_axis=false) |> display

texsampler = Makie.sampler(:viridis, rand(length(positions)))
coords = meta(positions, color=texsampler, normals=mesh_normals)
texture_mesh = GeometryBasics.Mesh(coords, triangles)

scren = mesh(texture_mesh, show_axis=false) |> display

glsl"""

out vec2 frag_uv;
out vec3 frag_normal;
out vec3 frag_position;
out vec4 frag_color;
out vec3 frag_lightdir;

uniform mat4 projectionMatrix;
uniform mat4 viewMatrix;
uniform mat4 modelMatrix;

vec3 tovec3(vec2 v){return vec3(v, 0.0);}
vec3 tovec3(vec3 v){return v;}

vec4 tovec4(vec3 v){return vec4(v, 1.0);}
vec4 tovec4(vec4 v){return v;}

void main(){
    // get_* gets the global inputs (uniform, sampler, vertex array)
    // those functions will get inserted by the shader creation pipeline
    vec3 vertex_position = tovec3(get_position());
    frag_normal = $(normals);
    frag_uv = get_texturecoordinates();

    vec3 lightpos = vec3(20);
    vec4 position_world = modelMatrix * vec4(vertex_position, 1);
    frag_lightdir = normalize(lightpos - position_world.xyz);
    // direction to camera
    frag_position = -position_world.xyz;
    frag_uv = vec2(1.0 - frag_uv.y, frag_uv.x);
    frag_color = tovec4(get_color());
    // screen space coordinates of the vertex
    gl_Position = projectionMatrix * viewMatrix * position_world;
}
"""