text:
# Array info
for (m, f) in ((:Base, :size), (:Base, :axes), (:Base, :firstindex), (:Base, :lastindex))
    @eval begin
        @inline $m.$f(A::AbstractDimArray, dims::AllDims) = $m.$f(A, dimnum(A, dims))
    end
end

# Reducing methods

# With a function arg version
for (m, f) in ((:Base, :sum), (:Base, :prod), (:Base, :maximum), (:Base, :minimum),
                     (:Base, :extrema), (:Statistics, :mean))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; dims=:, kw...) = $_f(A::AbstractDimArray, dims; kw...)
        $m.$f(f, A::AbstractDimArray; dims=:, kw...) = $_f(f, A::AbstractDimArray, dims; kw...)
        # Local dispatch methods
        # - Return a reduced DimArray
        $_f(A::AbstractDimArray, dims; kw...) =
            rebuild(A, $m.$f(parent(A); dims=dimnum(A, _astuple(dims)), kw...), reducedims(A, dims))
        $_f(f, A::AbstractDimArray, dims; kw...) =
            rebuild(A, $m.$f(f, parent(A); dims=dimnum(A, _astuple(dims)), kw...), reducedims(A, dims))
        # - Return a scalar
        $_f(A::AbstractDimArray, dims::Colon; kw...) = $m.$f(parent(A); dims, kw...)
        $_f(f, A::AbstractDimArray, dims::Colon; kw...) = $m.$f(f, parent(A); dims, kw...)
    end
end
# With no function arg version
for (m, f) in ((:Statistics, :std), (:Statistics, :var))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; corrected::Bool=true, mean=nothing, dims=:) =
            $_f(A, corrected, mean, dims)
        # Local dispatch methods - Returns a reduced array
        $_f(A::AbstractDimArray, corrected, mean, dims) =
            rebuild(A, $m.$f(parent(A); corrected=corrected, mean=mean, dims=dimnum(A, _astuple(dims))), reducedims(A, dims))
        # - Returns a scalar
        $_f(A::AbstractDimArray, corrected, mean, dims::Colon) =
            $m.$f(parent(A); corrected=corrected, mean=mean, dims=:)
    end
end
for (m, f) in ((:Statistics, :median), (:Base, :any), (:Base, :all))
    _f = Symbol('_', f)
    @eval begin
        # Base methods
        $m.$f(A::AbstractDimArray; dims=:) = $_f(A, dims)
        # Local dispatch methods - Returns a reduced array
        $_f(A::AbstractDimArray, dims) =
            rebuild(A, $m.$f(parent(A); dims=dimnum(A, _astuple(dims))), reducedims(A, dims))
        # - Returns a scalar
        $_f(A::AbstractDimArray, dims::Colon) = $m.$f(parent(A); dims=:)
    end
end

# These are not exported but it makes a lot of things easier using them
function Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractDimArray, dims)
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, _astuple(dims))), reducedims(A, dims))
end
function Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractDimArray, dims::Colon)
    Base._mapreduce_dim(f, op, nt, parent(A), dims)
end
function Base._mapreduce_dim(f, op, nt, A::AbstractDimArray, dims)
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
end
@static if VERSION >= v"1.6"
    function Base._mapreduce_dim(f, op, nt::Base._InitialValue, A::AbstractDimArray, dims)
        rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
    end
    function Base._mapreduce_dim(f, op, nt::Base._InitialValue, A::AbstractDimArray, dims::Colon)
        Base._mapreduce_dim(f, op, nt, parent(A), dims)
    end
end

# TODO: Unfortunately Base/accumulate.jl kw methods all force dims to be Integer.
# accumulate wont work unless that is relaxed, or we copy half of the file here.
# Base._accumulate!(op, B, A, dims::AllDims, init::Union{Nothing, Some}) =
    # Base._accumulate!(op, B, A, dimnum(A, dims), init)

# Dimension dropping

function Base.dropdims(A::AbstractDimArray; dims)
    dims = DD.dims(A, dims)
    data = Base.dropdims(parent(A); dims=dimnum(A, dims))
    rebuildsliced(A, data, _dropinds(A, dims))
end

@inline _dropinds(A, dims::DimTuple) = dims2indices(A, map(d -> basetypeof(d)(1), dims))
@inline _dropinds(A, dim::Dimension) = dims2indices(A, basetypeof(dim)(1))

# Function application

function Base.map(f, As::AbstractDimArray...)
    comparedims(As...)
    newdata = map(f, map(parent, As)...)
    rebuild(first(As); data=newdata)
end

function Base.mapslices(f, A::AbstractDimArray; dims=1, kw...)
    dimnums = dimnum(A, _astuple(dims))
    data = mapslices(f, parent(A); dims=dimnums, kw...)
    rebuild(A, data)
end

# This is copied from base as we can't efficiently wrap this function
# through the kw with a rebuild in the generator. Doing it this way
# also makes it faster to use a dim than an integer.
function Base.eachslice(A::AbstractDimArray; dims=1, kw...)
    if dims isa Tuple && length(dims) != 1
        throw(ArgumentError("only single dimensions are supported"))
    end
    dim = first(dimnum(A, dims))
    dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
    return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
end

# Duplicated dims

for fname in (:cor, :cov)
    @eval function Statistics.$fname(A::AbstractDimArray{<:Any,2}; dims=1, kw...)
        newdata = Statistics.$fname(parent(A); dims=dimnum(A, dims), kw...)
        removed_idx = dimnum(A, dims)
        newrefdims = $dims(A)[removed_idx]
        newdim = $dims(A)[3 - removed_idx]
        rebuild(A, newdata, (newdim, newdim), (newrefdims,))
    end
end

# Rotations

struct _Rot90 end
struct _Rot180 end
struct _Rot270 end
struct _Rot360 end

Base.rotl90(A::AbstractDimMatrix) = rebuild(A, rotl90(parent(A)), _rot(_Rot90(), dims(A)))
Base.rotl90(A::AbstractDimMatrix, k::Integer) =
    rebuild(A, rotl90(parent(A), k), _rot(_rottype(k), dims(A)))

Base.rotr90(A::AbstractDimMatrix) = rebuild(A, rotr90(parent(A)), _rot(_Rot270(), dims(A)))
Base.rotr90(A::AbstractDimMatrix, k::Integer) =
    rebuild(A, rotr90(parent(A), k), _rot(_rottype(-k), dims(A)))

Base.rot180(A::AbstractDimMatrix) = rebuild(A, rot180(parent(A)), _rot(_Rot180(), dims(A)))

# Not type stable - but we have to lose type stability somewhere when
# dims are being swapped, by an Int value, so it may as well be here
function _rottype(k)
    k = mod(k, 4)
    if k == 1
        _Rot90()
    elseif k == 2
        _Rot180()
    elseif k == 3
        _Rot270()
    else
        _Rot360()
    end
end

_rot(::_Rot90, (dima, dimb)) = (flip(Relation, dimb), dima)
_rot(::_Rot180, dims) = map(d -> flip(Relation, d), dims)
_rot(::_Rot270, (dima, dimb)) = (dimb, flip(Relation, dima))
_rot(::_Rot360, dims) = dims

# Dimension reordering

for (pkg, fname) in [(:Base, :permutedims), (:Base, :adjoint),
                     (:Base, :transpose), (:LinearAlgebra, :Transpose)]
    @eval begin
        @inline $pkg.$fname(A::AbstractDimArray{<:Any,2}) =
            rebuild(A, $pkg.$fname(parent(A)), reverse(dims(A)))
        @inline $pkg.$fname(A::AbstractDimArray{<:Any,1}) =
            rebuild(A, $pkg.$fname(parent(A)), (AnonDim(Base.OneTo(1)), dims(A)...))
    end
end
for fname in [:permutedims, :PermutedDimsArray]
    @eval begin
        @inline function Base.$fname(A::AbstractDimArray, perm)
            rebuild(A, $fname(parent(A), dimnum(A, Tuple(perm))), sortdims(dims(A), Tuple(perm)))
        end
    end
end

# Concatenation
function Base._cat(catdims::Union{Int,Base.Dims}, Xin::AbstractDimArray...)
    Base._cat(dims(first(Xin), catdims), Xin...)
end
function Base._cat(catdims::Tuple, Xin::AbstractDimArray...)
    A1 = first(Xin)
    comparedims(map(x -> otherdims(x, catdims), Xin)...)
    newcatdims = map(catdims) do catdim
        if all(x -> hasdim(x, catdim), Xin)
            # We concatenate an existing dimension
            newcatdim = if mode(A1, catdim) isa NoIndex
                # Colon will be converted to array axis in `formatdims`
                rebuild(catdim; val=:)
            else
                # vcat the index for the catdim in each of Xin
                reduce(vcat, map(x -> dims(x, catdim), Xin))
            end
        else
            # Concatenate new dims
            if all(map(x -> hasdim(refdims(x), catdim), Xin))
                # vcat the refdims
                reduce(vcat, map(x -> refdims(x, catdim), Xin))
            else
                # Use the catdim as the new dimension
                catdim
            end
        end
    end
    inserted_dims = dims(newcatdims, dims(A1))
    appended_dims = otherdims(newcatdims, inserted_dims)
    updated_dims = setdims(dims(A1), inserted_dims)
    newdims = (updated_dims..., appended_dims...)

    inserted_dnums = dimnum(A1, inserted_dims)
    appended_dnums = ntuple(i -> i + length(dims(A1)), length(appended_dims))
    cat_dnums = (inserted_dnums..., appended_dnums...)

    newrefdims = otherdims(refdims(A1), newcatdims)
    newA = Base._cat(cat_dnums, map(data, Xin)...)
    rebuild(A1, newA, formatdims(newA, newdims), newrefdims)
end
function Base._cat(catdim::DimType, Xin::AbstractDimArray...)
    Base._cat(catdim(; mode=NoIndex()), Xin...)
end
function Base._cat(catdim::DimOrDimType, Xin::AbstractDimArray...)
    Base._cat((catdim,), Xin...)
end
# function Base._cat(catdim::DimOrDimType, Xin::AbstractDimArray...)
#     A1 = first(Xin)
#     comparedims(map(x -> otherdims(x, catdim), Xin)...)
#     newrefdims = refdims(A1)
#     if all(A -> hasdim(A, catdim), Xin)
#         # We concatenate an existing dimension
#         if mode(A1, catdim) isa NoIndex
#             # Make new dims for NoIndex
#             dnum = dimnum(A1, catdim)
#             newdims = setdims(dims(A1), catdim)
#         else
#             # Use existing dims
#             dnum = dimnum(A1, catdim)
#             # cat the catdim, ignore others
#             newcatdim = reduce(vcat, map(x -> dims(x, catdim), Xin))
#             newdims = setdims(dims(A1), newcatdim)
#         end
#     else
#         # Concatenate a new dimension
#         if all(map(x -> hasdim(refdims(x), catdim), Xin))
#             # vcat the refdims
#             newcatdim = reduce(vcat, map(x -> refdims(x, catdim), Xin))
#             newdims = (dims(A1)..., newcatdim)
#             dnum = ndims(A1) + 1
#             newrefdims = otherdims(refdims(A1), catdim)
#         else
#             newdims = (dims(A1)..., catdim)
#             dnum = ndims(A1) + 1
#         end
#     end
#     newA = Base._cat(dnum, map(parent, Xin)...)
#     rebuild(A1; data=newA, dims=formatdims(newA, newdims), refdims=newrefdims)
# end

function Base.vcat(dims::Dimension...)
    newmode = _vcat_modes(mode(dims)...)
    newindex = _vcat_index(newmode, map(val, dims)...)
    rebuild(dims[1], newindex, newmode)
end

# IndexModes may need adjustment for `cat`
_vcat_modes(modes::IndexMode...) = first(modes)
function _vcat_modes(modes::AbstractSampled...)
    _vcat_modes(sampling(first(modes)), span(first(modes)), modes...)
end
function _vcat_modes(::Any, ::Regular, modes...)
    _step = step(first(modes))
    map(modes) do mode
        step(span(mode)) == _step || error("Step sizes $(step(span(mode))) and $_step do not match ")
    end
    first(modes)
end
function _vcat_modes(::Intervals, ::Irregular, modes...)
    allbounds = map(bounds ∘ span, modes)
    newbounds = minimum(map(first, allbounds)), maximum(map(last, allbounds))
    rebuild(modes[1]; span=Irregular(newbounds))
end
_vcat_modes(::Points, ::Irregular, modes...) =
    rebuild(first(modes); span=Irregular(nothing, nothing))

# Index vcat depends on mode: NoIndex is always Colon()
_vcat_index(mode::NoIndex, A...) = OneTo(sum(map(length, A)))
# TODO: handle vcat OffsetArrays?
# Otherwise just vcat. TODO: handle order breaking vcat?
_vcat_index(mode::IndexMode, A...) = vcat(A...)

Base.inv(A::AbstractDimArray{T,2}) where T =
    rebuild(A, inv(parent(A)), reverse(map(d -> flip(IndexOrder, d), dims(A))))

# Index breaking

# TODO: change the index and traits of the reduced dimension and return a DimArray.
Base.unique(A::AbstractDimArray; dims::Union{DimOrDimType,Colon}=:) = _unique(A, dims)
Base.unique(A::AbstractDimArray{<:Any,1}) = unique(parent(A))

_unique(A::AbstractDimArray, dims::DimOrDimType) =
    unique(parent(A); dims=dimnum(A, dims))
