text:
function _sqrt(a::Bool, n)
    if a
        return n > 0 ? 2.0 : 1.0
    else
        return 1.0
    end
end

function _napier(a::Bool, n)
    if a
        return n > 0 ? Float64(n) : 2.0
    else
        return n > 1 ? n - 1.0 : 1.0
    end
end

function _pi(a::Bool, n)
    if a
        return n > 0 ? 6.0 : 3.0
    else
        return (2.0 * n - 1.0) ^ 2.0 # exponentiation operator
    end
end

function calc(f::Function, expansions::Integer)
    a, b = true, false
    r = 0.0
    for i in expansions:-1:1
        r = f(b, i) / (f(a, i) + r)
    end
    return f(a, 0) + r
end

for (v, f) in (("√2", _sqrt), ("e", _napier), ("π", _pi))
    @printf("%3s = %f\n", v, calc(f, 1000))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Plot a real-time spectrogram
# see https://github.com/JuliaAudio/PortAudio.jl

using GR, PortAudio, SampledSignals

const N = 1024
const stream = PortAudioStream(1, 0, blocksize=N)
const buf = read(stream, N)
const fmin = 0Hz
const fmax = 10000Hz
const fs = Float32[float(f) for f in domain(fft(buf)[fmin..fmax])]

while true
    read!(stream, buf)
    plot(fs, abs.(fft(buf)[fmin..fmax]), xlim=(fs[1],fs[end]), ylim=(0,100))
end

ext: jl


text:
using ArgParse
using CGP
using YAML
include("clustering.jl")
include("datasets.jl")
include("seeds.jl")

settings = ArgParseSettings()
@add_arg_table settings begin
    "--seed"
    arg_type = Int
    default = 0
    "--expert"
    arg_type = String
    default = ""
    "--expfile"
    arg_type = String
    default = ""
    "--logfile"
    arg_type = String
    default = "stdp.log"
end

args = parse_args(settings)
srand(args["seed"])
Logging.configure(filename=args["logfile"], level=INFO)

CGP.Config.init("cfg/base.yaml")
CGP.Config.init("cfg/functions.yaml")
scfg = YAML.load_file("cfg/stdp.yaml")
fname = string(args["seed"])

function cluster_acc(c::Chromosome, X::Array{Float64}, Y::Array{Int64},
                     pname::String)
    n_cluster = length(unique(Y))
    nfunc = i->process(c, i)
    stdp_labels = stdp_cluster(
        X, Y, n_cluster, nfunc; seed=0, logfile=args["logfile"], problem=pname,
        fname=fname, train_epochs=scfg["train_epochs"],
        weight_mean=scfg["weight_mean"], weight_std=scfg["weight_std"],
        t_train=scfg["t_train"], t_blank=scfg["t_blank"], fr=scfg["fr"],
        pre_target=scfg["pre_target"], stdp_lr=scfg["stdp_lr"],
        stdp_mu=scfg["stdp_mu"], inhib_weight=scfg["inhib_weight"])
    acc = randindex(stdp_labels, Y)
    acc[1]
end

irisX, irisY = get_data("iris")

cluster_fit(c::Chromosome) = cluster_acc(c, irisX, irisY, "iris")

expert = nothing
if args["expert"] == "LIF"
    expert = to_chromo(lif_graph(-0.65, 0.3))
end
if args["expfile"] != ""
    experts = readdlm(args["expfile"], ',')
    s = mod(args["seed"], size(experts, 1))+1
    expert = PCGPChromo(experts[s, :], 5, 5)
end

maxfit, best = oneplus(PCGPChromo, 5, 5, cluster_fit; seed=args["seed"], expert=expert)

yeastX, yeastY = get_data("yeast")
yeast_acc = cluster_acc(best, yeastX, yeastY, "yeast")
Logging.info(@sprintf("F: %d %d %0.5f %0.5f %d %d",
seed, eval_count, maxfit, yeast_acc,
sum([n.active for n in best.nodes]),
length(best.nodes)))


text:
# from McCollum et al,
# "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior"
# Comp. Bio. and Chem., 30, pg. 39-49 (2006).

using DiffEqBase, DiffEqJump, DiffEqBiological
using Base.Test

using Plots; plotlyjs()
doplot = false
using BenchmarkTools
dobenchmark = true

dotestmean   = false
doprintmeans = false

# SSAs to test
#SSAalgs = (SortingDirect(),Direct()) #, DirectFW(), FRM(), FRMFW())
SSAalgs = (Direct(), SortingDirect(),)

# Nsims        = 80000
# tf           = 100000.0
Nsims        = 100
tf           = 500.0
u0           = [1,0,0,0]
expected_avg = 5.926553750000000e+02
reltol       = .01

# average number of proteins in a simulation
function runSSAs(jump_prob)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[3,end]
    end
    mean(Psamp)
end

# MODEL SETUP
N = 10
genenetwork = "@reaction_network gtype begin\n"
for i in 1:N
    genenetwork *= "\t 10.0, G$(2*i-1) --> G$(2*i-1) + M$(2*i-1)\n"
    genenetwork *= "\t 10.0, M$(2*i-1) --> M$(2*i-1) + P$(2*i-1)\n"
    genenetwork *= "\t 1.0,  M$(2*i-1) --> 0\n"
    genenetwork *= "\t 1.0,  P$(2*i-1) --> 0\n"

    genenetwork *= "\t 5.0, G$(2*i) --> G$(2*i) + M$(2*i)\n"
    genenetwork *= "\t 5.0, M$(2*i) --> M$(2*i) + P$(2*i)\n"
    genenetwork *= "\t 1.0,  M$(2*i) --> 0\n"
    genenetwork *= "\t 1.0,  P$(2*i) --> 0\n"

    genenetwork *= "\t 0.0001, G$(2*i) + P$(2*i-1) --> G$(2*i)_ind \n"
    genenetwork *= "\t 100., G$(2*i)_ind --> G$(2*i)_ind + M$(2*i)\n"
end
genenetwork *= "end"
rs = eval( parse(genenetwork) )
u0 = zeros(Int, length(rs.syms))
for i = 1:(2*N)
    u0[findfirst(rs.syms, Symbol("G$(i)"))] = 1
end

# TESTING:
prob = DiscreteProblem(u0, (0.0, tf))

