text:
module TestDatumStorage

using Mimi
using Test

@defcomp foo begin
    v = Variable(index = [time])
    function run_timestep(p, v, d, ts)
        v.v[ts] = gettime(ts)
    end
end

@defcomp bar begin
    region = Index()
    v = Variable(index = [time, region])
    function run_timestep(p, v, d, ts)
        # v.v[ts, 1:end] = gettime(ts)
        for d in d.region
            v.v[ts, d] = gettime(ts)
        end
    end
end

years = 2001:2010
regions = [:A, :B]
comp_first = 2003
comp_last = 2008

#------------------------------------------------------------------------------
# 1. Single dimension case, fixed timesteps
#------------------------------------------------------------------------------

m = Model()
set_dimension!(m, :time, years)
add_comp!(m, foo, first=comp_first, last=comp_last)

run(m)
v = m[:foo, :v]
@test length(v) == length(years) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years)
    if y < comp_first || y > comp_last
        @test ismissing(v[i])
    else
        @test v[i] == y
    end
end

#------------------------------------------------------------------------------
# 2. Multi-dimension case, fixed timesteps
#------------------------------------------------------------------------------

m2 = Model()

set_dimension!(m2, :time, years)
set_dimension!(m2, :region, regions)
add_comp!(m2, bar, first=comp_first, last=comp_last)

run(m2)
v2 = m2[:bar, :v]
@test size(v2) == (length(years), length(regions)) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years)
    if y < comp_first || y > comp_last
        [@test ismissing(v2[i, j]) for j in 1:length(regions)]
    else
        [@test v2[i, j]==y for j in 1:length(regions)]
    end
end

#------------------------------------------------------------------------------
# 3. Single dimension case, variable timesteps
#------------------------------------------------------------------------------

years_variable = [2000:2004..., 2005:5:2030...]
last = 2010

m = Model()
set_dimension!(m, :time, years_variable)
add_comp!(m, foo, first=comp_first, last=last)

run(m)
v = m[:foo, :v]
@test length(v) == length(years_variable) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years_variable)
    if y < comp_first || y > last
        @test ismissing(v[i])
    else
        @test v[i] == y
    end
end

#------------------------------------------------------------------------------
# 4. Multi-dimension case, variable timesteps
#------------------------------------------------------------------------------

m2 = Model()

set_dimension!(m2, :time, years_variable)
set_dimension!(m2, :region, regions)
add_comp!(m2, bar, first=comp_first, last=last)

run(m2)
v2 = m2[:bar, :v]
@test size(v2) == (length(years_variable), length(regions)) # Test that the array allocated for variable v is the full length of the time dimension

# Test that the missing values were filled in before/after the first/last values
for (i, y) in enumerate(years_variable)
    if y < comp_first || y > last
        [@test ismissing(v2[i, j]) for j in 1:length(regions)]
    else
        [@test v2[i, j]==y for j in 1:length(regions)]
    end
end

end # module


text:
# Testing:
# @assert 1==42 "Shouldn't this be so?"
#> ERROR: assertion failed: Shouldn't this be so?

using Test
# @test 1==3
#>
# Test Failed at REPL[5]:1
#   Expression: 1 == 3
#    Evaluated: 1 == 3
# ERROR: There was an error during testing.

# @test 1 ≈ 1.1
# Test Failed at REPL[58]:1
#   Expression: 1 ≈ 1.1
#    Evaluated: 1 ≈ 1.1
# ERROR: There was an error during testing
@test 1 ≈ 1.1 atol=0.2 #> Test Passed

# Debugging:
using InteractiveUtils # for @which
arr = [1, 2] #> 2-element Array{Int64,1}: 1 2
@which sort(arr)
#> sort(v::AbstractArray{T,1} where T) in Base.Sort at sort.jl:683
456 * 789 + (@show 2 + 3) #>  2 + 3 => 5     359789

# Benchmarking:
@time [x^2 for x in 1:1000] #> elapsed time: 3.911e-6 seconds (8064 bytes allocated) # 1000-element Array{Int64,1}: …
@timed [x^2 for x in 1:1000]
#> ([1,4,9,16,25,36,49,64,81,100  …  982081,984064,986049,988036,990025,992016,9940
# 09,996004,998001,1000000],3.911e-6,8064,0.0)
@elapsed [x^2 for x in 1:1000] #> 3.422e-6
@allocated [x^2 for x in 1:1000] #> 8064

# starting a task:
a = @async 1 + 2
# Task (done) @0x000000002d70faf0

ext: jl


text:
const ELEMENTS = [
    "H", "He",
    "Li", "Be", "B", "C", "N", "O", "F", "Ne",
    "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar",
    "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr",
    "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe",
    "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
    "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi"
    ]

const ANGULAR_MOMENTUM = [
    [(0,0,0)], #S
    [(1,0,0), (0,1,0), (0,0,1)], #P
    [(2,0,0), (1,1,0), (1,0,1), (0,2,0), (0,1,1), (0,0,2)], #D
    ]

const CODATA08_BOHR_TO_A = 0.5291772085936


text:
## Packages
using Parameters

################################### Create the Household instance #########################
Household = @with_kw (na = 8, # number of asset grid
    amax = 6.0, # asset max
    β = 0.96, # discount factor
    μ = 2.0, # inverse elasticity of substitution
    bc = 0, # borrowing constraint (must be weakly negative)
    lGrid = [0.2;0.9], # grid for labor process
    P_ϵ = [0.5 0.5; 0.5 0.5], # transition matrix
    nl = 2, # number of states for labor process
    aGrid = collect(LinRange(sqrt(bc), sqrt(amax), na).^2), # asset grid
    w = 1.0, # wage
    r = 0.02, # interest rate
    ζ = 50000)

################################### Useful Subroutines #########################
function ψ(x::Real,X::Vector,i::Real)
"""
    - This function constructs a linear basis, which is a function of a variable 'x'

    #### Fields

    - 'x': a variable
    - 'X': a vector of the domain of x
    - 'i': order of a basis function, ψ_{i}

    #### Returns

    - ψ_{i}(x): a piecewise function
"""
    if i == 1
        return (x >= X[i] && x <= X[i+1]) ? (X[i+1]-x)/(X[i+1]-X[i]) : 0.0
    elseif i == length(X)
        return (x >= X[i-1] && x <= X[i]) ? (x-X[i-1])/(X[i]-X[i-1]) : 0.0
    else
        if x >= X[i-1] && x <= X[i]
            return (x-X[i-1])/(X[i]-X[i-1]);
        elseif x >= X[i] && x <= X[i+1]
            return (X[i+1]-x)/(X[i+1]-X[i]);
        else
            return 0.0
        end
    end
end

function num_quad(f::Function, a::Real, b::Real, n::Int)
"""
    - This function performs numerical qaudrature using Newton-Cotes Method

    #### Fields

    - 'f': a function we want to integrate
    - 'a': lower bound
    - 'b': upper bound
    - 'n': number of intervals

    #### Returns

    - 'integral' : a numerical integration
"""
    h = (b-a)/(n-1)
    x = zeros(n)
    w = zeros(n)

    for i in 1:n
        x[i] = a + (i-1)*h
    end

    for i in 2:(n-1)
        w[i] = h
    end
    w[1] = h/2
    w[n] = h/2

    integral = sum([w[i]*f(x[i]) for i in 1:n])
    return integral
end

################################### Main algorithm #########################
function get_apol(hh,a,l,θ)
"""
    - This function represents capital policy function as a piecewise linear function: g^{a}(a,l;θ) = ∑θ_{i}ψ_{i}(a)

    #### Fields

    - 'hh': a household instance
    - 'a': asset
    - 'l': labor
    - 'θ': vector of coefficients

    #### Returns

    - apol: an estimated policy function, g^{a}(a,l;θ)

"""
    @unpack aGrid,na,nl,lGrid = hh
    l_index = findall(k->k==l,lGrid) # get index of l
    if first(l_index) == 1 # first(l_index) converts 1-element vector to a float
        apol = sum([θ[i]*ψ(a,aGrid,i) for i = 1:na])
    else
        apol = sum([θ[i+na]*ψ(a,aGrid,i) for i = 1:na])
    end
    return apol
end

function get_expected(hh,a,l,θ)
"""
    - This function calculates the expectation term in Euler equation

    #### Fields

    - 'hh': a household instance
    - 'a': asset
    - 'l': labor
    - 'θ': vector of coefficients

    #### Returns

    - expected: an expected value on the RHS of Euler

"""
    @unpack w,r,lGrid,ζ,μ = hh

    apol = get_apol(hh,a,l,θ)
    appol_L = get_apol(hh,first(apol),lGrid[1],θ) # first(apol) converts 1-element vector to a float
    appol_H = get_apol(hh,first(apol),lGrid[2],θ)

    expected = (w*lGrid[1]+(1+r)*first(apol)-first(appol_L))^(-μ) +
    (w*lGrid[2]+(1+r)*first(apol)-first(appol_H))^(-μ) +
    ζ*(min(first(apol),0)^2)

    return expected
end

function get_resid(hh,a,l,θ)
"""
    - This function represents the residual equation

    #### Fields

    - 'hh': a household instance
    - 'a': asset
    - 'l': labor
    - 'θ': vector of coefficients

    #### Returns

    - resid: a residual equation, R(a,l:θ)

"""
    @unpack w,r,μ,β = hh
    apol = get_apol(hh,a,l,θ)
    resid = (w*l+(1+r)*a-first(apol))^(-μ)-((β*(1+r))/2)*get_expected(hh,a,l,θ);

    return resid
end

function weighted_resid(hh,θ,i)
"""
    - This function represents the numerical integration of weighted residuals with Galerkin weights: ∫ψ_{i}(a)*R(a,l;θ)da

    #### Fields

    - 'hh': a household instance
    - 'θ': vector of coefficients
    - 'i': an integer indicates which weighted residual to compute

    #### Returns

    - w_resid: a weighted residual equation, ∫ψ_{i}(a)*R(a,l:θ)da
"""
    @unpack na,aGrid,lGrid = hh
    if i <= na
        l = lGrid[1];
        w_resid = num_quad(a->ψ(a,aGrid,i)*get_resid(hh,a,l,θ),aGrid[1],aGrid[end],100)
    else
        l = lGrid[2];
        w_resid = num_quad(a->ψ(a,aGrid,i-na)*get_resid(hh,a,l,θ),aGrid[1],aGrid[end],100)
    end
    return w_resid
end

function G(hh,θ)
"""
    - This function stacks up all weighted residuals into a matrix G

    #### Fields

    - 'hh': a household instance
    - 'θ': vector of coefficients

    #### Returns

    - 'G': a stacked system of weighted residuals
"""
    @unpack na,nl = hh
    return [weighted_resid(hh,θ,i) for i = 1:na*nl]
end


text:
abstract type SisoTf{T<:Number} end