_unique(A::AbstractDimArray, dims::Colon) = unique(parent(A); dims=:)

Base.diff(A::AbstractDimVector; dims=1) = _diff(A, dimnum(A, dims))
Base.diff(A::AbstractDimArray; dims) = _diff(A, dimnum(A, dims))

@inline function _diff(A::AbstractDimArray{<:Any,N}, dims::Integer) where {N}
    r = axes(A)
    # Copied from Base.diff
    r0 = ntuple(i -> i == dims ? UnitRange(1, last(r[i]) - 1) : UnitRange(r[i]), N)
    rebuildsliced(A, diff(parent(A); dims=dimnum(A, dims)), r0)
end

# Forward `replace` to parent objects
function Base._replace!(new::Base.Callable, res::AbstractDimArray, A::AbstractDimArray, count::Int)
    Base._replace!(new, parent(res), parent(A), count)
    return res
end
function Base._replace!(new::Base.Callable, res::AbstractArray, A::AbstractDimArray, count::Int)
    Base._replace!(new, res, parent(A), count)
    return res
end
function Base._replace!(new::Base.Callable, res::AbstractDimArray, A::AbstractArray, count::Int)
    Base._replace!(new, parent(res), A, count)
    return res
end


text:
using Printf

f = open( "target_a.sam", "r" )
n = countlines( f )  # 计算数据行数
seekstart( f )       # 返回文件起始位置

k_dict = Dict{String, String}()

for i = 1:n
    line = readline(f)

    if !(line[1]=='@')
        s = split( line, "\t" )
        chr_name = s[3]
        if s[7] == "="
            pos = parse(Int64, s[4])
            Interval_S = parse(Int64, s[9])

            if Interval_S > 0

                is_IN = false
                t_k=""
                for p = (pos-2):(pos+2)
                    for ii = (Interval_S-2):(Interval_S+2)

                        k_s = string(chr_name,":", p, "+", ii)
                        if haskey(k_dict, k_s)
                            is_IN = true
                            t_k = k_dict[k_s]
                            break
                        end
                    end
                    if is_IN == true
                        break
                    end

                end

                if is_IN == false
                    @printf("%s\t%s\n", string(chr_name, ":", pos, "+", Interval_S), line )
                    for p = (pos-2):(pos+2)
                        for ii = (Interval_S-2):(Interval_S+2)

                            k_s = string(chr_name,":", p, "+", ii)

                            k_dict[k_s] = string(chr_name, ":", pos, "+", Interval_S)     ## add region to dict
                        end
                    end

                else
                    @printf("%s\t%s\n", t_k, line )
                end
                #println(Interval_S)
            end
        end

    else
        continue
    end

end

#for (key, value) in k_dict
#    println(key," => ", value)
#end

ext: jl


text:
using Bokeh
using Compat
hold(true)
Bokeh.debug(true)

plotfile("glyphs_custom.html")

x = linspace(0, 2pi, 10)
y = 1 + 0.5sin(x)
# you can get a list of markers from https://github.com/bokeh/bokeh/blob/master/bokeh/models/markers.py
glyph_with_everything = Glyph(:DiamondCross,
							  linecolor="red",
							  linewidth=2,
							  linealpha=0.7,
							  dash=[4,1],
							  fillcolor="#88f",
							  fillalpha=0.5,
							  size=50)
plot([y y*-1], [Glyph(:Asterisk), glyph_with_everything])
palette = ["#000000", "#190000", "#320000", "#4b0000", "#640000", "#7d0000", "#960000", "#af0000", "#c80000", "#e10000", "#ff0000"]
low = minimum(y)
high= maximum(y)
y_ints = Int64[floor(9.99*(x-low)/(high-low)) + 1 for x in y]  # round to nearest value in range 1:11
y_colours = [palette[i] for i in y_ints]

data = @compat Dict{Symbol, Vector}(:x => x, :y => y, :y_colours => y_colours, :y_sizes => y)
fields = @compat Dict{Symbol, Symbol}(:fillcolor=>:y_colours, :size=>:y_sizes)
varying_glyph = Glyph(:Circle, size=30, linewidth=3, fields=fields)
plot([BokehDataSet(data, varying_glyph)])
showplot()

max_stars_repo_name: Laeeth/dbokeh

max_issues_repo_name: Laeeth/dbokeh

max_forks_repo_name: Laeeth/dbokeh


text:
ENV["GKSwstype"]="100"

using Test
using WaterWaves1D

include("./test_loadsave.jl")
include("./testtypes.jl")
include("./testmodels.jl")
include("./testsolvers.jl")
include("./testinit.jl")

include("./test_examples.jl")
include("./test_notebooks.jl")

param = ( ϵ  = 1/2, μ = 1)
paramX= ( N  = 2^8, L  = 10)
paramT= ( T  = 0.1, dt = 0.01)

@testset "Parameters" begin

    @test param.ϵ  == 0.5
    @test param.μ  == 1
    @test paramX.N  == 256
    @test paramX.L  == 10
    @test paramT.T  == 0.1
    @test paramT.dt == 0.01

end

@testset "Interpolate" begin
    f(x)=exp.(-x.^2);
    mesh=Mesh(paramX);x=mesh.x;f1=f(x);
    mesh2,f2 = interpolate(mesh,f1;n=2)
    @test f2≈f(mesh2.x)
    @test f2[1:2:end]≈f1
    x3 = [1 2.5 3]
    f3 = interpolate(mesh,f1,x3)
    @test f3≈f(x3)
    x4 = [1 2 3]
    f4 = interpolate(mesh,f1,x4;fast=true)
    @test f4≈f(x4)
end

### Preserved quantities
models=AbstractModel[];
#push!(models,WaterWaves(param;mesh=Mesh(paramX))); # will error if uncommented (because water waves generates non-equidistant meshes)
#push!(models,Matsuno(param;mesh=Mesh(paramX))); # will fail if uncommented (because Matsuno model does not preserve momentum and energy)

push!(models,Airy(param;mesh=Mesh(paramX)));
push!(models,SaintVenant(param;mesh=Mesh(paramX)));
push!(models,Boussinesq(param;mesh=Mesh(paramX)));
push!(models,WhithamBoussinesq(param;mesh=Mesh(paramX)));
push!(models,SerreGreenNaghdi(param;mesh=Mesh(paramX)));
push!(models,WhithamGreenNaghdi(param;mesh=Mesh(paramX)));
push!(models,NonHydrostatic(param;mesh=Mesh(paramX)));
push!(models,SquareRootDepth(param;mesh=Mesh(paramX)));
push!(models,IsobeKakinuma(param;mesh=Mesh(paramX)));
push!(models,WWn(param;mesh=Mesh(paramX)));

for model in models
    @testset "Preserved quantities : $(model.label)" begin
    pb=Problem(model,Init(x->exp.(-x.^2),x->(x.+1).*exp.(-(x.+1).^2)),paramT)
    solve!(pb;verbose=false)
    @test isapprox(mass(pb),mass(pb;T=0),rtol=1e-10)
    @test isapprox(momentum(pb),momentum(pb;T=0),rtol=1e-10)
    @test isapprox(energy(pb),energy(pb;T=0),rtol=1e-8)
    @test abs(massdiff(pb))<1e-10
    @test abs(momentumdiff(pb))<1e-10
    @test abs(energydiff(pb;rel=true))<1e-8
    end
end


text:
using JLD

function sample_from_prior()
    # this just does rejection sampling to
    # stay in the identified region of the
    # MA(2) model. I'm too lazy to reflect
    # the bad draws
    ok = false
    theta1 = 0.
    theta2 = 0.
    while !ok
        theta1 = 4. * rand() - 2.
        theta2 = 2. * rand() - 1.
        ok = (theta2 > -1. + theta1) & (theta2 > -1. - theta1)
    end
    return [theta1 theta2]

end

function ma2(theta, n)
    e = randn(n+2)
    y = e[3:end] + theta[1]*e[2:end-1] + theta[2]*e[1:end-2]
end

function lags(x,p)
    cols = size(x,2)
    result = zeros(size(x,1), cols*p)
    for i = 1:p
        result[:,i*cols-cols+1:i*cols] = [ones(i,cols); x[1:end-i,:]]
    end
    return result
end

# fits an AR(P) model to MA(2) data
function aux_stat(theta, n, P)
    y = ma2(theta,n)
    x = lags(y,P)
    y = y[P+1:end]
    x = [ones(n-P,1) x[P+1:end,:]]
    Z = x\y
    return Z
end

n = 100 # sample size
R = 1000000 # replications
P = 10 # order of AR fit to the data
thetas = zeros(R,2)
Zs = zeros(R,P+1)
for i = 1:R
    if mod(i,10000)== 0.
        println(i)
    end
    theta = sample_from_prior()
    thetas[i,:] = theta'
    Zs[i,:] = aux_stat(theta, n, P)
end

save("MA2data.jld","thetas", thetas, "Zs", Zs)


text:
type TransformedGrid
  grid::Vector{Float64}
  transformedGrid::Vector{Float64}
  dxm::Vector{Float64}
  dxp::Vector{Float64}
  dx::Vector{Float64}
end

function TransformedGrid(grid::Vector{Float64}, f::Function)
  n = length(grid)
  transformedGrid = f(grid)
  dxm = zeros(n)
  dxp = zeros(n)
  dx = zeros(n)

  for i = 2:n - 1
    dxm[i] = transformedGrid[i] - transformedGrid[i-1]
    dxp[i] = transformedGrid[i+1] - transformedGrid[i]
    dx[i] = dxm[i] + dxp[i]
  end

  return TransformedGrid(grid, transformedGrid, dxm, dxp, dx)
end

LogGrid(grid::Vector{Float64}) = TransformedGrid(grid, log)


text:
const DATASET_CHECKSUM_V1_0_0 = Dict(
    "pudl-eia860-eia923" => "2ab91241d96985102cd8917c7d55c7943bd3d22364b99b8d8dfb365b2fe43cb4",
    # "pudl-eia860-eia923-epacems" => "TODO"
    "pudl-ferc1" => "0051afa03da53d7655b586641f1cebe46e3df97cb7239f8f467daa8fc2ea28fa",
)
const DATASET_DESCRIPTION_V1_0_0 = Dict(
    "pudl-eia860-eia923" => """
        Data originally collected and published by the [US Energy Information Administration](https://www.eia.gov/)
        (US EIA). The data from [EIA Form 860](https://www.eia.gov/electricity/data/eia860/) covers the years 2011-2018.
        The [Form 923 data](https://www.eia.gov/electricity/data/eia923/) covers 2009-2018.
        A large majority of the data published in the original data sources has been included, but some parts,
        like fuel stocks on hand, and EIA 923 schedules 6, 7, & 8 have not yet been integrated.
        """,
    "pudl-eia860-eia923-epacems" => """
        This data package contains all of the same data as the `pudl-eia860-eia923` package,
        as well as the Hourly Emissions data from the US Environmental Protection Agency's (EPA's)
        [Continuous Emissions Monitoring System](https://www.epa.gov/emc/emc-continuous-emission-monitoring-systems)
        (CEMS) from 1995-2018. The EPA CEMS data covers thousands of power plants at hourly resolution for decades,
        and contains close to a billion records.
        """,
    "pudl-ferc1" => """
        Seven data tables from [FERC Form 1](https://www.ferc.gov/docs-filing/forms/form-1/data.asp)
        are included, primarily relating to individual power plants, and covering the years 1994-2018
        (the entire span of time for which FERC provides this data).

        These tables are the only ones which have been subjected to any cleaning or organization for
        programmatic use within PUDL. The complete, raw FERC Form 1 database contains 116 different
        tables with many thousands of columns of mostly financial data.
        We will archive a complete copy of the multi-year FERC Form 1 Database as a file-based SQLite
        database at Zenodo, independent of this data release.
        It can also be re-generated using the `catalystcoop.pudl` Python package and the original
        source data files archived as part of this data release.
        """,
)
const VERSION_DESCRIPTION_V1_0_0 = """
    Dataset: PUDL Data Release v1.0.0
    Website: https://zenodo.org/record/3653159
    Author: Selvans, Zane A., Gosnell, Christina M.
    Date of Publication: February 7, 2020
    License: https://creativecommons.org/licenses/by/4.0/legalcode

    This is the first data release from the
    [Public Utility Data Liberation (PUDL) project](https://catalyst.coop/pudl).
    It can be referenced & cited using https://doi.org/10.5281/zenodo.3653159

    For more information about the free and open source software used to generate this data release,
    see [Catalyst Cooperative's PUDL repository on Github](https://github.com/catalyst-cooperative/pudl),
    and the associated [documentation on Read TheDocs](https://catalystcoop-pudl.readthedocs.io/en/v0.3.1/).
    This data release was generated using v0.3.1 of the `catalystcoop.pudl` python package.
    """


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   International Geomagnetic Field Model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
#   [2] https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f
#   [3] https://www.mathworks.com/matlabcentral/fileexchange/34388-international-geomagnetic-reference-field--igrf--model
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export igrf12
export igrf12syn

################################################################################
#                                  Functions
################################################################################