# plotting one full trajectory
if doplot
    plothand = plot(reuse=false)
    for alg in SSAalgs
        jump_prob = JumpProblem(prob, alg, rs, save_positions=(false,false))
        sol = solve(jump_prob, SSAStepper(), saveat=(tf/1000.))
        plot!(plothand, sol.t, sol[[3,6],:]', lab=string(alg))
    end
    display(plothand)
end

# # test the means
# if dotestmean
#     means = zeros(Float64,length(SSAalgs))
#     for (i,alg) in enumerate(SSAalgs)
#         jump_prob = JumpProblem(prob, alg, majumps, save_positions=(false,false))
#         means[i]  = runSSAs(jump_prob)
#         relerr = abs(means[i] - expected_avg) / expected_avg
#         if doprintmeans
#             println("Mean from method: ", typeof(alg), " is = ", means[i], ", rel err = ", relerr)
#         end

#         if dobenchmark
#             @btime (runSSAs($jump_prob);)
#         end

#         @test abs(means[i] - expected_avg) < reltol*expected_avg
#     end
# end

# # benchmark performance
if dobenchmark
    # exact methods
    for alg in SSAalgs
        println("Solving with method: ", typeof(alg), ", using SSAStepper")
        jump_prob = JumpProblem(prob, alg, rs, save_positions=(false,false))
        @btime solve($jump_prob, SSAStepper(), saveat=(tf/1000.))
        #@btime (runSSAs($jump_prob);)
    end
    println()
end

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
module StructuralEquationModels

using LinearAlgebra, Optim,
    NLSolversBase, Statistics, SparseArrays, Symbolics,
    NLopt, FiniteDiff, ForwardDiff, PrettyTables,
    Distributions, StenoGraphs, LazyArtifacts, DelimitedFiles,
    DataFrames

import DataFrames: DataFrame
export *, ==, @StenoGraph, AbstractEdge, AbstractNode, DirectedEdge, Edge, EdgeModifier, MetaEdge, MetaNode, ModifiedEdge, ModifiedNode, Modifier, ModifyingNode, Node, NodeModifier, NodeOrEdgeModifier, SimpleNode, StenoGraphs, UndirectedEdge, convert, promote_rule, show, unarrow, unmeta, ←, →, ↔, ⇐, ⇒, ⇔

# type hierarchy
include("types.jl")
# fitted objects
include("frontend/fit/SemFit.jl")
# specification of models
include("frontend/specification/ParameterTable.jl")
include("frontend/specification/EnsembleParameterTable.jl")
include("frontend/specification/RAMMatrices.jl")
include("frontend/specification/StenoGraphs.jl")
include("frontend/fit/summary.jl")
# pretty printing
include("frontend/pretty_printing.jl")
# observed
include("observed/common.jl")
include("observed/missing.jl")
include("observed/EM.jl")
# constructor
include("frontend/specification/Sem.jl")
# imply
include("imply/RAM/symbolic.jl")
include("imply/RAM/generic.jl")
include("imply/empty.jl")
# loss
include("loss/ML/ML.jl")
include("loss/ML/FIML.jl")
include("loss/regularization/ridge.jl")
include("loss/WLS/WLS.jl")
include("loss/constant/constant.jl")
# diff
include("diff/optim.jl")
include("diff/NLopt.jl")
include("diff/Empty.jl")
# optimizer
include("optimizer/optim.jl")
include("optimizer/NLopt.jl")
# helper functions
include("additional_functions/helper.jl")
include("additional_functions/parameters.jl")
include("additional_functions/start_val/start_val.jl")
include("additional_functions/start_val/start_fabin3.jl")
include("additional_functions/start_val/start_partable.jl")
include("additional_functions/start_val/start_simple.jl")
include("additional_functions/artifacts.jl")
include("additional_functions/simulation.jl")
# identifier
include("additional_functions/identifier.jl")
# fit measures
include("frontend/fit/fitmeasures/AIC.jl")
include("frontend/fit/fitmeasures/BIC.jl")
include("frontend/fit/fitmeasures/chi2.jl")
include("frontend/fit/fitmeasures/df.jl")
include("frontend/fit/fitmeasures/minus2ll.jl")
include("frontend/fit/fitmeasures/n_par.jl")
include("frontend/fit/fitmeasures/n_obs.jl")
include("frontend/fit/fitmeasures/p.jl")
include("frontend/fit/fitmeasures/RMSEA.jl")
include("frontend/fit/fitmeasures/n_man.jl")
include("frontend/fit/fitmeasures/fit_measures.jl")
# standard errors
include("frontend/fit/standard_errors/hessian.jl")
include("frontend/fit/standard_errors/bootstrap.jl")

export  AbstractSem,
            AbstractSemSingle, AbstractSemCollection, Sem, SemFiniteDiff, SemForwardDiff, SemEnsemble,
        SemImply,
            RAMSymbolic, RAM, ImplyEmpty,
        start_val,
            start_fabin3, start_simple, start_parameter_table,
        SemLoss,
            SemLossFunction, SemML, SemFIML, em_mvn, SemLasso, SemRidge,
            SemConstant, SemWLS,
        SemDiff,
            SemDiffEmpty, SemDiffOptim, SemDiffNLopt, NLoptConstraint,
        SemObs,
            SemObsCommon, SemObsMissing,
        sem_fit,
        SemFit,
            minimum, solution,
        sem_summary,
        objective, objective!, gradient, gradient!, hessian, hessian!, objective_gradient!,
        ParameterTable,
            EnsembleParameterTable, update_partable!, update_estimate!, update_start!,
            Fixed, fixed, Start, start, Label, label,
            get_identifier_indices,
        RAMMatrices,
            RAMMatrices!,
        identifier,
        fit_measures,
            AIC, BIC, χ², df, fit_measures, minus2ll, n_par, n_obs, p_value, RMSEA, n_man,
            EmMVNModel,
        se_hessian, se_bootstrap,
        example_data,
        @StenoGraph, →, ←, ↔, ⇔
end


text:
module merge_sqlite

import Base
import SQLite
import DataFrames

const all_parameters = Base.Dict([("ID", "ID"),
                    ("LAT", "Latitude"),
                    ("LON", "Longitude"),
                    ("NN", "CloudCover2D"),
                    ("FI", "GeopotentialHeight"),
                    ("DD", "WindDirection"),
                    ("FF", "WindSpeed"),
                    ("GG", "WindGust"),
                    ("TT", "T2m"),
                    ("RH", "RH2m"),
                    ("PS", "MSLP"),
                    ("PSS", "PSurface"),
                    ("PE", "Precip12H"),
                    ("PE1", "Precip1H"),
                    ("PE3", "Precip3H"),
                    ("PE6", "Precip6H"),
                    ("PE24", "Precip24H"),
                    ("QQ", "Q2m"),
                    ("VI", "Visibility"),
                    ("TD", "TD2m"),
                    ("TX", "T2m_MaximaPast6H"),
                    ("TM", "T2m_MinimaPast6H"),
                    ("GM", "WindSpeed_Maxima1H"),
                    ("GX", "Max_Windgust1H"),
                    ("WX", "Unknown1"),
                    ("GW", "Unknown2"),
                    ("TIME", "Time")])

function merge(cmd_message)
    """Merges multiple SQLite files to a single SQLite file
    Assumes that the cmd_message comes in the following format:
    ("merge_sqlite", "db_prefix", "db_postfix", "indir", "query_string", "sqlitefile")
    """
    @info "Starting merge_sqlite"

    db_prefix     = cmd_message[2]
    db_postfix    = cmd_message[3]
    indir         = cmd_message[4]
    query         = cmd_message[5]
    sqlfile       = cmd_message[6]

    parameters = split_query_to_parameters(query)

    db_files = find_db_files(db_prefix, db_postfix, indir)

    merge_sqlite_files(db_files, parameters, query, sqlfile)

    @info "Finished"
end

function split_query_to_parameters(query::String)
    splitted_query = Base.split(query," ")
    parameters = String.(Base.split(splitted_query[2],",")) # Convert substring to array of strings

    if parameters[1] == "*"
        parameters = [k for (k,v) in all_parameters]
    end
    return parameters
end

function find_db_files(db_prefix::String, db_postfix::String, indir::String)

    files = readdir(indir, join=false)
    db_files = [indir*x for x in files if (startswith(x, db_prefix) && endswith(x, db_postfix))]

    return db_files
end

function merge_sqlite_files(db_files::Array, parameters::Array, query::String, sqlfile::String)

    db_target = SQLite.DB(sqlfile)
    make_target_table(db_target, parameters)

    get_and_put_data(db_files, parameters, query, db_target)

    return
end

function get_and_put_data(db_files::Array, parameters::Array, query::String, db_target::SQLite.DB)

    println(query)
    for f in db_files
        db_source = SQLite.DB(f)

        # TODO: If query returns too much data, StackOverflowError occurs. Instead one can possibly iterate over parameters
        # or limit the query. Until then:
        try
            data = SQLite.DBInterface.execute(db_source, query) |> DataFrames.DataFrame
        catch e
            if isa(e, StackOverflowError)
                @warn "LoadError: StackOverflowError. Try to reduce the data size returned by the query"
            end
        end

    end
end

function make_target_table(db::SQLite.DB, parameters::Array)
    """Makes the SQL Table if it does not exist"""

        id_present   = false
        time_present = false

        create_table_string = "CREATE TABLE IF NOT EXISTS vfld ("

        for par in parameters
            if par == "ID"
                parstring = "ID INT DEFAULT NULL, "
                id_present = true
            elseif par == "TIME"
                parstring = "TIME INT DEFAULT NULL, "
                time_present = true
            else
                parstring = String(par) * " REAL DEFAULT NULL, "
            end
            create_table_string = create_table_string * parstring
        end

        if id_present && time_present
            create_table_string = create_table_string * "PRIMARY KEY (ID, TIME));"
        elseif id_present
            create_table_string = create_table_string * "PRIMARY KEY (ID));"
        elseif time_present
            create_table_string = create_table_string * "PRIMARY KEY (TIME));"
        else
            create_table_string = create_table_string[1:length(create_table_string)-2] * ");"
        end

        SQLite.execute(db, create_table_string)
end

function inject_data(db::SQLite.DB, dataTable)
    """Inject data into SQLite Table"""
   try
       SQLite.load!(dataTable, db, "vfld")
   catch e
       if isa(e, LoadError)
         println("SQLite.SQLiteException, no such savepoint: Error caught but possibly data loss")
       end
   end
end

end


text:
using ImageFiltering
using Statistics
using StatsBase

function badpixremoval_clump(frame::AbstractMatrix{T}; max_niter=15, kwargs...) where T
    # create bad pixel map by sigma clipping outliers
    bad_pixels = badpixel_clip(frame; kwargs...)
    bad_pixels_cumul = copy(bad_pixels)
    total = sum(bad_pixels)

    # iterate replacing bad pixels with sigma clipping
    out = copy(frame)
    it = 0
    while total > 0 && it < max_niter
        it += 1
        out = sigma_filter!(out, bad_pixels)
        bad_pixels = badpixel_clip(out)
        bad_pixels_cumul .|= bad_pixels
        total = sum(bad_pixels)
    end
    return out, bad_pixels_cumul
end

function badpixel_clip(frame; sigma=4, box_size=9)
    med = mapwindow(median, frame, (box_size, box_size); border="reflect")
    std = mapwindow(mad, frame, (box_size, box_size); border="reflect")

    return @. abs(frame - med) > sigma * std
end

function sigma_filter!(frame, bad_pixels; box_size=9)
    bpm = copy(bad_pixels)

    min_neigh = sum(3:2:box_size)
    half_box = box_size ÷ 2
    while sum(bpm) > 0
        pad_bpm = BorderArray(bpm, Pad(:reflect, half_box, half_box))
        pad_im = BorderArray(frame, Pad(:reflect, half_box, half_box))
        for idx in findall(bpm)
            ax1 = idx.I[1] - half_box:idx.I[1] + half_box
            ax2 = idx.I[2] - half_box:idx.I[2] + half_box
            block_bpm = @view pad_bpm[ax1, ax2]

            # ensure enough "good" neighbors for imputation
            if sum(block_bpm) < min_neigh
                block_im = @view pad_im[ax1, ax2]
                frame[idx] = median(block_im)
                bpm[idx] = false
            end
        end
    end
    return frame
end

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
# Autogenerated wrapper script for rocminfo_jll for x86_64-linux-gnu-cxx11
export rocm_agent_enumerator, rocminfo

using hsa_rocr_jll
JLLWrappers.@generate_wrapper_header("rocminfo")
JLLWrappers.@declare_executable_product(rocm_agent_enumerator)
JLLWrappers.@declare_executable_product(rocminfo)
function __init__()
    JLLWrappers.@generate_init_header(hsa_rocr_jll)
    JLLWrappers.@init_executable_product(
        rocm_agent_enumerator,
        "bin/rocm_agent_enumerator",
    )

    JLLWrappers.@init_executable_product(
        rocminfo,
        "bin/rocminfo",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %%
using Printf
using Plots
default(size = (400, 250), legend=false, titlefontsize=12)

# %%
function laplacian!(v::Vector{Float64}, u::Vector{Float64}, dx::Float64)::Nothing
    v[1] = (u[end] + u[2] - 2u[1])/dx^2
    v[2:end-1] = (u[1:end-2] + u[3:end] - 2u[2:end-1])/dx^2
    v[end] = (u[end-1] + u[1] - 2u[1end])/dx^2
    return
end

# %%
n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
f(x) = sin(x) - cos(2x)
u = f.(x)
plot(x, u; label="", title="y = u(x)")

# %%
v = similar(u)
laplacian!(v, u, dx)
plot(x, v; label="", title="y = u''(x)")

# %%
function laplacian!(v::Vector{Float64}, u::Vector{Float64}, dx::Float64)::Nothing
    v[1] = (u[end] + u[2] - 2u[1])/dx^2
    v[2:end-1] = (u[1:end-2] + u[3:end] - 2u[2:end-1])/dx^2
    v[end] = (u[end-1] + u[1] - 2u[1end])/dx^2
    return
end

function heateq(u0::Vector{Float64}, dx::Float64, tmax::Float64, N::Int=200)::Tuple{Vector{Float64}, Matrix{Float64}}
    t = 0:dx:tmax
    dt = t[2] - t[1]
    u = Matrix{Float64}(undef, length(u0), length(t)+1)
    u[:, 1] = u0
    v = Vector{Float64}(undef, length(u0))
    for i in 2:length(t)+1
        u[:, i] = u[:, i-1]
        for _ in 1:N
            laplacian!(v, u[:, i], dx)
            u[:, i] += v*dt/N
        end
    end
    t, u
end

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
f(x) = sin(x) - cos(2x)
u0 = f.(x)
tmax = 1.0

@time t, u = heateq(u0, dx, tmax)
@time t, u = heateq(u0, dx, tmax)
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot(x, u[:, i]; ylim=extrema(u0), title)
end

# %%
function laplacian2!(v::Vector{Float64}, u::Vector{Float64}, dx::Float64)::Nothing
    v[1] = (u[end] + u[2] - 2u[1])/dx^2
    @. @views v[2:end-1] = (u[1:end-2] + u[3:end] - 2u[2:end-1])/dx^2
    v[end] = (u[end-1] + u[1] - 2u[1end])/dx^2
    return
end

function heateq2(u0::Vector{Float64}, dx::Float64, tmax::Float64, N::Int=200)::Tuple{Vector{Float64}, Matrix{Float64}}
    t = 0:dx:tmax
    dt = t[2] - t[1]
    u = Matrix{Float64}(undef, length(u0), length(t)+1)
    u[:, 1] = u0
    v = Vector{Float64}(undef, length(u0))
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian2!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
f(x) = sin(x) - cos(2x)
u0 = f.(x)
tmax = 1.0

t, u = heateq2(u0, dx, tmax)

# %%
function laplacian3!(v::Vector{Float64}, u, dx::Float64)::Nothing
    v[1] = (u[end] + u[2] - 2u[1])/dx^2
    @. @views v[2:end-1] = (u[1:end-2] + u[3:end] - 2u[2:end-1])/dx^2
    v[end] = (u[end-1] + u[1] - 2u[1end])/dx^2
    return
end

function heateq3(u0::Vector{Float64}, dx::Float64, tmax::Float64, N::Int=200)::Tuple{Vector{Float64}, Matrix{Float64}}
    t = 0:dx:tmax
    dt = t[2] - t[1]
    u = Matrix{Float64}(undef, length(u0), length(t)+1)
    u[:, 1] = u0
    v = Vector{Float64}(undef, length(u0))
    N = 100
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian3!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
f(x) = sin(x) - cos(2x)
u0 = f.(x)
tmax = 1.0

@time t, u = heateq3(u0, dx, tmax)
@time t, u = heateq3(u0, dx, tmax)
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot(x, u[:, i]; ylim=extrema(u0), title)
end

# %%
function laplacian4!(v, u, dx)
    v[begin] = (u[end] + u[begin+1] - 2u[begin])/dx^2
    @. @views v[begin+1:end-1] = (u[begin:end-2] + u[begin+2:end] - 2u[begin+1:end-1])/dx^2
    v[end] = (u[end-1] + u[begin] - 2u[1end])/dx^2
    return
end

function heateq4(u0, dx, tmax, N=200)
    t = 0:dx:tmax
    dt = step(t)
    u = similar(u0, length(u0), length(t)+1)
    u[:, 1] = u0
    v = similar(u0)
    N = 100
    for i in 2:length(t)+1
        @. @views u[:, i] = u[:, i-1]
        for _ in 1:N
            @views laplacian4!(v, u[:, i], dx)
            @. @views u[:, i] += v*dt/N
        end
    end
    t, u
end

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
f(x) = sin(x) - cos(2x)
u0 = f.(x)
tmax = 1.0

@time t, u = heateq4(u0, dx, tmax)
@time t, u = heateq4(u0, dx, tmax)
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot(x, u[:, i]; ylim=extrema(u0), title)
end

# %%
n = 200
x = range(-Float32(π), π; length=n+1)[1:end-1]
dx = step(x)
@show typeof(x) typeof(dx)
f(x) = sin(x) - cos(2x)
u0 = f.(x)
@show typeof(u0)
tmax = Float32(1)

t, u = heateq4(u0, dx, tmax)
@show typeof(t) typeof(u)
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot(x, u[:, i]; ylim=extrema(u0), title)
end

# %%
using StaticArrays

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
m = 200
noise1, noise2 = 0.3randn(m), 0.3randn(m)
f(x) = SVector{m}((1 .+ noise1)*sin(x) .- (1 .+ noise2)*cos(2x))
u0 = f.(x)

plot()
for k in 1:m
    plot!(x, (p -> p[k]).(u0); lw=0.1, color=:blue)
end
plot!()

# %%
tmax = 1.0
@time t, u = heateq4(u0, dx, tmax)
@time t, u = heateq4(u0, dx, tmax)

ylim = (minimum(minimum.(u0)), maximum(maximum.(u0)))
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot()
    for k in 1:m
        plot!(x, (p -> p[k]).(u[:, i]); ylim, lw=0.1, color=:blue)
    end
    title!(title)
end

# %%
# laplacian! と heateq で型を書くのをやめたもの (最適化がまったくされていない場合)

function laplacian1!(v, u, dx)
    v[1] = (u[end] + u[2] - 2u[1])/dx^2
    v[2:end-1] = (u[1:end-2] + u[3:end] - 2u[2:end-1])/dx^2
    v[end] = (u[end-1] + u[1] - 2u[1end])/dx^2
    return
end

function heateq1(u0, dx, tmax, N=200)
    t = 0:dx:tmax
    dt = t[2] - t[1]
    u = similar(u0, length(u0), length(t)+1)
    u[:, 1] = u0
    v = similar(u0)
    for i in 2:length(t)+1
        u[:, i] = u[:, i-1]
        for _ in 1:N
            laplacian1!(v, u[:, i], dx)
            u[:, i] += v*dt/N
        end
    end
    t, u
end

using StaticArrays

n = 200
x = range(-π, π; length=n+1)[1:end-1]
dx = step(x)
m = 200
noise1, noise2 = 0.3randn(m), 0.3randn(m)
f(x) = SVector{m}((1 .+ noise1)*sin(x) .- (1 .+ noise2)*cos(2x))
u0 = f.(x)
tmax = 1.0

@time t, u = heateq1(u0, dx, tmax)
@time t, u = heateq1(u0, dx, tmax)

ylim = (minimum(minimum.(u0)), maximum(maximum.(u0)))
@gif for i in 1:length(t)
    title = @sprintf("t = %4.2f", t[i])
    plot()
    for k in 1:m
        plot!(x, (p -> p[k]).(u[:, i]); ylim, lw=0.1, color=:blue)
    end
    title!(title)
end

# %%

max_stars_repo_path: 0004/1D heat equation example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0004/1D heat equation example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0004/1D heat equation example.jl

max_forks_repo_name: genkuroki/public


text:
Base.@kwdef struct NestedIntervalDisplayProperties{T₁ <: AbstractContext, T₂ <: NTuple} <: AbstractDisplayProperties
    id::String
    caption::String
    col::ImVec4 = ImVec4(0, 0, 0, 1)
    createwindow::Bool = true
    layout::RectangularLayout
    plotcontext::T₁
    padding::T₂ = (0, 0, 0 ,0)
end

function get_id(p::NestedIntervalDisplayProperties)
    p.id
end

function set_id!(p::NestedIntervalDisplayProperties, id::String)
    p.id = id
end

function get_caption(p::NestedIntervalDisplayProperties)
    p.caption
end

function set_caption!(p::NestedIntervalDisplayProperties, caption::String)
    p.caption = caption
end

function get_col(p::NestedIntervalDisplayProperties)
    p.col
end

function set_col!(p::NestedIntervalDisplayProperties, col::ImVec4)
    p.col = col
end

function is_new_window(p::NestedIntervalDisplayProperties)
    p.createwindow
end

function get_layout(p::NestedIntervalDisplayProperties)
    p.layout
end

function set_layout!(p::NestedIntervalDisplayProperties, layout::RectangularLayout)
    p.layout = layout
end

function get_padding(p::NestedIntervalDisplayProperties)
    p.padding
end

function set_padding!(p::NestedIntervalDisplayProperties, padding::NTuple)
    p.padding = padding
end

function get_plotcontext(p::NestedIntervalDisplayProperties)
   p.plotcontext
end

function set_plotcontext!(p::NestedIntervalDisplayProperties, plotcontext::AbstractContext)
   p.plotcontext = plotcontext
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule soapysdr_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("soapysdr")
JLLWrappers.@generate_main_file("soapysdr", UUID("343a40d9-ed99-5d34-8b56-649aaa4ecee6"))
end  # module soapysdr_jll


text:
# In the .mat file, the following code is used:
#  [chosen_juice, licking_result] = [1,1]   # animal chose A
#                                 = [2,1]   # animal chose B
#                                 = [0,1]   # animal made wrong choice during forced choice trial
#                                 = [0,0]   # animal made no licking responses
function mat_trialresult(nA, nB, choice, licked)
    leftA = nA < 0 || nB > 0
    nA, nB = abs(nA), abs(nB)
    licked == 0 && return TrialResult(nA, nB, leftA, missing)
    if choice == 0
        nA == 0 && nB == 0 && error("double-forced unexpected")
        nA == 0 || nB == 0 || error("expected this to be a forced-choice")
        choiceA = nA == 0
    else
        choiceA = choice == 1
    end
    return TrialResult(nA, nB, leftA, choiceA)
end

"""
    parsemat(filename)
    parsemat(readfunction, filename)

Parse the trial data in `filename`. Optionally pass a `readfunction` to extract additional elements;
`readfunction` should have syntax

```julia
function readfunction(data, args...)
    moredata = data["matlab_variable_name"]
    # maybe do some processing/validation
    return args..., moredata
end
```

`data` is the result of `MAT.matread(filename)`.

See also: [`positive_cells`](@ref).
"""
parsemat(filename::AbstractString) = parsemat(returntail, filename)

function parsemat(f, filename::AbstractString)
    data = matread(filename)
    celldata = data["celldata"]
    goodtrials = data["goodTrials"]
    psyphydata = data["psyphydata"]

    # parse celldata
    cts = OrderedDict{Int,CellsTrial{Float64}}()   # not all trial indexes are guaranteed to be present
    i0 = i = 1
    idx = Int(celldata[i,2])
    while i <= size(celldata, 1)
        i += 1
        idxnew = i <= size(celldata, 1) ? Int(celldata[i,2]) : idx+1
        if idxnew != idx
            rng = i0:i-1
            cts[idx] = CellsTrial{Float64}(celldata[rng,1]*ms, celldata[rng,3:end])
            i0 = i
            idx = idxnew
        end
    end

    # parse goodtrials (here we assume all trials are listed)
    tts = OrderedDict{Int,TrialResult}()
    for r in eachrow(goodtrials)
        tts[Int(r[1])] = mat_trialresult(r[2:end]...)
    end

    # parse psyphydata
    ets = OrderedDict{Int,EventTiming}()
    j = 0
    while j < size(psyphydata,1)
        trial_start = offer_on = offer_off = go = choice = trial_end = Tms(NaN)
        i = Int(psyphydata[j+1,2])
        while j < size(psyphydata,1) && psyphydata[j+1,2] == i
            j += 1
            t, code = psyphydata[j,1]*ms, psyphydata[j,3]
            if code == 21
                trial_start = t
            elseif code == 30
                offer_on = t
            elseif code == 31
                offer_off = t
            elseif code == 35
                go = t
            elseif code == 39
                choice = t
            elseif code == 43 || code == 44
                @assert t == choice
            elseif code == 60
                trial_end = t
            else
                error("code $code not recognized")
            end
        end
        ets[i] = EventTiming(trial_start, offer_on, offer_off, go, choice, trial_end)
    end

    keys(cts) == keys(tts) == keys(ets) || throw(DimensionMismatch("trial indices must agree"))

    return f(data, cts, tts, ets)
end

returntail(data, args...) = args

"""
    positive_cells(data, args...)

Extract the vector of `idx_positive_cells` from the .mat file.
"""
function positive_cells(data, args...)
    celldata = data["celldata"]
    poscells = data["idx_positive_cells"]
    size(poscells, 2) == 1 || throw(DimensionMismatch("expected idx_positive_cells to be a column vector"))
    poscells = vec(poscells)
    size(celldata, 2) == length(poscells) + 2 ||
        throw(DimensionMismatch("$filename: idx_positive_cells must match columns in celldata"))

    poscells = convert(AbstractVector{Bool}, poscells)

    return args..., poscells
end


text:
@register_symbolic Base.binomial(n,k)

@register_symbolic Base.sign(x)::Int
derivative(::typeof(sign), args::NTuple{1,Any}, ::Val{1}) = 0

@register_symbolic Base.signbit(x)::Bool
derivative(::typeof(signbit), args::NTuple{1,Any}, ::Val{1}) = 0
derivative(::typeof(abs), args::NTuple{1,Any}, ::Val{1}) = IfElse.ifelse(signbit(args[1]),-one(args[1]),one(args[1]))

function derivative(::typeof(min), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x < y, one(x), zero(x))
end
function derivative(::typeof(min), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x < y, zero(y), one(y))
end
function derivative(::typeof(max), args::NTuple{2,Any}, ::Val{1})
    x, y = args
    IfElse.ifelse(x > y, one(x), zero(x))
end
function derivative(::typeof(max), args::NTuple{2,Any}, ::Val{2})
    x, y = args
    IfElse.ifelse(x > y, zero(y), one(y))
end

@register_symbolic Base.ceil(x)
@register_symbolic Base.floor(x)
@register_symbolic Base.factorial(x)

function derivative(::Union{typeof(ceil),typeof(floor),typeof(factorial)}, args::NTuple{1,Any}, ::Val{1})
    zero(args[1])
end

@register_symbolic Base.rand(x)
@register_symbolic Base.randn(x)

@register_symbolic Distributions.pdf(dist,x)
@register_symbolic Distributions.logpdf(dist,x)
@register_symbolic Distributions.cdf(dist,x)
@register_symbolic Distributions.logcdf(dist,x)
@register_symbolic Distributions.quantile(dist,x)

@register_symbolic Distributions.Uniform(mu,sigma) false
@register_symbolic Distributions.Normal(mu,sigma) false

@register_symbolic ∈(x::Real, y::AbstractArray)::Bool
@register_symbolic ∪(x, y)
@register_symbolic ∩(x, y)
@register_symbolic ∨(x, y)
@register_symbolic ∧(x, y)
@register_symbolic ⊆(x, y)

LinearAlgebra.norm(x::Num, p::Real) = abs(x)


text:
#Regression Constants
const RANGEMAX = 62000.0 #ft
const ranges = [499, 800, 2000, 3038, 5316, 6450, 7200, 7950, 8725, 10633, 13671, 16709, 19747,
                22785, 25823, 28862, 31900, 34938, 37976, 41014, 48608, 60760]

const thetas = linspace(-pi,pi,41)
const psis   = linspace(-pi,pi,41)
const sos    = [100, 200, 300, 400, 500, 600, 700, 800]
const sis    = [0, 100, 200, 300, 400, 500, 600, 700, 800]
const taus   = [0, 1, 5, 10, 20, 40, 60, 80, 100]
const pas    = [1, 2, 3, 4, 5]
const pasTrue = [0, 1.5, -1.5, 3.0, -3.0]
const NSTATES = length(ranges)*length(thetas)*length(psis)*length(sos)*length(sis)*length(taus)*length(pas)
const ACTIONS = deg2rad([0.0 1.5 -1.5 3.0 -3.0])

#Deep RL Constants
const RangeMax = 3000.0 # meters
const RANGEMIN = 0.0 # meters
const PsiDim   = 41
const PsiMin   = -pi #[rad]
const PsiMax   = pi  #[rad]

const psis_drl      = linspace(PsiMin,PsiMax,PsiDim)
const sos_drl       = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
const sis_drl       = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

const Actions= [-20.0,-10.0,0.0,10.0,20.0,-6.0]

const STATE_DIM = 5
const ACTION_DIM = 6

const Thetamin = -pi #[rad]
const Thetamax = pi #[rad]
const Bearingmin = -pi  # [rad]
const Bearingmax = pi  # [rad]
const Speedmin = 10.0  # [m/s]
const Speedmax = 20.0  # [m/s]

const Rangedim = 67
const Thetadim = 41
const Bearingdim = 37 #37 #3
const Speeddim = 5 #5 #2
const Ranges = [0., 10., 20., 30., 40., 50., 60., 80., 100., 120.,
                140., 160., 180., 200., 230., 250., 270., 290., 310., 330.,
                350., 370., 390., 410., 430., 450., 470., 490., 510., 530.,
                570., 600., 630., 660., 690., 720., 750., 780., 810., 840.,
                870., 900., 930., 960., 1000., 1040., 1080., 1120., 1160., 1200.,
                1240., 1280., 1320., 1360., 1400., 1450., 1500., 1550., 1600., 1650.,
                1700., 1750., 1800., 1850., 1900., 1950., 2000.]
const Thetas = collect(linspace(Thetamin, Thetamax,Thetadim))
const Bearings = collect(linspace(Bearingmin, Bearingmax, Bearingdim))
const Speeds = collect(linspace(Speedmin, Speedmax, Speeddim))
const NSTATES_drl = Rangedim*Thetadim*Bearingdim*Speeddim^2+1

#Compatible with X and Y grids rather than range and theta
const Xdim = 71 #51 #8
const Ydim = 71 #51 #8
const Xs = collect(linspace(-RangeMax,RangeMax,Xdim))
const Ys = collect(linspace(-RangeMax,RangeMax,Ydim))
const NSTATES_drl_xandy = Xdim*Ydim*Bearingdim*Speeddim^2+1


text:
using FFTW

function split_and_zeropad(a, b)
    a_digits = digits(a)
    b_digits = digits(b)

    return vcat(a_digits, zeros(Int, length(b_digits))),
           vcat(b_digits, zeros(Int, length(a_digits)))
end

# TODO: there is probably a better way to do the carrying operations
#       check biginteger source code.
function vector_to_int(v)
    sum = 0
    for i = 1:length(v)
        sum += 10^(i-1)*v[i]
    end

    return sum
end

function carry!(v::Vector{Int})
    for i = 1:length(v)-1
        digs = digits(v[i])
        v[i+1] += vector_to_int(digs[2:end])
        v[i] = digs[1]
    end

    return v
end

function schonhage_strassen(a, b)
    va, vb = split_and_zeropad(a, b)

    vf = round.(Int, real(ifft(fft(va).*fft(vb))))

    carry!(vf)

    return vector_to_int(vf)
end

max_stars_repo_name: leios/simuleios

max_issues_repo_name: leios/simuleios

max_forks_repo_name: leios/simuleios


text:
@inline function getVelocityDelta(joint::Translational0, body1::AbstractBody, body2::Body{T}, v::SVector{3,T}) where T
    Δv = v # in body1 frame
    return Δv
end

@inline function getPositionDelta(joint::Translational0, body1::AbstractBody, body2::Body{T}, x::SVector{3,T}) where T
    Δx = x # in body1 frame
    return Δx
end

@inline function setForce!(joint::Translational0, body1::Body, body2::Body{T}, F::SVector{3,T}, No) where T
    F1 = vrotate(-F, body1.q[No])
    F2 = -F1

    body1.F[No] = F1
    body2.F[No] = F2
    return
end

@inline function setForce!(joint::Translational0, body1::Origin, body2::Body{T}, F::SVector{3,T}, No) where T
    body2.F[No] = F
    return
end

@inline function minimalCoordinates(joint::Translational0, body1::AbstractBody, body2::AbstractBody, No)
    body2.x[No]
end

@inline g(joint::Translational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = g(joint)

@inline ∂g∂posa(joint::Translational0, body1::AbstractBody, body2::AbstractBody, No) = ∂g∂posa(joint)
@inline ∂g∂posb(joint::Translational0, body1::AbstractBody, body2::AbstractBody, No) = ∂g∂posb(joint)
@inline ∂g∂vela(joint::Translational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = ∂g∂vela(joint)
@inline ∂g∂velb(joint::Translational0, body1::AbstractBody, body2::AbstractBody, Δt, No) = ∂g∂velb(joint)


text:
# Implementations of combinatorial functions from "Counting and Generating Terms in the Binary
# Lambda Calculus".

import Base

export grygiel_lescanne, grygiel_lescanne!, unrank, unrank!, terms

# TODO:
# - make k start from 0
# - adapt getindex of TermsIterator

########################
# COUNTING AND UNRANKING
########################

const Table{T} = Dict{Tuple{Int, Int}, T}

"""
     grygiel_lescanne(m, n[, ::Type{T}]) where {T}

Calculate the number of De Bruijn terms of binary size `n`, with at most `m` free variables.  This
corresponds to the series Sₘₙ from [this paper](https://arxiv.org/pdf/1511.05334v1.pdf).  The type
`T` can be used to specify the integer type used during the recursive construction.
"""
function grygiel_lescanne(m::Integer, n::Integer, ::Type{T} = Int) where {T<:Integer}
    return grygiel_lescanne!(m, n, Table{T}())
end

"""
    grygiel_lescanne!(m, n, table)

Calculate the number of de Bruijn terms of binary size `n`, with at most `m` free variables.  This
corresponds to the series Sₘₙ from [their paper](https://arxiv.org/pdf/1511.05334v1.pdf).

As a side effect, fill `table` with recursive calls whose entries correspond to `grygiel_lescanne(n, m)`.  This
method should be used instead of repeatedly calling `grygiel_lescanne`.
"""
function grygiel_lescanne!(m::Integer, n::Integer, table::Table{T}) where {T<:Integer}
    @assert(m >= 0)
    @assert(n >= 0)

    # Sₘ₀ = Sₘ₁ = 0
    # Sₘₙ = [m ≥ n - 1] + Sₘ₊₁,ₙ₋₂ + ∑_{k=0}^{n-2} Sₘ,k Sₘ,ₙ₋₂₋ₖ
    get!(table, (m, n)) do
        if n <= 1
            return 0
        else
            Int(m >= n - 1) +
                grygiel_lescanne!(m + 1, n - 2, table) +
                sum(grygiel_lescanne!(m, k, table) * grygiel_lescanne!(m, n - 2 - k, table) for k = 0:(n-2))
        end
    end
end

"""
    unrank{T}(m::Integer, n::Integer, k::Integer[, ::Type{T}]) -> Term

Calculate the `k`-th De Bruijn term of binary size `n`, with at most `m` free variables (indices
start at 1).  The type `T` can be used to specify the integer type used during the recursive
construction.
"""
function unrank(m::Integer, n::Integer, k::Integer, ::Type{T} = Int) where {T<:Integer}
    unrank!(m, n, k, Table{T}())
end

"""
    unrank!{T}(m::Integer, n::Integer, k::Integer, table::Table{T}) -> Term

Calculate the `k`-th De Bruijn term of binary size `n`, with at most `m` free variables,
using tabulated values for Sₙₘ.

This function should be used for repeated unrankings, so that the series Sₙₘ does
not have to be recomputed every time.
"""
function unrank!(m::Integer, n::Integer, k::Integer, table::Table{T}) where {T<:Integer}
    @assert(m >= 0)
    @assert(n >= 0)
    @assert(k >= 1)

    t = grygiel_lescanne!(m, n, table)
    @assert(t >= k)

    if m >= n - 1 && k == t
        return Var(n - 1)
    elseif k <= grygiel_lescanne!(m + 1, n - 2, table)
        return Abs(unrank!(m + 1, n - 2, k, table))
    else
        return unrank_app!(m, n - 2, 0, k - grygiel_lescanne!(m + 1, n - 2, table), table)
    end
end

function unrank_app!(m::Integer, n::Integer, j::Integer, r::Integer,
                     table::Table{T}) where {T<:Integer}
    tmnj = grygiel_lescanne!(m, n - j, table)
    tmjtmnj = grygiel_lescanne!(m, j, table) * tmnj

    if r <= tmjtmnj
        dv, rm = divrem(r - 1, tmnj)
        return App(unrank!(m, j, dv + 1, table), unrank!(m, n - j, rm + 1, table))
    else
        return unrank_app!(m, n, j + 1, r - tmjtmnj, table)
    end
end

#################
# ITERATING TERMS
#################

"""
    terms(m::Integer, n::Integer[, ::Type{T}])

An iterator for all de Bruijn terms of size `n` with at most `m` free variables.
"""
function terms(m, n, ::Type{T} = Int) where {T<:Integer}
    @assert(m >= 0)
    @assert(n >= 0)

    table = Table{T}()
    grygiel_lescanne!(m, n, table)
    TermsIterator{T}(m, n, table)
end

struct TermsIterator{T<:Integer}
    m::Int
    n::Int
    table::Table{T}
end

function Base.iterate(iter::TermsIterator, state = 1)
    if state ≤ iter.table[(iter.m, iter.n)]
        return (unrank!(iter.m, iter.n, state, iter.table), state + 1)
    else
        return nothing
    end
end

Base.IteratorSize(::Type{<:TermsIterator}) = Base.HasLength()
Base.length(iter::TermsIterator) = iter.table[(iter.m, iter.n)]
Base.IteratorEltype(::Type{<:TermsIterator}) = Base.HasEltype()
Base.eltype(::Type{<:TermsIterator}) = Term

Base.getindex(iter::TermsIterator, i::Integer) = unrank!(iter.m, iter.n, i, iter.table)
Base.firstindex(iter::TermsIterator) = 1
Base.lastindex(iter::TermsIterator) = grygiel_lescanne!(iter.m, iter.n, iter.table)


text:
@model function test0(x)
    λ ~ Gamma(2.0, inv(3.0))
    m ~ Normal(0, sqrt(1 / λ))
    x ~ Normal(m, sqrt(1 / λ))
end

@testdependencies(test0(1.4), λ, m, x)

@model function test1(x)
    s ~ Gamma(1.0, 1.0)
    x .~ Normal(0.0, s)
end

@testdependencies(test1([-0.5, 0.5]), s, x)

@model function test2(x)
    s ~ Gamma(1.0, 1.0)
    for i in eachindex(x)
        x[i] ~ Normal(float(i), s)
    end
end

@testdependencies(test2([-0.5, 0.5]), s, x[1], x[2])

@model function test3(w)
    p ~ Beta(1, 1)
    x ~ Bernoulli(p)
    y ~ Bernoulli(p)
    z ~ Bernoulli(p)
    w ~ MvNormal([x, y, z], 1.2)
end

@testdependencies(test3([1, 0, 1]), p, x, y, z, w)

@model function test4(x)
    μ ~ MvNormal(fill(0, 2), 2.0)
    z = Vector{Int}(undef, length(x))
    for i in 1:length(x)
        z[i] ~ Categorical([0.5, 0.5])
        x[i] ~ Normal(μ[z[i]], 0.1)
    end
end

@testdependencies(test4([1, 1, -1]), μ, z[1], z[2], z[3], x[1], x[2], x[3])

@model function test5(x)
    μ ~ MvNormal(fill(0, 2), 2.0)
    z = Vector{Int}(undef, length(x))
    z .~ Categorical.(fill([0.5, 0.5], length(x)))
    for i in 1:length(x)
        x[i] ~ Normal(μ[z[i]], 0.1)
    end
end

# there's a bug in DPPL occuring here: https://github.com/TuringLang/DynamicPPL.jl/issues/129
@testdependencies(test5([1, 1, -1]), μ, z, x[1], x[2], x[3])

@model function test6(x, y)
    λ ~ Gamma(2.0, inv(3.0))
    m ~ Normal(0, sqrt(1 / λ))
    D = Normal(m, sqrt(1 / λ))
    x ~ D
    y ~ D
end

@testdependencies(test6(1.4, 1.2), λ, m, x, y)

@model function test7(x)
    s ~ Gamma(1.0, 1.0)
    x[1] ~ Normal(0.0, s)
    for i in 2:length(x)
        x[i] ~ Normal(x[i-1], s)
    end
end

@testdependencies(test7([0.0, 0.1, -0.2]), s, x[1], x[2], x[3])

@model function test8()
    s ~ Gamma(1.0, 1.0)
    0.1 ~ Normal(0.0, s)
end

@testdependencies(test8(), s)

@model function test9(x)
    s ~ Gamma(1.0, 1.0)
    state = zeros(length(x) + 1)
    state[1] ~ Normal(0.0, s)
    for i in 1:length(x)
        state[i+1] ~ Normal(state[i], s)
        x[i] ~ Normal(state[i+1], s)
    end
end

@testdependencies(test9([0.0, 0.1, -0.2]), s, x[1], x[2], x[3], state[1], state[2], state[3], state[4])

@model function test10(x)
    s ~ Gamma(1.0, 1.0)
    state = zeros(length(x) + 1)
    state[1] = 42.0
    for i in 1:length(x)
        state[i+1] ~ Normal(state[i], s)
        x[i] ~ Normal(state[i+1], s)
    end
end

@testdependencies(test10([0.0, 0.1, -0.2]), s, x[1], x[2], x[3], state[2], state[3], state[4])

@model function test11(x)
    rpm = DirichletProcess(1.0)
    n = zeros(Int, length(x))
    z = zeros(Int, length(x))
    for i in eachindex(x)
        z[i] ~ ChineseRestaurantProcess(rpm, n)
        n[z[i]] += 1
    end
    K = findlast(!iszero, n)
    m ~ MvNormal(fill(0.0, K), 1.0)
    x ~ MvNormal(m[z], 1.0)
end

@testdependencies(test11([0.1, 0.05, 1.0]), m, x, z[1], z[2], z[3])


text:
module SputnikUtilities

using GroupedErrors
using JuliaDB, IndexedTables
using Parameters

import IndexedTables: AbstractIndexedTable

export Data2Select, SelectedData
export SelectPredicate, SelectValues
export Analysis, process

include("selectdata.jl")
include("process.jl")
include(joinpath("techniques", "statplotsrecipe.jl"))
include(joinpath("techniques", "groupederrors.jl"))
include(joinpath("techniques", "table.jl"))

end # module


text:
"""
    struct Alphabet{T}

A basic struct for storing alphabets. An alphabet consists of the symbols of a common type `T`.

# Example
```julia-repl
julia> Alphabet{Char}()
Empty alphabet of Char

julia> Alphabet(["a", "b", "c"])
Alphabet of String:
    1.  "a"
    2.  "b"
    3.  "c"
```
"""
struct Alphabet{T}
    letters::Vector{T}
    inversions::Vector{Int}

    function Alphabet(letters::AbstractVector{T}, inversions::AbstractVector{<:Integer}) where T
        @assert length(unique(letters)) == length(letters) "Non-unique set of letters"
        @assert !(T<:Integer) "Alphabets of integers are not supported."
        @assert length(letters) == length(inversions)
        @assert all(i-> 0 ≤ i ≤ length(letters), inversions)
        return new{T}(letters, inversions)
    end
end

Alphabet{T}() where T = Alphabet(T[])
Alphabet(init::AbstractVector{T}) where T = Alphabet(init, fill(0, length(init)))

letters(A::Alphabet) = A.letters

Base.length(A::Alphabet) = length(letters(A))
Base.isempty(A::Alphabet) = isempty(letters(A))
Base.:(==)(A::Alphabet, B::Alphabet) =
    letters(A) == letters(B) && A.inversions == B.inversions
Base.hash(A::Alphabet{T}, h::UInt) where T =
    hash(letters(A), hash(A.inversions, hash(h, hash(Alphabet{T}))))

Base.show(io::IO, A::Alphabet) = print(io, Alphabet, " ", letters(A))

hasinverse(i::Integer, A::Alphabet) = A.inversions[i] > 0
hasinverse(l::T, A::Alphabet{T}) where T = hasinverse(A[l], A)

function Base.show(io::IO, ::MIME"text/plain", A::Alphabet{T}) where T
    if isempty(A)
        print(io, "Empty alphabet of $(T)")
    else
        print(io, "Alphabet of $(T):")
        for (i, l) in pairs(letters(A))
            print(io, "\n\t$(i).\t")
            show(io, l)
            if hasinverse(i, A)
                print(io, " = (")
                show(io, A[-i])
                print(io, ")⁻¹")
            end
        end
    end
end

"""
    push!(A::Alphabet{T}, symbols::T...) where T

Push one or more elements of type `T` at the end of the alphabet `A`.

# Example
```julia-repl
julia> A = Alphabet{String}()
Empty alphabet of String

julia> push!(A, "a", "b")
Alphabet of String:
    1.  "a"
    2.  "b"
```
"""
function Base.push!(A::Alphabet{T}, symbols::T...) where T
    for s in symbols
        if findfirst(symbol -> symbol == s, letters(A)) !== nothing
            error("Symbol $(s) already in the alphabet.")
        end
        push!(A.letters, s)
        push!(A.inversions, 0)
    end
    A
end

"""
    set_inversion!(A::Alphabet{T}, x::T, y::T) where T

Set the inversion of `x` to `y` (and vice versa).

# Example
```julia-repl
julia> A = Alphabet(["a", "b", "c"])
Alphabet of String:
    1. "a"
    2. "b"
    3. "c"

julia> set_inversion!(A, "a", "c")
Alphabet of String:
    1. "a" = ("c")⁻¹
    2. "b"
    3. "c" = ("a")⁻¹

julia> set_inversion!(A, "a", "b")
Alphabet of String:
    1. "a" = ("b")⁻¹
    2. "b" = ("a")⁻¹
    3. "c"
```
"""
function set_inversion!(A::Alphabet{T}, x::T, y::T) where T
    if (ix = findfirst(symbol -> symbol == x, letters(A))) === nothing
        error("Element $(x) not found in the alphabet.")
    end
    if (iy = findfirst(symbol -> symbol == y, letters(A))) === nothing
        error("Element $(y) not found in the alphabet.")
    end

    if A.inversions[ix] > 0
        A.inversions[A.inversions[ix]] = 0
    end
    if A.inversions[iy] > 0
        A.inversions[A.inversions[iy]] = 0
    end

    A.inversions[ix] = iy
    A.inversions[iy] = ix
    A
end

"""
    getindex(A::Alphabet{T}, x::T) where T

Return the position of the symbol `x` in the alphabet `A`.

# Example
```julia-repl
julia> A = Alphabet(["a", "b", "c"])
Alphabet of String:
    1. "a"
    2. "b"
    3. "c"

julia> A["c"]
3
```
"""
function Base.getindex(A::Alphabet{T}, x::T) where T
    if (index = findfirst(symbol -> symbol == x, letters(A))) === nothing
        throw(DomainError("Element '$(x)' not found in the alphabet"))
    end
    return index
end

"""
    getindex(A::Alphabet{T}, p::Integer) where T

Return the symbol that holds the `p`th position in the alphabet `A`. If `p < 0`, then the inversion of the `|p|`th symbol is returned.

# Example
```julia-repl
julia> A = Alphabet(["a", "b", "c"])
Alphabet of String:
    1. "a"
    2. "b"
    3. "c"

julia> set_inversion!(A, "a", "c")
Alphabet of String:
    1. "a" = ("c")⁻¹
    2. "b"
    3. "c" = ("a")⁻¹

julia> A["a"]
1

julia> A[-A["a"]]
"c"
```
"""
Base.@propagate_inbounds function Base.getindex(A::Alphabet, p::Integer)
    @boundscheck checkbounds(letters(A), abs(p))
    if p > 0
        return @inbounds letters(A)[p]
    elseif p < 0 && hasinverse(-p, A)
        return @inbounds letters(A)[A.inversions[-p]]
    end

    throw(DomainError("Inversion of $(letters(A)[-p]) is not defined"))
end

function Base.inv(A::Alphabet, i::Integer)
    hasinverse(i, A) && return A.inversions[i]

    throw(DomainError(A[i], "is not invertible over $A"))
end

"""
    inv(A::Alphabet, w::AbstractWord)
Return the inverse of a word `w` in the context of alphabet `A`.
"""
function Base.inv(A::Alphabet, w::AbstractWord)
    res = similar(w)
    n = length(w)
    for (i, l) in zip(eachindex(res), Iterators.reverse(w))
        res[i] = inv(A, l)
    end
    return res
end

Base.inv(A::Alphabet{T}, a::T) where {T} = A[-A[a]]

function _print_syllable(io, symbol, pow)
    str = string(symbol)
    if length(str) > 3 && endswith(str, "^-1")
        print(io, first(str, length(str)-3), "^-", pow)
    else
        if pow == 1
            print(io, str)
        else
            print(io, str, "^", pow)
        end
    end
end

function print_repr(io::IO, w::AbstractWord, A::Alphabet, sep="*")
    if isone(w)
        print(io, w)
    else
        first_syllable = true
        idx = 1
        pow = 1
        while idx < length(w)
            if w[idx] == w[idx+1]
                pow += 1
                idx += 1
            else
                first_syllable || print(io, sep)
                _print_syllable(io, A[w[idx]], pow)
                first_syllable = false
                pow = 1
                idx += 1
            end
        end
        @assert idx == length(w)
        first_syllable || print(io, sep)
        _print_syllable(io, A[w[idx]], pow)
    end
end


text:
# include("../src/LinearRegression.jl")

@testset "from glm" begin

    fdf = DataFrame([[0.1,0.3,0.5,0.6,0.7,0.9],[0.086,0.269,0.446,0.538,0.626,0.782]], [:Carb, :OptDen])

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf)
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2
    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.ADJR2, 0.998808343507198)
    @test isapprox(lm1.AIC, -55.43694668654871) # using the SAS formula rather than the Julia-Statsmodel-GLM
    @test isapprox(lm1.stderrors,  [0.007833679251299831, 0.013534536322659505])
    @test isapprox(lm1.t_values, [0.6492114525712505, 64.74442074669666])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.ci_low, [-0.016664066127247305, 0.8387078171615459])
    @test isapprox(lm1.ci_up, [0.02683549469867456, 0.9138636114098853])
    @test isapprox(lm1.t_statistic, 2.7764451051977934)
    @test isapprox(lm1.VIF, [0.,  1.])

end

@testset "from glm regresspredict" begin

    fdf = DataFrame([[0.1,0.3,0.5,0.6,0.7,0.9],[0.086,0.269,0.446,0.538,0.626,0.782]], [:Carb, :OptDen])

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=["all"])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2

    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.ADJR2, 0.998808343507198)
    @test isapprox(lm1.AIC, -55.43694668654871) # using the SAS formula rather than the Julia-Statsmodel-GLM
    @test isapprox(lm1.stderrors,  [0.007833679251299831, 0.013534536322659505])
    @test isapprox(lm1.t_values, [0.6492114525712505, 64.74442074669666])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.ci_low, [-0.016664066127247305, 0.8387078171615459])
    @test isapprox(lm1.ci_up, [0.02683549469867456, 0.9138636114098853])
    @test isapprox(lm1.t_statistic, 2.7764451051977934)
    @test isapprox(lm1.VIF, [0.,  1.])

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=["none"])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=["r2", "P_values"])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2
    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])

    df = DataFrame(y=[1., 3., 3., 2., 2., 1.], x1=[1., 2., 3., 1., 2., 3.], x2=[1., 1., 1., -1., -1., -1.])
    lm2 = regress(@formula(y ~ 1 + x1 + x2), df, req_stats=["r2"])
    @test [1.5, 0.25, 0.3333333333333333] == lm2.coefs
    @test 0.22916666666666663 == lm2.R2