+(f1::SisoTf, f2::SisoTf) = +(promote(f1,f2)...)
-(f1::SisoTf, f2::SisoTf) = -(promote(f1,f2)...)

*(f1::SisoTf, f2::SisoTf) = *(promote(f1,f2)...)
/(f1::SisoTf, f2::SisoTf) = /(promote(f1,f2)...)


text:
module GridapMakie

using Gridap
using Gridap.Helpers
using Gridap.Arrays
using Gridap.CellData
using Gridap.Geometry
using Gridap.ReferenceFEs
using Gridap.Visualization
using FillArrays

import Makie
import GeometryBasics

include("conversions.jl")

include("recipes.jl")

end #module


text:
@adjoint function scatter_mul!(ys::CuArray{T}, us::CuArray{T}, xs::CuArray) where {T<:AbstractFloat}
    ys_ = copy(ys)
    scatter_mul!(ys_, us, xs)
    ys_, function (Δ)
        Δy = zero(ys) .+ Δ
        scatter_mul!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = gather(ys, xs) .* gather(zero(Δ)+Δ, xs)
        @inbounds for ind = CartesianIndices(xs)
            ind = Tuple(ind)
            inds = filter(x -> x != ind, rev_xs[xs[ind...]])
            for i = 1:size(us, 1)
                multiplier = one(T)
                for j = inds
                    multiplier *= us[i, j...]
                end
                Δu[i, ind...] *= multiplier
            end
        end
        (Δy, Δu, nothing)
    end
end

@adjoint function scatter_div!(ys::CuArray{T}, us::CuArray{T}, xs::CuArray) where {T<:AbstractFloat}
    ys_ = copy(ys)
    scatter_div!(ys_, us, xs)
    ys_, function (Δ)
        Δy = zero(ys) .+ Δ
        scatter_div!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = - gather(ys, xs)
        Δu .*= gather(zero(Δ)+Δ, xs)
        Δu ./= us.^2
        @inbounds for ind = CartesianIndices(xs)
            ind = Tuple(ind)
            inds = filter(x -> x != ind, rev_xs[xs[ind...]])
            for i = 1:size(us, 1)
                denom = one(T)
                for j = inds
                    denom *= us[i, j...]
                end
                Δu[i, ind...] /= denom
            end
        end
        (Δy, Δu, nothing)
    end
end

@adjoint function scatter_max!(ys::CuArray{T}, us::CuArray{T}, xs::CuArray) where {T<:AbstractFloat}
    max = copy(ys)
    scatter_max!(max, us, xs)
    max, function (Δ)
       Δy = numerical_cmp(ys, max) .* Δ
       Δu = gather(max, xs)
       Δu = numerical_cmp(us, Δu)
       Δu .*= gather(zero(Δ)+Δ, xs)
       (Δy, Δu, nothing)
    end
end

@adjoint function scatter_min!(ys::CuArray{T}, us::CuArray{T}, xs::CuArray) where {T<:AbstractFloat}
    min = copy(ys)
    scatter_min!(min, us, xs)
    min, function (Δ)
       Δy = numerical_cmp(ys, min) .* Δ
       Δu = gather(min, xs)
       Δu = numerical_cmp(us, Δu)
       Δu .*= gather(zero(Δ)+Δ, xs)
       (Δy, Δu, nothing)
    end
end


text:
"""
    module PlainLinearAlgebra

Some linear-algebra-related functions.

### Functions: `diagonal`, `idmat`
"""
module PlainLinearAlgebra

using LinearAlgebra: I, checksquare, diagind
export idmat, diagonal, superdiagonal, subdiagonal, offdiagonal
export altdiagonal

### diagonal

# See Julia issue #30250
"""
    diagonal(M::AbstractMatrix, k::Integer=0)

An iterator over the `k`th diagonal of `M`.

`collect(diagonal(M))` is equivalent to `diag(m)`.

# Examples
```jldoctest
julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> dA = diagonal(A, 1)
2-element view(::Array{Int64,1}, 4:4:8) with eltype Int64:
 2
 6

julia> dA[1] = 1
1

julia> A
3×3 Array{Int64,2}:
 1  1  3
 4  5  6
 7  8  9
```
"""
diagonal(A::AbstractMatrix, k::Integer=0) = view(A, diagind(A, k))
superdiagonal(A::AbstractMatrix) = diagonal(A, 1)
subdiagonal(A::AbstractMatrix) = diagonal(A, -1)
offdiagonal(A, k) = diagonal(A, k)

### idmat

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

Return an `n` x `n` identity matrix of type `Matrix{T}`.
"""
idmat(::Type{T}, n::Integer) where T = Matrix{T}(I, n, n)
idmat(n::Integer) = idmat(Float64, n)
idmat(::Type{T}, m::AbstractMatrix) where T = idmat(T, checksquare(m))
idmat(m::AbstractMatrix) = idmat(eltype(m), m)

include("fill.jl")
include("altdiagonal.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "IVerilog"
version = v"12.0.0"

# Collection of sources required to complete build
# This is a patched master including https://github.com/steveicarus/iverilog/pull/511/
# to allow for cross compilation
# TODO use upstream + stable once available
sources = [
   GitSource("https://github.com/sjkelly/iverilog.git", "1f09b041f1060116840b901c21938ad31d79bb98")
]

dependencies = [
    HostBuildDependency("Bison_jll"),
    Dependency("Readline_jll"; compat="8.1.1"),
    HostBuildDependency("gperf_jll"),
    Dependency("Zlib_jll"; compat="~1.2.11"),
]

# Bash recipe for building across all platforms
script = raw"""
cd iverilog
export CPPFLAGS="-I${includedir}"
sh ./autoconf.sh
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-zlib=${prefix}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("iverilog", :iverilog),
    ExecutableProduct("iverilog-vpi", :iverilog_vpi)
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
myamp=RHD2164("PortA1")

myfpga=FPGA(1,myamp)

mt=Task_NoTask()

mys=SaveNone()

(myrhd,sss,fpgas)=makeRHD([myfpga],sav=mys)

handles=makegui(myrhd,sss,mt,fpgas)


text:
# FWI on the 2D Overthrust model using spectral projected gradient descent
# Author: mlouboutin3@gatech.edu
# Date: February 2021
#

using DrWatson
@quickactivate :TimeProbeSeismic

# Load Sigsbee model
~isfile(datadir("models", "seam_model.bin")) && run(`curl -L https://www.dropbox.com/s/nyazel8g6ah4cld/seam_model.bin\?dl\=0 --create-dirs -o $(datadir("models", "seam_model.bin"))`)
~isfile(datadir("data", "seam_obn_data.bin"))&& run(`curl -L https://www.dropbox.com/s/ap2y4ny5n98kd5c/seam_obn_data.bin\?dl\=0 --create-dirs -o $(datadir("data", "seam_obn_data.bin"))`)
# Read model and set up background
vp, rho, n, d, o = deserialize(datadir("models", "bin.jld"))

m = vp.^(-2)
m0 = Float32.(imfilter(m, Kernel.gaussian(15)))
rho0 = Float32.(imfilter(rho[:, :], Kernel.gaussian(15)))
dm = m0 - m

n = size(m0)
d = (12.5, 10.)

# Setup info and model structure
model0 = Model(n, d, o, m0, rho=rho0; nb=40)

# Load data
dobs, q = deserialize(datadir("data", "seam_obn_data.bin"))
dobs = dobs[[1, 22]]
q = q[[1, 22]]

# Info structure for linear operators
ntComp = get_computational_nt(q.geometry, dobs.geometry, model0)    # no. of computational time steps
info = Info(prod(model0.n), dobs.nsrc, ntComp)

#################################################################################################
# Write shots as segy files to disk
opt = Options(space_order=16, isic=true)

# Setup operators
F0 = judiModeling(info, model0, q.geometry, dobs.geometry; options=opt)
J = judiJacobian(F0, q)

# data
# Nonlinear modeling
d0 = F0*q
residual = d0 - dobs

# Probe
ge = Array{Any}(undef, 6, 2)
for ps=1:6
    ge[ps, 1] = judiJacobian(F0[1], q[1], 2^ps, dobs[1])'*residual[1]
    ge[ps, 2] = judiJacobian(F0[2], q[2], 2^ps, dobs[2])'*residual[2]
end