"""
    function igrf12(date::Number, r::Number, λ::Number, Ω::Number, T; show_warns = true)

**IGRF v12 Model**

Compute the geomagnetic field vector [nT] at the date `date` [Year A.D.] and
position (`r`, `λ`, `Ω`).

The position representation is defined by `T`. If `T` is `Val{:geocentric}`,
then the input must be **geocentric** coordinates:

1. Distance from the Earth center `r` [m];
1. Geocentric latitude `λ` (-π/2, +π/2) \\[rad]; and
2. Geocentric longitude `Ω` (-π, +π) \\[rad].

If `T` is `Val{:geodetic}`, then the input must be **geodetic** coordinates:

1. Altitude above the reference ellipsoid `h` (WGS-84) \\[m];
2. Geodetic latitude `λ` (-π/2, +π/2) \\[rad]; and
3. Geodetic longitude `Ω` (-π, +π) \\[rad].

If `T` is omitted, then it defaults to `Val{:geocentric}`.

Notice that the output vector will be represented in the same reference system
selected by the parameter `T` (geocentric or geodetic). The Y-axis of the output
reference system always points East. In case of **geocentric coordinates**, the
Z-axis points toward the center of Earth and the X-axis completes a right-handed
coordinate system. In case of **geodetic coordinates**, the X-axis is tangent to
the ellipsoid at the selected location and points toward North, whereas the
Z-axis completes a right-hand coordinate system.

# Keywords

* `show_warns`: Show warnings about the data (**Default** = `true`).

# Remarks

The `date` must be greater or equal to 1900 and less than or equal 2025. Notice
that a warning message is printed for dates greater than 2020.

# Disclaimer

This function is an independent implementation of the IGRF model. It contains a
more readable code than the original one in FORTRAN, because it uses features
available in Julia language.

"""
igrf12(date::Number, r::Number, λ::Number, Ω::Number; show_warns = true) =
    igrf12(date, r, λ, Ω, Val{:geocentric}; show_warns = show_warns)

function igrf12(date::Number,
                r::Number,
                λ::Number,
                Ω::Number,
                ::Type{Val{:geocentric}};
                show_warns::Bool = true)

    # Input verification
    # ==================

    # Check the data, since this model is valid for years between 1900 and 2025.
    ( (date < 1900) || (date > 2025) ) &&
    error("This IGRF version will not work for years outside the interval [1900, 2025).")

    # Check if the latitude and longitude are valid.
    ( (λ < -pi/2) || (λ > pi/2) ) &&
    error("The latitude must be between -π/2 and +π/2 rad.")

    ( (Ω < -pi) || (Ω > pi) ) &&
    error("The longitude must be between -π and +π rad.")

    # Warn the user that for dates after the year 2020 the accuracy maybe
    # reduced.
    show_warns && (date > 2020) &&
    @warn("The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.")

    # Input variables conversion
    # ==========================

    # Convert latitude / longitude to colatitude and east-longitude.
    θ = pi/2 - λ
    ϕ = (Ω >= 0) ? Ω : 2pi + Ω

    # The input variable `r` is in [m], but all the algorithm requires it to be
    # in [km].
    r /= 1000

    # Preliminary setup
    # =================

    # Compute the epoch that will be used to obtain the coefficients. This is
    # necessary because the IGRF provides coefficients every 5 years. Between
    # two epochs, those coefficients must be interpolated.
    idx   = (date < 2020) ? floor(Int, (date-1900)*0.2+1) : 24
    epoch = 1900 + (idx-1)*5

    # Compute the fraction of time from the epoch of the coefficient selected by
    # `idx`.
    Δt = date - epoch

    # Compute the maximum spherical harmonic degree for the selected date.
    n_max = (epoch < 1995) ? 10 : 13

    # Compute the Schmidt quasi-normalized associated Legendre functions and
    # their first order derivative, neglecting the phase term.
    dP, P = dlegendre(Val{:schmidt}, θ, n_max, n_max, false)

    # Parameters and auxiliary variables
    # ==================================

    # Auxiliary variables to select the IGRF coefficients.
    H = H_igrf12
    G = G_igrf12

    # Reference radius [km].
    a = 6371.2

    # Auxiliary variables to decrease the computational burden.
    sin_ϕ,  cos_ϕ  = sincos(1ϕ)
    sin_2ϕ, cos_2ϕ = sincos(2ϕ)
    ratio   = a/r
    fact    = ratio

    # Initialization of variables
    # ===========================

    dVr = 0.0   # Derivative of the Geomagnetic potential w.r.t. r.
    dVθ = 0.0   # Derivative of the Geomagnetic potential w.r.t. θ.
    dVϕ = 0.0   # Derivative of the Geomagnetic potential w.r.t. ϕ.
    ΔG  = 0.0   # Auxiliary variable to interpolate the G coefficients.
    ΔH  = 0.0   # Auxiliary variable to interpolate the H coefficients.
    kg  = 1     # Index to obtain the values of the matrix `G`.
    kh  = 1     # Index to obtain the values of the matrix `H`.

    # Geomagnetic potential
    # =====================

    @inbounds for n = 1:n_max
        aux_dVr = 0.0
        aux_dVθ = 0.0
        aux_dVϕ = 0.0

        # Compute the contributions when `m = 0`
        # ======================================

        # Get the coefficients in the epoch and interpolate to the desired
        # time.
        Gnm_e0 = G[kg,idx+2]

        if date < 2015
            Gnm_e1 = G[kg,idx+3]
            ΔG     = (Gnm_e1-Gnm_e0)/5
        else
            ΔG     = G[kg,27]
        end

        Gnm  = Gnm_e0 + ΔG*Δt
        kg  += 1

        aux_dVr += -(n+1)/r*Gnm*P[n+1,1]
        aux_dVθ += Gnm*dP[n+1,1]

        # Sine and cosine with m = 1
        # ==========================
        #
        # This values will be used to update recursively `sin(m*ϕ)` and
        # `cos(m*ϕ)`, reducing the computational burden.
        sin_mϕ   = +sin_ϕ    # sin( 1*λ_gc)
        sin_m_1ϕ = 0.0       # sin( 0*λ_gc)
        sin_m_2ϕ = -sin_ϕ    # sin(-1*λ_gc)
        cos_mϕ   = +cos_ϕ    # cos( 1*λ_gc)
        cos_m_1ϕ = 1.0       # cos( 0*λ_gc)
        cos_m_2ϕ = +cos_ϕ    # cos(-2*λ_gc)

        # Other auxiliary variables that depend only on `n`
        # =================================================

        fact_dVr = (n+1)/r

        # Compute the contributions when `m ∈ [1,n]`
        # ==========================================

        for m = 1:n
            # Compute recursively `sin(m*ϕ)` and `cos(m*ϕ)`.
            sin_mϕ = 2cos_ϕ*sin_m_1ϕ-sin_m_2ϕ
            cos_mϕ = 2cos_ϕ*cos_m_1ϕ-cos_m_2ϕ

            # Compute the coefficients `G_nm` and `H_nm`
            # ==========================================

            # Get the coefficients in the epoch and interpolate to the
            # desired time.
            Gnm_e0 = G[kg,idx+2]
            Hnm_e0 = H[kh,idx+2]

            if date < 2015
                Gnm_e1 = G[kg,idx+3]
                Hnm_e1 = H[kh,idx+3]
                ΔG     = (Gnm_e1-Gnm_e0)/5
                ΔH     = (Hnm_e1-Hnm_e0)/5
            else
                ΔG     = G[kg,27]
                ΔH     = H[kh,27]
            end

            Gnm    = Gnm_e0 + ΔG*Δt
            Hnm    = Hnm_e0 + ΔH*Δt
            kg    += 1
            kh    += 1

            GcHs_nm = Gnm*cos_mϕ + Hnm*sin_mϕ
            GsHc_nm = Gnm*sin_mϕ - Hnm*cos_mϕ

            # Compute the contributions for `m`
            # =================================

            aux_dVr += -fact_dVr*GcHs_nm*P[n+1,m+1]
            aux_dVθ += GcHs_nm*dP[n+1,m+1]
            aux_dVϕ += (θ == 0) ? -m*GsHc_nm*dP[n+1,m+1] : -m*GsHc_nm*P[n+1,m+1]

            # Update the values for the next step
            # ===================================

            sin_m_2ϕ = sin_m_1ϕ
            sin_m_1ϕ = sin_mϕ
            cos_m_2ϕ = cos_m_1ϕ
            cos_m_1ϕ = cos_mϕ
        end

        # Perform final computations related to the summation in `n`
        # ==========================================================

        # fact = (a/r)^(n+1)
        fact    *= ratio

        # aux_<> *= (a/r)^(n+1)
        aux_dVr *= fact
        aux_dVϕ *= fact
        aux_dVθ *= fact

        dVr += aux_dVr
        dVϕ += aux_dVϕ
        dVθ += aux_dVθ
    end

    dVr *= a
    dVϕ *= a
    dVθ *= a

    # Compute the Geomagnetic field vector in the geocentric reference frame
    # ======================================================================

    x = +1/r*dVθ
    y = (θ == 0) ? -1/r*dVϕ : -1/(r*sin(θ))*dVϕ
    z = dVr

    B_gc = SVector{3,Float64}(x,y,z)
end

function igrf12(date::Number,
                h::Number,
                λ::Number,
                Ω::Number,
                ::Type{Val{:geodetic}};
                show_warns = true)

    # TODO: This method has a small error (≈ 0.01 nT) compared with the
    # `igrf12syn`.  However, the result is exactly the same as the MATLAB
    # function in [3]. Hence, this does not seem to be an error in the
    # conversion from geodetic to geocentric coordinates. This is probably
    # caused by a numerical error. Further verification is necessary.

    # Convert the geodetic coordinates to geocentric coordinates.
    (λ_gc, r) = GeodetictoGeocentric(λ, h)

    # Compute the geomagnetic field in geocentric coordinates.
    B_gc = igrf12(date, r, λ_gc, Ω, Val{:geocentric}; show_warns = show_warns)

    # Convert to geodetic coordinates.
    D_gd_gc = create_rotation_matrix(λ_gc - λ,:Y)
    B_gd    = D_gd_gc*B_gc
end

"""
    function igrf12syn(isv::Int, date::Number, itype::Int, alt::Number, colat::Number, elong::Number; show_warns = true)

This is a Julia implementation of the official IGRF source code, which was
written in Fortran [2]. The input and output variables are exactly the same as
the ones described in the function `igrf12syn` in [2].

# Args

* `isv`: `0` if main-field values are required, `1` if secular variation values
         are required.
* `date`: Year A.D.
* `itype`: `1` if geodetic (spheroid), `2` if geocentric (sphere).
* `alt`: Height above sea level [km] if `itype = 1`, or distance from the center of
         Earth [km] if `itype = 2` (must be > 3485 km).
* `colat`: Colatitude (0 - 180) [˚].
* `elong`: East-Longitude (0 - 360) [˚].

# Keywords

* `show_warns`: Show warnings about the data (**Default** = `true`).

# Returns

* The north component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The east component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The vertical component [nT] if `isv = 0`, or [nT/year] if `isv = 1`.
* The total intensity if `isv = 0`, or rubbish if `isv = 1`.

# Remarks

* The `date` must be greater or equal to 1900 and less than or equal 2025.
Notice that a warning message is printed for dates grated than 2020.

"""
function igrf12syn(isv::Int,
                   date::Number,
                   itype::Int,
                   alt::Number,
                   colat::Number,
                   elong::Number;
                   show_warns = true)

    # Check the data, since this model is valid for years between 1900 and 2025.
    ( (date < 1900) || (date > 2025) ) &&
    error("This IGRF version will not work for years outside the interval [1900, 2025).")

    # Warn the user that for dates after the year 2020 the accuracy maybe reduced.
    show_warns && (date > 2020) &&
    @warn("The magnetic field computed with this IGRF version may be of reduced accuracy for years greater than 2020.")

    # Declaration of variables
    gh          = gh_igrf12
    fn::Int     = 0
    gn::Int     = 0
    kmx::Int    = 0
    ll::Int     = 0
    nc::Int     = 0
    nmx::Int    = 0
    x::Float64  = 0.0
    y::Float64  = 0.0
    z::Float64  = 0.0
    t::Float64  = 0.0
    tc::Float64 = 0.0
    cl          = zeros(13)
    sl          = zeros(13)
    p           = zeros(105)
    q           = zeros(105)

    if date < 2015
        t  = 0.2*(date - 1900)
        ll = floor(Int,t)
        t  = t - ll

        # SH models before 1995.0 are only to degree 10.
        if date < 1995
            nmx = 10

            # nc = nmx*(nmx+2)
            nc  = 120

            ll  = nc*ll

            # kmx = (nmx+1)*(nmx+2)/2
            kmx = 66
        else
            nmx = 13

            # nc  = nmx*(nmx+2)
            nc  = 195

            ll  = floor(Int,0.2*(date - 1995))

            # 19 is the number of SH models that extend to degree 10.
            ll  = 120*19 + nc*ll

            # kmx = (nmx+1)*(nmx+2)/2
            kmx = 105
        end

        tc = 1 - t

        if isv == 1
            t  = +0.2
            tc = -0.2
        end
    else
        t  = date - 2015
        tc = 1.0

        if isv == 1
            t  = 1.0
            tc = 0.0
        end

        ll  = 3060
        nmx = 13

        # nc  = nmx*(nmx+2)
        nc = 195

        # kmx = (nmx+1)*(nmx+2)/2
        kmx = 105
    end

    r::Float64  = alt
    ct          = cos(colat*pi/180)
    st          = sin(colat*pi/180)
    cl[1]       = cos(elong*pi/180)
    sl[1]       = sin(elong*pi/180)
    cd::Float64 = 1.0
    sd::Float64 = 0.0
    l::Int      = 1
    m::Int      = 1
    n::Int      = 0

    if itype != 2
        # Conversion between geodetic and geocentric coordinates using WGS84
        # spheroid.
        a2    = 40680631.6
        b2    = 40408296.0
        one   = a2*st^2
        two   = b2*ct^2
        three = one + two
        rho   = sqrt(three)
        r     = sqrt(alt*(alt + 2*rho) + (a2*one + b2*two)/three)
        cd    = (alt + rho)/r
        sd    = (a2 - b2)/rho*ct*st/r
        one   = ct
        ct    = ct*cd -  st*sd
        st    = st*cd + one*sd
    end

    ratio = 6371.2/r
    rr    = ratio^2

    # Computation of Schmidt quasi-normal coefficients p and x(=q).
    # =============================================================

    p[1] = 1.0
    p[3] = st
    q[1] = 0.0
    q[3] = ct

    for k = 2:kmx
        # There is no need to check bounds here. The code guarantees that
        # everything is inside the right bounds. This increased the source code
        # performance by 13%.
        @inbounds begin
            if n < m
                m  = 0
                n  = n+1
                rr = rr*ratio
                fn = n
                gn = n-1
            end

            fm = m

            if (m == n)
                if k != 3
                    one   = sqrt(1 - 0.5/fm)
                    j     = k - n - 1
                    p[k]  = one*st*p[j]
                    q[k]  = one*(st*q[j] + ct*p[j])
                    cl[m] = cl[m-1]*cl[1] - sl[m-1]*sl[1]
                    sl[m] = sl[m-1]*cl[1] + cl[m-1]*sl[1]
                end
            else
                gmm   = m^2
                one   = sqrt(fn^2 - gmm)
                two   = sqrt(gn^2 - gmm)/one
                three = (fn + gn)/one
                i     = k - n
                j     = i - n + 1
                p[k]  = three*ct*p[i] - two*p[j]
                q[k]  = three*(ct*q[i] - st*p[i]) - two*q[j]
            end

            # Synthesis of x, y, and z in geocentric coordinates.
            lm  = ll + l
            one = (tc*gh[lm] + t*gh[lm+nc])*rr

            if m != 0
                two   = (tc*gh[lm+1] + t*gh[lm+nc+1])*rr
                three = one*cl[m] + two*sl[m]
                x     = x + three*q[k]
                z     = z - (fn + 1)*three*p[k]

                if st != 0
                    y = y + (one*sl[m] - two*cl[m])*fm*p[k]/st
                else
                    y = y + (one*sl[m] - two*cl[m])*q[k]*ct
                end

                l = l + 2
            else
                x = x + one*q[k]
                z = z - (fn + 1)*one*p[k]
                l = l + 1
            end

            m = m + 1
        end
    end

    # Conversion to coordinate system specified by itype.
    # ===================================================
    one   = x
    x     = x*cd +   z*sd
    z     = z*cd - one*sd
    f     = sqrt(x^2 + y^2 + z^2)

    (x,y,z,f)