end

@testset "predictions statistics" begin
    t_carb = [0.1, 0.3, 0.5, 0.6, 0.7, 0.9]
    t_optden = [0.086, 0.269, 0.446, 0.538, 0.626, 0.782]
    t_leverage = [0.5918367346938774, 0.2816326530612245, 0.1673469387755102, 0.1836734693877552, 0.24897959183673485, 0.5265306122448984]
    t_predicted = [0.09271428571428536, 0.26797142857142836, 0.44322857142857136, 0.5308571428571429, 0.6184857142857143, 0.7937428571428573]
    t_residuals = [-0.006714285714285367, 0.0010285714285716563, 0.002771428571428647, 0.0071428571428571175, 0.007514285714285696, -0.011742857142857277]
    t_stdp = [0.0066535244611503905, 0.00458978457544483, 0.0035380151243903845, 0.003706585424647827, 0.004315515434962329, 0.006275706318490394]
    t_stdi = [0.01091189203370195, 0.009791124677432761, 0.009344386069745653, 0.009409504530540022, 0.009665592246181276, 0.010685714285717253]
    t_stdr = [0.00552545131594831, 0.00733033952495041, 0.007891923021648551, 0.007814168189246363, 0.0074950868260910365, 0.005951093194035971]
    t_student = [-1.2151560714878948, 0.1403170242074994, 0.3511727830880084, 0.9140905301586566, 1.0025615297914663, -1.9732268946192373]
    t_rstudent = [-1.3249515797718379, 0.12181828539462737, 0.3089239861926457, 0.8900235681358146, 1.003419757326401, -10.478921731163984]
    t_lcli = [0.06241801648886679, 0.24078690838638886, 0.4172843964641476, 0.5047321700609886, 0.5916497280049665, 0.7640745580187355]
    t_ucli = [0.12301055493970393, 0.2951559487564679, 0.4691727463929951, 0.5569821156532972, 0.6453217005664621, 0.8234111562669791]
    t_lclp = [0.07424114029181057, 0.2552281436530222, 0.43340546665434193, 0.520566011897882, 0.6065039225799076, 0.7763187030532258]
    t_uclp = [0.11118743113676015, 0.2807147134898345, 0.4530516762028008, 0.5411482738164038, 0.630467505991521, 0.8111670112324888]
    t_press = [-0.016449999999999146, 0.001431818181818499, 0.0033284313725491102, 0.00874999999999997, 0.010005434782608673, -0.02480172413793134]
    t_cooksd = [1.0705381016035724, 0.0038594654616162225, 0.012392684477580572, 0.09400066844914513, 0.16661116000566892, 2.1649894168822432]

    fdf = DataFrame([[0.1,0.3,0.5,0.6,0.7,0.9],[0.086,0.269,0.446,0.538,0.626,0.782]], [:Carb, :OptDen])
    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf)
    results = predict_and_stats(lm1, fdf, α= 0.05, req_stats=["all"])
    @test t_leverage == results.leverage
    @test t_predicted == results.predicted
    @test t_residuals == results.residuals
    @test t_stdp == results.stdp
    @test t_stdi == results.stdi
    @test t_stdr == results.stdr
    @test t_student == results.student
    @test t_rstudent == results.rstudent
    @test t_lcli == results.lcli
    @test t_ucli == results.ucli
    @test t_lclp == results.lclp
    @test t_uclp == results.uclp
    @test t_press == results.press
    @test t_cooksd == results.cooksd

    results = predict_and_stats(lm1, fdf, α= 0.05, req_stats= ["none"])
    @test t_predicted == results.predicted
    @test_throws ArgumentError("column name :leverage not found in the data frame") t_leverage == results.leverage