for i=1:2
    clip = maximum(ge[ps, i])/10
    figure()
    for ps=1:6
        subplot(2,3,ps)
        imshow(ge[ps, i]', cmap="seismic", vmin=-clip, vmax=clip, aspect="auto")
        title("ps=$(2^ps)")
    end
    tight_layout()
end

# Probing vectors
figure();imshow(dobs.data[1]*dobs.data[1]', cmap="seismic", vmin=-10, vmax=10)
title(L"$d_{obs} d_{obs}^\top$")

figure()
for ps=1:9
    subplot(3,3,ps)
    Q = qr_data(dobs.data[1]*dobs.data[1]', 2^ps)
    imshow(Q, vmin=-.1, vmax=.1, cmap="seismic", aspect="auto")
    title("Probing vectors ps=$(2^ps)")
end
tight_layout()

max_stars_repo_name: slimgroup/TimeProbeSeismic

max_issues_repo_name: slimgroup/TimeProbeSeismic

max_forks_repo_name: slimgroup/TimeProbeSeismic


text:
@testset "enumtype" begin

@testset "instances" begin
    @test instances(GridType) == (PRIM, DUAL)
    @test instances(Sign) == (NEG, POS)
end

@testset "numel" begin
    @test numel(GridType) == 2
    @test numel(Sign) == 2
end

@testset "integers" begin
    @test Int.(SVec(instances(GridType))) == nPD == [nPR, nDL]
    @test Int.(SVec(instances(Sign))) == nNP == [nN, nP]
end

@testset "next and alter" begin
    @test alter(PRIM) == DUAL
    @test alter(DUAL) == PRIM

    @test alter(NEG) == POS
    @test alter(POS) == NEG
end

end  # @testset "enumtype"


text:
using VLConstraintBasedModelGenerationUtilities

# setup path to sequence file -
path_to_protein_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/P-MZ373340.fasta"

# load -
result = build_protein_table(path_to_protein_sequence_file) |> check


text:
using JSON

demandsDict = Dict{}()
# get ground trueth demands, indexed by 0
demandsDict[0] = iniDemand("../../data_original/Anaheim_trips.txt")
# get initial demands, indexed by 1
demandsDict[1] = iniDemand("../../data_original/Anaheim_trips.txt", 1)

numNodes, numLinks, numODpairs, capacity, free_flow_time, ta_data = paraNetwork("Anaheim")
numRoutes, odPairRoute, linkRoute, link_label_dict, link_label_dict_, link_length_dict, OD_pair_route_dict = furInfo()

start_node = ta_data.start_node
end_node = ta_data.end_node

#load OD pair labels
odPairLabel = readstring("../../temp_files/od_pair_label_dict_Anaheim_refined.json")
odPairLabel = JSON.parse(odPairLabel)

odPairLabel_ = readstring("../../temp_files/od_pair_label_dict__Anaheim_refined.json")
odPairLabel_ = JSON.parse(odPairLabel_)

#load node-link incidence
nodeLink = readstring("../../temp_files/node_link_incidence_Anaheim.json");
nodeLink = JSON.parse(nodeLink);

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
using Distributions
function discrete_normal(n,mu,sigma,width)

# creates equally spaced approximation to normal distribution
# n is number of points
# mu is mean
# sigma is standard deviation
# width is the multiple of stand deviation for the width of the grid
# f is the error in the approximation
# x gives the location of the points
# p is probabilities

    x = collect(range(mu-width*sigma,mu+width*sigma,length=n))
    if n==2
        p = 0.5.*ones(n)
    elseif n>2
        p  = zeros(n)
        p[1] = cdf(Normal(mu,sigma), x[1] + 0.5*(x[2]-x[1]))
        for i = 2:n-1
            p[i] = cdf(Normal(mu,sigma), x[i] + 0.5*(x[i+1]-x[i])) - cdf(Normal(mu,sigma), x[i] - 0.5*(x[i]-x[i-1]))
        end
        p[n] = 1 - sum(p[1:n-1])
    end

    Ex = x'*p
    SDx = sqrt((x.^2)'*p .- Ex.^2)

    f = SDx.-sigma

    return f,x,p

end


text:
using Statistics

using SIMD
using BenchmarkTools
using ChipSort
using JLD

const SUITE = BenchmarkGroup()

function randa(T, K...)
    data = reshape(valloc(T, div(64, sizeof(T)), prod(K)), K...)
    data .= rand(T, K...)
    data
end

T=UInt32
V = 2^4
J = 2^4
kk = 2 .^ (0:16)

for K in kk
    len = V*J*K
    SUITE[len] = BenchmarkGroup(["size-$len"])
    SUITE[len][:JuliaStd] = @benchmarkable sort!(data) setup=(data = randa($T, $len)) seconds=4
    if K >= 2^2
        SUITE[len][:ChipSortM] = @benchmarkable chipsort_medium!(data, Val($V), Val($J), Val($K)) setup=(data = randa($T, $len)) seconds=4
    end

    if len <= 2^18
        SUITE[len][:CombSort] = @benchmarkable combsort!(data) setup=(data = randa($T, $len)) seconds=4
    end
    SUITE[len][:CombSortMix] = @benchmarkable chipsort_serial!(data) setup=(data = randa($T, $len)) seconds=4

    if len > 2^10
        SUITE[len][:ChipSortL] = @benchmarkable chipsort_large(data, Val($V), Val(2^4)) setup=(data = randa($T, $len)) seconds=4
    end
    if K <= 2^3
        SUITE[len][:InsertionSort] = @benchmarkable insertion_sort!(data) setup=(data = randa($T, $len)) seconds=4
    end
end

bmk = run(SUITE)
save("chip-medium-bench.jld", "bmk", bmk, "T", string(T), "V", V, "J", J)


text:
"""
    DVFA.jl

Description:
    Includes all of the structures and logic for running a Dual-Vigilance Fuzzy ART (DVFA) module.

Authors:
    MATLAB implementation: Leonardo Enzo Brito da Silva
    Julia port: Sasha Petrenko <sap625@mst.edu>

References:
[1] L. E. Brito da Silva, I. Elnabarawy and D. C. Wunsch II, "Dual
Vigilance Fuzzy ART," Neural Networks Letters. To appear.
[2] G. Carpenter, S. Grossberg, and D. Rosen, "Fuzzy ART: Fast
stable learning and categorization of analog patterns by an adaptive
resonance system," Neural Networks, vol. 4, no. 6, pp. 759-771, 1991.
"""

"""
    opts_DVFA(;kwargs)

Dual Vigilance Fuzzy ART options struct.

# Keyword Arguments
- `rho_lb::Float`: lower-bound vigilance value, [0, 1], default 0.55.
- `rho_ub::Float`: upper-bound vigilance value, [0, 1], default 0.75.
- `alpha::Float`: choice parameter, alpha > 0, default 1e-3.
- `beta::Float`: learning parameter, (0, 1], default 1.0.
- `display::Bool`: display flag, default true.
- `max_epoch::Int`: maximum number of epochs during training, default 1.
"""
@with_kw mutable struct opts_DVFA <: ARTOpts @deftype Float
    # Lower-bound vigilance parameter: [0, 1]
    rho_lb = 0.55; @assert rho_lb >= 0.0 && rho_lb <= 1.0
    # Upper bound vigilance parameter: [0, 1]
    rho_ub = 0.75; @assert rho_ub >= 0.0 && rho_ub <= 1.0
    # Choice parameter: alpha > 0
    alpha = 1e-3; @assert alpha > 0.0
    # Learning parameter: (0, 1]
    beta = 1.0; @assert beta > 0.0 && beta <= 1.0
    # Display flag
    display::Bool = true
    # Maximum number of epochs during training
    max_epochs::Int = 1
end # opts_DVFA

"""
    DVFA <: ART

Dual Vigilance Fuzzy ARTMAP module struct.

For module options, see [`AdaptiveResonance.opts_DVFA`](@ref).

# Option Parameters
- `opts::opts_DVFA`: DVFA options struct.
- `config::DataConfig`: data configuration struct.

# Working Parameters
- `threshold_ub::Float`: operating upper bound module threshold value, a function of the upper bound vigilance parameter.
- `threshold_lb::Float`: operating lower bound module threshold value, a function of the lower bound vigilance parameter.
- `labels::IntegerVector`: incremental list of labels corresponding to each F2 node, self-prescribed or supervised.
- `W::RealMatrix`: category weight matrix.
- `T::RealVector`: activation values for every weight for a given sample.
- `M::RealVector`: match values for every weight for a given sample.
- `n_categories::Int`: number of category weights (F2 nodes).
- `n_clusters::Int`: number of labeled clusters, may be lower than `n_categories`
- `epoch::Int`: current training epoch.

# References:
1. L. E. Brito da Silva, I. Elnabarawy and D. C. Wunsch II, "Dual Vigilance Fuzzy ART," Neural Networks Letters. To appear.
2. G. Carpenter, S. Grossberg, and D. Rosen, "Fuzzy ART: Fast stable learning and categorization of analog patterns by an adaptive resonance system," Neural Networks, vol. 4, no. 6, pp. 759-771, 1991.
"""
mutable struct DVFA <: ART
    # Get parameters
    opts::opts_DVFA
    config::DataConfig

    # Working variables
    threshold_ub::Float
    threshold_lb::Float
    labels::IntegerVector
    W::RealMatrix
    T::RealVector
    M::RealVector
    n_categories::Int
    n_clusters::Int
    epoch::Int
end # DVFA

"""
    DVFA()

Implements a DVFA learner with default options.

# Examples
```julia-repl
julia> DVFA()
DVFA
    opts: opts_DDVFA
    ...
```
"""
function DVFA()
    opts = opts_DVFA()
    DVFA(opts)
end # DVFA()

"""
    DVFA(;kwargs...)

Implements a DVFA learner with keyword arguments.

# Examples
```julia-repl
julia> DVFA(rho=0.7)
DVFA
    opts: opts_DDVFA
    ...
```
"""
function DVFA(;kwargs...)
    opts = opts_DVFA(;kwargs...)
    DVFA(opts)
end # DVFA(;kwargs...)

"""
    DVFA(opts::opts_DVFA)

Implements a DVFA learner with specified options.

# Examples
```julia-repl
julia> my_opts = opts_DVFA()
julia> DVFA(my_opts)
DVFA
    opts: opts_DVFA
    ...
```
"""
function DVFA(opts::opts_DVFA)
    DVFA(
        opts,                           # opts
        DataConfig(),                   # config
        0.0,                            # threshold_ub
        0.0,                            # threshold_lb
        Array{Int}(undef, 0),           # labels
        Array{Float}(undef, 0, 0),      # W
        Array{Float}(undef, 0),         # M
        Array{Float}(undef, 0),         # T
        0,                              # n_categories
        0,                              # n_clusters
        0                               # epoch
    )
end # DDVFA(opts::opts_DDVFA)

"""
    set_threshold!(art::DVFA)

Configure the threshold values of the DVFA module.
"""
function set_threshold!(art::DVFA)
    # DVFA thresholds
    art.threshold_ub = art.opts.rho_ub * art.config.dim
    art.threshold_lb = art.opts.rho_lb * art.config.dim
end # set_threshold!(art::DVFA)

# Incremental DVFA training method
function train!(art::DVFA, x::RealVector ; y::Integer=0, preprocessed::Bool=false)
    # Flag for if training in supervised mode
    supervised = !iszero(y)

    # Run the sequential initialization procedure
    sample = init_train!(x, art, preprocessed)

    # Initialization
    if isempty(art.W)
        # Set the threshold
        set_threshold!(art)
        # Set the first label as either 1 or the first provided label
        y_hat = supervised ? y : 1
        # Create a new category and cluster
        art.W = ones(art.config.dim_comp, 1)
        art.n_categories = 1
        art.n_clusters = 1
        push!(art.labels, y_hat)
        return y_hat
    end

    # If label is new, break to make new category
    if supervised && !(y in art.labels)
        y_hat = y
        # Update sample labels
        push!(art.labels, y)
        # Fast commit the sample
        art.W = hcat(art.W, sample)
        art.n_categories += 1
        art.n_clusters += 1
        return y_hat
    end

    # Compute the activation and match for all categories
    activation_match!(art, sample)
    # Sort activation function values in descending order
    index = sortperm(art.T, rev=true)

    # Default to mismatch
    mismatch_flag = true
    # Loop over all categories
    for j = 1:art.n_categories
        # Best matching unit
        bmu = index[j]
        # Vigilance test upper bound
        if art.M[bmu] >= art.threshold_ub
            # If supervised and the label differs, trigger mismatch
            if supervised && (art.labels[bmu] != y)
                break
            end
            # Learn the sample
            learn!(art, sample, bmu)
            # Update sample label for output
            # y_hat = supervised ? y : art.labels[bmu]
            y_hat = art.labels[bmu]
            # No mismatch
            mismatch_flag = false
            break
        # Vigilance test lower bound
        elseif art.M[bmu] >= art.threshold_lb
            # If supervised and the label differs, trigger mismatch
            if supervised && (art.labels[bmu] != y)
                break
            end
            # Update sample labels
            y_hat = supervised ? y : art.labels[bmu]
            push!(art.labels, y_hat)
            # Fast commit the sample
            art.W = hcat(art.W, sample)
            art.n_categories += 1
            # No mismatch
            mismatch_flag = false
            break
        end
    end

    # If there was no resonant category, make a new one
    if mismatch_flag
        # Create a new category-to-cluster label
        y_hat = supervised ? y : art.n_clusters + 1
        push!(art.labels, y_hat)
        # Fast commit the sample
        art.W = hcat(art.W, sample)
        # Increment the number of categories and clusters
        art.n_categories += 1
        art.n_clusters += 1
    end

    return y_hat
end # train!(art::DVFA, x::RealVector ; y::Integer=0, preprocessed::Bool=false)

# Incremental DVFA classify method
function classify(art::DVFA, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)
    # Preprocess the data
    sample = init_classify!(x, art, preprocessed)

    # Compute activation and match functions
    activation_match!(art, sample)
    # Sort activation function values in descending order
    index = sortperm(art.T, rev=true)
    mismatch_flag = true
    for jx in 1:art.n_categories
        bmu = index[jx]
        # Vigilance check - pass
        if art.M[bmu] >= art.threshold_ub
            # Current winner
            y_hat = art.labels[bmu]
            mismatch_flag = false
            break
        end
    end

    # If we did not find a resonant category
    if mismatch_flag
        # Create new weight vector
        @debug "Mismatch"
        # Report either the best matching unit or the mismatch label -1
        y_hat = get_bmu ? art.labels[index[1]] : -1
    end

    return y_hat
end # classify(art::DVFA, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)

"""
    activation_match!(art::DVFA, x::RealVector)

Compute and store the activation and match values for the DVFA module.
"""
function activation_match!(art::DVFA, x::RealVector)
    art.T = zeros(art.n_categories)
    art.M = zeros(art.n_categories)
    for jx = 1:art.n_categories
        numerator = norm(element_min(x, art.W[:, jx]), 1)
        art.T[jx] = numerator/(art.opts.alpha + norm(art.W[:, jx], 1))
        art.M[jx] = numerator
    end
end # activation_match!(art::DVFA, x::RealVector)

"""
    learn(art::DVFA, x::RealVector, W::RealVector)

Return the modified weight of the DVFA module conditioned by sample x.
"""
function learn(art::DVFA, x::RealVector, W::RealVector)
    # Update W
    return art.opts.beta .* element_min(x, W) .+ W .* (1 - art.opts.beta)
end # learn(art::DVFA, x::RealVector, W::RealVector)

"""
    learn!(art::DVFA, x::RealVector, index::Integer)

In place learning function.
"""
function learn!(art::DVFA, x::RealVector, index::Integer)
    # Update W
    art.W[:, index] = learn(art, x, art.W[:, index])
end # learn!(art::DVFA, x::RealVector, index::Integer)

"""
    stopping_conditions(art::DVFA)

Stopping conditions for a DVFA module.
"""
function stopping_conditions(art::DVFA)
    return art.epoch >= art.opts.max_epochs
end # stopping_conditions(art::DVFA)


text:
using Calphad
using Documenter

DocMeta.setdocmeta!(Calphad, :DocTestSetup, :(using Calphad); recursive=true)

makedocs(;
    modules=[Calphad],
    authors="Brandon Bocklund",
    repo="https://github.com/bocklund/Calphad.jl/blob/{commit}{path}#{line}",
    sitename="Calphad.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bocklund.github.io/Calphad.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bocklund/Calphad.jl",
    devbranch="main",
)


text:
# ---
# title: 1156. Swap For Longest Repeated Character Substring
# id: problem1156
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/swap-for-longest-repeated-character-substring/description/>
# hidden: true
# ---
#
# Given a string `text`, we are allowed to swap two of the characters in the
# string. Find the length of the longest substring with repeated characters.
#
#
#
# **Example 1:**
#
#
#
#     Input: text = "ababa"
#     Output: 3
#     Explanation: We can swap the first 'b' with the last 'a', or the last 'b' with the first 'a'. Then, the longest repeated character substring is "aaa", which its length is 3.
#
#
# **Example 2:**
#
#
#
#     Input: text = "aaabaaa"
#     Output: 6
#     Explanation: Swap 'b' with the last 'a' (or the first 'a'), and we get longest repeated character substring "aaaaaa", which its length is 6.
#
#
# **Example 3:**
#
#
#
#     Input: text = "aaabbaaa"
#     Output: 4
#
#
# **Example 4:**
#
#
#
#     Input: text = "aaaaa"
#     Output: 5
#     Explanation: No need to swap, longest repeated character substring is "aaaaa", length is 5.
#
#
# **Example 5:**
#
#
#
#     Input: text = "abcdef"
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= text.length <= 20000`
#   * `text` consist of lowercase English characters only.
#
#
## @lc code=start
using LeetCode

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


text:
vcov(m, ::Nothing) = nothing

"""
    criticalvalue(vce::CovarianceEstimator, level::Real, T::Int)

Return the critical value for `vce` at the given `level` when the sample has `T` periods.
"""
criticalvalue(::CovarianceEstimator, level::Real, T::Int) = norminvccdf((1.0-level)/2.0)

"""
    pvalue(vce::CovarianceEstimator, ts::Real, T::Int)

Return the p-value for test statistic `ts` associated with `vce`
when the sample has `T` periods.
"""
pvalue(::CovarianceEstimator, ts::Real, T::Int) = 2.0*normccdf(abs(ts))

"""
    SimpleVCE <: CovarianceEstimator

Simple variance-covariance estimator.
"""
struct SimpleVCE <: CovarianceEstimator end

function vcov(m::OLS, ::SimpleVCE)
    T, K = size(m.X)
    return kron(m.invXX, m.resid'm.resid)./(T-K)
end

show(io::IO, ::MIME"text/plain", ::SimpleVCE) =
    print(io, "Simple variance-covariance estimator")

"""
    HRVCE <: CovarianceEstimator

Eicker-Huber-White heteroskedasticity-robust variance-covariance estimator.
"""
struct HRVCE <: CovarianceEstimator end

function vcov(m::OLS, ::HRVCE)
    T, K = size(m.X)
    return T/(T-K).*kron_fastr(m.invXX, kron_fastr(m.invXX, m.score'm.score)')'
end

function vcov(m::Ridge, ::HRVCE)
    T = size(m.C, 1)
    S = m.score'm.score
    return T/(m.dof_res-m.dof_adj) * m.invCCP * S * m.invCCP
end

show(io::IO, ::MIME"text/plain", ::HRVCE) =
    print(io, "Heteroskedasticity-robust variance-covariance estimator")

"""
    LongRunVariance

Supertype for all long-run variance estimators.
"""
abstract type LongRunVariance end

"""
    HARVCE{LR<:LongRunVariance} <: CovarianceEstimator

Heteroskedasticity-autocorrelation-robust variance-covariance estimator
with a long-run variance estimator of type `LR`.

# Fields
- `lr::LR`: the long-run variance estimator.
- `bw::Function`: bandwidth as a function of the number of periods in sample.
- `cv::Function`: critical value as a function of the level of the test and the bandwidth.
- `pv::Function`: p-value as a function of the test statistic and the bandwidth.
"""
struct HARVCE{LR<:LongRunVariance} <: CovarianceEstimator
    lr::LR
    bw::Function
    cv::Function
    pv::Function
end

"""
    EqualWeightedCosine <: LongRunVariance

Equal-weighted cosine transform. The alias [`EWC`](@ref) is provided for convenience.

# Reference
Lazarus, Eben, Daniel J. Lewis, James H. Stock, and Mark W. Watson. 2018.
"HAR Inference: Recommendations for Practice."
Journal of Business & Economic Statistics 36 (4): 541-559.
"""
struct EqualWeightedCosine <: LongRunVariance end

"""
    EWC <: LongRunVariance

An alias for [`EqualWeightedCosine`](@ref).
"""
const EWC = EqualWeightedCosine

HARVCE(lr::EqualWeightedCosine; bw=_ewcbw, cv=_ewccv, pv=_ewcpv) = HARVCE(lr, bw, cv, pv)

_ewcbw(T::Integer) = round(Int, 0.4*T^(2/3))
_ewccv(level::Real, bw::Int) = tdistinvccdf(bw, (1.0-level)/2.0)
_ewcpv(ts::Real, bw::Int) = 2.0*tdistccdf(bw, abs(ts))

function lrv(vce::HARVCE{EqualWeightedCosine}, sc::AbstractMatrix)
    T, K = size(sc)
    TF = eltype(sc)
    bw = vce.bw(T)
    Z = zeros(Complex{TF}, 4*T, K)
    inds = iseven.(1:2*T)
    copyto!(view(view(Z, 1:2*T, :), inds, :), sc)
    fft!(Z, 1)
    rffts = real.(view(Z, 2:bw+1, :))
    return 2.0/bw * (rffts'rffts)
end

show(io::IO, ::MIME"text/plain", ::EqualWeightedCosine) =
    print(io, "Equal-weighted cosine transform")

function vcov(m::OLS, vce::HARVCE)
    T, K = size(m.X)
    return T/(T-K).*kron_fastr(m.invXX, kron_fastr(m.invXX, lrv(vce, m.score))')'
end

function vcov(m::Ridge, vce::HARVCE)
    T = size(m.C, 1)
    return T/(m.dof_res-m.dof_adj) * m.invCCP * lrv(vce, m.score) * m.invCCP
end

criticalvalue(vce::HARVCE, level::Real, T::Int) = vce.cv(level, vce.bw(T))
pvalue(vce::HARVCE, ts::Real, T::Int) = vce.pv(ts, vce.bw(T))

function show(io::IO, ::MIME"text/plain", vce::HARVCE)
    println(io, "Heteroskedasticity-autocorrelation-robust variance-covariance estimator:")
    print(io, "  Long-run variance estimator: ", sprint(show, MIME("text/plain"), vce.lr))
end


text:
# ---
# title: 1038. Binary Search Tree to Greater Sum Tree
# id: problem1038
# author: Indigo
# date: 2021-07-31
# difficulty: Medium
# categories: Binary Search Tree
# link: <https://leetcode.com/problems/binary-search-tree-to-greater-sum-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a Binary Search Tree (BST), convert it to a Greater Tree
# such that every key of the original BST is changed to the original key plus
# sum of all keys greater than the original key in BST.
#
# As a reminder, a _binary search tree_ is a tree that satisfies these
# constraints:
#
#   * The left subtree of a node contains only nodes with keys  **less than**  the node's key.
#   * The right subtree of a node contains only nodes with keys  **greater than**  the node's key.
#   * Both the left and right subtrees must also be binary search trees.
#
# **Note:** This question is the same as 538:
# <https://leetcode.comhttps://leetcode.com/problems/convert-bst-to-greater-tree/>
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/02/tree.png)
#
#
#
#     Input: root = [4,1,6,0,2,5,7,null,null,null,3,null,null,null,8]
#     Output: [30,36,21,36,35,26,15,null,null,null,33,null,null,null,8]
#
#
# **Example 2:**
#
#
#
#     Input: root = [0,null,1]
#     Output: [1,null,1]
#
#
# **Example 3:**
#
#
#
#     Input: root = [1,0,2]
#     Output: [3,3,2]
#
#
# **Example 4:**
#
#
#
#     Input: root = [3,2,4,1]
#     Output: [7,9,4,10]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[1, 100]`.
#   * `0 <= Node.val <= 100`
#   * All the values in the tree are **unique**.
#   * `root` is guaranteed to be a valid binary search tree.
#
#
## @lc code=start
using LeetCode

bst_to_gst(::Nothing) = nothing
function bst_to_gst(root::TreeNode{Int})::TreeNode{Int}
    s = 0
    rev_first_ord(::Nothing) = nothing
    function rev_first_ord(node::TreeNode{Int})
        rev_first_ord(node.right)
        s += node.val
        node.val = s
        return rev_first_ord(node.left)
    end
    !isnothing(root) && rev_first_ord(root)
    return root
end
## @lc code=end


text:
# ---
# title: 1115. Print FooBar Alternately
# id: problem1115
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/print-foobar-alternately/description/>
# hidden: true
# ---
#
# Suppose you are given the following code:
#
#
#
#     class FooBar {
#       public void foo() {
#         for (int i = 0; i < n; i++) {
#           print("foo");
#         }
#       }
#
#       public void bar() {
#         for (int i = 0; i < n; i++) {
#           print("bar");
#         }
#       }
#     }
#
#
# The same instance of `FooBar` will be passed to two different threads. Thread
# A will call `foo()` while thread B will call `bar()`. Modify the given program
# to output "foobar" _n_ times.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: "foobar"
#     Explanation: There are two threads being fired asynchronously. One of them calls foo(), while the other calls bar(). "foobar" is being output 1 time.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: "foobarfoobar"
#     Explanation: "foobar" is being output 2 times.
#
#
#
## @lc code=start
using LeetCode

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


text:
using Mustache
using Test

@testset " delimiters " begin

	## The equals sign (used on both sides) should permit delimiter changes.
tpl = """{{=<% %>=}}(<%text%>)"""

	@test Mustache.render(tpl, Dict{Any,Any}("text"=>"Hey!")) == """(Hey!)"""

	## Characters with special meaning regexen should be valid delimiters.
tpl = """({{=[ ]=}}[text])"""

	@test Mustache.render(tpl, Dict{Any,Any}("text"=>"It worked!")) == """(It worked!)"""

	## Delimiters set outside sections should persist.
tpl = """[
{{#section}}
  {{data}}
  |data|
{{/section}}

{{= | | =}}
|#section|
  {{data}}
  |data|
|/section|
]
"""

	@test Mustache.render(tpl, Dict{Any,Any}("section"=>true,"data"=>"I got interpolated.")) == """[
  I got interpolated.
  |data|

  {{data}}
  I got interpolated.
]
"""

	## Delimiters set outside inverted sections should persist.
tpl = """[
{{^section}}
  {{data}}
  |data|
{{/section}}

{{= | | =}}
|^section|
  {{data}}
  |data|
|/section|
]
"""

	@test Mustache.render(tpl, Dict{Any,Any}("section"=>false,"data"=>"I got interpolated.")) == """[
  I got interpolated.
  |data|

  {{data}}
  I got interpolated.
]
"""

	## Delimiters set in a parent template should not affect a partial.
tpl = """[ {{>include}} ]
{{= | | =}}
[ |>include| ]
"""

	@test Mustache.render(tpl, Dict{Any,Any}("value"=>"yes","include"=>".{{value}}.")) == """[ .yes. ]
[ .yes. ]
"""

	## Delimiters set in a partial should not affect the parent template.
tpl = """[ {{>include}} ]
[ .{{value}}.  .|value|. ]
"""

	@test Mustache.render(tpl, Dict{Any,Any}("value"=>"yes","include"=>".{{value}}. {{= | | =}} .|value|.")) == """[ .yes.  .yes. ]
[ .yes.  .|value|. ]
"""

	## Surrounding whitespace should be left untouched.
tpl = """| {{=@ @=}} |"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """|  |"""

	## Whitespace should be left untouched.
tpl = """ | {{=@ @=}}
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """ |
"""

	## Standalone lines should be removed from the template.
tpl = """Begin.
{{=@ @=}}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## Indented standalone lines should be removed from the template.
tpl = """Begin.
  {{=@ @=}}
End.
"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """Begin.
End.
"""

	## "\r\n" should be considered a newline for standalone tags.
tpl = """|
{{= @ @ =}}
|"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """|
|"""

	## Standalone tags should not require a newline to precede them.
tpl = """  {{=@ @=}}
="""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """="""

	## Standalone tags should not require a newline to follow them.
tpl = """=
  {{=@ @=}}"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """=
"""

	## Superfluous in-tag whitespace should be ignored.
tpl = """|{{= @   @ =}}|"""

	@test Mustache.render(tpl, Dict{Any,Any}()) == """||"""
end


text:
@testset "ch1" begin
    function create_toy_data(func, sample_size, std)
        x = collect(range(0, stop = 1.0, length = sample_size))
        noise = rand(Normal(0.0, std), sample_size)
        return x, func(x) + noise
    end

    function sinusoidal(x)
        return sin.(2 * pi * x)
    end

    function rmse(a::Array{Float64,1}, b::Array{Float64,1})
        return sum((b - a) .^ 2) / size(a)[1]
    end

    x_train, y_train = create_toy_data(sinusoidal, 10, 0.25)
    x_test = collect(range(0, stop = 1.0, length = 100))
    y_test = sinusoidal(x_test)

    @testset "ch1_1" begin
        feature = PolynomialFeature(9)
        X_train = transform(feature, x_train)
        X_test = transform(feature, x_test)
        #x = collect(reshape(range(0, stop=1.0, length=10), 10));
        model = LinearRegressor([0], 0)
        #Phi = collect(reshape(transpose(X_train), size(X_train)[2], size(X_train)[1]));
        fitting(model, X_train, y_train)
        #tmp = collect(reshape(transpose(X_test), size(X_test)[2], size(X_test)[1]));
        y, y_std = predict(model, X_test, true)
    end

    @testset "ch1_2" begin

        training_errors = []
        test_errors = []

        for i = 0:10
            feature = PolynomialFeature(i)
            X_train = transform(feature, x_train)
            X_test = transform(feature, x_test)

            model = LinearRegressor([0], 0)
            fitting(model, X_train, y_train)
            y_trained = predict(model, X_train, false)
            push!(training_errors, rmse(predict(model, X_train, false), y_train))
            push!(
                test_errors,
                rmse(
                    predict(model, X_test, false),
                    y_test + rand(Normal(0.0, 0.25), size(y_test)[1]),
                ),
            )
        end
    end
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
################################################################################
### Basic extension
################################################################################

# Sites
convert_arguments(::PointBased, s::S) where {T, D, S<:AbstractSite{T, D}} =
    ([s |> point |> Point{D, Float32}],)
convert_arguments(::PointBased, s::Vector{S}) where {T, D, S<:AbstractSite{T, D}} =
    (s .|> point .|> Point{D, Float32},)

# Bonds

# These shouldn't really exist (single bond from two sites)
for T in (LineSegments, Lines)
    @eval function convert_arguments(P::Type{<:$T}, s1::S, s2::S)  where {
            T, D, S<:AbstractSite{T, D}
        }
        p1 = point(s1)
        p2 = point(s2)
        convert_arguments(P, Point{D, Float32}[p1, p2])
    end
end

# Maybe also wireframe?
function convert_arguments(
        P::Type{<:LineSegments}, bonds::Vector{B}, sites::Vector{S}
    ) where {
        T, D, S <: AbstractSite{T, D}, B <: AbstractBond
    }
    ps = [point(sites[i]) for b in bonds for i in [from(b), to(b)]]
    convert_arguments(P, Point{D, Float32}.(ps))
end

function convert_arguments(
        P::Type{<:MeshScatter}, bonds::Vector{B}, sites::Vector{S}
    ) where {
        T, D, S <: AbstractSite{T, D}, B <: AbstractBond
    }
    convert_arguments(P, sites[bonds .|> from] .|> point .|> Point{D, Float32})
end

# Lattice
convert_arguments(P::Type{<:Scatter}, l::AbstractLattice) =
    convert_arguments(P, sites(l))
convert_arguments(P::Type{<:MeshScatter}, l::AbstractLattice) =
    convert_arguments(P, sites(l))
convert_arguments(P::Type{<:LineSegments}, l::AbstractLattice) =
    convert_arguments(P, bonds(l), sites(l))

# ReciprocalPath
function convert_arguments(::PointBased, p::P) where {
        D, P <: AbstractReciprocalPath{<: AbstractReciprocalPoint{D}}
    }
    (Point{D, Float32}.(point.(points(p))),)
end

# Brillouin Zone
function convert_arguments(P::Type{<:Scatter}, bz::BZ) where {
        D, PT <: AbstractReciprocalPoint{D},
        R <: AbstractReciprocalUnitcell{PT},
        BZ <: AbstractBrillouinZone{R}
    }
    convert_arguments(P, Point{D, Float32}[c for c in corners(bz)])
end

# TODO mabye duplicate points for sharp edges?
function BZ_to_mesh(bz::BZ) where {
        D, PT <: AbstractReciprocalPoint{D},
        R <: AbstractReciprocalUnitcell{PT},
        BZ <: AbstractBrillouinZone{R}
    }
    # I assume that every face is flat and rotates around a central point
    points = Point{D, Float32}.(corners(bz))
    _faces = GLTriangleFace[]
    for f in LatPhysReciprocal.faces(bz)
        center = mean(points[f])
        push!(points, center)
        c = length(points)
        append!(_faces, GLTriangleFace.(c, f[1:end-1], f[2:end]))
        push!(_faces, GLTriangleFace(c, f[end], f[1]))
    end
    normal_mesh(points, _faces)
end
function convert_arguments(P::Type{<: MakieCore.Mesh}, bz::AbstractBrillouinZone)
    convert_arguments(P, BZ_to_mesh(bz))
end
function convert_arguments(P::Type{<: Wireframe}, bz::AbstractBrillouinZone)
    convert_arguments(P, BZ_to_mesh(bz))
end

# Only really works with one face
function convert_arguments(P::Type{<:Lines}, bz::AbstractBrillouinZone)
    points = Point{D, Float32}.(corners(bz))
    face_points = [points[i] for f in LatPhysReciprocal.faces(bz) for i in [f..., f[1]]]
    convert_arguments(P, face_points)
end

function convert_arguments(P::Type{<:LineSegments}, bz::BZ) where {
        D, PT <: AbstractReciprocalPoint{D},
        R <: AbstractReciprocalUnitcell{PT},
        BZ <: AbstractBrillouinZone{R}
    }
    points = Point{D, Float32}.(corners(bz))
    face_points = Point{D, Float32}[]
    for f in LatPhysReciprocal.faces(bz)
        ps = [points[i] for i in [f..., f[1]] for _ in 1:2][2:end-1]
        append!(face_points, ps)
    end
    convert_arguments(P, face_points)
end

# Bandstructure (doesn't really make sense for anything but scatter, whatever)
function convert_arguments(P::Type{<:Scatter}, bs::AbstractBandstructure)
    k_point_indices = zeros(Int64, length(energies(bs))+1)
    points = Point2f0[]
    for s in 1:length(energies(bs))
        k_point_indices[s+1] = length(energies(bs)[s][1]) + k_point_indices[s] - 1
        xvals = range(k_point_indices[s], stop=k_point_indices[s+1], length=length(energies(bs)[s][1]))
        for band in energies(bs)[s]
            append!(points, Point2f0.(xvals, band))
        end
    end
    convert_arguments(P, points)
end

function convert_arguments(P::Type{<:LineSegments}, bs::AbstractBandstructure)
    k_point_indices = zeros(Int64, length(energies(bs))+1)
    points = Point2f0[]
    for s in 1:length(energies(bs))
        k_point_indices[s+1] = length(energies(bs)[s][1]) + k_point_indices[s] - 1
        xvals = range(k_point_indices[s], stop=k_point_indices[s+1], length=length(energies(bs)[s][1]))
        for band in energies(bs)[s]
            append!(points, [p for p in Point2f0.(xvals, band) for _ in 1:2][2:end-1])
        end
    end
    convert_arguments(P, points)
end

# Energy Manifold
function convert_arguments(P::Type{<:Scatter}, em::EM) where {
        LS, D, S <: AbstractSite{LS,D},
        L, UC <: AbstractUnitcell{S},
        H <: AbstractHamiltonian{L,UC},
        EM <: AbstractEnergyManifold{H}
    }
    convert_arguments(P, Point{D, Float32}.(kpoints(em)))
end
function convert_arguments(P::Type{<:Contour}, em::AbstractEnergyManifold)
    convert_arguments(P,
        range(-2pi, 2pi, length=100),
        range(-2pi, 2pi, length=100),
        (x, y) -> minimum(eigvals(matrixAtK(em.h, [x, y])))
    )
end
# Why doesn't this work?
# function convert_arguments(P::Type{<:Contour}, em::EM) where {
#         LS, S <: AbstractSite{LS,3},
#         L, UC <: AbstractUnitcell{S},
#         H <: AbstractHamiltonian{L,UC},
#         EM <: AbstractEnergyManifold{H}
#     }
#     xs = range(-2pi, 2pi, length=100)
#     ys = range(-2pi, 2pi, length=100)
#     zs = range(-2pi, 2pi, length=100)
#     values = [
#         minimum(eigvals(matrixAtK(em.h, [x, y, z])))
#         for x in xs, y in ys, z in zs
#     ]
#     convert_arguments(Contour3d, values)
# end

################################################################################
### Full extensions
################################################################################

# TODO: this probably exists in Makie...
const_lift(o::Node) = o
const_lift(x) = Node(x)
#Attributes(p::AbstractPlot; kwargs...) = merge(Attributes(; kwargs...), Attributes(p))

# 3D Lattice
########################################

# TODO: How do I get all default attributes in here and pass them?
function default_theme(scene::SceneLike, ::Type{<: Plot(LT)}) where {
        T, S <: AbstractSite{T, 3}, LT <: AbstractLattice{S}
    }
    Attributes(
        site_color = :orange,
        site_size = 0.1,
        site_marker = Sphere(Point3f0(0), 1f0),
        site_matcap = nothing,
        bond_color = :white,
        bond_width = 0.05,
        bond_matcap = nothing
    )
end

function plot!(p::Plot(LT)) where {
        T, S <: AbstractSite{T, 3}, LT <: AbstractLattice{S}
    }
    lattice = const_lift(p[1])

    bond_width  = pop!(p.attributes, :bond_width)
    bond_attributes = Attributes(
        color  = pop!(p.attributes, :bond_color),
        matcap = pop!(p.attributes, :bond_matcap)
    )

    site_attributes = Attributes(
        color       = pop!(p.attributes, :site_color),
        markersize  = pop!(p.attributes, :site_size),
        marker      = pop!(p.attributes, :site_marker),
        matcap      = pop!(p.attributes, :site_matcap)
    )
    merge!(site_attributes, p.attributes)
    meshscatter!(p, lattice; site_attributes...)

    # Need to do some explicit construction to make cylinder work
    _sites = map(sites, lattice)
    _bonds = map(bonds, lattice)

    dir = map(_sites, _bonds) do s, bs
        dir = Point{3, Float32}[]
        for b in bs
            p1 = Point{3, Float32}(point(s[from(b)]))
            p2 = Point{3, Float32}(point(s[to(b)]))
            push!(dir, p2 - p1)
        end
        dir
    end

    marker = map(bond_width) do w
        Cylinder(Point3f0(0), Point3f0(0, 0, 1), Float32(w))
    end

    merge!(bond_attributes, Attributes(
        marker = marker,
        rotations = map(vs -> normalize.(vs), dir),
        markersize = map(vs -> norm.(vs), dir)
    ))
    meshscatter!(p, _bonds, _sites; bond_attributes...)
end

# 2D Lattice
########################################

function default_theme(scene::SceneLike, ::Type{<: Plot(L)}) where {
        LS,S<:AbstractSite{LS,2},L<:AbstractLattice{S}
    }
    scatter_defaults = default_theme(scene, Scatter)
    line_defaults = default_theme(scene, Lines)
    Attributes(
        site_color = :orange,
        site_size = scatter_defaults[:markersize][],
        site_marker = scatter_defaults[:marker][],
        bond_color = :grey,
        bond_width = line_defaults[:linewidth][]
    )
end

function plot!(p::Plot(L)) where {
        LS,S<:AbstractSite{LS,2},L<:AbstractLattice{S}
    }
    bond_attributes = Attributes(
        linewidth  = pop!(p.attributes, :bond_width),
        color      = pop!(p.attributes, :bond_color)
    )

    site_attributes = Attributes(
        color       = pop!(p.attributes, :site_color),
        markersize  = pop!(p.attributes, :site_size),
        marker      = pop!(p.attributes, :site_marker),
    )

    merge!(site_attributes, p.attributes)
    scatter!(p, p[1]; site_attributes...)

    merge!(bond_attributes, p.attributes)
    linesegments!(p, p[1]; bond_attributes...)
end

# Reciprocal Path
########################################

# Not gonna do special 3D paths now...
# see lattice recipe if you want to add it
function default_theme(scene::SceneLike, ::Type{<: Plot(P)}) where {P<:AbstractReciprocalPath}
    Attributes(
        markercolor = :gray65,
        linecolor = :black
    )
end

function plot!(p::Plot(P)) where {P <: AbstractReciprocalPath}
    stripped = Attributes(p)
    mc = pop!(stripped, :markercolor)
    lc = pop!(stripped, :linecolor)
    lines!(p, p[1], color = lc; stripped...)
    scatter!(p, p[1], color = mc; stripped...)
end

# Reciprocal Path
########################################

# TODO: use spheres as corners?
function default_theme(scene::SceneLike, ::Type{<: Plot(BZ)}) where {BZ<:AbstractBrillouinZone}
    Attributes(
        BZ_corners = true
    )
end

function plot!(p::Plot(BZ)) where {BZ <: AbstractBrillouinZone}
    linesegments!(p, p[1]; Attributes(p)...)
    scatter!(p, p[1]; Attributes(p)..., visible=p[:BZ_corners])
end

# Bandstructure
########################################

# Not using lines here to make this very... interactive?
# Using scatter allows us to collect all points here, which means
# the number of bands in a bandstructure can also change without breaking the plot
function default_theme(scene::SceneLike, ::Type{<: Plot(BS)}) where {BS <: AbstractBandstructure}
    Attributes()
end

function plot!(p::Plot(BS)) where {BS <: AbstractBandstructure}
    linesegments!(p, const_lift(p[1]); Attributes(p)...)
end


text:
using NS2D
using Plots
using FFTW
using DifferentialEquations

function plot_vorticity(ω, p)

    u, v = compute_velocity(ω, p)

    ω = real(ifft(ω))

    x, y = mesh(p)

    fig = heatmap(
        x, x, ω,
        aspect_ratio=:equal,
        c=cgrad(:RdBu),
        xlims=(-p.Lx/2, p.Lx/2),
        ylims=(-p.Ly/2, p.Ly/2),
        clims=(-maximum(abs.(ω)), maximum(abs.(ω)))
    )

    return fig
end

function plot_velocity(ω, p)

    #
    u, v = compute_velocity(ω, p)

    #
    velocity_magnitude = .√(u.^2 + v.^2)

    x, y = mesh(p)

    fig = heatmap(
        x, y, velocity_magnitude,
        aspect_ratio=:equal,
        c=cgrad(:viridis),
        xlims=(-p.Lx/2, p.Lx/2),
        ylims=(-p.Ly/2, p.Ly/2),
        clims=(0, maximum(velocity_magnitude)),
    )

    return fig
end

function main(; T=1)

    # --> Simulation parameters (default).
    p = SimParams()

    # --> Isotropic enstrophy spectrum
    #     (unit energy for all Fourier modes)
    spectrum(x) = one(x)
    ω = isotropic_turbulence(spectrum, p)

    # -->
    prob = UnforcedProblem(p, ω)

    # --> Callback.
    data = SavedValues(Float64, Float64)
    cb = SavingCallback((u, t, integrator) -> maximum(abs.(u)), data)

    # --> Keyword arguments for DifferentialEquations.jl
    kwargs = Dict(
        :saveat => 0:0.1:T,
        :callback => cb
    )

    # --> Run the simulation.
    sol = simulate(prob, T ; kwargs...)

    # --> Plot the final condition.
    ω = sol[end]
    fig1 = plot_vorticity(ω, p)
    title!("Vorticity")

    fig2 = plot_velocity(ω, p)
    title!("Kinetic energy")

    fig = plot(fig1, fig2, layout=(1, 2), size=(800, 400))
    display(fig)

    return sol, data
end


text:
struct TimeInvariantKeldyshTimeGF{T, scalar} <: AbstractTimeGF{T, scalar}
  grid::KeldyshTimeGrid
  gtr::AntiHermitianToeplitzStorage{T,scalar}
  les::AntiHermitianToeplitzStorage{T,scalar}
  ξ::GFSignEnum
end

norbitals(G::TimeInvariantKeldyshTimeGF) = G.gtr.norb

function TimeInvariantKeldyshTimeGF(::Type{T}, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  nt = grid.nt

  gtr = AntiHermitianToeplitzStorage(T, nt, norb, scalar)
  les = AntiHermitianToeplitzStorage(T, nt, norb, scalar)

  TimeInvariantKeldyshTimeGF(grid, gtr, les, ξ)
end
TimeInvariantKeldyshTimeGF(grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = TimeInvariantKeldyshTimeGF(ComplexF64, grid, norb, ξ, scalar)

@inline function Base.getindex(G::TimeInvariantKeldyshTimeGF, k, l, t1::TimeGridPoint, t2::TimeGridPoint, greater=true)
  greater = t1 == t2 ? greater : heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx

  return greater ? G.gtr[k,l,i,j] : G.les[k,l,i,j]
end

function Base.setindex!(G::TimeInvariantKeldyshTimeGF, v, k, l, t1::TimeGridPoint, t2::TimeGridPoint)
  greater = heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx

  return greater ? G.gtr[k,l,i,j] = v : G.les[k,l,i,j] = v
end

function Base.similar(G::T, ξ::GFSignEnum) where T <: TimeInvariantKeldyshTimeGF
  T(G.grid, similar(G.gtr), similar(G.les), ξ)
end

function Base.zero(G::T, ξ::GFSignEnum) where T <: TimeInvariantKeldyshTimeGF
  T(G.grid, zero(G.gtr), zero(G.les), ξ)
end

function Base.similar(G::T) where T <: TimeInvariantKeldyshTimeGF
  similar(G, G.ξ)
end

function Base.zero(G::T) where T <: TimeInvariantKeldyshTimeGF
  zero(G, G.ξ)
end

function TimeDomain(G::TimeInvariantKeldyshTimeGF)
  grid = G.grid
  tplus = grid[forward_branch]
  tminus = reverse(grid[backward_branch])

  nt = length(tplus)

  points = Tuple{TimeGridPoint, TimeGridPoint}[]

  for i in 1:nt
    push!(points, (tminus[1], tminus[i])) # greater
    push!(points, (tplus[1], tminus[i])) # lesser
  end

  return TimeDomain(points)
end

function TimeInvariantKeldyshTimeGF(f::Function, ::Type{T}, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  G = TimeInvariantKeldyshTimeGF(T, grid, norb, ξ, scalar)
  D = TimeDomain(G)

  for (t1, t2) in D.points
    G[t1, t2] = f(t1, t2)
  end

  return G
end

TimeInvariantKeldyshTimeGF(f::Function, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = TimeInvariantKeldyshTimeGF(f, ComplexF64, grid, norb, ξ, scalar)

function TimeInvariantKeldyshTimeGF(dos::AbstractDOS, β, grid::KeldyshTimeGrid)
  TimeInvariantKeldyshTimeGF(grid, 1, fermionic, true) do t1, t2
    Keldysh.dos2gf(dos, β, t1.bpoint, t2.bpoint)
  end
end


text:
module AISIS10024

using ..AISIS10016

export h411, h42

function h411(Cw, Fy, Wn, ASDorLRFD)

    if ASDorLRFD==0
        StrengthFactor=1/1.67
    elseif ASDorLRFD==1
        StrengthFactor=0.90
    else
        StrengthFactor=1.0   #to just get nominal strength
    end

    #calculate bimoment yield strength
    Bn=Cw.*Fy./Wn

    # #calculate nominal bimoment strength including local buckling
    # #for now use DSM local buckling flexural curve
    # Bn, notused = AISIS10016.f321(By, Bcrℓ, ASDorLRFD)

    eBn = Bn * StrengthFactor

    return Bn, eBn

end

function h42(Mxbar,Mybar,Bbar,Mybar_freeflange, Maxℓo,Mayℓo,Ba, Mayℓo_freeflange)

    ActionMx = abs.(Mxbar./(Maxℓo))
    ActionMy = abs.(Mybar./(Mayℓo))
    ActionB = abs.(Bbar./(Ba))
    ActionMy_freeflange = abs.(Mybar_freeflange./(Mayℓo_freeflange))

    Interaction = ActionMx .+ ActionMy .+ ActionB .+ ActionMy_freeflange

    return ActionMx, ActionMy, ActionB, ActionMy_freeflange, Interaction

end

end #module


text:
# The methods below are specialized on the mortar type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::P4estMesh{3}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
  # TODO: Taal compare performance of different types
  fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2), 4)
                    for _ in 1:Threads.nthreads()]

  fstar_tmp_threaded = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                         for _ in 1:Threads.nthreads()]
  u_threaded          = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                         for _ in 1:Threads.nthreads()]

  (; fstar_threaded, fstar_tmp_threaded, u_threaded)
end

function prolong2interfaces!(cache, u,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack interfaces = cache

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for interface in eachinterface(dg, cache)
    primary_element   = interfaces.element_ids[1, interface]
    secondary_element = interfaces.element_ids[2, interface]

    primary_indices   = interfaces.node_indices[1, interface]
    secondary_indices = interfaces.node_indices[2, interface]

    # Use Tuple `node_indices` and `evaluate_index` to copy values
    # from the correct face and in the correct orientation
    for j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        interfaces.u[1, v, i, j, interface] = u[v, evaluate_index(primary_indices, size_, 1, i, j),
                                                   evaluate_index(primary_indices, size_, 2, i, j),
                                                   evaluate_index(primary_indices, size_, 3, i, j),
                                                   primary_element]

        interfaces.u[2, v, i, j, interface] = u[v, evaluate_index(secondary_indices, size_, 1, i, j),
                                                   evaluate_index(secondary_indices, size_, 2, i, j),
                                                   evaluate_index(secondary_indices, size_, 3, i, j),
                                                   secondary_element]
      end
    end
  end

  return nothing
end

function calc_interface_flux!(surface_flux_values,
                              mesh::P4estMesh{3},
                              nonconservative_terms::Val{false},
                              equations, surface_integral, dg::DG, cache)
  @unpack surface_flux = surface_integral
  @unpack u, element_ids, node_indices = cache.interfaces

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for interface in eachinterface(dg, cache)
    # Get neighboring elements
    primary_element   = element_ids[1, interface]
    secondary_element = element_ids[2, interface]

    primary_indices   = node_indices[1, interface]
    secondary_indices = node_indices[2, interface]

    primary_direction   = indices2direction(primary_indices)
    secondary_direction = indices2direction(secondary_indices)

    # Use Tuple `node_indices` and `evaluate_index` to access node indices
    # at the correct face and in the correct orientation to get normal vectors
    for j in eachnode(dg), i in eachnode(dg)
      u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, j, interface)

      normal_vector = get_normal_vector(primary_direction, cache,
                                        evaluate_index(primary_indices, size_, 1, i, j),
                                        evaluate_index(primary_indices, size_, 2, i, j),
                                        evaluate_index(primary_indices, size_, 3, i, j),
                                        primary_element)

      flux_ = surface_flux(u_ll, u_rr, normal_vector, equations)

      # Use Tuple `node_indices` and `evaluate_index_surface` to copy flux
      # to left and right element storage in the correct orientation
      for v in eachvariable(equations)
        surf_i = evaluate_index_surface(primary_indices, size_, 1, i, j)
        surf_j = evaluate_index_surface(primary_indices, size_, 2, i, j)
        surface_flux_values[v, surf_i, surf_j, primary_direction, primary_element] = flux_[v]

        surf_i = evaluate_index_surface(secondary_indices, size_, 1, i, j)
        surf_j = evaluate_index_surface(secondary_indices, size_, 2, i, j)
        surface_flux_values[v, surf_i, surf_j, secondary_direction, secondary_element] = -flux_[v]
      end
    end
  end

  return nothing
end

function prolong2boundaries!(cache, u,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for boundary in eachboundary(dg, cache)
    element       = boundaries.element_ids[boundary]
    node_indices  = boundaries.node_indices[boundary]

    # Use Tuple `node_indices` and `evaluate_index` to copy values
    # from the correct face and in the correct orientation
    for j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        boundaries.u[v, i, j, boundary] = u[v, evaluate_index(node_indices, size_, 1, i, j),
                                               evaluate_index(node_indices, size_, 2, i, j),
                                               evaluate_index(node_indices, size_, 3, i, j),
                                               element]
      end
    end
  end

  return nothing
end

function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  @unpack surface_flux_values, node_coordinates = cache.elements
  @unpack surface_flux = surface_integral

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for local_index in eachindex(boundary_indexing)
    # Use the local index to get the global boundary index from the pre-sorted list
    boundary = boundary_indexing[local_index]

    element       = boundaries.element_ids[boundary]
    node_indices  = boundaries.node_indices[boundary]
    direction     = indices2direction(node_indices)

    # Use Tuple `node_indices` and `evaluate_index` to access node indices
    # at the correct face and in the correct orientation to get normal vectors
    for j in eachnode(dg), i in eachnode(dg)
      node_i = evaluate_index(node_indices, size_, 1, i, j)
      node_j = evaluate_index(node_indices, size_, 2, i, j)
      node_k = evaluate_index(node_indices, size_, 3, i, j)

      # Extract solution data from boundary container
      u_inner = get_node_vars(boundaries.u, equations, dg, i, j, boundary)

      # Outward-pointing normal vector
      normal_vector = get_normal_vector(direction, cache, node_i, node_j, node_k, element)

      # Coordinates at boundary node
      x = get_node_coords(node_coordinates, equations, dg, node_i, node_j, node_k, element)

      flux_ = boundary_condition(u_inner, normal_vector, x, t, surface_flux, equations)

      # Use Tuple `node_indices` and `evaluate_index_surface` to copy flux
      # to left and right element storage in the correct orientation
      for v in eachvariable(equations)
        surf_i = evaluate_index_surface(node_indices, size_, 1, i, j)
        surf_j = evaluate_index_surface(node_indices, size_, 2, i, j)
        surface_flux_values[v, surf_i, surf_j, direction, element] = flux_[v]
      end
    end
  end
end

function prolong2mortars!(cache, u,
                          mesh::P4estMesh{3}, equations,
                          mortar_l2::LobattoLegendreMortarL2,
                          surface_integral, dg::DGSEM)
  # temporary buffer for projections
  @unpack fstar_tmp_threaded = cache
  @unpack element_ids, node_indices = cache.mortars

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for mortar in eachmortar(dg, cache)
    fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

    small_indices = node_indices[1, mortar]
    large_indices = node_indices[2, mortar]

    # Copy solution small to small
    for pos in 1:4
      for j in eachnode(dg), i in eachnode(dg)
        for v in eachvariable(equations)
          # Use Tuple `node_indices` and `evaluate_index` to copy values
          # from the correct face and in the correct orientation
          cache.mortars.u[1, v, pos, i, j, mortar] = u[v, evaluate_index(small_indices, size_, 1, i, j),
                                                          evaluate_index(small_indices, size_, 2, i, j),
                                                          evaluate_index(small_indices, size_, 3, i, j),
                                                          element_ids[pos, mortar]]
        end
      end
    end

    # Buffer to copy solution values of the large element in the correct orientation
    # before interpolating
    u_buffer = cache.u_threaded[Threads.threadid()]

    # Copy solution of large element face to buffer in the correct orientation
    for j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        # Use Tuple `node_indices` and `evaluate_index` to copy values
        # from the correct face and in the correct orientation
        u_buffer[v, i, j] = u[v, evaluate_index(large_indices, size_, 1, i, j),
                                 evaluate_index(large_indices, size_, 2, i, j),
                                 evaluate_index(large_indices, size_, 3, i, j),
                                 element_ids[5, mortar]]
      end
    end

    # Interpolate large element face data from buffer to small face locations
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, :, mortar),
                            mortar_l2.forward_lower,
                            mortar_l2.forward_lower,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, :, mortar),
                            mortar_l2.forward_upper,
                            mortar_l2.forward_lower,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 3, :, :, mortar),
                            mortar_l2.forward_lower,
                            mortar_l2.forward_upper,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 4, :, :, mortar),
                            mortar_l2.forward_upper,
                            mortar_l2.forward_upper,
                            u_buffer,
                            fstar_tmp)
  end

  return nothing
end

function calc_mortar_flux!(surface_flux_values,
                           mesh::P4estMesh{3},
                           nonconservative_terms::Val{false}, equations,
                           mortar_l2::LobattoLegendreMortarL2,
                           surface_integral, dg::DG, cache)
  @unpack u, element_ids, node_indices = cache.mortars
  @unpack fstar_threaded, fstar_tmp_threaded = cache
  @unpack surface_flux = surface_integral

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  @threaded for mortar in eachmortar(dg, cache)
    # Choose thread-specific pre-allocated container
    fstar = fstar_threaded[Threads.threadid()]
    fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

    small_indices = node_indices[1, mortar]
    small_direction = indices2direction(small_indices)

    # Use Tuple `node_indices` and `evaluate_index` to access node indices
    # at the correct face and in the correct orientation to get normal vectors
    for pos in 1:4
      for j in eachnode(dg), i in eachnode(dg)
        u_ll, u_rr = get_surface_node_vars(u, equations, dg, pos, i, j, mortar)

        normal_vector = get_normal_vector(small_direction, cache,
                                          evaluate_index(small_indices, size_, 1, i, j),
                                          evaluate_index(small_indices, size_, 2, i, j),
                                          evaluate_index(small_indices, size_, 3, i, j),
                                          element_ids[pos, mortar])

        flux_ = surface_flux(u_ll, u_rr, normal_vector, equations)

        # Copy flux to buffer
        set_node_vars!(fstar, flux_, equations, dg, i, j, pos)
      end
    end

    # Buffer to interpolate flux values of the large element to before copying
    # in the correct orientation
    u_buffer = cache.u_threaded[Threads.threadid()]

    mortar_fluxes_to_elements!(surface_flux_values,
                               mesh, equations, mortar_l2, dg, cache,
                               mortar, fstar, u_buffer, fstar_tmp)
  end

  return nothing
end

@inline function mortar_fluxes_to_elements!(surface_flux_values,
                                            mesh::P4estMesh{3}, equations,
                                            mortar_l2::LobattoLegendreMortarL2,
                                            dg::DGSEM, cache, mortar, fstar, u_buffer, fstar_tmp)
  @unpack element_ids, node_indices = cache.mortars

  small_indices  = node_indices[1, mortar]
  large_indices  = node_indices[2, mortar]

  small_direction = indices2direction(small_indices)
  large_direction = indices2direction(large_indices)

  size_ = (nnodes(dg), nnodes(dg), nnodes(dg))

  # Copy solution small to small
  for pos in 1:4
    for j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        # Use Tuple `node_indices` and `evaluate_index_surface` to copy flux
        # to left and right element storage in the correct orientation
        surface_index1 = evaluate_index_surface(small_indices, size_, 1, i, j)
        surface_index2 = evaluate_index_surface(small_indices, size_, 2, i, j)
        surface_flux_values[v, surface_index1, surface_index2, small_direction,
                            element_ids[pos, mortar]] = fstar[v, i, j, pos]
      end
    end
  end

  large_element = element_ids[5, mortar]

  # Project small fluxes to large element.
  multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_lower, mortar_l2.reverse_lower,
    view(fstar, .., 1),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_upper, mortar_l2.reverse_lower,
    view(fstar, .., 2),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_lower, mortar_l2.reverse_upper,
    view(fstar, .., 3),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_upper, mortar_l2.reverse_upper,
    view(fstar, .., 4),
    fstar_tmp)

  # The flux is calculated in the outward direction of the small elements,
  # so the sign must be switched to get the flux in outward direction
  # of the large element.
  # The contravariant vectors of the large element (and therefore the normal vectors
  # of the large element as well) are four times as large as the contravariant vectors
  # of the small elements. Therefore, the flux need to be scaled by a factor of 4
  # to obtain the flux of the large element.
  u_buffer .*= -4

  # Copy interpolated flux values from buffer to large element face in the correct orientation
  for j in eachnode(dg), i in eachnode(dg)
    for v in eachvariable(equations)
      # Use Tuple `node_indices` and `evaluate_index_surface` to copy flux
      # to surface flux storage in the correct orientation
      surface_index1 = evaluate_index_surface(large_indices, size_, 1, i, j)
      surface_index2 = evaluate_index_surface(large_indices, size_, 2, i, j)
      surface_flux_values[v, surface_index1, surface_index2,
                          large_direction, large_element] = u_buffer[v, i, j]
    end
  end

  return nothing
end

function calc_surface_integral!(du, u,
                                mesh::P4estMesh{3},
                                equations,
                                surface_integral::SurfaceIntegralWeakForm,
                                dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  # Note that all fluxes have been computed with outward-pointing normal vectors
  @threaded for element in eachelement(dg, cache)
    for m in eachnode(dg), l in eachnode(dg)
      for v in eachvariable(equations)
        # surface at -x
        du[v, 1,          l, m, element] += (surface_flux_values[v, l, m, 1, element]
                                             * boundary_interpolation[1, 1])
        # surface at +x
        du[v, nnodes(dg), l, m, element] += (surface_flux_values[v, l, m, 2, element]
                                             * boundary_interpolation[nnodes(dg), 2])
        # surface at -y
        du[v, l, 1,          m, element] += (surface_flux_values[v, l, m, 3, element]
                                             * boundary_interpolation[1, 1])
        # surface at +y
        du[v, l, nnodes(dg), m, element] += (surface_flux_values[v, l, m, 4, element]
                                             * boundary_interpolation[nnodes(dg), 2])
        # surface at -z
        du[v, l, m, 1,          element] += (surface_flux_values[v, l, m, 5, element]
                                             * boundary_interpolation[1,          1])
        # surface at +z
        du[v, l, m, nnodes(dg), element] += (surface_flux_values[v, l, m, 6, element]
                                             * boundary_interpolation[nnodes(dg), 2])
      end
    end
  end

  return nothing
end


text:
immutable Weibull <: ContinuousUnivariateDistribution
    shape::Float64
    scale::Float64
    function Weibull(sh::Real, sc::Real)
    	zero(sh) < sh && zero(sc) < sc || error("Both shape and scale must be positive")
    	new(float64(sh), float64(sc))
    end
end

Weibull(sh::Real) = Weibull(sh, 1.0)

@_jl_dist_2p Weibull weibull

function cdf(d::Weibull, x::Real)
	if 0.0 < x
		return 1.0 - exp(-((x / d.scale)^d.shape))
	else
		0.0
	end
end

function entropy(d::Weibull)
    k, l = d.shape, d.scale
    return ((k - 1.0) / k) * -digamma(1.0) + log(l / k) + 1.0
end

insupport(::Weibull, x::Real) = zero(x) <= x < Inf
insupport(::Type{Weibull}, x::Real) = zero(x) <= x < Inf

function kurtosis(d::Weibull)
    λ, k = d.scale, d.shape
    μ = mean(d)
    σ = std(d)
    γ = skewness(d)
    den = λ^4 * gamma(1.0 + 4.0 / k) -
          4.0 * γ * σ^3 * μ -
          6.0 * μ^2 * σ^2 - μ^4
    num = σ^4
    return den / num - 3.0
end

mean(d::Weibull) = d.scale * gamma(1.0 + 1.0 / d.shape)

median(d::Weibull) = d.scale * log(2.0)^(1.0 / d.shape)

function modes(d::Weibull)
    if d.shape <= 1.0
        return [0.0]
    else
        return [d.scale * ((d.shape - 1.0) / d.shape)^(1.0 / d.shape)]
    end
end

function skewness(d::Weibull)
    tmp = gamma(1.0 + 3.0 / d.shape) * d.scale^3
    tmp -= 3.0 * mean(d) * var(d)
    tmp -= mean(d)^3
    return tmp / std(d)^3
end

var(d::Weibull) = d.scale^2 * gamma(1.0 + 2.0 / d.shape) - mean(d)^2

function mode(d::Weibull)
    inv_k = 1.0 / d.shape
    d.shape > 1.0 ? d.scale * (1.0 - inv_k) ^ inv_k : 0.0
end

modes(d::Weibull) = [mode(d)]


text:
"""
    @check_args(D, cond)

macro for argument checking.

source: https://github.com/JuliaStats/Distributions.jl/blob/v0.25.32/src/utils.jl#L3
"""
macro check_args(D, cond)
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(D)), ": the condition ", $(string(cond)), " is not satisfied.")))
        end
    end
end


text:
##### Alias Table #####

immutable AliasTable <: Sampleable{Univariate,Discrete}
    accept::Vector{Float64}
    alias::Vector{Int}
    isampler::RandIntSampler
end
numcategories(s::AliasTable) = length(s.accept)

function make_alias_table!(a::AbstractVector{Float64})
    # input probabilities via a, which is then
    # overriden by acceptance probabilities
    # and used internally by the AliasTable instance
    #
    n = length(a)
    for i=1:n
        @inbounds a[i] *= n
    end
    alias = Array(Int,n)
    larges = Int[]
    smalls = Int[]

    for i = 1:n
        acci = a[i]
        if acci > 1.0
            push!(larges,i)
        elseif acci < 1.0
            push!(smalls,i)
        end
    end

    while !isempty(larges) && !isempty(smalls)
        s = pop!(smalls)
        l = pop!(larges)
        @inbounds alias[s] = l
        @inbounds a[l] = (a[l] - 1.0) + a[s]
        if a[l] > 1
            push!(larges,l)
        else
            push!(smalls,l)
        end
    end

    # this loop should be redundant, except for rounding
    for s in smalls
        @inbounds a[s] = 1.0
    end
    AliasTable(a, alias, RandIntSampler(n))
end

function AliasTable{T<:Real}(probs::AbstractVector{T})
    n = length(probs)
    n > 0 || error("The input probability vector is empty.")
    a = Array(Float64, n)
    for i = 1:n
        @inbounds a[i] = probs[i]
    end
    make_alias_table!(a)
end

rand(s::AliasTable) =
    (i = rand(s.isampler); u = rand(); @inbounds r = u < s.accept[i] ? i : s.alias[i]; r)

show(io::IO, s::AliasTable) = @printf(io, "AliasTable with %d entries", numcategories(s))