end


text:
import BasisFunctions: interpolation_grid
export interpolation_grid
@ap_interface interpolation_grid
@ap_sampling_property_to_platform interpolation_grid

interpolation_grid(ss::SamplingStyle, platform::Platform, plt_par, smpl_par; options...) =
    interpolation_grid(dictionary(platform, plt_par))

export oversampling_grid
@ap_interface oversampling_grid
@ap_sampling_property_to_platform oversampling_grid

oversampling_grid(ss::SamplingStyle, platform::Platform, plt_par, smpl_par; options...) =
    oversampling_grid(dictionary(platform, plt_par), smpl_par)
include("oversampling_grid.jl")

export platform_grid
@ap_interface platform_grid
@ap_sampling_property_to_platform platform_grid

# there is no default, signal a hopefully informative error
platform_grid(ss::SamplingStyle, platform::Platform, plt_par, smpl_par; grid = nothing, options...) =
    grid == nothing ? error("This platform does not implement a `platform_grid`.") : grid

export sampling_grid
@ap_interface sampling_grid
@ap_sampling_property sampling_grid

"Return the sampling grid that would be used for an approximation problem."
sampling_grid(ss::InterpolationStyle, ap::ApproximationProblem; options...) =
    interpolation_grid(ap; options...)
sampling_grid(ss::GridStyle, ap::ApproximationProblem; grid = nothing, options...) =
    grid == nothing ? platform_grid(ap; options...) : grid
sampling_grid(ss::OversamplingStyle, ap::ApproximationProblem; options...) =
    oversampling_grid(ap; options...)
sampling_grid(ss::WeightedDiscreteStyle, ap::ApproximationProblem; options...) =
    sampling_grid(unweighted(ss), ap; options...)

to_product(grids::AbstractGrid...) = productgrid(grids...)
# function sampling_grid(ss::ProductSamplingStyle, ap::ApproximationProblem; options...)
#     @assert nfactors(ss)==nfactors(platform(ap))==length(platformparameter(ap))==length(samplingparameter(ap))
#     productgrid( map((st,pl,p_par,s_par)->sampling_grid(pl, p_par, s_par; samplingstyle=st, options...),
#         factors(ss), factors(platform(ap)), productparameter(ap), samplingparameter(ap))...)
# end

# sampling_grid(ss::ProductSamplingStyle, ap::ApproximationProblem; options...) =
#     sampling_grid(ss, platform(ap), productparameter(ap), samplingparameter(ap); options...)
#
# function sampling_grid(ss::ProductSamplingStyle, platform::ProductPlatform, plt_par, smpl_par; options...)
#     @assert nfactors(ss)==nfactors(platform)==length(plt_par)==length(smpl_par)
#     productgrid(map(
#         (pl,pl_par,s_par, st) -> sampling_grid(pl, pl_par, s_par; samplingstyle=st, options...),
#         factors(platform), plt_par, smpl_par, factors(ss))...)
# end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random, OrdinaryDiffEq

a = 0.2
b = 0.2
c = 5.7
amu = 0.02
ω1 = 1. + amu
ω2 = 1. - amu
k1 = 0.01
k2 = k1

y_init = [0.0, 0.2, 0., 0.11, 0.19, 0.1]

croesslers = Systems.coupled_roessler(y_init; ω1, ω2, a, b, c, k1, k2)
diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9, maxiters = Inf)
Δt = 0.01
T = 1000
t = 0:Δt:T
L = length(t)÷10

ks = [0.03, 0.04, 0.05] # estimated from the Lyapunov exponents plot

bbox2 = Dict(:boxstyle => "round,pad=0.3", :facecolor=>"white", :alpha => 0.9)
axs_for_labels = []

fig = figure(figsize = (figx, 1.5figy))
for (i, k) in enumerate(ks)
    M = length(ks)
    @show k
    set_parameter!(croesslers, 6, k)
    set_parameter!(croesslers, 7, k)
    X = trajectory(croesslers, T; Ttr = 5000, Δt, diffeq...)

    row = (i-1)*3
    ax1 = fig.add_subplot(M, 3, (1, 2) .+ row)
    ax2 = fig.add_subplot(M, 3, 3+row)

    ax1.plot(t[1:L], X[1:L, 1])
    ax1.plot(t[1:L], X[1:L, 4]; lw = 3.0, color = COLORS[2])
    ax2.plot(X[:, 1], X[:, 4]; lw = 1.0, color = COLORS[1])
    push!(axs_for_labels, ax1)

    λs = lyapunovspectrum(croesslers, 50000, 4; Ttr = 1000, diffeq...)
    # ax1.text(0.02, 0.9, "\$k = $(k)\$\n\$ \\lambda_1 = $(round(λs[1];digits=2))\$\n\$\\lambda_2 = $(round(λs[2];digits=2))\$",
    #     transform = ax1.transAxes, va = :top, bbox = bbox2, fontsize = 28
    # )
    ax1.text(0.02, 0.5, "\$k = $(k),\\, \\lambda_1 = $(round(λs[1];digits=2))\\,\\lambda_2 = $(round(λs[2];digits=2))\$",
        transform = ax1.transAxes, va = :center, bbox = bbox2, fontsize = 28
    )
    @show λs[1:4]
    ax1.set_xlim(0, t[1:L][end])
    ax1.set_ylim(-12,12)
    ax2.set_xlim(-12,12)
    ax2.set_ylim(-12,12)
    ax2.set_yticklabels([])
    ax1.set_ylabel("\$x_1, y_1\$",labelpad = -20)
    ax2.set_ylabel("\$y_1\$")
    if i < M
        ax1.set_xticklabels([])
        ax2.set_xticklabels([])
    else
        ax1.set_xlabel("\$t\$"; labelpad = -20)
        ax2.set_xticks([-10, 0, 10])
        ax2.set_xticklabels(["-10", "", "10"])
        ax2.set_xlabel("\$x_1\$", labelpad = -20)
    end