end


text:
function make_input_mms(opts0, degree; dg=false, operator="SBPOmega")
  opts = copy(opts0)
  opts["order"] = degree
  opts["IC_name"] = "ICp$degree"
  opts["BC1_name"] = string("p", degree, "BC")
  opts["SRCname"] = "SRCp$degree"
  opts["operator_type"] = operator

  if dg
    opts["use_DG"] = true
    opts["Flux_name"] = "LFFlux"
  else
    opts["use_DG"] = false
  end

  read_input(opts)  # get default values

  return opts
end

"""
  Test that the residual is zero for polynomials up to a specified degree

  **Inputs**

   * opts: options dictionary
   * maxdegree: maximum degree polynomials to test
"""
function test_mms_poly(opts::Dict, maxdegree::Integer)

  mesh, sbp, eqn, opts = createObjects(opts)

  for d=0:maxdegree
    println("testing degree ", d, " polynomials")

    # get IC
    AdvectionEquationMod.ICDict["ICp$d"](mesh, sbp, eqn, opts, eqn.q_vec)
    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

    # set BCs
    opts["BC1_name"] = string("p", d, "BC")
    AdvectionEquationMod.getBCFunctors(mesh, sbp, eqn, opts)

    # set source
    opts["SRCname"] = string("SRCp", d)
    AdvectionEquationMod.getSRCFunctors(mesh, sbp, eqn, opts)

    # compute residual
    fill!(eqn.res, 0.0)
    AdvectionEquationMod.evalResidual(mesh, sbp, eqn, opts)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    @test isapprox( calcNorm(eqn, eqn.res_vec), 0.0) atol=1e-12
  end

  return nothing
end

"""
  Test manufactured solutions, including cases that should be exact, for both
  gamma and omega operators.
"""
function test_mms()
  @testset "----- Testing using manufactured polynomials -----" begin

    # test 2d DG
    println("  ----- Testing 2d -----")
    opts0 = read_input_file("input_vals_mms.jl")
    operator_names = ["SBPOmega", "SBPGamma", "SBPDiagonalE", "SBPDiagonalE2",
                      "SBPOmega3"]
    max_degree =     [         4,          4,              4,              2,
                               4,]

    for i=1:length(operator_names)
      for degree=1:max_degree[i]
        println(  "  ----- testing degree $degree $(operator_names[i])")
        opts = make_input_mms(opts0, degree, dg=true, operator=operator_names[i])
        test_mms_poly(opts, degree)
      end
    end

    println("  ----- Testing 3d -----")
    opts0 = read_input_file("input_vals_3d.jl")
    operator_names = ["SBPOmega", "SBPGamma", "SBPDiagonalE"]
    max_degree =     [         4,          4,              4]

    for i=1:length(operator_names)
      for degree=1:max_degree[i]
        println(  "  ----- testing degree $degree $(operator_names[i])")
        opts = make_input_mms(opts0, degree, dg=true, operator=operator_names[i])
        test_mms_poly(opts, degree)
      end
    end

  end  # end facts block

  return nothing
end  # end function

#test_mms()
add_func1!(AdvectionTests, test_mms, [TAG_MMS, TAG_SHORTTEST, TAG_TMP])


text:
# Lazy concatenation of AbstractVector's.
# Similar to Iterators.Flatten and some code has been reused from julia/base/iterators.jl

# This fixes type inference issues in v1.0
for op in  (:hcat, :vcat)
    @eval call(::ApplyLayout{typeof($op)}, _) = $op
end

const Vcat{T,N,I<:Tuple} = ApplyArray{T,N,typeof(vcat),I}

Vcat(A...) = ApplyArray(vcat, A...)
Vcat{T}(A...) where T = ApplyArray{T}(vcat, A...)
Vcat() = Vcat{Any}()

Vcat(A::AbstractVector...) = ApplyVector(vcat, A...)
Vcat{T}(A::AbstractVector...) where T = ApplyVector{T}(vcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(vcat)})
    isempty(A.args) && return A
    m = size(A.args[1],2)
    for k=2:length(A.args)
        size(A.args[k],2) == m || throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(vcat)}) = promote_type(map(eltype,A.args)...)