end
add_identifiers!(fig, axs_for_labels)
fig.tight_layout(;pad=0.33)
wsave(plotsdir("9", "coupled_roessler"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module AnnealingAPI

using CorrelationFunctions: Directional

"""
Abstract type for a rollback token returned by `update_corrfns!`.
"""
abstract type AbstractRollbackToken end

"""
    SimpleRollbackToken(value, index)

Create a simple rollback token which just holds an index and an old
value at that index.
"""
struct SimpleRollbackToken{T, N} <: AbstractRollbackToken
    value :: T
    index :: CartesianIndex{N}
end

"""
    update_corrfns!(tracker, value, index)

This function is equivalent to writing `tracker[index] = value` with
exception that it also returns a rollback handle which can fastly
bring the tracker to the previous state by calling `rollback!`.

See also: [`rollback!`](@ref).
"""
function update_corrfns! end

function update_corrfns!(array :: AbstractArray{T, N},
                         val,
                         index :: CartesianIndex{N}) where {T, N}
    token = SimpleRollbackToken(array[index], index)
    array[index] = val
    return token
end

"""
    rollback!(tracker, token)

Bring the system to the state before `update_corrfns!` was
called. `token` must be an object returned by `update_corrfns!`.

See also: [`update_corrfns!`](@ref).
"""
function rollback! end

rollback!(array :: AbstractArray, rollback :: AbstractRollbackToken) =
    array[rollback.index] = 1 - array[rollback.index]

rollback!(array    :: AbstractArray{T, N},
          rollback :: SimpleRollbackToken{T, N}) where {T, N} =
    array[rollback.index] = rollback.value

"""
    Generic type for correlation function tracker.

See also: [`L2Tracker`](@ref), [`S2Tracker`](@ref),
[`SSTracker`](@ref).
"""
abstract type AbstractTracker{T} end

"""
    L2Tracker(phase)

Descriptor for line-segment correlation function for the phase
`phase`.

See also: [`S2Tracker`](@ref), [`SSTracker`](@ref),
[`AbstractTracker`](@ref).
"""
struct L2Tracker{T} <: AbstractTracker{T}
    phase :: T
end

"""
    S2Tracker(phase)

Descriptor for two-point correlation function for the phase
`phase`.

See also: [`L2Tracker`](@ref), [`SSTracker`](@ref),
[`AbstractTracker`](@ref).
"""
struct S2Tracker{T} <: AbstractTracker{T}
    phase :: T
end

"""
    SSTracker(phase)

Descriptor for surface-surface correlation function for the phase
`phase`.

See also: [`S2Tracker`](@ref), [`L2Tracker`](@ref),
[`AbstractTracker`](@ref).
"""
struct SSTracker{T} <: AbstractTracker{T}
    phase :: T
end

"""
    SVTracker(phase)

Descriptor for surface-void correlation function for the phase
`phase`.

See also: [`S2Tracker`](@ref), [`L2Tracker`](@ref),
[`SSTracker`](@ref), [`AbstractTracker`](@ref).
"""
struct SVTracker{T} <: AbstractTracker{T}
    phase :: T
end

(tracked :: L2Tracker{T})(system :: AbstractArray{T}; kwargs...) where T =
    Directional.l2(system, tracked.phase; kwargs...)
(tracked :: S2Tracker{T})(system :: AbstractArray{T}; kwargs...) where T =
    Directional.s2(system, tracked.phase; kwargs...)
(tracked :: SSTracker{T})(system :: AbstractArray{T}; kwargs...) where T =
    Directional.surfsurf(system, tracked.phase; kwargs...)
(tracked :: SVTracker{T})(system :: AbstractArray{T}; kwargs...) where T =
    Directional.surfvoid(system, tracked.phase; kwargs...)

"""
    tracked_data(x :: AbstractArray)

Return an iterator over correlation function descriptors which are
tracked by the tracker `x`.
"""
function tracked_data end

"""
    tracked_length(x :: AbstractArray)

Return maximal tracked correlation length
"""
function tracked_length end

"""
    tracked_directions(x :: AbstractArray)

Return directions along which correlation functions are tracked.
"""
function tracked_directions end

export
    AbstractRollbackToken,
    update_corrfns!,
    rollback!,

    AbstractTracker,
    S2Tracker,
    L2Tracker,
    SSTracker,
    SVTracker,

    tracked_data,
    tracked_length,
    tracked_directions

end # module


text:
const WEBIO_CORE_PACKAGE_PATH = normpath(joinpath(
    @__DIR__, "..", "packages", "webio",
))
const JUPYTER_LAB_PROVIDER_PATH = normpath(joinpath(
    @__DIR__, "..", "packages", "jupyter-lab-provider",
))
const JUPYTER_NOTEBOOK_PROVIDER_PATH = normpath(joinpath(
    @__DIR__, "..", "packages", "jupyter-notebook-provider",
))

const CONFIG_BEGIN_MARKER = "###JULIA-WEBIO-CONFIG-BEGIN"
const CONFIG_END_MARKER = "###JULIA-WEBIO-CONFIG-END"

"""
    install_notebook_config()

Install necessary configuration for the `jlstaticserve` notebook (server)
extension. This function only configures the notebook extension, not the browser
nbextension or Jupyter Lab extension.
* Adds the path to `./deps` to Python's `sys.path` so that we can load the
    `jlstaticserve.py` extension. This is done in `jupyter_notebook_config.py`
    because there's no way to add to `sys.path` from the JSON config file.
* Adds `jlstaticserve` to the list of extensions loaded in the notebook server.
    This is done in `jupyter_notebook_config.json` because that file has higher
    precedence when both the `.py` and `.json` files exist (IPyWidgets, for
    example, writes to the JSON file, so if we only wrote to the `.py` file,
    that directive would take precedence and the `jlstaticserve` extension would
    **not** be loaded).
"""
function install_jupyter_serverextension()
    config_dir = jupyter_config_dir()
    mkpath(config_dir)
    config_file_py = joinpath(config_dir, "jupyter_notebook_config.py")
    config_py = isfile(config_file_py) ? read(config_file_py, String) : ""

    # Remove previous config
    config_regex = Regex(
        "\n?" * CONFIG_BEGIN_MARKER * ".*" * CONFIG_END_MARKER * "\n?",
        "s",
    )
    config_py = replace(config_py, config_regex => "")

    # We do a repr to make sure that all the necessary characters are escaped
    # and the whole path is wrapped in quotes.
    deps_dir = dirname(@__FILE__)
    deps_dir_esc = repr(deps_dir)
    error_msg = strip("""
        Directory $deps_dir could not be found; WebIO will not work as expected.
        """)
    config_py *= """

    $CONFIG_BEGIN_MARKER
    # Add the path to WebIO/deps so that we can load the `jlstaticserve` extension.
    import os, sys, warnings
    webio_deps_dir = $deps_dir_esc
    if os.path.isfile(os.path.join(webio_deps_dir, "jlstaticserve.py")):
        sys.path.append(webio_deps_dir)
    else:
        warning_msg = (
            'Directory %s could not be found; WebIO.jl will not work as expected. '
            + 'Make sure WebIO.jl is installed correctly (try running '
            + 'Pkg.add("WebIO") and Pkg.build("WebIO") from the Julia console to '
            + 'make sure it is).'
        ) % webio_deps_dir
        warnings.warn(warning_msg)
    $CONFIG_END_MARKER
    """

    config_file_json = joinpath(config_dir, "jupyter_notebook_config.json")
    config_json = isfile(config_file_json) ? read(config_file_json, String) : "{}"
    config_json = try
        JSON.parse(config_json)
    catch exc
        @error "Unable to parse Jupyter notebook config file ($config_file_json). Please fix it and rebuild WebIO." exception=exc
        rethrow()
    end

    # Magic to safely access nested JSON objects and set them to empty objects
    # if they don't exist yet.
    app_config = get!(config_json, "NotebookApp", Dict())
    extensions_config = get!(app_config, "nbserver_extensions", Dict())
    extensions_config["jlstaticserve"] = true

    # Defer writing until both files until the end to avoid inconsistent state.
    write(config_file_py, config_py)
    write(config_file_json, json(config_json, 4))
    return nothing
end

"""
    find_jupyter_cmd()

Find the most likely candidate for the `jupyter` executable.
This will locate `jupyter` by searching the `PATH` environment variable and,
if not found, tries to return Conda.jl's jupyter.
If both of these approaches fail, an error is thrown.
"""
function find_jupyter_cmd()::Cmd
    jupyter = Sys.which("jupyter")
    if jupyter !== nothing
        return `$jupyter`
    end
    # Try to find Conda.jl's Jupyter.
    # This is a heuristic - it might be different. I'm not sure that there's
    # an easy way to load Conda.jl (I tried adding it to extras but no dice).
    conda_root = joinpath(first(Base.DEPOT_PATH), "conda", "3")
    jupyter = joinpath(conda_root, "bin", "jupyter")
    if isfile(jupyter)
        return `$jupyter`
    end
    @error(
        "Unable to find a jupyter executable after searching PATH and Conda.jl.",
        PATH=ENV["PATH"], conda_root,
    )
    error("Could not find the Jupyter executable.")
end

"""
    install_jupyter_labextension([jupyter])

Install the Jupyter Lab extension for WebIO using the specified `jupyter`
executable.
The executable defaults to the first one found in the `PATH` or Jupyter
installed via Conda.jl.

The [IJulia provider documentation](https://juliagizmos.github.io/WebIO.jl/latest/providers/ijulia/)
provides some more information (and caveats) about the relationship between
Jupyter Lab and WebIO.
"""
function install_jupyter_labextension(jupyter::Cmd=find_jupyter_cmd())
    install_jupyter_serverextension()
    run(`$jupyter labextension link --no-build $WEBIO_CORE_PACKAGE_PATH`)
    run(`$jupyter labextension install --no-build $JUPYTER_LAB_PROVIDER_PATH`)
    run(`$jupyter lab build`)
end

"""
    install_webio_nbextension()

Install the Jupyter Notebook extension (nbextension) for WebIO.
This copies the nbextension code to the appropriate place and wr9tes the
appropriate configuration files.
"""
function install_jupyter_nbextension()
    install_jupyter_serverextension()
    extensions_dir = jupyter_nbextensions_dir()
    mkpath(extensions_dir)
    extension_dir = joinpath(extensions_dir, "webio")

    # Copy the nbextension files.
    @info "Copying WebIO nbextension files to $(extension_dir)."
    cp(
        joinpath(@__DIR__, "../packages/jupyter-notebook-provider/dist"),
        extension_dir,
        ; force=true
    )

    # Enable the notebook extension.
    config_dir = jupyter_nbconfig_dir()
    mkpath(config_dir)
    config_file = joinpath(config_dir, "notebook.json")
    config_data = try
        isfile(config_file) ? JSON.parse(read(config_file, String)) : Dict()
    catch exc
        @error(
            "Error parsing Jupyter config file $config_file - fix it and build "
                * "again or delete it to enable WebIO.",
            exception=exc,
        )
        error("Unable to parse Jupyter config file.")
    end
    config_data["load_extensions"] = get(config_data, "load_extensions", Dict())
    config_data["load_extensions"]["webio/main"] = true
    open(config_file, "w") do io
        JSON.print(JSON.Writer.PrettyContext(io, 4), config_data)
    end
end

### BEGIN BORROWED CODE ###
# This code is adapted from IJulia.jl (MIT license).
# https://github.com/JuliaLang/IJulia.jl/blob/68748c2b88b4b394e2802e911eb154b7008ebdd2/deps/kspec.jl
@static if Sys.iswindows()
    function appdata() # return %APPDATA%
        path = zeros(UInt16, 300)
        CSIDL_APPDATA = 0x001a
        result = ccall((:SHGetFolderPathW,:shell32), stdcall, Cint,
            (Ptr{Cvoid},Cint,Ptr{Cvoid},Cint,Ptr{UInt16}),C_NULL,CSIDL_APPDATA,C_NULL,0,path)
        return result == 0 ? transcode(String, resize!(path, findfirst(iszero, path)-1)) : get(ENV, "APPDATA", "")
    end
    function default_jupyter_data_dir()
        APPDATA = appdata()
        return !isempty(APPDATA) ? joinpath(APPDATA, "jupyter") : joinpath(get(ENV, "JUPYTER_CONFIG_DIR", joinpath(homedir(), ".jupyter")), "data")
    end
elseif Sys.isapple()
    default_jupyter_data_dir() = joinpath(homedir(), "Library/Jupyter")
else
    function default_jupyter_data_dir()
        xdg_data_home = get(ENV, "XDG_DATA_HOME", "")
        data_home = !isempty(xdg_data_home) ? xdg_data_home : joinpath(homedir(), ".local", "share")
        joinpath(data_home, "jupyter")
    end
end

function jupyter_data_dir()
    env_data_dir = get(ENV, "JUPYTER_DATA_DIR", "")
    !isempty(env_data_dir) ? env_data_dir : default_jupyter_data_dir()
end
### END BORROWED CODE ###

function jupyter_config_dir()
    env_config_dir = get(ENV, "JUPYTER_CONFIG_DIR", "")
    !isempty(env_config_dir) ? env_config_dir : joinpath(homedir(), ".jupyter")
end
jupyter_nbextensions_dir() = joinpath(jupyter_data_dir(), "nbextensions")
jupyter_nbconfig_dir() = joinpath(jupyter_config_dir(), "nbconfig")


text:
using Pkg
Pkg.instantiate()

max_stars_repo_name: JonPallbo/Laby

max_issues_repo_name: JonPallbo/Laby

max_forks_repo_name: JonPallbo/Laby


text:
"""
Created on Fri May 22 2020

@author: Yoann Pradat

Functions for scaling matrices in NMF iterations.
"""

function scale_col_W!(W, H, p=1)
    #### col norm
    norms_W = mapslices(x -> norm(x, p), W, dims=1)

    #### one or more columns of W may be null vectors
    indices = (1:size(W,2))[vec(norms_W .> 0)]

    W[:, indices] .= W[:, indices] .* repeat(norms_W[:, indices] .^ -1, size(W, 1), 1)
    H[indices, :] .= H[indices, :] .* repeat(transpose(norms_W[:, indices]), 1, size(H, 2))

    nothing
end

function scale_row_H!(W, H, p=1)
    #### col norm
    norms_H = mapslices(x -> norm(x, p), H, dims=2)

    #### one or more columns of W may be null vectors
    indices = (1:size(H,1))[vec(norms_H .> 0)]

    W[:, indices] .= W[:, indices] .* repeat(transpose(norms_H[indices, :]), size(W, 1), 1)
    H[indices, :] .= H[indices, :] .* repeat(norms_H[indices, :] .^ -1, 1, size(H, 2))

    nothing
end


text:
# Sort for dicts
import Base: sort, sort!

function sort!(d::OrderedDict; byvalue::Bool=false, args...)
    if d.ndel > 0
        rehash!(d)
    end

    if byvalue
        p = sortperm(d.vals; args...)
    else
        p = sortperm(d.keys; args...)
    end
    d.keys = d.keys[p]
    d.vals = d.vals[p]
    rehash!(d)
    return d
end

sort(d::OrderedDict; args...) = sort!(copy(d); args...)
sort(d::Dict; args...) = sort!(OrderedDict(d); args...)


text:
module ProxSDP

import Arpack
import KrylovKit
import MathOptInterface
import TimerOutputs
import TimerOutputs: @timeit

import Printf
import SparseArrays
import LinearAlgebra
import LinearAlgebra: BlasInt
import Random
import Logging

include("structs.jl")
include("options.jl")
include("util.jl")
include("printing.jl")
include("scaling.jl")
include("equilibration.jl")
include("pdhg.jl")
include("residuals.jl")
include("eigsolver.jl")
include("prox_operators.jl")

include("MOI_wrapper.jl")

end


text:
# Author: Yoshimitsu Okada (https://www.researchgate.net/profile/Yoshimitsu_Okada)
# Translated by Pengcheng Shi (shipengcheng1230@gmail.com) 06/2018

# An example wrapper for DC3D in julia as below:
# function dc3d_fortran(x::T, y::T, z::T, α::T, dep::T, dip::T, al1::T, al2::T, aw1::T, aw2::T,
#     disl1::T, disl2::T, disl3::T) where {T <: AbstractFloat}
#
#     # initial return values
#     # `RefValue{T}` may be also viable other than `Array{T, 1}`
#     ux = Array{Float64}(1)
#     uy = Array{Float64}(1)
#     uz = Array{Float64}(1)
#     uxx = Array{Float64}(1)
#     uyx = Array{Float64}(1)
#     uzx = Array{Float64}(1)
#     uxy = Array{Float64}(1)
#     uyy = Array{Float64}(1)
#     uzy = Array{Float64}(1)
#     uxz = Array{Float64}(1)
#     uyz = Array{Float64}(1)
#     uzz = Array{Float64}(1)
#     iret = Array{Int64}(1)
#
#     # call okada's code which is renamed as "__dc3d__" (see binding rename shown below)
#     # input args tuple must be syntactically written instead of a variable assigned
#     # macros could be used to simplify this in the future
#     ccall((:__dc3d__, "dc3d.so"), Void,
#         (
#             Ref{Float64},
#             Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
#             Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
#             Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64},
#             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
#             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
#             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
#             Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}}, Ptr{Array{Float64,1}},
#             Ref{Int64},
#         ),
#         α, x, y, z, dep, dip, al1, al2, aw1, aw2, disl1, disl2, disl3,
#         ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz,
#         iret,
#     )
#
#     # results valid iff iret[1] == 0
#     return (
#         iret[1],
#         ux[1], uy[1], uz[1],
#         uxx[1], uyx[1], uzx[1],
#         uxy[1], uyy[1], uzy[1],
#         uxz[1], uyz[1], uzz[1]
#     )
# end
# ```
#
#
# The corresponding fortran module is:
# ```fortran
# MODULE okada
#   USE, INTRINSIC :: iso_c_binding
#   IMPLICIT NONE
# CONTAINS
#
#   SUBROUTINE dc3d_wrapper(&
#        & alpha, &
#        & x, y, z, &
#        & depth, dip, &
#        & al1, al2, &
#        & aw1, aw2, &
#        & disl1, disl2, disl3, &
#        & ux, uy, uz, &
#        & uxx, uyx, uzx, &
#        & uxy, uyy, uzy, &
#        & uxz, uyz, uzz, &
#        & iret) BIND(C, NAME='__dc3d__')
#
#     REAL*8 :: &
#          & alpha, &
#          & x, y, z, &
#          & depth, dip, &
#          & al1, al2, &
#          & aw1, aw2, &
#          & disl1, disl2, disl3, &
#          & ux, uy, uz, &
#          & uxx, uyx, uzx, &
#          & uxy, uyy, uzy, &
#          & uxz, uyz, uzz
#
#     INTEGER*8 :: iret
#
#     CALL dc3d(&
#          & alpha, &
#          & x, y, z, &
#          & depth, dip, &
#          & al1, al2, &
#          & aw1, aw2, &
#          & disl1, disl2, disl3, &
#          & ux, uy, uz, &
#          & uxx, uyx, uzx, &
#          & uxy, uyy, uzy, &
#          & uxz, uyz, uzz, &
#          & iret)
#
#   END SUBROUTINE dc3d_wrapper
#
# END MODULE okada
# ```
#
#
# A sample of makefile is as below:
# ```make
# # Build Okada's code for calculating deformation due to a fault model
# #
# CC = gfortran
# CFLAGS = -fPIC -w -O3
# LDFLAGS = -shared
#
# SRCS = dc3d.f okada.f90
# OBJS = \$(SRCS:.c=.o)
#
# TARGET = dc3d.so
#
# \$(TARGET): \$(OBJS)
#     \$(CC) \$(CFLAGS) \$(LDFLAGS) -o \$(TARGET) \$(OBJS)
#

export dc3d

dc3d_cache(T) = ntuple(x -> zeros(T, 12), Val(10))

"""
    dc3d(x::T, y::T, z::T, α::T, dep::T, dip::T,
        al::Union{A, SubArray}, aw::Union{A, SubArray}, disl::A
        ) where {T <: Number, A <: AbstractVecOrMat{T}}

Calculate displacements and gradient of displacements due to a rectangular dislocation
    in an elastic isotropic halfspace.

Please see [dc3d](http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html) for details.
    Also this fault coordinate system is widely used in this package.

## Arguments
- `x`, `y`, `z`: observational position, where ``z ≤ 0``
- `α`: elastic constant
- `dep`: depth of fault origin
- `dip`: dip angle in degree
- `al`: a vector of 2 numbers, indicating along strike (x-axis) spanning
- `aw`: a vector of 2 numbers, indicating along downdip (y-z plane) spanning
- `disl`: a vector of 3 numbers, indicating dislocation in along-strike,
    along-downdip and tensile respectively.

## Output
A vector of 12 numbers, each is ``u_{x}``, ``u_{y}``, ``u_{z}``, ``u_{x,x}``,
    ``u_{y,x}``, ``u_{z,x}``, ``u_{x,y}``, ``u_{y,y}``, ``u_{z,y}``, ``u_{x,z}``
    ``u_{y,z}``, ``u_{z,z}``.
"""
function dc3d(x::T, y::T, z::T, α::T, dep::T, dip::T, al::A, aw::A, disl::AbstractVector{<:Real},
    cache::U=dc3d_cache(T)) where {T<:Real, A<:AbstractVector, U}
    dc3d(x, y, z, α, dep, dip, al[1], al[2], aw[1], aw[2], disl[1], disl[2], disl[3], cache)
end

function dc3d(x::T, y::T, z::T, α::T, dep::T, dip::T, al1::T, al2::T, aw1::T, aw2::T, disl1::T, disl2::T, disl3::T,
    cache::U=dc3d_cache(T)) where {T<:Real, U}

    u, du, dua, dub, duc, zerosol, xi, et, kxi, ket = cache
    fill!(u, zero(T))
    fill!(kxi, zero(T))
    fill!(ket, zero(T))

    z > zero(T) && return zerosol

    xi[1] = x - al1
    xi[2] = x - al2
    d = dep + z
    sc1 = shared_constants_1(α, dip)
    sd, cd = sc1[6], sc1[7]
    p = y * cd + d * sd
    q = y * sd - d * cd
    et[1] = p - aw1
    et[2] = p - aw2

    if q ≈ zero(T) && ((xi[1] * xi[2] ≤ zero(T) && et[1] * et[2] ≈ zero(T)) ||	(et[1] * et[2] ≤ zero(T) && xi[1] * xi[2] ≈ zero(T)))
        return zerosol
    end

    r12 = hypot(xi[1], et[2], q)
    r21 = hypot(xi[2], et[1], q)
    r22 = hypot(xi[2], et[2], q)

    (xi[1] ≤ zero(T) && (r21 + xi[2]) ≈ zero(T)) && (kxi[1] = one(T))
    (xi[1] ≤ zero(T) && (r22 + xi[2]) ≈ zero(T)) && (kxi[2] = one(T))
    (et[1] ≤ zero(T) && (r12 + et[2]) ≈ zero(T)) && (ket[1] = one(T))
    (et[1] ≤ zero(T) && (r22 + et[2]) ≈ zero(T)) && (ket[2] = one(T))

    for k in 1:2, j in 1:2
        sc2 = shared_constants_2(xi[j], et[k], q, sd, cd, kxi[k], ket[j])
        ua(dua, sc1, sc2, xi[j], et[k], q, disl1, disl2, disl3)
        for i = 1:3:10
            du[i] = -dua[i]
            du[i + 1] = -dua[i + 1] * cd + dua[i + 2] * sd
            du[i + 2] = -dua[i + 1] * sd - dua[i + 2] * cd
            i < 10 && continue
            du[i] = -du[i]
            du[i + 1] = -du[i + 1]
            du[i + 2] = -du[i + 2]
        end
        for i = 1:12
            if j + k ≠ 3  u[i] += du[i] end
            if j + k == 3  u[i] -= du[i] end
        end
    end

    d = dep - z
    p = y * cd + d * sd
    q = y * sd - d * cd
    et[1] = p - aw1
    et[2] = p - aw2

    if q ≈ zero(T) && ((xi[1] * xi[2] ≤ zero(T) && et[1] * et[2] ≈ zero(T)) ||	(et[1] * et[2] ≤ zero(T) && xi[1] * xi[2] ≈ zero(T)))
        return zerosol
    end

    fill!(kxi, zero(T))
    fill!(ket, zero(T))
    r12 = hypot(xi[1], et[2], q)
    r21 = hypot(xi[2], et[1], q)
    r22 = hypot(xi[2], et[2], q)

    (xi[1] ≤ zero(T) && (r21 + xi[2]) ≈ zero(T)) && (kxi[1] = one(T))
    (xi[1] ≤ zero(T) && (r22 + xi[2]) ≈ zero(T)) && (kxi[2] = one(T))
    (et[1] ≤ zero(T) && (r12 + et[2]) ≈ zero(T)) && (ket[1] = one(T))
    (et[1] ≤ zero(T) && (r22 + et[2]) ≈ zero(T)) && (ket[2] = one(T))

    for k in 1:2, j in 1:2
        sc2 = shared_constants_2(xi[j], et[k], q, sd, cd, kxi[k], ket[j])
        ua(dua, sc1, sc2, xi[j], et[k], q, disl1, disl2, disl3)
        ub(dub, sc1, sc2, xi[j], et[k], q, disl1, disl2, disl3)
        uc(duc, sc1, sc2, xi[j], et[k], q, z, disl1, disl2, disl3)
        for i = 1:3:10
            du[i] = dua[i] + dub[i] + z * duc[i]
            du[i + 1] = (dua[i + 1] + dub[i + 1] + z * duc[i + 1]) * cd - (dua[i + 2] + dub[i + 2] + z * duc[i + 2]) * sd
            du[i + 2] = (dua[i + 1] + dub[i + 1] - z * duc[i + 1]) * sd + (dua[i + 2] + dub[i + 2] - z * duc[i + 2]) * cd
            i < 10 && continue
            du[10] += duc[1]
            du[11] += duc[2] * cd - duc[3] * sd
            du[12] -= duc[2] * sd + duc[3] * cd
        end
        for i = 1:12
            j + k ≠ 3 && (u[i] += du[i])
            j + k == 3 && (u[i] -= du[i])
        end
    end
    return u
end

@inline function ua(u::A, sc1::B1, sc2::B2, xi::T, et::T, q::T, disl1::T, disl2::T, disl3::T,
    ) where {T <: Number, A <: AbstractArray{T}, B1 <: NTuple{12,T}, B2 <: NTuple{24,T}}
    alp1, alp2, alp3, alp4, alp5, sd, cd, sdsd, cdcd, sdcd, s2d, sc2d = sc1
    xi2, et2, q2, r, r2, r3, r5, y, d, tt, alx, ale, x11, y11, x32, y32, ey, ez, fy, fz, gy, gz, hy, hz = sc2

    fill!(u, zero(T))
    xy = xi * y11
    qx = q * x11
    qy = q * y11

    coeff = disl1 / 2π
    u[1]  += coeff * (tt / 2 + alp2 * xi * qy)
    u[2]  += coeff * (alp2 * q / r)
    u[3]  += coeff * (alp1 * ale - alp2 * q * qy)
    u[4]  += coeff * (-alp1 * qy - alp2 * xi2 * q * y32)
    u[5]  += coeff * (-alp2 * xi * q / r3)
    u[6]  += coeff * (alp1 * xy + alp2 * xi * q2 * y32)
    u[7]  += coeff * (alp1 * xy * sd + alp2 * xi * fy + d / 2 * x11)
    u[8]  += coeff * (alp2 * ey)
    u[9]  += coeff * (alp1 * (cd / r + qy * sd) - alp2 * q * fy)
    u[10] += coeff * (alp1 * xy * cd + alp2 * xi * fz + y / 2 * x11)
    u[11] += coeff * (alp2 * ez)
    u[12] += coeff * (-alp1 * (sd / r - qy * cd) - alp2 * q * fz)

    coeff = disl2 / 2π
    u[1]  += coeff * (alp2 * q / r)
    u[2]  += coeff * (tt / 2 + alp2 * et * qx)
    u[3]  += coeff * (alp1 * alx - alp2 * q * qx)
    u[4]  += coeff * (-alp2 * xi * q / r3)
    u[5]  += coeff * (-qy / 2 - alp2 * et * q / r3)
    u[6]  += coeff * (alp1 / r + alp2 * q2 / r3)
    u[7]  += coeff * (alp2 * ey)
    u[8]  += coeff * (alp1 * d * x11 + xy / 2 * sd + alp2 * et * gy)
    u[9]  += coeff * (alp1 * y * x11 - alp2 * q * gy)
    u[10] += coeff * (alp2 * ez)
    u[11] += coeff * (alp1 * y * x11 + xy / 2 * cd + alp2 * et * gz)
    u[12] += coeff * (-alp1 * d * x11 - alp2 * q * gz)

    coeff = disl3 / 2π
    u[1]  += coeff * (-alp1 * ale - alp2 * q * qy)
    u[2]  += coeff * (-alp1 * alx - alp2 * q * qx)
    u[3]  += coeff * (tt / 2 - alp2 * (et * qx + xi * qy))
    u[4]  += coeff * (-alp1 * xy + alp2 * xi * q2 * y32)
    u[5]  += coeff * (-alp1 / r + alp2 * q2 / r3)
    u[6]  += coeff * (-alp1 * qy - alp2 * q * q2 * y32)
    u[7]  += coeff * (-alp1 * (cd / r + qy * sd) - alp2 * q * fy)
    u[8]  += coeff * (-alp1 * y * x11 - alp2 * q * gy)
    u[9]  += coeff * (alp1 * (d * x11 + xy * sd) + alp2 * q * hy)
    u[10] += coeff * (alp1 * (sd / r - qy * cd) - alp2 * q * fz)
    u[11] += coeff * (alp1 * d * x11 - alp2 * q * gz)
    u[12] += coeff * (alp1 * (y * x11 + xy * cd) + alp2 * q * hz)
end

@inline function ub(u::A, sc1::B1, sc2::B2, xi::T, et::T, q::T, disl1::T, disl2::T, disl3::T,
    ) where {T <: Number,A <: AbstractArray{T},B1 <: NTuple{12,T},B2 <: NTuple{24,T}}
    alp1, alp2, alp3, alp4, alp5, sd, cd, sdsd, cdcd, sdcd, s2d, sc2d = sc1
    xi2, et2, q2, r, r2, r3, r5, y, d, tt, alx, ale, x11, y11, x32, y32, ey, ez, fy, fz, gy, gz, hy, hz = sc2

    fill!(u, zero(T))
    rd = r + d
    d11 = one(T) / (r * rd)
    aj2 = xi * y / rd * d11
    aj5 = -(d + y * y / rd) * d11
    if (cd ≉ zero(T))
        if (xi ≈ zero(T))
            ai4 = zero(T)
        else
            x = sqrt(xi2 + q2)
            ai4 = 1 / cdcd * (xi / rd * sdcd + 2 * atan((et * (x + q * cd) + x * (r + x) * sd) / (xi * (r + x) * cd)))
        end
        ai3 = (y * cd / rd - ale + sd * log(rd)) / cdcd
        ak1 = xi * (d11 - y11 * sd) / cd
        ak3 = (q * y11 - y * d11) / cd
        aj3 = (ak1 - aj2 * sd) / cd
        aj6 = (ak3 - aj5 * sd) / cd
    else
        rd2 = rd * rd
        ai3 = (et / rd + y * q / rd2 - ale) / 2
        ai4 = xi * y / rd2 / 2
        ak1 = xi * q / rd * d11
        ak3 = sd / rd * (xi2 * d11 - 1)
        aj3 = -xi / rd2 * (q2 * d11 - 1 / 2)
        aj6 = -y / rd2 * (xi2 * d11 - 1 / 2)
    end

    xy = xi * y11
    ai1 = -xi / rd * cd - ai4 * sd
    ai2 = log(rd) + ai3 * sd
    ak2 = 1 / r + ak3 * sd
    ak4 = xy * cd - ak1 * sd
    aj1 = aj5 * cd - aj6 * sd
    aj4 = -xy - aj2 * cd + aj3 * sd

    qx = q * x11
    qy = q * y11

    coeff = disl1 / 2π
    u[1]  += coeff * (-xi * qy - tt - alp3 * ai1 * sd)
    u[2]  += coeff * (-q / r + alp3 * y / rd * sd)
    u[3]  += coeff * (q * qy - alp3 * ai2 * sd)
    u[4]  += coeff * (xi2 * q * y32 - alp3 * aj1 * sd)
    u[5]  += coeff * (xi * q / r3 - alp3 * aj2 * sd)
    u[6]  += coeff * (-xi * q2 * y32 - alp3 * aj3 * sd)
    u[7]  += coeff * (-xi * fy - d * x11 + alp3 * (xy + aj4) * sd)
    u[8]  += coeff * (-ey + alp3 * (1 / r + aj5) * sd)
    u[9]  += coeff * (q * fy - alp3 * (qy - aj6) * sd)
    u[10] += coeff * (-xi * fz - y * x11 + alp3 * ak1 * sd)
    u[11] += coeff * (-ez + alp3 * y * d11 * sd)
    u[12] += coeff * (q * fz + alp3 * ak2 * sd)

    coeff = disl2 / 2π
    u[1]  += coeff * (-q / r + alp3 * ai3 * sdcd)
    u[2]  += coeff * (-et * qx - tt - alp3 * xi / rd * sdcd)
    u[3]  += coeff * (q * qx + alp3 * ai4 * sdcd)
    u[4]  += coeff * (xi * q / r3 + alp3 * aj4 * sdcd)
    u[5]  += coeff * (et * q / r3 + qy + alp3 * aj5 * sdcd)
    u[6]  += coeff * (-q2 / r3 + alp3 * aj6 * sdcd)
    u[7]  += coeff * (-ey + alp3 * aj1 * sdcd)
    u[8]  += coeff * (-et * gy - xy * sd + alp3 * aj2 * sdcd)
    u[9]  += coeff * (q * gy + alp3 * aj3 * sdcd)
    u[10] += coeff * (-ez - alp3 * ak3 * sdcd)
    u[11] += coeff * (-et * gz - xy * cd - alp3 * xi * d11 * sdcd)
    u[12] += coeff * (q * gz - alp3 * ak4 * sdcd)

    coeff = disl3 / 2π
    u[1]  += coeff * (q * qy - alp3 * ai3 * sdsd)
    u[2]  += coeff * (q * qx + alp3 * xi / rd * sdsd)
    u[3]  += coeff * (et * qx + xi * qy - tt - alp3 * ai4 * sdsd)
    u[4]  += coeff * (-xi * q2 * y32 - alp3 * aj4 * sdsd)
    u[5]  += coeff * (-q2 / r3 - alp3 * aj5 * sdsd)
    u[6]  += coeff * (q * q2 * y32 - alp3 * aj6 * sdsd)
    u[7]  += coeff * (q * fy - alp3 * aj1 * sdsd)
    u[8]  += coeff * (q * gy - alp3 * aj2 * sdsd)
    u[9]  += coeff * (-q * hy - alp3 * aj3 * sdsd)
    u[10] += coeff * (q * fz + alp3 * ak3 * sdsd)
    u[11] += coeff * (q * gz + alp3 * xi * d11 * sdsd)
    u[12] += coeff * (-q * hz + alp3 * ak4 * sdsd)
end

@inline function uc(u::A, sc1::B1, sc2::B2, xi::T, et::T, q::T, z::T, disl1::T, disl2::T, disl3::T,
    ) where {T <: Number,A <: AbstractArray{T},B1 <: NTuple{12,T},B2 <: NTuple{24,T}}
    alp1, alp2, alp3, alp4, alp5, sd, cd, sdsd, cdcd, sdcd, s2d, sc2d = sc1
    xi2, et2, q2, r, r2, r3, r5, y, d, tt, alx, ale, x11, y11, x32, y32, ey, ez, fy, fz, gy, gz, hy, hz = sc2

    fill!(u, zero(T))
    c = d + z
    x53 = (8 * r2 + 9 * r * xi + 3 * xi2) * x11 * x11 * x11 / r2
    y53 = (8 * r2 + 9 * r * et + 3 * et2) * y11 * y11 * y11 / r2
    h = q * cd - z
    z32 = sd / r3 - h * y32
    z53 = 3 * sd / r5 - h * y53
    y0 = y11 - xi2 * y32
    z0 = z32 - xi2 * z53
    ppy = cd / r3 + q * y32 * sd
    ppz = sd / r3 - q * y32 * cd
    qq = z * y32 + z32 + z0
    qqy = 3 * c * d / r5 - qq * sd
    qqz = 3 * c * y / r5 - qq * cd + q * y32
    xy = xi * y11
    # qx = q * x11
    qy = q * y11
    qr = 3 * q / r5
    # cqx = c * q * x53
    cdr = (c + d) / r3
    yy0 = y / r3 - y0 * cd

    coeff = disl1 / 2π
    u[1]  += coeff * (alp4 * xy * cd - alp5 * xi * q * z32)
    u[2]  += coeff * (alp4 * (cd / r + 2 * qy * sd) - alp5 * c * q / r3)
    u[3]  += coeff * (alp4 * qy * cd - alp5 * (c * et / r3 - z * y11 + xi2 * z32))
    u[4]  += coeff * (alp4 * y0 * cd - alp5 * q * z0)
    u[5]  += coeff * (-alp4 * xi * (cd / r3 + 2 * q * y32 * sd) + alp5 * c * xi * qr)
    u[6]  += coeff * (-alp4 * xi * q * y32 * cd + alp5 * xi * (3 * c * et / r5 - qq))
    u[7]  += coeff * (-alp4 * xi * ppy * cd - alp5 * xi * qqy)
    u[8]  += coeff * (alp4 * 2 * (d / r3 - y0 * sd) * sd - y / r3 * cd - alp5 * (cdr * sd - et / r3 - c * y * qr))
    u[9]  += coeff * (-alp4 * q / r3 + yy0 * sd + alp5 * (cdr * cd + c * d * qr - (y0 * cd + q * z0) * sd))
    u[10] += coeff * (alp4 * xi * ppz * cd - alp5 * xi * qqz)
    u[11] += coeff * (alp4 * 2 * (y / r3 - y0 * cd) * sd + d / r3 * cd - alp5 * (cdr * cd + c * d * qr))
    u[12] += coeff * (yy0 * cd - alp5 * (cdr * sd - c * y * qr - y0 * sdsd + q * z0 * cd))

    coeff = disl2 / 2π
    u[1]  += coeff * (alp4 * cd / r - qy * sd - alp5 * c * q / r3)
    u[2]  += coeff * (alp4 * y * x11 - alp5 * c * et * q * x32)
    u[3]  += coeff * (-d * x11 - xy * sd - alp5 * c * (x11 - q2 * x32))
    u[4]  += coeff * (-alp4 * xi / r3 * cd + alp5 * c * xi * qr + xi * q * y32 * sd)
    u[5]  += coeff * (-alp4 * y / r3 + alp5 * c * et * qr)
    u[6]  += coeff * (d / r3 - y0 * sd + alp5 * c / r3 * (1 - 3 * q2 / r2))
    u[7]  += coeff * (-alp4 * et / r3 + y0 * sdsd - alp5 * (cdr * sd - c * y * qr))
    u[8]  += coeff * (alp4 * (x11 - y * y * x32) - alp5 * c * ((d + 2 * q * cd) * x32 - y * et * q * x53))
    u[9]  += coeff * (xi * ppy * sd + y * d * x32 + alp5 * c * ((y + 2 * q * sd) * x32 - y * q2 * x53))
    u[10] += coeff * (-q / r3 + y0 * sdcd - alp5 * (cdr * cd + c * d * qr))
    u[11] += coeff * (alp4 * y * d * x32 - alp5 * c * ((y - 2 * q * sd) * x32 + d * et * q * x53))
    u[12] += coeff * (-xi * ppz * sd + x11 - d * d * x32 - alp5 * c * ((d - 2 * q * cd) * x32 - d * q2 * x53))

    coeff = disl3 / 2π
    u[1]  += coeff * (-alp4 * (sd / r + qy * cd) - alp5 * (z * y11 - q2 * z32))
    u[2]  += coeff * (alp4 * 2 * xy * sd + d * x11 - alp5 * c * (x11 - q2 * x32))
    u[3]  += coeff * (alp4 * (y * x11 + xy * cd) + alp5 * q * (c * et * x32 + xi * z32))
    u[4]  += coeff * (alp4 * xi / r3 * sd + xi * q * y32 * cd + alp5 * xi * (3 * c * et / r5 - 2 * z32 - z0))
    u[5]  += coeff * (alp4 * 2 * y0 * sd - d / r3 + alp5 * c / r3 * (1 - 3 * q2 / r2))
    u[6]  += coeff * (-alp4 * yy0 - alp5 * (c * et * qr - q * z0))
    u[7]  += coeff * (alp4 * (q / r3 + y0 * sdcd) + alp5 * (z / r3 * cd + c * d * qr - q * z0 * sd))
    u[8]  += coeff * (-alp4 * 2 * xi * ppy * sd - y * d * x32 + alp5 * c * ((y + 2 * q * sd) * x32 - y * q2 * x53))
    u[9]  += coeff * (-alp4 * (xi * ppy * cd - x11 + y * y * x32) + alp5 * (c * ((d + 2 * q * cd) * x32 - y * et * q * x53) + xi * qqy))
    u[10] += coeff * (-et / r3 + y0 * cdcd - alp5 * (z / r3 * sd - c * y * qr - y0 * sdsd + q * z0 * cd))
    u[11] += coeff * (alp4 * 2 * xi * ppz * sd - x11 + d * d * x32 - alp5 * c * ((d - 2 * q * cd) * x32 - d * q2 * x53))
    u[12] += coeff * (alp4 * (xi * ppz * cd + y * d * x32) + alp5 * (c * ((y - 2 * q * sd) * x32 + d * et * q * x53) + xi * qqz))
end

@inline function shared_constants_1(α::T, dip::T) where {T <: Number}
    alp1 = (1 - α) / 2
    alp2 = α / 2
    alp3 = (1 - α) / α
    alp4 = 1 - α
    alp5 = α
    sd = sind(dip)
    cd = cosd(dip)
    sdsd = sd^2
    cdcd = cd^2
    sdcd = sd * cd
    s2d = 2 * sdcd
    sc2d = cdcd - sdsd
    (alp1, alp2, alp3, alp4, alp5, sd, cd, sdsd, cdcd, sdcd, s2d, sc2d)
end

@inline function shared_constants_2(xi::T, et::T, q::T, sd::T, cd::T, kxi::T, ket::T) where {T <: Number}
    xi2 = xi * xi
    et2 = et * et
    q2 = q * q
    r2 = xi2 + et2 + q2
    r = sqrt(r2)

    r3 = r * r2
    r5 = r3 * r2
    y = et * cd + q * sd
    d = et * sd - q * cd

    if q ≈ zero(T)
        tt = zero(T)
    else
        tt = atan(xi * et / (q * r))
    end

    if kxi ≈ one(T)
        alx = -log(r - xi)
        x11 = zero(T)
        x32 = zero(T)
    else
        rxi = r + xi
        alx = log(rxi)
        x11 = 1 / (r * rxi)
        x32 = (r + rxi) * x11 * x11 / r
    end

    if ket ≈ one(T)
        ale = -log(r - et)
        y11 = zero(T)
        y32 = zero(T)
    else
        ret = r + et
        ale = log(ret)
        y11 = 1 / (r * ret)
        y32 = (r + ret) * y11 * y11 / r
    end

    ey = sd / r - y * q / r3
    ez = cd / r + d * q / r3
    fy = d / r3 + xi2 * y32 * sd
    fz = y / r3 + xi2 * y32 * cd
    gy = 2 * x11 * sd - y * q * x32
    gz = 2 * x11 * cd + d * q * x32
    hy = d * q * x32 + xi * q * y32 * sd
    hz = y * q * x32 + xi * q * y32 * cd

    (xi2, et2, q2, r, r2, r3, r5, y, d, tt, alx, ale, x11, y11, x32, y32, ey, ez, fy, fz, gy, gz, hy, hz)
end


text:
using ReservoirComputing
using GaussianProcesses

#model parameters
const approx_res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sigma = 0.1
const beta = 0.1
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 1
const out_size = 1
extended_states = true
const h_steps = 2

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

esn = ESN(approx_res_size,
    train,
    degree,
    radius,
    activation = activation,
    sigma = sigma,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size+out_size)
#test predict
output = ESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)
#test h steps predict
output = ESNpredict_h_steps(esn, predict_len, h_steps, test, W_out)
@test size(output) == (out_size, predict_len)

#test esgp
mean = MeanZero()
kernel = Lin(1.0)

gp = ESGPtrain(esn, mean, kernel, lognoise = -2.0, optimize = false)

esgp_output, sigma_esgp = ESGPpredict(esn, predict_len, gp)
@test size(output) == (out_size, predict_len)
esgp_output, sigma_esgp = ESGPpredict_h_steps(esn, predict_len, h_steps, test, gp)
@test size(output) == (out_size, predict_len)

#model parameters
const first_activation = tanh
const second_activation = identity
const first_lambda = 0.45
const second_lambda = 0.3

esn = dafESN(approx_res_size,
            train,
            degree,
            radius,
            first_lambda,
            second_lambda,
            first_activation = first_activation,
            second_activation = second_activation,
            sigma = sigma,
            alpha = alpha,
            nla_type = nla_type,
            extended_states = extended_states)

W_out = ESNtrain(esn, beta)
@test size(W_out) == (out_size, esn.res_size+out_size)
#test predict
output = dafESNpredict(esn, predict_len, W_out)
@test size(output) == (out_size, predict_len)
#test h steps predict
output = dafESNpredict_h_steps(esn, predict_len, h_steps, test, W_out)
@test size(output) == (out_size, predict_len)


text:
@with_kw mutable struct Domain{FT}
    # Something CUDA
    max_streams::IndexT
    streams::Nothing

    # Elem-centered
    matElemlist::VD{IndexT} # material indexset
    nodelist::VD{IndexT}    # elemToNode connectivity

    lxim::VD{IndexT}         # element connectivity through face g
    lxip::VD{IndexT}
    letam::VD{IndexT}
    letap::VD{IndexT}
    lzetam::VD{IndexT}
    lzetap::VD{IndexT}

    elemBC::VD{Int}         # elem face symm/free-surf flag g

    e::VD{FT}             # energy g
    d_e::VD{FT}           # change in energy g

    p::VD{FT}             # pressure g

    q::VD{FT}             # q g
    ql::VD{FT}            # linear term for q g
    qq::VD{FT}            # quadratic term for q g

    v::VD{FT}             # relative volume g

    volo::VD{FT}          # reference volume g
    delv::VD{FT}          # m_vnew - m_v g
    vdov::VD{FT}          # volume derivative over volume g

    arealg::VD{FT}        # char length of an element g

    ss::VD{FT}            # "sound speed" g

    elemMass::VD{FT}      # mass g

    # TODO: Not making distinction between pointers to arrays down here

    vnew::VD{FT}          # new relative volume -- temporary g

    delv_xi::VD{FT}       # velocity gradient -- temporary g
    delv_eta::VD{FT}
    delv_zeta::VD{FT}

    delx_xi::VD{FT}       # coordinate gradient -- temporary g
    delx_eta::VD{FT}
    delx_zeta::VD{FT}

    dxx::VD{FT}           # principal strains -- temporary g
    dyy::VD{FT}
    dzz::VD{FT}

    # Node-centered g

    x::VD{FT}             # coordinates g
    y::VD{FT}
    z::VD{FT}

    xd::VD{FT}            # velocities g
    yd::VD{FT}
    zd::VD{FT}

    xdd::VD{FT}           # accelerations g
    ydd::VD{FT}
    zdd::VD{FT}

    fx::VD{FT}            # forces g
    fy::VD{FT}
    fz::VD{FT}

    dfx::VD{FT}          # AD of the forces g
    dfy::VD{FT}
    dfz::VD{FT}

    nodalMass::VD{FT}     # mass g
    h_nodalMass::Vector{FT}     # mass - host g

    # device pointers for comms g
    # TODO: Not sure how to store the pointers or if even necessary
    # Real_t *d_delv_xi       # velocity gradient -- temporary g
    # Real_t *d_delv_eta
    # Real_t *d_delv_zeta

    # Real_t *d_x             # coordinates g
    # Real_t *d_y
    # Real_t *d_z

    # Real_t *d_xd            # velocities g
    # Real_t *d_yd
    # Real_t *d_zd

    # Real_t *d_fx            # forces g
    # Real_t *d_fy
    # Real_t *d_fz
    # Boundary nodesets

    symmX::VD{IndexT}      # symmetry plane nodesets
    symmY::VD{IndexT}
    symmZ::VD{IndexT}

    nodeElemCount::VD{Int}
    nodeElemStart::VD{Int}
    nodeElemCornerList::VD{IndexT}

    # Parameters

    dtfixed::FT                # fixed time increment g
    deltatimemultlb::FT
    deltatimemultub::FT
    stoptime::FT               # end time for simulation g
    dtmax::FT                  # maximum allowable time increment g
    cycle::Int                 # iteration count for simulation g

    dthydro_h::FT             # hydro time constraint g
    d_dthydro_h::FT           # AD change of the hydro time constraint g
    dtcourant_h::FT           # courant time constraint g
    d_dtcourant_h::FT         # AD of the courant time constraint g
    bad_q_h::IndexT              # flag to indicate Q error g
    bad_vol_h::IndexT            # flag to indicate volume error g

    # cuda Events to indicate completion of certain kernels
    # TODO Will check later how this works with KA
    # cudaEvent_t time_constraint_computed;
    time_h::FT                # current time g
    deltatime_h::FT           # variable time increment g

    u_cut::FT                 # velocity tolerance g
    hgcoef::FT                # hourglass control g
    qstop::FT                 # excessive q indicator g
    monoq_max_slope::FT
    monoq_limiter_mult::FT
    e_cut::FT                 # energy tolerance g
    p_cut::FT                 # pressure tolerance g
    ss4o3::FT
    q_cut::FT                 # q tolerance g
    v_cut::FT                 # relative volume tolerance g
    qlc_monoq::FT             # linear term coef for q g
    qqc_monoq::FT             # quadratic term coef for q g
    qqc::FT
    eosvmax::FT
    eosvmin::FT
    pmin::FT                  # pressure floor g
    emin::FT                  # energy floor g
    dvovmax::FT               # maximum allowable volume change g
    refdens::FT               # reference density g

    m_colLoc::IndexT
    m_rowLoc::IndexT
    m_planeLoc::IndexT
    m_tp::IndexT
    sizeX::IndexT
    sizeY::IndexT
    sizeZ::IndexT
    maxPlaneSize::IndexT
    maxEdgeSize::IndexT

    numElem::IndexT
    # padded_numElem::IndexT

    numNode::IndexT
    # padded_numNode::IndexT

    numSymmX::IndexT
    numSymmY::IndexT
    numSymmZ::IndexT

    octantCorner::IndexT

    # Region information
    numReg::Int                     # number of regions (def:11)
    balance::Int                    # Load balance between regions of a domain (def: 1)
    cost::Int                       # imbalance cost (def: 1)
    regElemSize::Vector{Int}        # Size of region sets
    regCSR::VD{Int}                 # records the begining and end of each region
    regReps::VD{Int}                # records the rep number per region
    regNumList::VD{IndexT}         # Region number per domain element
    regElemlist::VD{IndexT}        # region indexset
    regSorted::VD{IndexT}          # keeps index of sorted regions

    # MPI-Related additional data

    # TODO I think we can handle this differently

    # IndexT m_numRanks;
    # IndexT& numRanks() { return m_numRanks ; }

    # # Used in setup
    # IndexT m_rowMin, m_rowMax;
    # IndexT m_colMin, m_colMax;
    # IndexT m_planeMin, m_planeMax ;

    # # Communication Work space
    # Real_t *commDataSend ;
    # Real_t *commDataRecv ;

    # Real_t *d_commDataSend ;
    # Real_t *d_commDataRecv ;

    # # Maximum number of block neighbors
    # MPI_Request recvRequest[26] ; # 6 faces + 12 edges + 8 corners
    # MPI_Request sendRequest[26] ; # 6 faces + 12 edges + 8 corners
end


text:
import Base: ∈,
             isempty

export Complement

"""
    Complement{N<:Real, S<:LazySet{N}}

Type that represents the complement of a convex set.

### Fields

- `X` -- convex set

### Notes

Since `X` is assumed to be closed, unless `X` is empty or the universe, its
complement is open (i.e., not closed).
Since `X` is assumed to be closed, unless `X` is empty, the universe, or a
half-space, its complement is not convex.

The complement of the complement is the original set again.

### Examples

```jldoctest
julia> B = BallInf(zeros(2), 1.); C = Complement(B)
Complement{Float64,BallInf{Float64}}(BallInf{Float64}([0.0, 0.0], 1.0))

julia> Complement(C)
BallInf{Float64}([0.0, 0.0], 1.0)
```
"""
struct Complement{N<:Real, S<:LazySet{N}}
    X::S
end

# the complement of the complement is the original set again
Complement(C::Complement) = C.X

"""
    dim(C::Complement)

Return the dimension of the complement of a convex set.

### Input

- `C` -- complement of a convex set

### Output

The dimension of the complement of a convex set.
"""
function dim(C::Complement)::Int
    return dim(C.X)
end

"""
    ∈(x::AbstractVector{N}, C::Complement{N})::Bool where {N<:Real}

Check whether a given point is contained in the complement of a convex set.

### Input

- `x` -- point/vector
- `C` -- complement of a convex set

### Output

`true` iff the vector is contained in the complement.

### Algorithm

```math
    x ∈ X^C ⟺ x ∉ X
```
"""
function ∈(x::AbstractVector{N}, C::Complement{N})::Bool where {N<:Real}
    @assert length(x) == dim(C)
    return x ∉ C.X
end

"""
    isempty(C::Complement)::Bool

Return if the complement of a convex set is empty or not.

### Input

- `C` -- complement of a convex set

### Output

`false` unless the original set is universal.

### Algorithm

We use the `isuniversal` method.
"""
function isempty(C::Complement)::Bool
    return isuniversal(C.X)
end


text:
module contactForceLaw_2_boxes

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "DryWood"

zDim = 1.0
boxDim = 0.3

@assembly ThreeDFiles begin
  world = Modia3D.Object3D(visualizeFrame=false)

  boxMoving     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(boxDim,boxDim,boxDim,rsmall=0.0) , "DryWood", vmat1; contactMaterial = cmat); r=[1.0, 0.0, boxDim/2], fixed=false ) # ,R=ModiaMath.rot2(-pi/3) ) #

  box     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(5.0,2.0,zDim,rsmall=0.0) , "DryWood", vmat1; contactMaterial = cmat); r=[0.0, 0.0, -zDim/2], fixed=true )