@inline eltype(A::Applied{<:Any,typeof(vcat),Tuple{}}) = Any
@inline ndims(A::Applied{<:Any,typeof(vcat),I}) where I = max(1,maximum(map(ndims,A.args)))
@inline ndims(A::Applied{<:Any,typeof(vcat),Tuple{}}) = 1
@inline axes(f::Vcat{<:Any,1,Tuple{}}) = (OneTo(0),)
@inline axes(f::Vcat{<:Any,1}) = tuple(oneto(+(map(length,f.args)...)))
@inline axes(f::Vcat{<:Any,2}) = (oneto(+(map(a -> size(a,1), f.args)...)), axes(f.args[1],2))
Base.IndexStyle(::Type{<:Vcat{T,1}}) where T = Base.IndexLinear()
Base.IndexStyle(::Type{<:Vcat{T,2}}) where T = Base.IndexCartesian()

function ==(a::Vcat{T,N}, b::Vcat{T,N}) where {N,T}
    a_args = arguments(a)
    b_args = arguments(b)
    if length(a_args) ≠ length(b_args) || any(map(size,a_args) .≠ map(size,b_args))
        return Base.invoke(==, NTuple{2,AbstractArray}, a, b)
    end
    all(a_args .== b_args)
end

@propagate_inbounds @inline vcat_getindex(f, idx...) =
    vcat_getindex_recursive(f, idx, f.args...)

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer}, A, args...)
    k, = idx
    T = eltype(f)
    n = length(A)
    k ≤ n && return convert(T, A[k])::T
    vcat_getindex_recursive(f, (k - n, ), args...)
end

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer,Integer}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 1)
    k ≤ n && return convert(T, A[k, j])::T
    vcat_getindex_recursive(f, (k - n, j), args...)
end

@propagate_inbounds @inline function vcat_getindex_recursive(
        f, idx::Tuple{Integer,Union{Colon,AbstractVector}}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 1)
    k ≤ n && return convert(AbstractVector{T}, A[k, j])
    vcat_getindex_recursive(f, (k - n, j), args...)
end

@inline vcat_getindex_recursive(f, idx) = throw(BoundsError(f, idx))

@propagate_inbounds @inline getindex(f::Vcat{<:Any,1}, k::Integer) = vcat_getindex(f, k)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::Integer) = vcat_getindex(f, k, j)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::AbstractVector) = vcat_getindex(f, k, j)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::Colon) = vcat_getindex(f, k, j)

getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)

# since its mutable we need to make a copy
copy(f::Vcat) = Vcat(map(copy, f.args)...)
map(::typeof(copy), f::Vcat) = Vcat(map.(copy, f.args)...)

@propagate_inbounds @inline function vcat_setindex_recursive!(
        f::Vcat{T,1} where T, v, idx::NTuple{1}, A, args...)
    k, = idx
    n = length(A)
    k ≤ n && return setindex!(A, v, idx...)
    vcat_setindex_recursive!(f, v, (k - n, ), args...)
end

@propagate_inbounds @inline function vcat_setindex_recursive!(
        f::Vcat{T,2} where T, v, idx::NTuple{2}, A, args...)
    k, j = idx
    n = size(A, 1)
    k ≤ n && return setindex!(A, v, idx...)
    vcat_setindex_recursive!(f, v, (k - n, j), args...)
end

@inline vcat_setindex_recursive!(f, v, idx) = throw(BoundsError(f, idx))

@propagate_inbounds @inline function setindex!(
        f::Vcat{T,N}, v, idx::Vararg{Integer,N}) where {T,N}
    vcat_setindex_recursive!(f, v, idx, f.args...)
end

reverse(f::Vcat{<:Any,1}) = Vcat((reverse(itr) for itr in reverse(f.args))...)

const Hcat{T,I<:Tuple} = ApplyArray{T,2,typeof(hcat),I}

Hcat(A...) = ApplyArray(hcat, A...)
Hcat{T}(A...) where T = ApplyArray{T}(hcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(hcat)})
    isempty(A.args) && return A
    m = size(A.args[1],1)
    for k=2:length(A.args)
        size(A.args[k],1) == m || throw(ArgumentError("number of rows of each array must match (got $(map(x->size(x,1), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(hcat)}) = promote_type(map(eltype,A.args)...)
ndims(::Applied{<:Any,typeof(hcat)}) = 2
size(f::Applied{<:Any,typeof(hcat)}) = (size(f.args[1],1), +(map(a -> size(a,2), f.args)...))
Base.IndexStyle(::Type{<:Hcat}) where T = Base.IndexCartesian()

@inline hcat_getindex(f, k, j::Integer) =
    hcat_getindex_recursive(f, (k, j), f.args...)

@inline function hcat_getindex_recursive(
        f, idx::Tuple{Integer,Integer}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return convert(T, A[k, j])::T
    hcat_getindex_recursive(f, (k, j - n), args...)
end

@inline function hcat_getindex_recursive(
        f, idx::Tuple{Union{Colon,AbstractVector},Integer}, A, args...)
    kr, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return convert(AbstractVector{T}, A[kr, j])
    hcat_getindex_recursive(f, (kr, j - n), args...)
end

@inline hcat_getindex_recursive(f, idx) = throw(BoundsError(f, idx))

getindex(f::Hcat, k::Integer, j::Integer) = hcat_getindex(f, k, j)
getindex(f::Hcat, k::AbstractVector, j::Integer) = hcat_getindex(f, k, j)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)

# since its mutable we need to make a copy
copy(f::Hcat) = Hcat(map(copy, f.args)...)

@inline function hcat_setindex_recursive!(
        f, v, idx::NTuple{2}, A, args...)
    k, j = idx
    T = eltype(f)
    n = size(A, 2)
    j ≤ n && return setindex!(A, v, k, j)
    hcat_setindex_recursive!(f, v, (k, j - n), args...)
end

@inline hcat_setindex_recursive!(f, v, idx) = throw(BoundsError(f, idx))

function setindex!(f::Hcat{T}, v, k::Integer, j::Integer) where T
    hcat_setindex_recursive!(f, v, (k, j), f.args...)
end

## copyto!
# based on Base/array.jl, Base/abstractarray.jl
_copyto!(_, LAY::ApplyLayout{typeof(vcat)}, dest::AbstractArray{<:Any,N}, V::AbstractArray{<:Any,N}) where N =
    vcat_copyto!(dest, arguments(LAY, V)...)
function vcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest,1)
    nrows == sum(a->size(a, 1), arrays) || throw(DimensionMismatch("sum of rows each matrix must equal $nrows"))
    ncols = size(dest, 2)
    for a in arrays
        if size(a, 2) != ncols
            throw(DimensionMismatch("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    pos = 1
    for a in arrays
        p1 = pos+size(a,1)-1
        copyto!(view(dest,pos:p1, :), a)
        pos = p1+1
    end
    return dest
end

function vcat_copyto!(arr::AbstractVector, arrays...)
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    i = 0
    @inbounds for a in arrays
        m = length(a)
        copyto!(view(arr,i+1:i+m), a)
        i += m
    end
    arr
end

function vcat_copyto!(arr::Vector{T}, arrays::Vector{T}...) where T
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    ptr = pointer(arr)
    if isbitstype(T)
        elsz = Core.sizeof(T)
    elseif isbitsunion(T)
        elsz = bitsunionsize(T)
        selptr = convert(Ptr{UInt8}, ptr) + n * elsz
    else
        elsz = Core.sizeof(Ptr{Cvoid})
    end
    t = @_gc_preserve_begin arr
    for a in arrays
        na = length(a)
        nba = na * elsz
        if isbitstype(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
        elseif isbitsunion(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
            # copy selector bytes
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  selptr, convert(Ptr{UInt8}, pointer(a)) + nba, na)
            selptr += na
        else
            ccall(:jl_array_ptr_copy, Cvoid, (Any, Ptr{Cvoid}, Any, Ptr{Cvoid}, Int),
                  arr, ptr, a, pointer(a), na)
        end
        ptr += nba
    end
    @_gc_preserve_end t
    return arr
end

# special case for adjoints of hcat. This is useful for catching fast paths
# for vector case, e.g., _fast_blockbradcast_copyto! in BlockArrays.jl
function vcat_copyto!(dest::AbstractMatrix, arrays::Adjoint{<:Any,<:AbstractVector}...)
    hcat_copyto!(dest', map(adjoint, arrays)...)
    dest
end

_copyto!(_, LAY::ApplyLayout{typeof(hcat)}, dest::AbstractMatrix, H::AbstractMatrix) =
    hcat_copyto!(dest, arguments(LAY,H)...)
function hcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest, 1)
    ncols = 0
    dense = true
    for a in arrays
        dense &= isa(a,Array)
        nd = ndims(a)
        ncols += (nd==2 ? size(a,2) : 1)
    end

    nrows == size(first(arrays),1) || throw(DimensionMismatch("Destination rows must match"))
    ncols == size(dest,2) || throw(DimensionMismatch("Destination columns must match"))

    pos = 1
    if dense
        for a in arrays
            n = length(a)
            copyto!(dest, pos, a, 1, n)
            pos += n
        end
    else
        for a in arrays
            p1 = pos+(isa(a,AbstractMatrix) ? size(a, 2) : 1)-1
            copyto!(view(dest,:, pos:p1), a)
            pos = p1+1
        end
    end
    return dest
end

function hcat_copyto!(dest::AbstractMatrix, arrays::AbstractVector...)
    height = size(dest, 1)
    for j = 1:length(arrays)
        if length(arrays[j]) != height
            throw(DimensionMismatch("vectors must have same lengths"))
        end
    end
    for j=1:length(arrays)
        copyto!(view(dest,:,j), arrays[j])
    end

    dest
end

#####
# adjoint/transpose
#####

for adj in (:adjoint, :transpose)
    @eval begin
        $adj(A::Hcat{T}) where T = Vcat{T}(map($adj,A.args)...)
        $adj(A::Vcat{T,2}) where T = Hcat{T}(map($adj,A.args)...)
    end
end

_vec(a) = a
_vec(a::AbstractArray) = vec(a)
_vec(a::Adjoint{<:Number,<:AbstractVector}) = _vec(parent(a))
vec(A::Hcat) = Vcat(map(_vec,A.args)...)

copy(f::Adjoint{<:Any,<:Union{Vcat,Hcat}}) = copy(parent(f))'
copy(f::Transpose{<:Any,<:Union{Vcat,Hcat}}) = transpose(copy(parent(f)))

_permutedims(a) = a
_permutedims(a::AbstractArray) = permutedims(a)

permutedims(A::Hcat{T}) where T = Vcat{T}(map(_permutedims,A.args)...)
permutedims(A::Vcat{T}) where T = Hcat{T}(map(_permutedims,A.args)...)

transposelayout(::ApplyLayout{typeof(vcat)}) = ApplyLayout{typeof(hcat)}()
transposelayout(::ApplyLayout{typeof(hcat)}) = ApplyLayout{typeof(vcat)}()

arguments(::ApplyLayout{typeof(vcat)}, A::Adjoint) = map(adjoint, arguments(ApplyLayout{typeof(hcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(hcat)}, A::Adjoint) = map(adjoint, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(vcat)}, A::Transpose) = map(transpose, arguments(ApplyLayout{typeof(hcat)}(), parent(A)))
arguments(::ApplyLayout{typeof(hcat)}, A::Transpose) = map(transpose, arguments(ApplyLayout{typeof(vcat)}(), parent(A)))

#####
# broadcasting
#
# We want broadcasting for numbers with concaenations to pass through
# to take advantage of special implementations of the sub-components
######

BroadcastStyle(::Type{<:Vcat{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(::Type{<:Hcat{<:Any}}) where N = LazyArrayStyle{2}()

broadcasted(::LazyArrayStyle, op, A::Vcat) =
    Vcat(broadcast(x -> broadcast(op, x), A.args)...)

for Cat in (:Vcat, :Hcat)
    @eval begin
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Number) =
            $Cat(broadcast((x,y) -> broadcast(op, x, y), A.args, c)...)
        broadcasted(::LazyArrayStyle, op, c::Number, A::$Cat) =
            $Cat(broadcast((x,y) -> broadcast(op, x, y), c, A.args)...)
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Ref) =
            $Cat(broadcast((x,y) -> broadcast(op, x, Ref(y)), A.args, c)...)
        broadcasted(::LazyArrayStyle, op, c::Ref, A::$Cat) =
            $Cat(broadcast((x,y) -> broadcast(op, Ref(x), y), c, A.args)...)
    end
end

# determine indices of components of a vcat
_vcat_axes(::Tuple{}) = (1,)
_vcat_axes(a::Tuple{<:AbstractUnitRange}) = (first(a),)
_vcat_axes(::Tuple{}, b, c...) = tuple(1, broadcast(x -> broadcast(+, 1, x), _vcat_axes(b, c...))...)
_vcat_axes(a::Tuple{<:AbstractUnitRange}, b, c...) = tuple(first(a), broadcast((α,x) -> broadcast(+, α, x), last(first(a)),
                                                            _vcat_axes(b, c...))...)

_vcat_getindex_eval(y) = ()
_vcat_getindex_eval(y, a, b...) = tuple(y[a], _vcat_getindex_eval(y, b...)...)

# let it stay lazy
layout_broadcasted(::ApplyLayout{typeof(vcat)}, ::AbstractLazyLayout, op, A::AbstractVector, B::AbstractVector) =
    Broadcasted{LazyArrayStyle{1}}(op, (A, B))
layout_broadcasted(::AbstractLazyLayout, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    Broadcasted{LazyArrayStyle{1}}(op, (A, B))

function layout_broadcasted(::ApplyLayout{typeof(vcat)}, _, op, A::AbstractVector, B::AbstractVector)
    kr = _vcat_axes(map(axes,A.args)...)  # determine how to break up B
    B_arrays = _vcat_getindex_eval(B,kr...)    # evaluate B at same chunks as A
    ApplyVector(vcat, broadcast((a,b) -> broadcast(op,a,b), A.args, B_arrays)...)
end

function layout_broadcasted(_, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector)
    kr = _vcat_axes(axes.(B.args)...)
    A_arrays = _vcat_getindex_eval(A,kr...)
    Vcat(broadcast((a,b) -> broadcast(op,a,b), A_arrays, B.args)...)
end

layout_broadcasted(::ApplyLayout{typeof(vcat)}, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    Broadcasted{LazyArrayStyle{1}}(op, (A, B))

layout_broadcasted(::ApplyLayout{typeof(vcat)}, Blay::CachedLayout, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(UnknownLayout(), Blay, op, A, B)
layout_broadcasted(Alay::CachedLayout, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(Alay, UnknownLayout(), op, A, B)

broadcasted(::LazyArrayStyle, op, a::Vcat{<:Any,1}, b::AbstractVector) = layout_broadcasted(op, a, b)
broadcasted(::LazyArrayStyle, op, a::AbstractVector, b::Vcat{<:Any,1}) = layout_broadcasted(op, a, b)
broadcasted(::LazyArrayStyle{1}, op, a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, op, a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), op, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(\), a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), \, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(/), a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), /, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Vcat{<:Any,1}, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Zeros{<:Any,1}, b::Vcat{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)

# Cannot broadcast Vcat's in a lazy way so stick to BroadcastArray
broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::Vcat{<:Any,1}) = layout_broadcasted(op, A, B)

# ambiguities
broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::CachedVector) = layout_broadcasted(op, A, B)
broadcasted(::LazyArrayStyle, op, A::CachedVector, B::Vcat{<:Any,1}) = layout_broadcasted(op, A, B)

function +(A::Vcat, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::Vcat, B::AbstractArray)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::AbstractArray, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end

####
# Cumsum
####

sum(V::Vcat) = mapreduce(sum, +, V.args)
all(V::Vcat) = all(all.(V.args))
any(V::Vcat) = any(any.(V.args))
all(f::Function, V::Vcat) = all(all.(f, V.args))
any(f::Function, V::Vcat) = any(any.(f, V.args))

_dotplus(a,b) = broadcast(+, a, b)

@inline _cumsum(x::Number) = x
@inline _cumsum(x) = cumsum(x)
_cumsum_last(x::AbstractVector{T}) where T = isempty(x) ? zero(T) : last(x)
_cumsum_last(x) = last(x)

_tuple_cumsum() = ()
_tuple_cumsum(a) = (a,)
_tuple_cumsum(a, b...) = (a, broadcast(+,a,_tuple_cumsum(b...))...)
function _vcat_cumsum(x...)
    cs = map(_cumsum,x)
    cslasts = tuple(0,_tuple_cumsum(map(_cumsum_last,most(cs))...)...)
    map((a,b) -> broadcast(+,a,b), cslasts, cs)
end

@inline cumsum(V::Vcat{<:Any,1}) = ApplyVector(vcat,_vcat_cumsum(V.args...)...)

_vcat_diff(x::Number) = ()
_vcat_diff(x) = (diff(x),)

_vcat_diff(a::Number, b, c...) = (first(b)-a, _vcat_diff(b,c...)...)
_vcat_diff(a, b, c...) = (diff(a), first(b)-last(a), _vcat_diff(b,c...)...)
@inline diff(V::Vcat{T,1}) where T = ApplyVector{T}(vcat,_vcat_diff(V.args...)...)

####
# maximum/minimum
####

for op in (:maximum, :minimum)
    @eval $op(V::Vcat) = $op($op.(V.args))
end

function in(x, V::Vcat)
    for a in V.args
        in(x, a) && return true
    end
    false
end

_fill!(a, x) = fill!(a,x)
function _fill!(a::Number, x)
    a == x || throw(ArgumentError("Cannot set $a to $x"))
    a
end

function fill!(V::Union{Vcat,Hcat}, x)
    for a in V.args
        _fill!(a, x)
    end
    V
end

###
# *
###

function materialize!(M::MatMulVecAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 function materialize!(M::MatMulMatAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 ####
 # col/rowsupport
 ####

most(a) = reverse(tail(reverse(a)))
colsupport(M::Vcat, j) = first(colsupport(first(M.args),j)):(size(Vcat(most(M.args)...),1)+last(colsupport(last(M.args),j)))

function rowsupport(V::Vcat, k::Integer)
    ξ = k
    for A in arguments(V)
        n = size(A,1)
        ξ ≤ n && return rowsupport(A, ξ)
        ξ -= n
    end
    return 1:0
end

function colsupport(H::Hcat, j::Integer)
    ξ = j
    for A in arguments(H)
        n = size(A,2)
        ξ ≤ n && return colsupport(A, ξ)
        ξ -= n
    end
    return 1:0
end

rowsupport(M::Hcat, k) = first(rowsupport(first(M.args),k)):(size(Hcat(most(M.args)...),2)+last(rowsupport(last(M.args),k)))

###
# padded
####

struct PaddedLayout{L} <: MemoryLayout end
applylayout(::Type{typeof(vcat)}, ::A, ::ZerosLayout) where A = PaddedLayout{A}()
applylayout(::Type{typeof(vcat)}, ::ScalarLayout, ::ScalarLayout, ::ZerosLayout) = PaddedLayout{ApplyLayout{typeof(vcat)}}()
applylayout(::Type{typeof(vcat)}, ::A, ::PaddedLayout) where A = PaddedLayout{ApplyLayout{typeof(vcat)}}()
applylayout(::Type{typeof(vcat)}, ::ScalarLayout, ::ScalarLayout, ::PaddedLayout) = PaddedLayout{ApplyLayout{typeof(vcat)}}()
cachedlayout(::A, ::ZerosLayout) where A = PaddedLayout{A}()
sublayout(::PaddedLayout{Lay}, sl::Type{<:Tuple{Slice,Integer}}) where Lay =
    PaddedLayout{typeof(sublayout(Lay(), sl))}()

paddeddata(A::CachedArray{<:Any,N,<:Any,<:Zeros}) where N = cacheddata(A)
_vcat_paddeddata(A, B::Zeros) = A
_vcat_paddeddata(A, B) = Vcat(A, paddeddata(B))
_vcat_paddeddata(A, B, C...) = Vcat(A, _vcat_paddeddata(B, C...))
paddeddata(A::Vcat) = _vcat_paddeddata(A.args...)

colsupport(::PaddedLayout, A, j) = colsupport(paddeddata(A),j)

function _vcat_resizedata!(::PaddedLayout, B, m...)
    Base.checkbounds(paddeddata(B), m...)
    B
end

resizedata!(B::Vcat, m...) = _vcat_resizedata!(MemoryLayout(B), B, m...)

function ==(A::CachedVector{<:Any,<:Any,<:Zeros}, B::CachedVector{<:Any,<:Any,<:Zeros})
    length(A) == length(B) || return false
    n = max(A.datasize[1], B.datasize[1])
    resizedata!(A,n); resizedata!(B,n)
    view(A.data,OneTo(n)) == view(B.data,OneTo(n))
end

# special copyto! since `similar` of a padded returns a cached
for Typ in (:Number, :AbstractVector)
    @eval begin
        function _copyto!(::PaddedLayout, ::PaddedLayout, dest::CachedVector{<:$Typ}, src::AbstractVector)
            length(src) ≤ length(dest)  || throw(BoundsError())
            a = paddeddata(src)
            n = length(a)
            resizedata!(dest, n) # make sure we are padded enough
            copyto!(view(dest.data,OneTo(n)), a)
            dest
        end
        function _copyto!(::PaddedLayout, ::PaddedLayout, dest::SubArray{<:Any,1,<:CachedVector{<:$Typ}}, src::AbstractVector)
            length(src) ≤ length(dest)  || throw(BoundsError())
            a = paddeddata(src)
            n = length(a)
            k = first(parentindices(dest)[1])
            resizedata!(parent(dest), k+n-1) # make sure we are padded enough
            copyto!(view(parent(dest).data,k:k+n-1), a)
            dest
        end

        _copyto!(::PaddedLayout, ::PaddedLayout, dest::CachedVector{<:$Typ}, src::CachedVector) =
            _padded_copyto!(dest, src)
    end
end

function _padded_copyto!(dest::CachedVector, src::CachedVector)
    length(src) ≤ length(dest)  || throw(BoundsError())
    n = src.datasize[1]
    resizedata!(dest, n)
    copyto!(view(dest.data,OneTo(n)), view(src.data,OneTo(n)))
    dest
end

_copyto!(::PaddedLayout, ::PaddedLayout, dest::CachedVector, src::CachedVector) =
    _padded_copyto!(dest, src)

function _copyto!(::PaddedLayout, ::ZerosLayout, dest::AbstractVector, src::AbstractVector)
    zero!(paddeddata(dest))
    dest
end

######
# Special Vcat broadcasts
#
# We use Vcat for infinite padded vectors, so we need to special case
# two arrays. This may be generalisable in the future
######

function broadcasted(::LazyArrayStyle{1}, op, A::Vcat{<:Any,1,<:Tuple{AbstractVector,AbstractFill}},
                                              B::Vcat{<:Any,1,<:Tuple{AbstractVector,AbstractFill}})
    (Ahead, Atail) = A.args
    (Bhead, Btail) = B.args
    T = promote_type(eltype(A), eltype(B))

    if length(Ahead) ≥ length(Bhead)
        M,m = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(view(Ahead,1:m), Bhead)
        view(Chead,m+1:M) .= op.(view(Ahead,m+1:M),Btail[1:M-m])

        Ctail = op.(Atail, Btail[M-m+1:end])
    else
        m,M = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(Ahead, view(Bhead,1:m))
        view(Chead,m+1:M) .= op.(Atail[1:M-m],view(Bhead,m+1:M))

        Ctail = op.(Atail[M-m+1:end], Btail)
    end

    Vcat(Chead, Ctail)
end

function broadcasted(::LazyArrayStyle{1}, op, A::Vcat{<:Any,1,<:Tuple{Number,AbstractFill}},
                                              B::Vcat{<:Any,1,<:Tuple{Number,AbstractFill}})
    (Ahead, Atail) = A.args
    (Bhead, Btail) = B.args
    Vcat(op.(Ahead,Bhead), op.(Atail,Btail))
end

function broadcasted(::LazyArrayStyle{1}, op, A::Vcat{<:Any,1,<:Tuple{SVector{M},AbstractFill}},
                                              B::Vcat{<:Any,1,<:Tuple{SVector{M},AbstractFill}}) where M
    (Ahead, Atail) = A.args
    (Bhead, Btail) = B.args
    Vcat(op.(Ahead,Bhead), op.(Atail,Btail))
end

# special case handle broadcasting with padded and cached arrays
function layout_broadcasted(::PaddedLayout, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector)
    a,b = paddeddata(A),paddeddata(B)
    n,m = length(a),length(b)
    dat = if n ≤ m
        [broadcast(op, a, view(b,1:n)); broadcast(op, zero(eltype(A)), @view(b[n+1:end]))]
    else
        [broadcast(op, view(a,1:m), b); broadcast(op, @view(a[m+1:end]), zero(eltype(B)))]
    end
    CachedArray(dat, broadcast(op, Zeros{eltype(A)}(length(A)), Zeros{eltype(B)}(length(B))))
end

function layout_broadcasted(_, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector)
    b = paddeddata(B)
    m = length(b)
    zB = Zeros{eltype(B)}(size(B)...)
    CachedArray(convert(Array,broadcast(op, view(A,1:m), b)), broadcast(op, A, zB))
end

function layout_broadcasted(::PaddedLayout, _, op, A::AbstractVector, B::AbstractVector)
    a = paddeddata(A)
    n = length(a)
    zA = Zeros{eltype(A)}(size(A)...)
    CachedArray(convert(Array,broadcast(op, a, view(B,1:n))), broadcast(op, zA, B))
end

function layout_broadcasted(::PaddedLayout, ::CachedLayout, op, A::AbstractVector, B::AbstractVector)
    a = paddeddata(A)
    n = length(a)
    resizedata!(B,n)
    Bdata = cacheddata(B)
    b = view(Bdata,1:n)
    zA1 = Zeros{eltype(A)}(size(Bdata,1)-n)
    zA = Zeros{eltype(A)}(size(A)...)
    CachedArray([broadcast(op, a, b); broadcast(op, zA1, @view(Bdata[n+1:end]))], broadcast(op, zA, B.array))
end

function layout_broadcasted(::CachedLayout, ::PaddedLayout, op, A::AbstractVector, B::AbstractVector)
    b = paddeddata(B)
    n = length(b)
    resizedata!(A,n)
    Adata = cacheddata(A)
    a = view(Adata,1:n)
    zB1 = Zeros{eltype(B)}(size(Adata,1)-n)
    zB = Zeros{eltype(B)}(size(B)...)
    CachedArray([broadcast(op, a, b); broadcast(op, @view(Adata[n+1:end]), zB1)], broadcast(op, A.array, zB))
end

layout_broadcasted(::ApplyLayout{typeof(vcat)}, lay::PaddedLayout, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(UnknownLayout(), lay, op, A, B)
layout_broadcasted(lay::PaddedLayout, ::ApplyLayout{typeof(vcat)}, op, A::AbstractVector, B::AbstractVector) =
    layout_broadcasted(lay, UnknownLayout(), op, A, B)

###
# Dot/Axpy
###

struct Axpy{StyleX,StyleY,T,XTyp,YTyp}
    α::T
    X::XTyp
    Y::YTyp
end

Axpy(α::T, X::XTyp, Y::YTyp) where {T,XTyp,YTyp} = Axpy{typeof(MemoryLayout(XTyp)), typeof(MemoryLayout(YTyp)), T, XTyp, YTyp}(α, X, Y)
materialize!(d::Axpy{<:Any,<:Any,<:Number,<:AbstractArray,<:AbstractArray}) = Base.invoke(BLAS.axpy!, Tuple{Number,AbstractArray,AbstractArray}, d.α, d.X, d.Y)
function materialize!(d::Axpy{<:PaddedLayout,<:PaddedLayout,U,<:AbstractVector{T},<:AbstractVector{V}}) where {U,T,V}
    x = paddeddata(d.X)
    resizedata!(d.Y, length(x))
    y = paddeddata(d.Y)
    BLAS.axpy!(d.α, x, view(y,1:length(x)))
    y
end
axpy!(α, X, Y) = materialize!(Axpy(α,X,Y))
BLAS.axpy!(α, X::LazyArray, Y::AbstractArray) = axpy!(α, X, Y)
BLAS.axpy!(α, X::SubArray{<:Any,N,<:LazyArray}, Y::AbstractArray) where N = axpy!(α, X, Y)

###
# l/rmul!
###

function materialize!(M::Lmul{ScalarLayout,<:PaddedLayout})
    lmul!(M.A, paddeddata(M.B))
    M.B
end

function materialize!(M::Rmul{<:PaddedLayout,ScalarLayout})
    rmul!(paddeddata(M.A), M.B)
    M.A
end

###
# norm
###

for Cat in (:Vcat, :Hcat)
    for (op,p) in ((:norm1,1), (:norm2,2), (:normInf,Inf))
        @eval $op(a::$Cat) = $op(norm.(a.args,$p))
    end
    @eval normp(a::$Cat, p) = norm(norm.(a.args, p), p)
end

_norm2(::PaddedLayout, a) = norm(paddeddata(a),2)
_norm1(::PaddedLayout, a) = norm(paddeddata(a),1)
_normInf(::PaddedLayout, a) = norm(paddeddata(a),Inf)
_normp(::PaddedLayout, a, p) = norm(paddeddata(a),p)

function copy(D::Dot{<:PaddedLayout, <:PaddedLayout})
    a = paddeddata(D.A)
    b = paddeddata(D.B)
    m = min(length(a), length(b))
    convert(eltype(D), dot(view(a, 1:m), view(b, 1:m)))
end

function copy(D::Dot{<:PaddedLayout})
    a = paddeddata(D.A)
    m = length(a)
    convert(eltype(D), dot(a, view(D.B, 1:m)))
end

function copy(D::Dot{<:Any, <:PaddedLayout})
    b = paddeddata(D.B)
    m = length(b)
    convert(eltype(D), dot(view(D.A, 1:m), b))
end

###
# subarrays
###

sublayout(::ApplyLayout{typeof(vcat)}, _) = ApplyLayout{typeof(vcat)}()
sublayout(::ApplyLayout{typeof(hcat)}, _) = ApplyLayout{typeof(hcat)}()
# a row-slice of an Hcat is equivalent to a Vcat
sublayout(::ApplyLayout{typeof(hcat)}, ::Type{<:Tuple{Number,AbstractVector}}) = ApplyLayout{typeof(vcat)}()

_vcat_lastinds(sz) = _vcat_cumsum(sz...)
_vcat_firstinds(sz) = (1, (1 .+ most(_vcat_lastinds(sz)))...)

_argsindices(sz) = broadcast(:, _vcat_firstinds(sz), _vcat_lastinds(sz))

_view_vcat(a::Number, kr) = Fill(a,length(kr))
_view_vcat(a::Number, kr, jr) = Fill(a,length(kr), length(jr))
_view_vcat(a, kr...) = _viewifmutable(a, kr...)

function _vcat_sub_arguments(::ApplyLayout{typeof(vcat)}, A, V, kr)
    sz = size.(arguments(A),1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    map(_view_vcat, arguments(A), skr2)
end

function _vcat_sub_arguments(::ApplyLayout{typeof(vcat)}, A, V, kr, jr)
    sz = size.(arguments(A),1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(arguments(A), skr2, Ref(jr))
end

_vcat_sub_arguments(LAY::ApplyLayout{typeof(vcat)}, A, V) = _vcat_sub_arguments(LAY, A, V, parentindices(V)...)

function _vcat_sub_arguments(L::ApplyLayout{typeof(hcat)}, A, V)
    A = parent(V)
    args = arguments(L, A)
    k,jr = parentindices(V)
    sz = size.(args,2)
    sjr = intersect.(_argsindices(sz), Ref(jr))
    sjr2 = broadcast((a,b) -> a .- b .+ 1, sjr, _vcat_firstinds(sz))
    _view_hcat.(args, k, sjr2)
end

_vcat_sub_arguments(::PaddedLayout, A, V) = _vcat_sub_arguments(ApplyLayout{typeof(vcat)}(), A, V)
_vcat_sub_arguments(::DualLayout{ML}, A, V) where ML = _vcat_sub_arguments(ML(), A, V)
_vcat_sub_arguments(A, V) = _vcat_sub_arguments(MemoryLayout(typeof(A)), A, V)
arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,1}) = _vcat_sub_arguments(parent(V), V)

function arguments(L::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2})
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(args, skr2, Ref(jr))
end

@inline _view_hcat(a::Number, kr, jr) = Fill(a,length(kr),length(jr))
@inline _view_hcat(a::Number, kr::Number, jr) = Fill(a,length(jr))
@inline _view_hcat(a, kr, jr) = _viewifmutable(a, kr, jr)
@inline _view_hcat(a::AbstractVector, kr, jr::Colon) = _viewifmutable(a, kr)

# equivalent to broadcast but written to maintain type stability
__view_hcat(::Tuple{}, _, ::Tuple{}) = ()
__view_hcat(::Tuple{}, _, ::Colon) = ()
@inline __view_hcat(args::Tuple, kr, jrs::Tuple) = (_view_hcat(args[1], kr, jrs[1]), __view_hcat(tail(args), kr, tail(jrs))...)
@inline __view_hcat(args::Tuple, kr, ::Colon) = (_view_hcat(args[1], kr, :), __view_hcat(tail(args), kr, :)...)

function arguments(L::ApplyLayout{typeof(hcat)}, V::SubArray)
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,2)
    sjr = intersect.(_argsindices(sz), Ref(jr))
    sjr2 = broadcast((a,b) -> a .- b .+ 1, sjr, _vcat_firstinds(sz))
    __view_hcat(args, kr, sjr2)
end

arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Slice,<:Any}}) =
    _viewifmutable.(arguments(parent(V)), Ref(:), Ref(parentindices(V)[2]))
arguments(::ApplyLayout{typeof(hcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) =
    __view_hcat(arguments(parent(V)), parentindices(V)[1], :)

function sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, _)
    ret = similar(V)
    n = 0
    _,jr = parentindices(V)
    for a in arguments(V)
        m = size(a,1)
        copyto!(view(ret,n+1:n+m,:), a)
        n += m
    end
    ret
end

sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractVector, _) = ApplyVector(V)

function sub_materialize(::ApplyLayout{typeof(hcat)}, V, _)
    ret = similar(V)
    n = 0
    kr,_ = parentindices(V)
    for a in arguments(V)
        m = size(a,2)
        copyto!(view(ret,:,n+1:n+m), a)
        n += m
    end
    ret
end

# temporarily allocate. In the future, we add a loop over arguments
materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))
materialize!(M::MatMulVecAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))

sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()
sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange,AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()

_lazy_getindex(dat, kr...) = view(dat, kr...)
_lazy_getindex(dat::Number, _) = dat

function sub_paddeddata(_, S::SubArray{<:Any,1,<:AbstractVector})
    dat = paddeddata(parent(S))
    (kr,) = parentindices(S)
    kr2 = kr ∩ axes(dat,1)
    _lazy_getindex(dat, kr2)
end

function sub_paddeddata(_, S::SubArray)
    P = parent(S)
    (kr,jr) = parentindices(S)
    resizedata!(P, 1, maximum(jr)) # ensure enough rows
    dat = paddeddata(P)
    kr2 = kr ∩ axes(dat,1)
    _lazy_getindex(dat, kr2, jr)
end

paddeddata(S::SubArray) = sub_paddeddata(MemoryLayout(parent(S)), S)

function _padded_sub_materialize(v::AbstractVector{T}) where T
    dat = paddeddata(v)
    Vcat(dat, Zeros{T}(length(v) - length(dat)))
end

sub_materialize(::PaddedLayout, v::AbstractVector{T}, _) where T =
    _padded_sub_materialize(v)

function sub_materialize(::PaddedLayout, v::AbstractMatrix{T}, _) where T
    dat = paddeddata(v)
    Vcat(dat, Zeros{T}(size(v,1) - size(dat,1), size(dat,2)))
end

## print

_replace_in_print_matrix(A::AbstractArray, k, j, s) = replace_in_print_matrix(A, k, j, s)
_replace_in_print_matrix(_, k, j, s) = s

function layout_replace_in_print_matrix(LAY::ApplyLayout{typeof(vcat)}, f::AbstractVecOrMat, k, j, s)
    κ = k
    for A in arguments(LAY, f)
        n = size(A,1)
        κ ≤ n && return _replace_in_print_matrix(A, κ, j, s)
        κ -= n
    end
    throw(BoundsError(f, (k,j)))
end

function layout_replace_in_print_matrix(::PaddedLayout, f::AbstractVecOrMat, k, j, s)
    data = paddeddata(f)
    k in axes(data,1) ? _replace_in_print_matrix(data, k, j, s) : Base.replace_with_centered_mark(s)
end

# searchsorted

_searchsortedfirst(a, x) = searchsortedfirst(a, x)
_searchsortedfirst(a::Number, x) = 1 + (x > a)
_searchsortedlast(a, x) = searchsortedlast(a, x)
_searchsortedlast(a::Number, x) = 0 + (x ≥ a)

searchsortedfirst(f::Vcat{<:Any,1}, x) =
    searchsortedfirst_recursive(0, x, arguments(f)...)

searchsortedlast(f::Vcat{<:Any,1}, x) =
    searchsortedlast_recursive(length(f), x, reverse(arguments(f))...)

@inline searchsortedfirst_recursive(n, x) = n + 1

@inline function searchsortedfirst_recursive(n, x, a, args...)
    m = length(a)
    r = _searchsortedfirst(a, x)
    r ≤ m && return n + r
    return searchsortedfirst_recursive(n + m, x, args...)
end

@inline searchsortedlast_recursive(n, x) = n

@inline function searchsortedlast_recursive(n, x, a, args...)
    n -= length(a)
    r = _searchsortedlast(a, x)
    r > 0 && return n + r
    return searchsortedlast_recursive(n, x, args...)
end

searchsorted(f::Vcat{<:Any,1}, x) = searchsortedfirst(f, x):searchsortedlast(f,x)

# avoid ambiguity in LazyBandedMatrices
copy(M::Mul{<:DiagonalLayout,<:PaddedLayout}) = copy(Lmul(M))

# Triangular columns

sublayout(::TriangularLayout{'U','N', ML}, INDS::Type{<:Tuple{KR,Integer}}) where {KR,ML} =
    sublayout(PaddedLayout{typeof(sublayout(ML(), INDS))}(), Tuple{KR})

sublayout(::TriangularLayout{'L','N', ML}, INDS::Type{<:Tuple{Integer,JR}}) where {JR,ML} =
    sublayout(PaddedLayout{typeof(sublayout(ML(), INDS))}(), Tuple{JR})

resizedata!(A::UpperTriangular, k::Integer, j::Integer) = resizedata!(parent(A), min(k,j), j)
resizedata!(A::LowerTriangular, k::Integer, j::Integer) = resizedata!(parent(A), k, min(k,j))

function sub_paddeddata(::TriangularLayout{'U','N'}, S::SubArray{<:Any,1,<:AbstractMatrix})
    P = parent(S)
    (kr,j) = parentindices(S)
    view(triangulardata(P), kr ∩ (1:j), j)
end

function sub_paddeddata(::TriangularLayout{'L','N'}, S::SubArray{<:Any,1,<:AbstractMatrix})
    P = parent(S)
    (k,jr) = parentindices(S)
    view(triangulardata(P), k, jr ∩ (1:k))
end


text:
# Plots.jl recipes

# plot config object, used for dispatch
mutable struct ParetoShapePlotConfig
    args
end
function _paretoshapeplot(::Val{:Plots}, args...; kw...)
    return RecipesBase.plot(ParetoShapePlotConfig(args); kw...)
end
function _paretoshapeplot!(::Val{:Plots}, args...; kw...)
    return RecipesBase.plot!(ParetoShapePlotConfig(args); kw...)
end
function paretoshapeplot!(plt::RecipesBase.AbstractPlot, args...; kw...)
    return RecipesBase.plot!(plt, ParetoShapePlotConfig(args); kw...)
end

# pre-process to make PSISResult if necessary
RecipesBase.@recipe function f(config::ParetoShapePlotConfig; showlines=false)
    showlines && RecipesBase.@series begin
        seriestype := :hline
        primary := false
        linestyle --> [:dot :dashdot :dash :solid]
        linealpha --> 0.7
        linecolor --> :grey
        y := [0 0.5 0.7 1]
    end
    xlabel --> "Parameter index"
    ylabel --> "Pareto shape"
    seriestype --> :scatter
    arg = first(config.args)
    k = as_array(missing_to_nan(arg isa PSISResult ? pareto_shape(arg) : arg))
    return (k,)
end

# plot PSISResult using paretoshapeplot if seriestype not specified
RecipesBase.@recipe function f(result::PSISResult)
    if haskey(plotattributes, :seriestype)
        k = as_array(missing_to_nan(pareto_shape(result)))
        return (k,)
    else
        return ParetoShapePlotConfig((result,))
    end
end