end

gravField = Modia3D.UniformGravityField(g=0.0, n=[0,0,-1])
threeD = ThreeDFiles(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=false, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

# Modia3D.visualizeAssembly!( threeD )

model = Modia3D.SimulationModel( threeD )
ModiaMath.print_ModelVariables(model)
result = ModiaMath.simulate!(model; stopTime=2.0, tolerance=1e-8,interval=0.001, log=true)

ModiaMath.plot(result, ["boxMoving.r","boxMoving.v"])

println("... success of contactForceLaw_2_boxes.jl!")

end


text:
function test_l2_regularizer(backend::Backend, T, eps)
  println("-- Testing L2 regularizer on $(typeof(backend)){$T}...")

  coef = rand()
  param = rand(T, 2,3,4,5) - convert(T, 0.5)
  param_blob = make_blob(backend, param)
  regu = L2Regu(coef)

  loss = forward(backend, regu, 1.0, param_blob)
  expected_loss = coef * vecnorm(param)^2
  @test -eps < loss - expected_loss < eps

  grad_blob = make_zero_blob(backend, T, size(param))
  backward(backend, regu, 1.0, param_blob, grad_blob)
  grad = zeros(T, size(param))
  copy!(grad, grad_blob)

  @test all(-eps .< grad - 2coef*param .< eps)
end

function test_l2_regularizer(backend::Backend)
  test_l2_regularizer(backend, Float32, 1e-5)
  test_l2_regularizer(backend, Float64, 1e-10)
end

if test_cpu
  test_l2_regularizer(backend_cpu)
end
if test_gpu
  test_l2_regularizer(backend_gpu)
end


text:
# Configuration for nearly neutral simuation
function dfe_disadv( x::Int64 )
  dfe_deleterious( x::Int64; alpha=dfe_disadv_alpha, beta=dfe_disadv_beta )
end

const nn_simtype = 1    # nearly neutral infinite alleles
const type_str = "D4"
const mu_list_flag = false
@everywhere const N_list     = [25,50,100,200,400,800,1600,3200,6400]     # popsize
const N_mu_list= [0.5,1.0,2.0,4.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 400000      # Generations after burn-in
const fix_minimum = 0.45  # Fraction of population needed to "fix"
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_disadv
dfe_str = "dfe_disadv alpha_disadv=0.2 beta_disadv=1.0"
dfe_disadv_alpha=0.2
dfe_disadv_beta=1.0

const use_poplist=false   # if true, save a list of populations as a check of correctness

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral