text:
# # Single Inverted Pendulum
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Single-Pendulum.ipynb)
#
# This is a classical inverted pendulum. A ball of mass $m$ is attached to a
# massless beam of length $L$.  The beam is actuated with a torque $T$ and
# we assume viscous friction exists with a coefficient of $c$.

module SinglePendulum  #jl

using NeuralNetworkAnalysis
using NeuralNetworkAnalysis: SingleEntryVector

# The following option determines whether the falsification settings should be
# used or not. The falsification settings are sufficient to show that the safety
# property is violated. Concretely, we start from an initial point and use a
# smaller time step.
const falsification = true;

# ## Model
#
# The governing equation of motion can be obtained as:
#
# ```math
# \ddot\theta = \dfrac{g}{L} \sin(\theta) + \dfrac{1}{m L^2} (T - c\dot\theta)
# ```
# where $θ$ is the angle that the link makes with the upward vertical axis.
# The state vector is $[θ, θ']$.

# model constants
const m = 0.5
const L = 0.5
const c = 0.0
const g = 1.0
const gL = g/L
const mL = 1/(m*L^2)

@taylorize function single_pendulum!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = gL * sin(x[1]) + mL*(x[3] - c*x[2])
    dx[3] = zero(x[3])
    return dx
end

controller = read_nnet(@modelpath("Single-Pendulum",
                                  "controller_single_pendulum.nnet"));

# ## Specification
#
# The initial set is
# ```math
# (\theta, \dot\theta) \in [1, 1.2] \times [0, 0.2].
# ```
# The safety specification is that $θ ∈ [0, 1]$ for all control steps
# $10 ≤ n_t ≤ 20$.

X0 = Hyperrectangle([1.1, 0.1], [0.1, 0.1])
if falsification
    ## choose a single point in the initial states (here: the top-most one)
    X0 = Singleton(high(X0))
end
U0 = ZeroSet(1)
ivp = @ivp(x' = single_pendulum!(x), dim: 3, x(0) ∈ X0 × U0)
vars_idx = Dict(:state_vars=>1:2, :control_vars=>3)

period = 0.05  # control period
if falsification
    k = 11
else
    k = 20
end
T = k * period  # time horizon

prob = ControlledPlant(ivp, controller, vars_idx, period)

unsafe_states = HalfSpace(SingleEntryVector(1, 3, -1.0), -1.0)

predicate = X -> overapproximate(X, Hyperrectangle) ⊆ unsafe_states
function predicate_sol(sol; silent::Bool=false)
    for F in sol
        for R in F
            t = tspan(R)
            if t.lo >= 0.5 && t.hi <= 1.0 && predicate(R)
                silent || println("violation for time range $t")
                return true
            end
        end
    end
    return false
end;

# ## Results

import DifferentialEquations

alg = TMJets(abstol=1e-7, orderT=4, orderQ=1)
alg_nn = Ai2()

function benchmark(; silent::Bool=false)
    ## We solve the controlled system:
    silent || println("flowpipe construction")
    res_sol = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg)
    sol = res_sol.value
    silent || print_timed(res_sol)

    ## Next we check the property for an overapproximated flowpipe:
    silent || println("property checking")
    solz = overapproximate(sol, Zonotope)
    res_pred = @timed predicate_sol(solz; silent=silent)
    silent || print_timed(res_pred)
    if res_pred.value
        silent || println("The property is violated.")
    else
        silent || println("The property may be satisfied.")
    end

    ## We also compute some simulations:
    silent || println("simulation")
    trajectories = falsification ? 1 : 10
    res_sim = @timed simulate(prob, T=T, trajectories=trajectories,
                              include_vertices=!falsification)
    sim = res_sim.value
    silent || print_timed(res_sim)

    return solz, sim
end;

benchmark(silent=true)  # warm-up
res = @timed benchmark()  # benchmark
sol, sim = res.value
println("total analysis time")
print_timed(res);

# Finally we plot the results:

using Plots
import DisplayAs

vars = (0, 1)
fig = plot(ylab="θ")
unsafe_states_projected = project(unsafe_states, [vars[2]])
time = Interval(0.5, 1.0)
unsafe_states_projected = cartesian_product(time, unsafe_states_projected)
plot!(fig, unsafe_states_projected, color=:red, alpha=:0.2, lab="unsafe states")
time0 = Singleton([0.0])
if !falsification
    plot!(fig, cartesian_product(time0, project(initial_state(prob), [vars[2]])),
          lab="X₀")
end
plot!(fig, sol, vars=vars, color=:yellow, lab="")
if falsification
    xlims!(0, T)
    ylims!(0.95, 1.22)
else
    xlims!(0, T)
    ylims!(0.55, 1.3)
end
lab_sim = falsification ? "simulation" : ""
plot_simulation!(fig, sim; vars=vars, color=:black, lab=lab_sim)
fig = DisplayAs.Text(DisplayAs.PNG(fig))
## savefig("Single-Pendulum.png")
fig

#-

end  #jl
nothing  #jl


text:
using HypothesisTests
using Test
using DelimitedFiles

@testset "Hotelling utility functions" begin
    HT = HypothesisTests

    @test HT.checkdims(rand(3, 2), rand(4, 2)) == (2, 3, 4)
    # Mismatched number of variables
    @test_throws DimensionMismatch HT.checkdims(rand(4, 3), rand(4, 6))
    # Empty observations
    @test_throws ArgumentError HT.checkdims(rand(0, 4), rand(2, 4))

    Sx = [ 3.0 -1.5 0.0
          -1.5  1.0 0.5
           0.0  0.5 1.0]
    Sy = fill(4.5, (3, 3))
    out = [3.5 0.5      1.5
           0.5 2.166667 1.833333
           1.5 1.833333 2.166667]
    @test HT.poolcov!(Sx, 2, Sy, 1) ≈ out atol=1e-6
    # Input is modified
    @test Sx ≈ out atol=1e-6

    # Positive semi-definite but not positive definite
    P = [ 1.0000  0.7426  0.1601 -0.7000 0.5500
          0.7426  1.0000 -0.2133 -0.5818 0.5000
          0.1601 -0.2133  1.0000 -0.1121 0.1000
         -0.7000 -0.5818 -0.1121  1.0000 0.4500
          0.5500  0.5000  0.1000  0.4500 1.0000]
    # Positive definite
    Q = [ 1.0 -0.5  0.0
         -0.5  1.0 -0.5
          0.0 -0.5  1.0]
    @test @inferred(HT.At_Binv_A(ones(5), P)) ≈ -0.8008792 atol=1e-6
    @test @inferred(HT.At_Binv_A(ones(3), Q)) ≈ 10.0
end

@testset "One sample Hotelling's T²" begin
    # Columns are calcium, iron, protein, vitamin A, vitamin C
    nutrient = readdlm(joinpath(@__DIR__, "data", "nutrient.txt"))[:,2:end]

    t = OneSampleHotellingT2Test(nutrient, [1000, 15, 60, 800, 75])
    @test nobs(t) == 737
    @test dof(t) == (5, 732)
    @test pvalue(t) ≈ 0.0 atol=eps()
    @test t.T² ≈ 1758.5413137 atol=1e-6
    @test t.F ≈ 349.7968048 atol=1e-6
    let out = sprint(show, t)
        @test occursin("reject h_0", out) && !occursin("fail to", out)
    end

    # Columns are survey answers: husbands' to questions 1-4, wives' to questions 1-4
    spouse = readdlm(joinpath(@__DIR__, "data", "spouse.txt"))

    # Paired
    p = OneSampleHotellingT2Test(spouse[:,1:4], spouse[:,5:end])
    @test nobs(p) == 30
    @test dof(p) == (4, 26)
    @test pvalue(p) ≈ 0.039369144 atol=1e-6
    @test p.T² ≈ 13.127840261 atol=1e-6
    @test p.F ≈ 2.942446955 atol=1e-6
    let out = sprint(show, p)
        @test occursin("reject h_0", out) && !occursin("fail to", out)
    end
end

@testset "Two sample Hotelling's T²" begin
    # Columns are type, length, left, right, bottom, top, diag
    swiss = readdlm(joinpath(@__DIR__, "data", "swiss3.txt"))
    genuine = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "real", 2:end])
    counterfeit = convert(Matrix{Float64}, swiss[view(swiss, :, 1) .== "fake", 2:end])

    eq = EqualCovHotellingT2Test(genuine, counterfeit)
    @test nobs(eq) == (100, 100)
    @test dof(eq) == (6, 193)
    @test pvalue(eq) ≈ 0.0 atol=eps()
    @test eq.T² ≈ 2412.4506855 atol=1e-6
    @test eq.F ≈ 391.9217023 atol=1e-6
    @test occursin("reject h_0", sprint(show, eq))

    un = UnequalCovHotellingT2Test(genuine, counterfeit)
    @test nobs(un) == (100, 100)
    @test dof(un) == (6, 193)
    @test pvalue(un) ≈ 0.0 atol=eps()
    @test un.T² ≈ 2412.4506855 atol=1e-6
    @test un.F ≈ 391.9217023 atol=1e-6
    @test occursin("reject h_0", sprint(show, un))
end

ext: jl


text:
import Pkg
Pkg.add(["DelimitedFiles"])
using DelimitedFiles

# Inclusion of all part to run the simulation
include("typs_and_structs.jl")
include("initialize.jl")
include("equilibriumcalculation.jl")

# Reading the input file
path = pwd()
len, maxtime, dt, g, γ = readinput(path * "/input/input.txt")
# len, maxtime, dt, g, γ = readinput()
println("Chain length : ", len, "\nRuntime : ", maxtime, "\nδt : ", dt, "\ngravity : ", g, "\nγ : ", γ)

lbemoments = moments1d_64(ones(len), zeros(len), zeros(len))
println(lbemoments.velocity)

equilibriumdist = calc_equilibrium_distribution(lbemoments, g)
println(equilibriumdist.f0)

# n = []
# for i in lbemoments.height
#     push!(n, i)
# end
# println(n)


text:
function hydraulic_variables(h, T, constants, soil_parameters)

  d_gamma = -0.0001425 .- 4.76e-7 .*(T.-273);
  rho_w = 1000 .- 7.3e-3.*(T.-(273+4)).^2+3.79e-5.*(T.-(273+4)).^3;  #water density [kg m^-3]
  D_a = 2.12e-5.*(T./273).^2; # air diffusivity
  rho_sv = 1e-3.*exp.(19.84 .- 4975.9 ./T); # saturated vapor pressure [kg/m^3]
  drho_vs = 4.9759./(T.^2).*exp.(19.84 .-4975.9 ./T); #saturated vapor pressure derivative [kg m^-3 K^-1]

  Se = 1. ./(1 .+(soil_parameters.alpha.*abs.(h)).^soil_parameters.n).^soil_parameters.m;
  Se[h.>0] .= 1;
  #Se(Se<1e-3) = 1e-3;
  theta_l = soil_parameters.theta_r.+(soil_parameters.theta_s-soil_parameters.theta_r).*Se;
  theta_a = soil_parameters.theta_s.-theta_l;
  Hr = exp.(h.*constants.M.*constants.g./(constants.Rg.*T));
  Hr[h.>0] .= 1;
  rho_vs = Hr.*rho_sv;
  theta_v = rho_vs.*theta_a./rho_w;

  Cm = -(soil_parameters.alpha.*abs.(h)).^soil_parameters.n.*soil_parameters.m.*soil_parameters.n.*
      (soil_parameters.theta_s.-soil_parameters.theta_r).*((soil_parameters.alpha.*abs.(h)).^soil_parameters.n.+1).^(.-soil_parameters.m.-1)./h;
  Cm[h.>=0] .= 1e-13;
  Cv = rho_sv./rho_w.*(soil_parameters.theta_s.-soil_parameters.theta_r).*(constants.g.*constants.M./(constants.Rg.*T).*
      Hr.*(1 .-Se).-Cm.*Hr./(soil_parameters.theta_s.-soil_parameters.theta_r));
  Cv[h.>0] .= 0;

  tau = theta_a.^(7/3)./soil_parameters.theta_s.^2;
  D = tau.*theta_a.*D_a;
  eta = 9.5 .+ 3 .*theta_l./soil_parameters.theta_s.-8.5.*exp.(-((1 .+2.6 ./sqrt.(soil_parameters.fc)).*(theta_l./soil_parameters.theta_s)).^4);

  #hydraulic parameters
  tau = theta_a.^(7/3)./soil_parameters.theta_s.^2;
  D = tau.*theta_a.*D_a;
  eta = 9.5 .+ 3 .* theta_l./soil_parameters.theta_s .- 8.5 .* exp.(-((1+2.6/sqrt(soil_parameters.fc)).*(theta_l./soil_parameters.theta_s)).^4);

  soil_hydraulicC = Cm+Cv;
  soil_hydraulicK = soil_parameters.Ks.*sqrt.(Se).*(1 .-(1 .-Se.^(1/soil_parameters.m)).^soil_parameters.m).^2;
  soil_hydraulicK_vh = D./(constants.Rg.*T.*rho_w).*rho_sv.*constants.M.*constants.g.*Hr;
  soil_hydraulicK_h = soil_hydraulicK+soil_hydraulicK_vh;

  soil_hydraulicK_lT = soil_hydraulicK.*h.*d_gamma./(soil_parameters.G_wT.*constants.gamma_o);
  soil_hydraulicK_vT = D.*eta.*Hr.*drho_vs./rho_w;
  soil_hydraulicK_T = soil_hydraulicK_lT + soil_hydraulicK_vT;

  soil_hydraulic=shydraulic(soil_hydraulicC,soil_hydraulicK,soil_hydraulicK_vh,soil_hydraulicK_h,
                    soil_hydraulicK_lT,soil_hydraulicK_vT,soil_hydraulicK_T)

  rho_vs = rho_vs[end];

  return theta_l, theta_v, soil_hydraulic, rho_vs

end


text:
module terra

using GDAL
using ArchGDAL; const AG = ArchGDAL

using CoordinateTransformations
using StaticArrays
using GeoFormatTypes
const GFT = GeoFormatTypes
# include("shp/GeoDataFrames.jl")

import Statistics
import Statistics: mean

include("GeoArray.jl")
include("st_bbox.jl")
include("st_affine.jl")
include("st_crs.jl")

include("st_dim.jl")
include("st_coordinates.jl")

include("tools_ratser.jl")
include("tools_Ipaper.jl")

include("utils/utils.jl")
include("raster/Raster.jl")

include("gdal_info.jl")
include("gdal_polygonize.jl")

include("st_plot.jl")
include("st_subset.jl")
include("st_read.jl")
include("st_write.jl")
include("st_convert.jl")

# include("hydro/hydro.jl")
# include("hydro/snap_pour_points.jl")

export GeoArray

export coords
export centercoords
export indices

export compose!

export epsg!
export crs!

export -,+,*,/

end


text:
module _XApEn
export XApEn
using Statistics: mean, std
    """
        XAp, Phi = XApEn(Sig)

    Returns the cross-approximate entropy estimates (`XAp`) and the average
    number of matched vectors (`Phi`) for m = [0,1,2], estimated for the data
    sequences contained in 'Sig' using the default parameters:
    embedding dimension = 2, time delay = 1,
    radius distance threshold= 0.2*SD(`Sig`), logarithm = natural

    **NOTE**: XApEn is direction-dependent. Thus, the first column of
    `Sig` is used as the template data sequence, and the second
    column is the matching sequence.``

        XAp, Phi = XApEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    Returns the cross-approximate entropy estimates (`XAp`) between the data
    sequences contained in 'Sig' using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer  [default: 2] \n
    `tau`   - Time Delay, a positive integer        [default: 1]  \n
    `r`     - Radius Distance Threshold, a positive scalar [default: 0.2*SD(`Sig`)] \n
    `Logx`  - Logarithm base, a positive scalar     [default: natural]  \n

    # See also `XSampEn`, `XFuzzEn`, `XMSEn`, `ApEn`, `SampEn`, `MSEn`

    # References:
        [1] Steven Pincus and Burton H. Singer,
            "Randomness and degrees of irregularity."
            Proceedings of the National Academy of Sciences
            93.5 (1996): 2083-2088.

        [2] Steven Pincus,
            "Assessing serial irregularity and its implications for health."
            Annals of the New York Academy of Sciences
            954.1 (2001): 245-267.

    """
    function XApEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2,
        tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    size(Sig,2) > size(Sig,1) ? Sig = transpose(Sig) : nothing

    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (r >=0) ? nothing : error("r:     must be a positive value")
    (Logx>0) ? nothing : error("Logx: must be a positive number > 0")

    S1 = Sig[:,1]; S2 = Sig[:,2]
    Counter = 1*(abs.(S1 .- transpose(S2)) .<= r)
    M = vcat(m*ones(Int,N-(m*tau)), repeat((m-1):-1:1,inner=tau))
    XAp = zeros(m+1)
    Phi = zeros(m+2)
    for n = 1:N - tau
        ix = findall(Counter[n, :] .== 1)
        for k = 1:M[n]
            ix = ix[ix .+ (k*tau) .<= N]
            isempty(ix) ? break : nothing
            p1 = repeat(transpose(S1[n:tau:n+(tau*k)]), outer=length(ix))
            p2 = S2[ix .+ transpose(collect(0:tau:(k*tau)))]
            ix = ix[findall(maximum(abs.(p1 - p2),dims=2) .<= r)]
            Counter[n, ix] .+= 1
        end
    end

    Phi[1] = log(Logx, N)/N
    #Phi[2] = mean(log.(Logx, sum(Counter.>0,dims=1)/N))
    Temp = sum(Counter.>0,dims=1); Temp = Temp[Temp.!=0]
    Phi[2] = mean(log.(Logx, Temp/N))
    XAp[1] = Phi[1] - Phi[2]
    for k = 0:m-1
        ai = sum(Counter.>(k+1),dims=1)/(N-(k+1)*tau)
        bi = sum(Counter.>k,dims=1)/(N-(k*tau))
        ai = ai[ai.!=0]
        bi = bi[bi.!=0]
        Phi[k+3] = sum(log.(Logx, ai))/(N-(k+1)*tau)
        XAp[k+2]= sum(log.(Logx, bi))/(N-(k*tau)) - Phi[k+3]
    end

    return XAp, Phi
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
abstract type NACABase{T<:AbstractFloat} <: AbstractAirfoil{T} end

include("NACA4.jl")

function NACA(::Type{T}, num::AbstractString, sharp_TE::Bool=false) where {T<:AbstractFloat}
    if length(num) == 4
        NACA4(T, num, sharp_TE)
    else
        throw(ArgumentError("unsupported airfoil NACA$num"))
    end
end

NACA(num, sharp_TE) = NACA(Float64, num, sharp_TE)

macro NACA_str(num)
    NACA(num, false)
end

macro NACA_str(num, s)
    if (s != "s")
        error("use NACA\"XXXX\"s for sharp trailing edge")
    end
    NACA(num, true)
end


text:
module ShiftRightTests

include(joinpath(dirname(@__FILE__), "..", "src", "AndersonMoore.jl"))
using .AndersonMoore

# test shiftRight! firmvalue example
function firmvalue()::Bool

neq=4::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[0.  -0.4  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  -0.4  0.  0.  0.  1.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.1e-16::Float64,atol=0.0::Float64)
end;

# test shiftRight! firmvalue3Leads2Lags example
function firmvalue3Leads2Lags()::Bool

neq=4::Int64;nlag=2::Int64;nlead=3::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  -0.4  0.  0.  0.  0.  0.  0.  0.  1.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

# test shiftRight! example7 example
function example7()::Bool

neq=4::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[0.  0.  0.  1.  0.  0.  1.  -1.  0.  0.  0.  0.;
0.  0.  0.  -0.66  0.  -1.1  0.  1.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  1.  -1.;
0.  0.  0.  0.  0.  0.  0.  -0.66  0.  -1.1  0.  1.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

# test shiftRight! oneEquationNoLead example
function oneEquationNoLead()::Bool

neq=1::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[2.  3.  0.]::Array{Float64,2}

shifted=[0.  2.  3.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

# test shiftRight! reliablePaperExmpl example
function reliablePaperExmpl()::Bool

neq=5::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[0.  0.  0.  -0.5  0.  0.  -1.  0.7  -0.5  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  0.  0.  0.  -0.5  0.  0.  -1.  0.7  -0.5  1.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

# test shiftRight! athan example
function athan()::Bool

neq=9::Int64;nlag=1::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[-0.293  -0.764  0.293  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  -2.  -1.  1.  1.  0.  0.  2.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  -0.911  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.293  -0.764  0.293  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -2.  -1.  1.  1.  0.  0.  2.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.911  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  -1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  1.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

# test shiftRight! habitmod example
function habitmod()::Bool

neq=12::Int64;nlag=4::Int64;nlead=1::Int64
qRows=(neq*nlead)::Int64;qCols=(neq*(nlag+nlead))::Int64

toShift=[0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  -1.  -1.  0.  0.  3.45419  0.  0.  -0.74118  0.00208517  0.494116  -0.1149  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.88879  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.52802  0.  0.11279  0.031833  0.024734  0.  0.  0.  0.  0.  0.  0.  -0.71526  0.  0.042076  -0.22829  -0.012267  0.  0.  0.  0.  0.  0.  0.  0.22215  0.  -0.14172  0.26679  -0.11959  0.  0.  0.  0.  0.  0.  0.  -0.42308  0.  -0.61076  -0.10946  0.17172  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  -0.00011277  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.46436  0.  0.16131  0.02744  0.17087  0.  0.  0.  0.  0.  0.  0.  -1.3425  0.  0.16794  -0.35818  -0.09837  0.  0.  0.  0.  0.  0.  0.  1.06382  0.  -0.02255  0.51594  0.08949  0.  0.  0.  0.  0.  0.  0.  -0.46672  0.  -0.11369  -1.10657  -0.26021  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  -0.00141  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
-0.23292  0.  0.08512  -0.04213  0.14149  0.  0.  0.  0.  0.  0.  0.  0.03313  0.  0.15223  0.01572  -0.34762  0.  0.  0.  0.  0.  0.  0.  0.70552  0.  -0.24697  0.21888  -0.07805  0.  0.  0.  0.  0.  0.  0.  -0.56316  0.  -0.08268  -0.16054  -0.68637  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  -0.00348  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.99853  0.  0.  0.  0.  0.  -0.00147  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]::Array{Float64,2}

shifted=[0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  -1.  -1.  0.  0.  3.45419  0.  0.  -0.74118  0.00208517  0.494116  -0.1149;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.88879  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.52802  0.  0.11279  0.031833  0.024734  0.  0.  0.  0.  0.  0.  0.  -0.71526  0.  0.042076  -0.22829  -0.012267  0.  0.  0.  0.  0.  0.  0.  0.22215  0.  -0.14172  0.26679  -0.11959  0.  0.  0.  0.  0.  0.  0.  -0.42308  0.  -0.61076  -0.10946  0.17172  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  -0.00011277;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.46436  0.  0.16131  0.02744  0.17087  0.  0.  0.  0.  0.  0.  0.  -1.3425  0.  0.16794  -0.35818  -0.09837  0.  0.  0.  0.  0.  0.  0.  1.06382  0.  -0.02255  0.51594  0.08949  0.  0.  0.  0.  0.  0.  0.  -0.46672  0.  -0.11369  -1.10657  -0.26021  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  -0.00141;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.23292  0.  0.08512  -0.04213  0.14149  0.  0.  0.  0.  0.  0.  0.  0.03313  0.  0.15223  0.01572  -0.34762  0.  0.  0.  0.  0.  0.  0.  0.70552  0.  -0.24697  0.21888  -0.07805  0.  0.  0.  0.  0.  0.  0.  -0.56316  0.  -0.08268  -0.16054  -0.68637  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  -0.00348;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -0.99853  0.  0.  0.  0.  0.  -0.00147  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.;
0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.]::Array{Float64,2}

shiftResult=shiftRight!(toShift,neq)
isapprox(shiftResult,shifted,rtol=0.0::Float64,atol=0.1e-16::Float64)
end;

end


text:
function replace_headers(s::AbstractString)

    subs =  [
        r"\[(https?://.+?) (.+?)\]"        => s"[\2](\1)",  # hrefs http and https

        r"\[source:(.+?) (.+?)\]"         => s"[\2](https://hirlam.org/trac/browser/\1)",  # Source links
        r"\[raw-attachment:(.+?) (.+?)\]" => s"[\2](\1)",  # Attachements
        r"^==\s*'''\s*(.+?)\s*'''\s*=="    => s"## \1",     # level 2 header (not bold in md)
        r"^=\s*'''\s*(.+)\s*'''\s*="       => s"# \1",      # level 1 headers (not bold in md)
        r"^===\s*(.+?)\s*==="              => s"### \1",    # level 3 headers
        r"^==\s*(.+?)\s*=="                => s"## \1",     # level 2 headers

        r"^=\s*(.+?)\s*="                  => s"# \1",      # level 1 headers
        r"'''\s*(.+?)\s*'''"              => s"**\1**",    # Bold
        r"''\s*(.+?)\s*''"                => s"*\1*",      # Italic
    ]

    return foldl(replace, subs, init = s)
end


text:
# MIRT/z-all.jl

include("utility/z-list.jl")
include("plot/z-list.jl")
include("../data/z-list.jl")

include("algorithm/z-list.jl")
include("fbp/z-list.jl")
include("io/z-list.jl")
include("mri/z-list.jl")
include("nufft/z-list.jl")

include("regularize/z-list.jl")
include("system/z-list.jl")


text:
export count_not_missing, count_missing, countmissing, countnotmissing

countmissing(args...) = count_missing(args...)
countnotmissing(args...) = count_not_missing(args...)

count_not_missing(x) = length(x) - count(ismissing, x)

count_missing(x) = count(ismissing, x)

function count_not_missing(::Type{S}, x::Vector{Union{T, Missing}}) where {S, T}
    @assert isbitstype(T)
    res = zero(S)
    @inbounds for i in 1:length(x)
        res += !ismissing(x[i])
    end
    res
end

count_not_missing(x::Vector{Union{T, Missing}}) where T = count_not_missing(Int, x)

count_missing(::Type{S}, x::Vector{Union{T, Missing}}) where {S, T} =
    length(x) - count_not_missing(S, x)

count_missing(x::Vector{Union{T, Missing}}) where T = count_missing(Int, x)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
push!(LOAD_PATH,"../src/")
using Documenter, PosDefManifold

makedocs(
   sitename="PosDefManifold",
   modules=[PosDefManifold],
   pages =  [
      "index.md",
      "introToRiemannianGeometry.md",
      "MainModule.md",
      "riemannianGeometry.md",
      "linearAlgebra.md",
      "statistics.md",
      "signalProcessing.md",
      "test.md",
   ]
)

deploydocs(
    repo = "github.com/Marco-Congedo/PosDefManifold.jl.git",
    target = "build",
    devurl = "dev",
)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DataDiskImage <: SwaggerModel
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun

    function DataDiskImage(;lun=nothing)
        o = new()
        validate_property(DataDiskImage, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        o
    end
end # type DataDiskImage

const _property_map_DataDiskImage = Dict{Symbol,Symbol}(Symbol("lun")=>Symbol("lun"))
const _property_types_DataDiskImage = Dict{Symbol,String}(Symbol("lun")=>"Int32")
Base.propertynames(::Type{ DataDiskImage }) = collect(keys(_property_map_DataDiskImage))
Swagger.property_type(::Type{ DataDiskImage }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DataDiskImage[name]))}
Swagger.field_name(::Type{ DataDiskImage }, property_name::Symbol) =  _property_map_DataDiskImage[property_name]

function check_required(o::DataDiskImage)
    true
end

function validate_property(::Type{ DataDiskImage }, name::Symbol, val)
end


text:
module SplitApplyPlot

using Tables: rows, columns, getcolumn, columnnames
using StructArrays: components, GroupPerm, StructArray
using Colors: RGB
using AbstractPlotting
using AbstractPlotting: automatic, Automatic, PlotFunc
import AbstractPlotting.MakieLayout: hidexdecorations!,
                                     hideydecorations!,
                                     hidedecorations!,
                                     linkaxes!,
                                     linkxaxes!,
                                     linkyaxes!
using PooledArrays: PooledArray
using KernelDensity: kde, pdf
using StatsBase: fit, histrange, Histogram, normalize, weights, AbstractWeights, sturges
using DataAPI: refarray
import FileIO

export hideinnerdecorations!, deleteemptyaxes!
export arguments, Entry, Entries, AxisEntries
export renamer, nonnumeric
export density, histogram, linear, visual, data, dims, mapping
export draw, draw!
export facet!

include("arguments.jl")
include("scales.jl")
include("entries.jl")
include("utils.jl")
include("facet.jl")
include("helpers.jl")
include("algebra/layer.jl")
include("algebra/layers.jl")
include("algebra/processing.jl")
include("transformations/grouping.jl")
include("transformations/visual.jl")
include("transformations/linear.jl")
include("transformations/density.jl")
include("transformations/histogram.jl")

end


text:
import Observables
import AbstractPlotting

export Torus
export update

"""
    Represents a torus of revolution.

fields: q, r, R, segments, color and observables.
"""
mutable struct Torus <: Sprite
    q::Biquaternion
    r::Float64
    R::Float64
    segments::Int
    color::AbstractPlotting.RGBAf0
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Torus(joint, scene)

Construct a torus of revolution with the given configuration `q` and `scene`, and the
optional arguments: smaller radius `r`, bigger radius `R`, `segments`, `color` and
`transparency`.
"""
function Torus(q::Biquaternion,
               scene::AbstractPlotting.Scene;
               r::Float64 = 0.025,
               R::Float64 = 1.0,
               segments::Int = 36,
               color::AbstractPlotting.RGBAf0 = AbstractPlotting.RGBAf0(0.1, 0.1, 0.1, 0.9),
               transparency::Bool = false)
    torus = constructtorus(q, r, R, segments = segments)
    colorarray = fill(color, segments, segments)
    observable = buildsurface(scene, torus, colorarray, transparency = transparency)
    Torus(q, r, R, segments, color, observable)
end

"""
    update(torus, q)

Update a Torus by changing its observable with the given `torus` and configuration `q`.
"""
function update(torus::Torus, q::Biquaternion)
    torus.q = q
    value = constructtorus(torus.q, torus.r, torus.R, segments = torus.segments)
    updatesurface(value, torus.observable)
end


text:
#################################################
# Classical ExpRK method caches
abstract type ExpRKCache <: OrdinaryDiffEqMutableCache end
abstract type ExpRKConstantCache <: OrdinaryDiffEqConstantCache end

# Precomputation of exponential-like operators
"""
    expRK_operators(alg,dt,A) -> ops

Compute operator(s) for an ExpRK algorithm. `dt` is the time step and `A` is
the matrix form of the linear operator (from either a linear problem or a
SplitODEProblem). All ExpRK methods that use caching operators should implement
this method.
"""
expRK_operators(alg::ExponentialAlgorithm, dt, A) =
  error("$alg does not support caching operators at the moment.")
expRK_operators(::LawsonEuler, dt, A) = exp(dt * A)
expRK_operators(::NorsettEuler, dt, A) = phi(dt * A, 1)[2]
function expRK_operators(::ETDRK2, dt, A)
  P = phi(dt * A, 2)
  return P[2], P[3] # ϕ1(hA), ϕ2(hA)
end
function expRK_operators(::ETDRK3, dt, A)
  Phalf = phi(dt/2 * A, 1)
  A21 = 0.5Phalf[2]
  P = phi(dt * A, 3)
  phi1, phi2, phi3 = P[2], P[3], P[4]
  A3 = phi1
  B1 = 4phi3 - 3phi2 + phi1
  B2 = -8phi3 + 4phi2
  B3 = 4phi3 - phi2
  return A21, A3, B1, B2, B3
end
function expRK_operators(::ETDRK4, dt, A)
  P = phi(dt * A, 3)
  Phalf = phi(dt/2 * A, 1)
  A21 = 0.5Phalf[2] # A32 = A21
  A41 = (dt/4 * A) * Phalf[2]^2
  A43 = Phalf[2]
  B1 = P[2] - 3P[3] + 4P[4]
  B2 = 2P[3] - 4P[4] # B3 = B2
  B4 = -P[3] + 4P[4]
  return A21, A41, A43, B1, B2, B4
end
function expRK_operators(::HochOst4, dt, A)
  P = phi(dt * A, 3)
  Phalf = phi(dt/2 * A, 3)
  A21 = 0.5Phalf[2]
  A31 = A21 - Phalf[3]
  A32 = Phalf[3]
  A41 = P[2] - 2P[3]
  A42 = P[3]
  A52 = 0.5Phalf[3] - P[4] + 0.25P[3] - 0.5Phalf[4]
  A54 = 0.25Phalf[3] - A52
  A51 = 0.5Phalf[2] - 2A52 - A54
  B1 = P[2] - 3P[3] + 4P[4]
  B4 = -P[3] + 4P[4]
  B5 = 4P[3] - 8P[4]
  return A21, A31, A32, A41, A42, A51, A52, A54, B1, B4, B5
end

# Unified constructor for constant caches
for (Alg, Cache) in [(:LawsonEuler, :LawsonEulerConstantCache),
                     (:NorsettEuler, :NorsettEulerConstantCache),
                     (:ETDRK2, :ETDRK2ConstantCache),
                     (:ETDRK3, :ETDRK3ConstantCache),
                     (:ETDRK4, :ETDRK4ConstantCache),
                     (:HochOst4, :HochOst4ConstantCache)]
  @eval struct $Cache{opType,FType} <: ExpRKConstantCache
    ops::opType # precomputed operators
    uf::FType   # derivative wrapper
  end

  @eval function alg_cache(alg::$Alg,u,rate_prototype,uEltypeNoUnits,
    uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
    if alg.krylov
      ops = nothing # no caching
    else
      isa(f, SplitFunction) || throw(ArgumentError("Caching can only be used with SplitFunction"))
      A = size(f.f1.f) == () ? convert(Number, f.f1.f) : convert(AbstractMatrix, f.f1.f)
      ops = expRK_operators(alg, dt, A)
    end
    if isa(f, SplitFunction) || DiffEqBase.has_jac(f)
      uf = nothing
    else
      uf = UDerivativeWrapper(f,t,p)
    end
    return $Cache(ops,uf)
  end
end

"""
    alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist)

Construct the non-standard caches (not uType or rateType) for ExpRK integrators.

`plist` is a list of integers each corresponding to the order of a `phiv(!)`
call in `perform_step!`.
"""
function alg_cache_expRK(alg::OrdinaryDiffEqExponentialAlgorithm,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist)
  n = length(u); T = eltype(u)
  # Allocate cache for ForwardDiff
  if isa(f, SplitFunction) || DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  # Allocate cache for the Jacobian
  if isa(f, SplitFunction)
    J = nothing
  elseif DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), n, n)
  end
  if alg.krylov
    ops = nothing # no caching
    # Build up caches used by Krylov phiv
    m = min(alg.m, n)
    Ks = KrylovSubspace{T}(n, m)
    phiv_cache = PhivCache{T}(m, maximum(plist))
    ws = [Matrix{T}(undef, n, plist[i] + 1) for i = 1:length(plist)]
    KsCache = (Ks, phiv_cache, ws)
  else
    KsCache = nothing
    # Precompute the operators
    A = size(f.f1.f) == () ? convert(Number, f.f1.f) : convert(AbstractMatrix, f.f1.f)
    ops = expRK_operators(alg, dt, A)
  end
  return uf, jac_config, J, ops, KsCache
end

@cache struct LawsonEulerCache{uType,rateType,JCType,FType,JType,expType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  G::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  exphA::expType
  KsCache::KsType
end

function alg_cache(alg::LawsonEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                            # uType caches
  rtmp, G, du1 = (zero(rate_prototype) for i = 1:3)             # rateType caches
  # other caches
  # This is different from other ExpRK integrators because LawsonEuler only
  # needs expv.
  n = length(u); T = eltype(u)
  # Allocate caches for ForwardDiff
  if isa(f, SplitFunction) || DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  # Allocate cache for the Jacobian
  if isa(f, SplitFunction)
    J = nothing
  elseif DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), n, n)
  end
  if alg.krylov
    exphA = nothing # no caching
    m = min(alg.m, n)
    Ks = KrylovSubspace{T}(n, m)
    expv_cache = ExpvCache{T}(m)
    KsCache = (Ks, expv_cache)
  else
    KsCache = nothing
    A = size(f.f1.f) == () ? convert(Number, f.f1.f) : convert(AbstractMatrix, f.f1.f)
    exphA = expRK_operators(alg, dt, A)
  end
  LawsonEulerCache(u,uprev,tmp,dz,rtmp,G,du1,jac_config,uf,J,exphA,KsCache)
end

@cache struct NorsettEulerCache{uType,rateType,JCType,FType,JType,expType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  G::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  phihA::expType
  KsCache::KsType
end

function alg_cache(alg::NorsettEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                # uType caches
  rtmp, G, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  plist = (1,)
  uf,jac_config,J,phihA,KsCache = alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist) # other caches
  NorsettEulerCache(u,uprev,tmp,dz,rtmp,G,du1,jac_config,uf,J,phihA,KsCache)
end

@cache struct ETDRK2Cache{uType,rateType,JCType,FType,JType,opType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  F2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  ops::opType
  KsCache::KsType
end

function alg_cache(alg::ETDRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                 # uType caches
  rtmp, F2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  plist = (2,2)
  uf,jac_config,J,ops,KsCache = alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist) # other caches
  ETDRK2Cache(u,uprev,tmp,dz,rtmp,F2,du1,jac_config,uf,J,ops,KsCache)
end

@cache struct ETDRK3Cache{uType,rateType,JCType,FType,JType,opType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  Au::rateType
  F2::rateType
  F3::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  ops::opType
  KsCache::KsType
end

function alg_cache(alg::ETDRK3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                         # uType caches
  rtmp, Au, F2, F3, du1 = (zero(rate_prototype) for i = 1:5) # rateType caches
  plist = (1,3,3,3)
  uf,jac_config,J,ops,KsCache = alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist) # other caches
  ETDRK3Cache(u,uprev,tmp,dz,rtmp,Au,F2,F3,du1,jac_config,uf,J,ops,KsCache)
end

@cache struct ETDRK4Cache{uType,rateType,JCType,FType,JType,opType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  Au::rateType
  F2::rateType
  F3::rateType
  F4::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  ops::opType
  KsCache::KsType
end

function alg_cache(alg::ETDRK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                             # uType caches
  rtmp, Au, F2, F3, F4, du1 = (zero(rate_prototype) for i = 1:6) # rateType caches
  plist = (1,1,3,3,3,3)
  uf,jac_config,J,ops,KsCache = alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist) # other caches
  ETDRK4Cache(u,uprev,tmp,dz,rtmp,Au,F2,F3,F4,du1,jac_config,uf,J,ops,KsCache)
end

@cache struct HochOst4Cache{uType,rateType,JCType,FType,JType,opType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  Au::rateType
  F2::rateType
  F3::rateType
  F4::rateType
  F5::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  ops::opType
  KsCache::KsType
end

function alg_cache(alg::HochOst4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                                        # uType caches
  rtmp, rtmp2, Au, F2, F3, F4, F5, du1 = (zero(rate_prototype) for i = 1:8) # rateType caches
  plist = (3,3,3,3,3,3,3,3,3)
  uf,jac_config,J,ops,KsCache = alg_cache_expRK(alg,u,uEltypeNoUnits,uprev,f,t,dt,p,du1,tmp,dz,plist) # other caches
  HochOst4Cache(u,uprev,tmp,dz,rtmp,rtmp2,Au,F2,F3,F4,F5,du1,jac_config,uf,J,ops,KsCache)
end

####################################
# EPIRK method caches
function _phiv_timestep_caches(u_prototype, maxiter::Int, p::Int)
  n = length(u_prototype); T = eltype(u_prototype)
  u = similar(u_prototype)                      # stores the current state
  W = Matrix{T}(undef, n, p+1)                  # stores the w vectors
  P = Matrix{T}(undef, n, p+2)                  # stores output from phiv!
  Ks = KrylovSubspace{T}(n, maxiter)            # stores output from arnoldi!
  phiv_cache = PhivCache{T}(maxiter, p+1)       # cache used by phiv! (need +1 for error estimation)
  return u, W, P, Ks, phiv_cache
end

# Unified constructor for constant caches
for (Alg, Cache) in [(:Exp4, :Exp4ConstantCache),
                     (:EPIRK4s3A, :EPIRK4s3AConstantCache),
                     (:EPIRK4s3B, :EPIRK4s3BConstantCache),
                     (:EPIRK5s3, :EPIRK5s3ConstantCache),
                     (:EXPRB53s3, :EXPRB53s3ConstantCache),
                     (:EPIRK5P1, :EPIRK5P1ConstantCache),
                     (:EPIRK5P2, :EPIRK5P2ConstantCache)]
  @eval struct $Cache{FType} <: ExpRKConstantCache
    uf::FType   # derivative wrapper
  end
  @eval function alg_cache(alg::$Alg,u,rate_prototype,uEltypeNoUnits,
    uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
    if DiffEqBase.has_jac(f)
      uf = nothing
    else
      uf = UDerivativeWrapper(f,t,p)
    end
    return $Cache(uf)
  end
end

@cache struct Exp4Cache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::Exp4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                    # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  # TODO: units
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 3)
  B = fill(zero(T), n, 2)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 1)
  Exp4Cache(u,uprev,tmp,dz,rtmp,rtmp2,du1,jac_config,uf,K,J,B,KsCache)
end

@cache struct EPIRK4s3ACache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EPIRK4s3A,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                    # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 2)
  B = fill(zero(T), n, 5)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 4)
  EPIRK4s3ACache(u,uprev,tmp,dz,rtmp,rtmp2,du1,jac_config,uf,K,J,B,KsCache)
end

@cache struct EPIRK4s3BCache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EPIRK4s3B,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                    # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 2)
  B = fill(zero(T), n, 5)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 4)
  EPIRK4s3BCache(u,uprev,tmp,dz,rtmp,rtmp2,du1,jac_config,uf,K,J,B,KsCache)
end

@cache struct EPIRK5s3Cache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  k::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EPIRK5s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz, k = (zero(u) for i = 1:3)                 # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  B = fill(zero(T), n, 5)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 4)
  EPIRK5s3Cache(u,uprev,tmp,dz,k,rtmp,rtmp2,du1,jac_config,uf,J,B,KsCache)
end

@cache struct EXPRB53s3Cache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EXPRB53s3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                    # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 2)
  B = fill(zero(T), n, 5)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 4)
  EXPRB53s3Cache(u,uprev,tmp,dz,rtmp,rtmp2,du1,jac_config,uf,K,J,B,KsCache)
end

@cache struct EPIRK5P1Cache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EPIRK5P1,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                    # uType caches
  rtmp, rtmp2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 3)
  B = fill(zero(T), n, 4)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 3)
  EPIRK5P1Cache(u,uprev,tmp,dz,rtmp,rtmp2,du1,jac_config,uf,K,J,B,KsCache)
end

@cache struct EPIRK5P2Cache{uType,rateType,JCType,FType,matType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  rtmp2::rateType
  dR::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  K::matType
  J::JType
  B::matType
  KsCache::KsType
end
function alg_cache(alg::EPIRK5P2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
  tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  tmp, dz = (zero(u) for i = 1:2)                        # uType caches
  rtmp, rtmp2, dR, du1 = (zero(rate_prototype) for i = 1:4) # rateType caches
  # Allocate jacobian and caches for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), length(u), length(u))
  end
  # Allocate matrices
  n = length(u); T = eltype(u)
  K = Matrix{T}(undef, n, 3)
  B = fill(zero(T), n, 4)
  # Allocate caches for phiv_timestep
  maxiter = min(alg.m, n)
  KsCache = _phiv_timestep_caches(u, maxiter, 3)
  EPIRK5P2Cache(u,uprev,tmp,dz,rtmp,rtmp2,dR,du1,jac_config,uf,K,J,B,KsCache)
end

####################################
# Adaptive exponential Rosenbrock method caches

## Constant caches
for (Alg, Cache) in [(:Exprb32, :Exprb32ConstantCache),
                     (:Exprb43, :Exprb43ConstantCache)]
  @eval struct $Cache{FType} <: ExpRKConstantCache
    uf::FType   # derivative wrapper
  end
  @eval function alg_cache(alg::$Alg,u,rate_prototype,uEltypeNoUnits,
    uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
    if DiffEqBase.has_jac(f)
      uf = nothing
    else
      uf = UDerivativeWrapper(f,t,p)
    end
    return $Cache(uf)
  end
end

## Mutable caches
"""
    alg_cache_exprb(alg,uEltypeNoUnits,uprev,f,t,p,du1,tmp,dz,plist)

Construct the non-standard caches (not uType or rateType) for Exprb integrators.

`plist` is a list of integers each corresponding to the order of a `phiv(!)`
call in `perform_step!`.
"""
function alg_cache_exprb(alg::OrdinaryDiffEqAdaptiveExponentialAlgorithm,u,uEltypeNoUnits,uprev,f,t,p,du1,tmp,dz,plist)
  if f isa SplitFunction
    error("Algorithm $alg cannnot be used for split problems. Consider reformat to a regular `ODEProblem`")
  end
  n = length(u); T = eltype(u)
  # Allocate cache for ForwardDiff
  if DiffEqBase.has_jac(f)
    uf = nothing
    jac_config = nothing
  else
    uf = UJacobianWrapper(f,t,p)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end
  # Allocate cache for the Jacobian
  if DiffEqBase.has_jac(f) && f.jac_prototype !== nothing
    J = deepcopy(f.jac_prototype)
  else
    J = fill(zero(uEltypeNoUnits), n, n)
  end
  # Build up caches used by Krylov phiv
  m = min(alg.m, n)
  Ks = KrylovSubspace{T}(n, m)
  phiv_cache = PhivCache{T}(m, maximum(plist))
  ws = [Matrix{T}(undef, n, plist[i] + 1) for i = 1:length(plist)]
  KsCache = (Ks, phiv_cache, ws)
  return uf, jac_config, J, KsCache
end

@cache struct Exprb32Cache{uType,rateType,JCType,FType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  utilde::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  F2::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  KsCache::KsType
end
function alg_cache(alg::Exprb32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  utilde, tmp, dz = (zero(u) for i = 1:3)         # uType caches
  rtmp, F2, du1 = (zero(rate_prototype) for i = 1:3) # rateType caches
  plist = (3, 3)
  uf,jac_config,J,KsCache = alg_cache_exprb(alg,u,uEltypeNoUnits,uprev,f,t,p,du1,tmp,dz,plist) # other caches
  Exprb32Cache(u,uprev,utilde,tmp,dz,rtmp,F2,du1,jac_config,uf,J,KsCache)
end

struct Exprb43Cache{uType,rateType,JCType,FType,JType,KsType} <: ExpRKCache
  u::uType
  uprev::uType
  utilde::uType
  tmp::uType
  dz::uType
  rtmp::rateType
  Au::rateType
  F2::rateType
  F3::rateType
  du1::rateType
  jac_config::JCType
  uf::FType
  J::JType
  KsCache::KsType
end
function alg_cache(alg::Exprb43,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  utilde, tmp, dz = (zero(u) for i = 1:3)                 # uType caches
  rtmp, Au, F2, F3, du1 = (zero(rate_prototype) for i = 1:5) # rateType caches
  plist = (1,4,4,4)
  uf,jac_config,J,KsCache = alg_cache_exprb(alg,u,uEltypeNoUnits,uprev,f,t,p,du1,tmp,dz,plist) # other caches
  Exprb43Cache(u,uprev,utilde,tmp,dz,rtmp,Au,F2,F3,du1,jac_config,uf,J,KsCache)
end

####################################
# Multistep exponential method caches

#=
  Fsal separately the linear and nonlinear part, as well as the nonlinear
  part in the previous time step.
=#
mutable struct ETD2Fsal{rateType}
  lin::rateType
  nl::rateType
  nlprev::rateType
end
ETD2Fsal(rate_prototype) = ETD2Fsal(zero(rate_prototype),zero(rate_prototype),zero(rate_prototype))
function recursivecopy!(dest::ETD2Fsal, src::ETD2Fsal)
  recursivecopy!(dest.lin, src.lin)
  recursivecopy!(dest.nl, src.nl)
  recursivecopy!(dest.nlprev, src.nlprev)
end

struct ETD2ConstantCache{expType} <: OrdinaryDiffEqConstantCache
  exphA::expType
  phihA::expType
  B1::expType # ϕ1(hA) + ϕ2(hA)
  B0::expType # -ϕ2(hA)
end

function alg_cache(alg::ETD2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  A = size(f.f1.f) == () ? convert(Number, f.f1.f) : convert(AbstractMatrix, f.f1.f)
  Phi = phi(dt*A, 2)
  ETD2ConstantCache(Phi[1], Phi[2], Phi[2] + Phi[3], -Phi[3])
end

@cache struct ETD2Cache{uType,rateType,expType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  utmp::uType
  rtmp1::rateType
  rtmp2::rateType
  exphA::expType
  phihA::expType
  B1::expType # ϕ1(hA) + ϕ2(hA)
  B0::expType # -ϕ2(hA)
end

function alg_cache(alg::ETD2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  A = size(f.f1.f) == () ? convert(Number, f.f1.f) : convert(AbstractMatrix, f.f1.f)
  Phi = phi(dt*A, 2)
  ETD2Cache(u,uprev,zero(u),zero(rate_prototype),zero(rate_prototype),Phi[1],Phi[2],Phi[2]+Phi[3],-Phi[3])
end


text:
using Surrogates

#=
#1D MOE
n = 30
lb = 0.0
ub = 5.0
x = Surrogates.sample(n,lb,ub,UniformSample())
f = x-> 2*x
y = f.(x)
#Standard definition
my_moe = MOE(x,y,lb,ub)
val = my_moe(3.0)

#Local surrogates redefinition
my_local_kind = [InverseDistanceStructure(p = 1.0),
                 KrigingStructure(p=1.0, theta=1.0)]
my_moe = MOE(x,y,lb,ub,k = 2,local_kind = my_local_kind)

#ND MOE
n = 30
lb = [0.0,0.0]
ub = [5.0,5.0]
x = sample(n,lb,ub,LatinHypercubeSample())
f = x -> x[1]*x[2]
y = f.(x)
my_moe_ND = MOE(x,y,lb,ub)
val = my_moe_ND((1.0,1.0))

#Local surr redefinition
my_locals = [InverseDistanceStructure(p = 1.0),
             RandomForestStructure(num_round=1)]
my_moe_redef = MOE(x,y,lb,ub,k = 2,local_kind = my_locals)
=#


text:
module CommunalHelperStationBlockTrack

using ..Ahorn, Maple

const switchStates = ["None", "On", "Off"]
const moveModes = ["None", "ForwardOneWay", "BackwardOneWay", "ForwardForce", "BackwardForce"]

@mapdef Entity "CommunalHelper/StationBlockTrack" StationBlockTrack(
    x::Integer,
    y::Integer,
    width::Integer=24,
    height::Integer=24,
    horizontal::Bool=false,
    trackSwitchState::String="None",
    moveMode::String="None",
    multiBlockTrack::Bool=false,
    indicator::Bool=true,
    indicatorColor::String="008080",
    indicatorIncomingColor::String="c92828",
)

const placements = Ahorn.PlacementDict(
    "Station Block Track ($orientation, $stateLabel) (Communal Helper)" => Ahorn.EntityPlacement(
        StationBlockTrack,
        "rectangle",
        Dict{String,Any}(
            "horizontal" => orientation == "Horizontal",
            "trackSwitchState" => state,
        ),
    ) for orientation in ["Vertical", "Horizontal"], (state, stateLabel) in zip(switchStates, ["No Switching", "Switch On", "Switch Off"])
)

Ahorn.minimumSize(entity::StationBlockTrack) = 24, 24

function Ahorn.resizable(entity::StationBlockTrack)
    horiz = Bool(get(entity.data, "horizontal", false))
    return (horiz, !horiz)
end

Ahorn.editingOptions(entity::StationBlockTrack) = Dict{String,Any}(
    "trackSwitchState" => switchStates,
    "moveMode" => moveModes,
)

# very weird rotation
function Ahorn.rotated(entity::StationBlockTrack, steps::Int)
    horiz = Bool(get(entity.data, "horizontal", false))
    trackSwitchState = get(entity.data, "trackSwitchState", "None")
    return StationBlockTrack(entity.x, entity.y, entity.height, entity.width, !horiz, trackSwitchState)
end

function Ahorn.selection(entity::StationBlockTrack)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))
    horiz = Bool(get(entity.data, "horizontal", false))

    return horiz ? Ahorn.Rectangle(x, y, width, 8) : Ahorn.Rectangle(x, y, 8, height)
end

const nodeSprite = "objects/CommunalHelper/stationBlock/tracks/outline/node"
const hTrack = "objects/CommunalHelper/stationBlock/tracks/outline/h"
const vTrack = "objects/CommunalHelper/stationBlock/tracks/outline/v"
const arrows = "objects/CommunalHelper/stationBlock/tracks/outline/arrows"

const noneColor = (255, 255, 255, 255) ./ 255
const onColor = (66, 167, 255, 255) ./ 255
const offColor = (255, 48, 131, 255) ./ 255

function getTrackColor(entity::StationBlockTrack)
    trackSwitchState = get(entity.data, "trackSwitchState", "None")
    if trackSwitchState == "On"
        return onColor
    elseif trackSwitchState == "Off"
        return offColor
    else
        return noneColor
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::StationBlockTrack)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))

    horiz = Bool(get(entity.data, "horizontal", false))

    color = getTrackColor(entity)

    moveMode = String(get(entity.data, "moveMode", "None"))
    backwardArrowCount = moveMode == "BackwardForce" ? 2 : (moveMode == "BackwardOneWay" ? 1 : 0)
    forwardArrowCount = moveMode == "ForwardForce" ? 2 : (moveMode == "ForwardOneWay" ? 1 : 0)

    if horiz
        tilesWidth = div(width, 8)

        for i in 2:tilesWidth-1
            Ahorn.drawImage(ctx, hTrack, x + (i - 1) * 8, y, tint=color)
        end

        Ahorn.drawImage(ctx, nodeSprite, x, y, 0, 0, 8, 8, tint=color)
        Ahorn.drawImage(ctx, nodeSprite, x + width - 8, y, 8, 0, 8, 8, tint=color)
        for i in 0:backwardArrowCount-1
            Ahorn.drawImage(ctx, arrows, x - 6 - 3 * i, y, 0, 8, 8, 8, tint=color)
        end
        for i in 0:forwardArrowCount-1
            Ahorn.drawImage(ctx, arrows, x + width - 2 + 3 * i, y, 8, 0, 8, 8, tint=color)
        end
    else
        tilesHeight = div(height, 8)

        for i in 2:tilesHeight-1
            Ahorn.drawImage(ctx, vTrack, x, y + (i - 1) * 8, tint=color)
        end

        Ahorn.drawImage(ctx, nodeSprite, x, y, 8, 8, 8, 8, tint=color)
        Ahorn.drawImage(ctx, nodeSprite, x, y + height - 8, 0, 8, 8, 8, tint=color)
        for i in 0:backwardArrowCount-1
            Ahorn.drawImage(ctx, arrows, x, y - 6 - i * 3, 0, 0, 8, 8, tint=color)
        end
        for i in 0:forwardArrowCount-1
            Ahorn.drawImage(ctx, arrows, x, y + height - 2 + i * 3, 8, 8, 8, 8, tint=color)
        end
    end
end

end


text:
#=
Main program for FWI
=#
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--generate_data"
            arg_type = Bool
            default = false
        "--version"
            arg_type = String
            default = "0000"
        "--verbose"
            arg_type = Bool
            default = false
    end

    return parse_args(s)
end

args = parse_commandline()
if !isdir("./$(args["version"])")
    mkdir("./$(args["version"])")
end

using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)
include("ops_imseq.jl")
include("../Ops/FWI/fwi_util.jl")
include("fwi_util_op.jl")
np = pyimport("numpy")

nz = 134
nx = 384
dz = 24. # meters
dx = 24.
nSteps = 2001
dt = 0.0025
f0 = 4.5
filter_para = [0, 0.1, 100.0, 200.0]
isAc = true
nPml = 32
nPad = 32 - mod((nz+2*nPml), 32)
nz_pad = nz + 2*nPml + nPad
nx_pad = nx + 2*nPml

# reflection
x_src = collect(4:8:384)
z_src = 2ones(Int64, size(x_src))
x_rec = collect(3:381)
z_rec = 2ones(Int64, size(x_rec))

# xwell
# z_src = collect(5:10:nz-5) #14->11srcs 10->15srcs
# x_src = 5ones(Int64, size(z_src))
# z_rec = collect(5:1:nz-5)
# x_rec = (nx-5) .* ones(Int64, size(z_rec))

para_fname = "./$(args["version"])/para_file.json"
survey_fname = "./$(args["version"])/survey_file.json"
data_dir_name = "./$(args["version"])/Data"
paraGen(nz_pad, nx_pad, dz, dx, nSteps, dt, f0, nPml, nPad, filter_para, isAc, para_fname, survey_fname, data_dir_name)
surveyGen(z_src, x_src, z_rec, x_rec, survey_fname)

tf_cp = constant(reshape(reinterpret(Float32,read("Mar_models/Model_Cp_true.bin")),(nz_pad, nx_pad)), dtype=Float64)
cs = zeros(nz_pad, nx_pad)
den = 1000.0 .* ones(nz_pad, nx_pad)
cp_pad_value = 3000.0

# tf_cp = constant(cp)
tf_cs = constant(cs)
tf_den = constant(den)

src = Matrix{Float64}(undef, 1, 2001)
# # src[1,:] = Float64.(reinterpret(Float32, read("../Ops/FWI/Src/params/ricker_10Hz.bin")))
src[1,:] = Float64.(reinterpret(Float32, read("../Ops/FWI/Src/params/Mar_source_2001.bin")))
# src = sourceGene(f0, nSteps, dt)
tf_stf = constant(repeat(src, outer=length(z_src)))
# tf_para_fname = tf.strings.join([para_fname])
tf_gpu_id0 = constant(0, dtype=Int32)
tf_gpu_id1 = constant(1, dtype=Int32)
nGpus = 2
tf_gpu_id_array = constant(collect(0:nGpus-1), dtype=Int32)
tf_shot_ids0 = constant(collect(Int32, 0:length(x_src)-1), dtype=Int32)
shot_id_points = Int32.(trunc.(collect(LinRange(0, length(z_src)-1, nGpus+1))))

function pad_cp(cp)
    tran_cp = cast(cp, Float64)
    return tf.pad(tran_cp, [nPml (nPml+nPad); nPml nPml], constant_values=3000.0)
end

# NOTE Generate Data
if args["generate_data"]
    println("Generate Test Data...")

    if !isdir("./$(args["version"])/Data")
        mkdir("./$(args["version"])/Data")
    end

    res = fwi_obs_op(tf_cp, tf_cs, tf_den, tf_stf, tf_gpu_id0, tf_shot_ids0, para_fname)

    config = tf.ConfigProto()
    config.intra_op_parallelism_threads = 24
    config.inter_op_parallelism_threads = 24
    sess = Session(config=config); init(sess);
    run(sess, res)
    error("Generate Data: Stop")
end

cp_init = reshape(reinterpret(Float32,read("Mar_models/Model_Cp_init_1D.bin")),(nz_pad, nx_pad))
tf_cp_inv = Variable(cp_init, dtype=Float64)

Mask = ones(nz_pad, nx_pad)
Mask[nPml+1:nPml+10,:] .= 0.0
tf_cp_inv_msk = tf_cp_inv .* constant(Mask) + constant(cp_init[1,1] .* (1. .- Mask))

# NOTE Compute FWI loss
# loss = constant(0.0)
# for i = 1:nGpus
#     global loss
#     tf_shot_ids = constant(collect(shot_id_points[i] : shot_id_points[i+1]), dtype=Int32)
#     loss += fwi_op(tf_cp_inv_msk, tf_cs, tf_den, tf_stf, tf_gpu_id_array[i], tf_shot_ids, para_fname)
# end
loss = fwi_op(tf_cp_inv_msk, tf_cs, tf_den, tf_stf, tf_gpu_id_array[1], tf_shot_ids0, para_fname)
gradCp = gradients(loss, tf_cp_inv)

if args["verbose"]
    sess = Session(); init(sess)
    println("Initial loss = ", run(sess, loss))
    g = gradients(loss, tfCtxInit.K)
    G = run(sess, g)
    pcolormesh(G); savefig("test.png"); close("all")
end

# Optimization
__cnt = 0
# invK = zeros(m,n)
function print_loss(l, Cp, gradCp)
    global __cnt, __l, __Cp, __gradCp
    if mod(__cnt,1)==0
        println("\niter=$__iter, eval=$__cnt, current loss=",l)
        # println("a=$a, b1=$b1, b2=$b2")
    end
    __cnt += 1
    __l = l
    __Cp = Cp
    __gradCp = gradCp
end

__iter = 0
function print_iter(rk)
    global __iter, __l
    if mod(__iter,1)==0
        println("\n************* ITER=$__iter *************\n")
    end
    __iter += 1
    open("./$(args["version"])/loss.txt", "a") do io
        writedlm(io, Any[__iter __l])
    end
    open("./$(args["version"])/Cp$__iter.txt", "w") do io
        writedlm(io, __Cp)
    end
    open("./$(args["version"])/gradCp$__iter.txt", "w") do io
        writedlm(io, __gradCp)
    end
end

config = tf.ConfigProto()
config.intra_op_parallelism_threads = 24
config.inter_op_parallelism_threads = 24
sess = Session(config=config); init(sess);
# cp_low_bd = 1500. .* ones(nz_pad, nx_pad)
# cp_high_bd = 5500. .* ones(nz_pad, nx_pad)
# cp_high_bd[nPml+1:nPml+10,:] .= 1500.0
opt = ScipyOptimizerInterface(loss, var_list=[tf_cp_inv], var_to_bounds=Dict(tf_cp_inv=> (1500.0, 5500.0)), method="L-BFGS-B",
    options=Dict("maxiter"=> 100, "ftol"=>1e-6, "gtol"=>1e-6))
@info "Optimization Starts..."
ScipyOptimizerMinimize(sess, opt, loss_callback=print_loss, step_callback=print_iter, fetches=[loss,tf_cp_inv,gradCp])

# adam = AdamOptimizer(learning_rate=50.0)
# op = minimize(adam, loss)
# sess = Session(); init(sess);
# for iter = 1:1000
#     _, misfit, cp, cpgrad = run(sess, [op, loss, tf_cp_inv, gradCp])

#     open("./$(args["version"])/Cp$iter.txt", "w") do io
#         writedlm(io, cp)
#     end

#     open("./$(args["version"])/loss.txt", "a") do io
#         writedlm(io, Any[iter misfit])
#     end

#     open("./$(args["version"])/gradCp$iter.txt", "w") do io
#         writedlm(io, cpgrad)
#     end
# end


text:
using Test
using SeisModels

@testset "1D models" begin
    @testset "Broadcasting" begin
        @test vp.((PREM, AK135), 1000) == (11.10541122721928, 11.103425443786982)
        @test vs.(PREM, 0:3) ≈ [3.6678, 3.667799890427708, 3.6677995617108317, 3.6677990138493715]
    end

    @testset "Find layer" begin
        @test_throws DomainError SeisModels.findlayer(PREM, -1)
        @test_throws DomainError SeisModels.findlayer(PREM, 1e6)
    end

    @testset "Evaluation" begin
        @testset "SteppedLayeredModel" begin
            let m = SteppedLayeredModel(1, 2, [0.25, 1], [2, 1], [0.5, 1], [4, 3],
                    true, [1, 2], [3, 4], [5, 6], [7, 8], [9, 10],
                    true, [100, 200], [300, 400])
                @test vp(m, 0) == 2
                @test vp(m, 0.26) == 1
                @test vs(m, 0.1) == 0.5
                @test vs(m, 1) == 1
                @test density(m, 0.1) == 4
                @test density(m, 0.9) == 3
                @test vph(m, 0.2) == 1
                @test vph(m, 0.25) == 2
                @test vpv(m, 0) == 3
                @test vpv(m, 1) == 4
                @test vsh(m, 0.24) == 5
                @test vsh(m, nextfloat(0.25)) == 6
                @test vsv(m, prevfloat(0.25)) == 7
                @test vsv(m, nextfloat(0.25)) == 8
                @test eta(m, prevfloat(0.25)) == 9
                @test eta(m, nextfloat(0.25)) == 10
                @test Qmu(m, 0) == Qμ(m, 0) == 100
                @test Qμ(m, 1) == 200
                @test Qκ(m, 0.2) == Qkappa(m, 0.2) == 300
                @test Qκ(m, 0.3) == 400
                @test mass(m, surface_radius(m)) == 4/3*π*(0.25e3^3*4e3 + (1e3^3 - 0.25e3^3)*3e3)
                @test mass(m) == mass(m, surface_radius(m))
                @test gravity(m, 0.75) ≈ 6.67428e-11*4/3*π*(0.25e3^3*4e3 + (0.75e3^3 - 0.25e3^3)*3e3)/0.75e3^2
                @test shear_modulus(m, 1) == 3e3*1e3^2/1e9
                @test bulk_modulus(m, 0.5) == (3e3*1e3^2/1e9 - 4/3*shear_modulus(m, 0.5))
                @test poissons_ratio(m, 0) == 0.4666666666666667
                # Evaluation via property symbol
                for (f, sym) in zip((vp, vs, density, vph, vpv, vsh, vsv, eta, Qμ,  Qκ),
                    (:vp, :vs, :density, :vph, :vpv, :vsh, :vsv, :eta, :Qμ,  :Qκ))
                    @test f(m, 0.9) == evaluate(m, sym, 0.9)
                end
            end
        end

        @testset "PREMPolyModels" begin
            # Polynomial evaluation
            @test SeisModels._evalpoly(0.25, hcat([0, 0, 0, 0], [1, 2, 3, 4]), 2) == 1.75
            @test SeisModels._evalpoly(0.25, hcat([1, 2, 3, 4], [0, 0, 0, 0]), 1) == 1.75
            # Velocities default to reference frequency
            for f in (vp, vs, vpv, vsv, vph, vsh)
                @test f(PREM, 5000) == f(PREM, 5000, freq=reffrequency(PREM))
            end
            # Value from separate implementation at https://github.com/andreww/prem4derg
            @test vs(PREM, 1000, freq=0.1) == 3.5274008549889566
            @test mass(PREM) == mass(PREM, surface_radius(PREM))
        end

        # TODO: Add tests for LinearLayeredModel and PREMPolyModel
    end

    # Keyword constructors
    @testset "Construction" begin
        let n = 5, r = [0; sort(10_000rand(n-1))], Vp = rand(n), Vs = 0.5rand(n),
                Rho = rand(n), Vpv = 1.01.*Vp, Vph = 0.99.*Vp,
                Vsv = 0.99.*Vs, Vsh = 1.01.*Vs, Eta = 1 .+ 0.1rand(n),
                qμ = rand(n), qκ = rand(n), R = last(r), fref = rand(),
                null = Array{Float64}(undef, 0, 0)
            for model_type in (LinearLayeredModel, SteppedLayeredModel)
                # Arrays wrong length
                @test_throws ArgumentError model_type(r=r, vp=[Vp; 1])
                # r must be provided
                @test_throws UndefKeywordError model_type(vp=Vp, vs=Vs)
                # r must increase
                @test_throws ArgumentError model_type(r=[3, 2, 1])
                # Equivalence of default and keyword constructors
                @test model_type(R, n, r, Vp, Vs, Rho, true, Vph, Vpv,
                    Vsh, Vsv, Eta, true, qμ, qκ) ==
                    model_type(r=r, vp=Vp, vs=Vs, density=Rho, vph=Vph,
                        vpv=Vpv, vsh=Vsh, vsv=Vsv, Qμ=qμ, Qκ=qκ, eta=Eta)
                @test model_type(R, n, r, Vp, Vs, [], false, [],
                                         [], [], [], [], true, qμ, qκ) ==
                    model_type(r=r, vp=Vp, vs=Vs, Qμ=qμ, Qκ=qκ)
                @test model_type(R, n, r, Vp, Vs, Rho, true,
                                         Vph, Vpv, Vsh, Vsv, Eta, false, [], []) ==
                    model_type(r=r, vp=Vp, vs=Vs, density=Rho,
                                       vph=Vph, vpv=Vpv, vsh=Vsh, vsv=Vsv, eta=Eta)
                # ASCII arguments are the same and take precedence
                @test model_type(r=r, Qmu=qμ, Qkappa=qκ) ==
                    model_type(r=r, Qμ=qμ, Qκ=qκ) ==
                    model_type(r=r, Qμ=Vs, Qκ=Vp, Qmu=qμ, Qkappa=qκ)
            end
            ## Model-specific tests
            # LinearLayeredModel
            # r must start at 0
            @test_throws ArgumentError LinearLayeredModel(r=r[2:end])

            # PREMPolyModel
            @test_throws ArgumentError PREMPolyModel(r=r, vp=[Vp; 1])
            @test_throws UndefKeywordError PREMPolyModel(vp=Vp, vs=Vs)
            # r must increase
            @test_throws ArgumentError PREMPolyModel(r=[3, 2, 1])
            # Equivalence of default and keyword constructors
            @test PREMPolyModel(R, n, r, Vp', Vs', Rho', true, Vph', Vpv',
                Vsh', Vsv', Eta', true, qμ', qκ', fref) ==
                PREMPolyModel(r=r, vp=Vp, vs=Vs, density=Rho, vph=Vph,
                    vpv=Vpv, vsh=Vsh, vsv=Vsv, Qμ=qμ, Qκ=qκ, eta=Eta, fref=fref)
            @test PREMPolyModel(R, n, r, Vp', Vs', null, false, null,
                                     null, null, null, null, true, qμ', qκ', NaN) ==
                PREMPolyModel(r=r, vp=Vp, vs=Vs, Qμ=qμ, Qκ=qκ)
            @test PREMPolyModel(R, n, r, Vp', Vs', Rho', true, Vph', Vpv',
                                Vsh', Vsv', Eta', false, null, null, fref) ==
                PREMPolyModel(r=r, vp=Vp, vs=Vs, density=Rho, vph=Vph,
                              vpv=Vpv, vsh=Vsh, vsv=Vsv, eta=Eta, fref=fref)
            # ASCII arguments are the same and take precedence
            @test PREMPolyModel(r=r, Qmu=qμ, Qkappa=qκ) ==
                PREMPolyModel(r=r, Qμ=qμ, Qκ=qκ) ==
                PREMPolyModel(r=r, Qμ=Vs, Qκ=Vp, Qmu=qμ, Qkappa=qκ)
            # Arrays are converted to the correct shape
            @test_throws ArgumentError PREMPolyModel(r=r, vp=rand(3, n+1))
            Vp′ = rand(3, n)
            @test PREMPolyModel(r=r, vp=Vp′) == PREMPolyModel(R, n, r, Vp′,
                null, null, false, null, null, null, null, null, false, null, null, NaN)
        end
    end

    # Test for equivalence of depth versions of functions
    @testset "Depth" begin
        let n = 5, radii = sort(rand(n)), R = radii[end],
                (Vp, Vs, rho, Vpv, Vph, Vsv, Vsh, Eta, Qkappa, Qmu) = rand.(repeat([n], 10))
            for model_type in (SteppedLayeredModel, LinearLayeredModel, PREMPolyModel)
                m = if model_type <: PREMPolyModel
                    model_type(R, n, radii, Vp', Vs', rho', true, Vpv', Vph', Vsv', Vsh', Eta',
                        true, Qkappa', Qmu', NaN)
                else
                    model_type(R, n, radii, Vp, Vs, rho, true, Vpv, Vph, Vsv, Vsh, Eta,
                        true, Qkappa, Qmu)
                end
                r = R*0.99*rand()
                z = R - r
                for func in (vp, vs, density, vpv, vph, vsv, vsh, eta, Qκ, Qμ,
                        pressure, bulk_modulus, shear_modulus,
                        youngs_modulus, poissons_ratio, mass, surface_mass, gravity)
                    @test func(m, z, depth=true) ≈ func(m, r)
                end
                @test evaluate(m, :vp, z, depth=true) ≈ evaluate(m, :vp, r)
            end
        end
    end
end


text:
"""
$(TYPEDEF)

Matrix epigraph of perspective function of real or complex matrix outer product
for a matrix (stacked column-wise) of `nrows` rows and `ncols` columns with
`nrows ≤ ncols`.

    $(FUNCTIONNAME){T, R}(nrows::Int, ncols::Int, use_dual::Bool = false)
"""
mutable struct MatrixEpiPerSquare{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d1::Int
    d2::Int
    is_complex::Bool
    rt2::T

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    hess_aux_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    U_idxs::UnitRange{Int}
    v_idx::Int
    W_idxs::UnitRange{Int}
    U::Hermitian{R, Matrix{R}}
    W::Matrix{R}
    Z::Hermitian{R, Matrix{R}}
    fact_Z
    Zi::Hermitian{R, Matrix{R}}
    ZiW::Matrix{R}
    ZiUZi::Hermitian{R, Matrix{R}}
    WtZiW::Hermitian{R, Matrix{R}}
    tempd2d2::Matrix{R}
    tempd1d1::Matrix{R}
    tempd1d1b::Matrix{R}
    tempd1d1c::Matrix{R}
    tempd1d1d::Matrix{R}
    tempd1d2::Matrix{R}
    ZiUZiW::Matrix{R}
    Hvv::T

    function MatrixEpiPerSquare{T, R}(
        d1::Int,
        d2::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert 1 <= d1 <= d2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.is_complex = (R <: Complex)
        cone.v_idx = 1 + svec_length(R, d1)
        cone.dim = cone.v_idx + vec_length(R, d1 * d2)
        cone.d1 = d1
        cone.d2 = d2
        cone.rt2 = sqrt(T(2))
        cone.U_idxs = 1:(cone.v_idx - 1)
        cone.W_idxs = (cone.v_idx + 1):cone.dim
        return cone
    end
end

reset_data(cone::MatrixEpiPerSquare) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated =
    cone.hess_aux_updated = false)

function setup_extra_data!(
    cone::MatrixEpiPerSquare{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    (d1, d2) = (cone.d1, cone.d2)
    cone.U = Hermitian(zeros(R, d1, d1), :U)
    cone.W = zeros(R, d1, d2)
    cone.Z = Hermitian(zeros(R, d1, d1), :U)
    cone.Zi = Hermitian(zeros(R, d1, d1), :U)
    cone.ZiW = zeros(R, d1, d2)
    cone.ZiUZi = Hermitian(zeros(R, d1, d1), :U)
    cone.WtZiW = Hermitian(zeros(R, d2, d2), :U)
    cone.tempd2d2 = zeros(R, d2, d2)
    cone.tempd1d1 = zeros(R, d1, d1)
    cone.tempd1d1b = zeros(R, d1, d1)
    cone.tempd1d1c = zeros(R, d1, d1)
    cone.tempd1d1d = zeros(R, d1, d1)
    cone.tempd1d2 = zeros(R, d1, d2)
    cone.ZiUZiW = zeros(R, d1, d2)
    return cone
end

get_nu(cone::MatrixEpiPerSquare) = cone.d1 + 1

function set_initial_point!(
    arr::AbstractVector,
    cone::MatrixEpiPerSquare{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    incr = (cone.is_complex ? 2 : 1)
    arr .= 0
    k = 1
    @inbounds for i in 1:cone.d1
        arr[k] = 1
        k += incr * i + 1
    end
    arr[cone.v_idx] = 1
    return arr
end

function update_feas(cone::MatrixEpiPerSquare{T}) where T
    @assert !cone.feas_updated
    v = cone.point[cone.v_idx]

    if v > eps(T)
        @views U = svec_to_smat!(cone.U.data, cone.point[cone.U_idxs], cone.rt2)
        @views W = vec_copyto!(cone.W, cone.point[cone.W_idxs])
        copyto!(cone.Z.data, U)
        mul!(cone.Z.data, cone.W, cone.W', -1, 2 * v)
        cone.fact_Z = cholesky!(cone.Z, check = false)
        cone.is_feas = isposdef(cone.fact_Z)
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::MatrixEpiPerSquare{T}) where T
    v = cone.dual_point[cone.v_idx]
    if v > eps(T)
        @views svec_to_smat!(cone.tempd1d1b, cone.dual_point[cone.U_idxs],
            cone.rt2)
        F = cholesky!(Hermitian(cone.tempd1d1b, :U), check = false)
        isposdef(F) || return false
        @views W = vec_copyto!(cone.tempd1d2, cone.dual_point[cone.W_idxs])
        LW = ldiv!(F.U', W)
        trLW = sum(abs2, LW)
        return (2 * v - trLW > eps(T))
    end
    return false
end

function update_grad(cone::MatrixEpiPerSquare)
    @assert cone.is_feas
    U = cone.U
    W = cone.W
    dim = cone.dim
    v = cone.point[cone.v_idx]

    Zi = cone.Zi
    inv_fact!(Zi.data, cone.fact_Z)
    @views smat_to_svec!(cone.grad[cone.U_idxs], Zi, cone.rt2)
    @views cone.grad[cone.U_idxs] .*= -2 * v
    cone.grad[cone.v_idx] = -2 * dot(Zi, U) + (cone.d1 - 1) / v
    ldiv!(cone.ZiW, cone.fact_Z, W)
    @views vec_copyto!(cone.grad[cone.W_idxs], cone.ZiW)
    @. @views cone.grad[cone.W_idxs] *= 2

    cone.grad_updated = true
    return cone.grad
end

function update_hess_aux(cone::MatrixEpiPerSquare)
    @assert cone.grad_updated
    ZiUZi = cone.ZiUZi

    ldiv!(ZiUZi.data, cone.fact_Z, cone.U)
    rdiv!(ZiUZi.data, cone.fact_Z)
    mul!(cone.WtZiW.data, cone.W', cone.ZiW) # TODO not used for hess prod
    v = cone.point[cone.v_idx]
    cone.Hvv = 4 * dot(ZiUZi, cone.U) - (cone.d1 - 1) / v / v
    mul!(cone.ZiUZiW, cone.ZiUZi, cone.W)

    cone.hess_aux_updated = true
    return cone.hess_aux_updated
end

function update_hess(cone::MatrixEpiPerSquare)
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    d1 = cone.d1
    d2 = cone.d2
    U_idxs = cone.U_idxs
    v_idx = cone.v_idx
    W_idxs = cone.W_idxs
    v = cone.point[v_idx]
    tempd2d2 = cone.tempd2d2
    ZiUZi = cone.ZiUZi
    tempd1d1 = cone.tempd1d1
    Zi = cone.Zi
    ZiW = cone.ZiW
    idx_incr = (cone.is_complex ? 2 : 1)

    # H_W_W part
    copyto!(tempd2d2, I)
    tempd2d2 .+= cone.WtZiW

    # TODO parallelize loops
    r_idx = v_idx + 1
    for i in 1:d2, j in 1:d1
        c_idx = r_idx
        @inbounds for k in i:d2
            ZiWjk = ZiW[j, k]
            tempd2d2ik = tempd2d2[i, k]
            lstart = (i == k ? j : 1)
            @inbounds for l in lstart:d1
                term1 = Zi[l, j] * tempd2d2ik
                term2 = ZiW[l, i] * ZiWjk
                spectral_kron_element!(H, r_idx, c_idx, term1, term2)
                c_idx += idx_incr
            end
        end
        r_idx += idx_incr
    end
    @views H[W_idxs, W_idxs] .*= 2

    # H_U_U part
    @views H_U_U = H[U_idxs, U_idxs]
    symm_kron!(H_U_U, Zi, cone.rt2)
    @. H_U_U *= 4 * abs2(v)

    # H_v_v part
    @views H[v_idx, v_idx] = cone.Hvv

    # H_U_W part
    # TODO parallelize loops
    # TODO use dispatch for complex part and clean up
    row_idx = 1
    for i in 1:d1, j in 1:i # U lower tri idxs
        col_idx = v_idx + 1
        for l in 1:d2, k in 1:d1 # W idxs
            @inbounds if cone.is_complex
                term1 = Zi[k, i] * ZiW[j, l]
                term2 = Zi[k, j] * ZiW[i, l]
                if i != j
                    term1 *= cone.rt2
                    term2 *= cone.rt2
                end
                H[row_idx, col_idx] = real(term1) + real(term2)
                H[row_idx, col_idx + 1] = imag(term1) + imag(term2)
                if i != j
                    H[row_idx + 1, col_idx] = imag(term2) - imag(term1)
                    H[row_idx + 1, col_idx + 1] = real(term1) - real(term2)
                end
            else
                term = Zi[i, k] * ZiW[j, l] + Zi[k, j] * ZiW[i, l]
                if i != j
                    term *= cone.rt2
                end
                H[row_idx, col_idx] = term
            end
            col_idx += idx_incr
        end
        row_idx += (i == j ? 1 : idx_incr)
    end
    @. @views H[U_idxs, W_idxs] *= -2 * v

    # H_v_W part
    @views vec_copyto!(H[v_idx, W_idxs], cone.ZiUZiW)
    @. @views H[v_idx, W_idxs] *= -4

    # H_U_v part
    copyto!(tempd1d1, ZiUZi)
    axpby!(-2, Zi, 4 * v, tempd1d1)
    @views smat_to_svec!(H[U_idxs, v_idx], tempd1d1, cone.rt2)

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::MatrixEpiPerSquare,
    )
    cone.hess_aux_updated || update_hess_aux(cone)
    U_idxs = cone.U_idxs
    v_idx = cone.v_idx
    W_idxs = cone.W_idxs
    v2 = 2 * cone.point[cone.v_idx]
    ZiUZi = cone.ZiUZi
    tempd1d1 = Hermitian(cone.tempd1d1, :U)
    temp_U = Hermitian(cone.tempd1d1b, :U)
    temp_U2 = Hermitian(cone.tempd1d1c, :U)
    temp_U3 = cone.tempd1d1d
    temp_W = cone.tempd1d2

    @inbounds for i in 1:size(arr, 2)
        @views svec_to_smat!(temp_U.data, arr[U_idxs, i], cone.rt2)
        @views vec_copyto!(temp_W, arr[W_idxs, i])
        v_arr = arr[v_idx, i]
        @views U_prod = prod[U_idxs, i]
        @views W_prod = prod[W_idxs, i]

        ldiv!(cone.fact_Z, temp_W)
        mul!(temp_U3, temp_W, cone.ZiW')
        copyto!(tempd1d1.data, temp_U)
        rdiv!(ldiv!(cone.fact_Z, tempd1d1.data), cone.fact_Z)
        @. temp_U2.data = temp_U3 + temp_U3' - v2 * tempd1d1.data

        copyto!(tempd1d1, temp_U2)
        axpy!(-2 * v_arr, ZiUZi.data, tempd1d1.data)
        vec_copyto!(W_prod, mul!(temp_W, tempd1d1, cone.W, 2, 2))

        copyto!(tempd1d1, ZiUZi)
        axpby!(-2, cone.Zi.data, 2 * v2, tempd1d1.data)
        prod[v_idx, i] = real(dot(tempd1d1, temp_U)) - 4 *
            real(dot(cone.U, temp_U3)) + cone.Hvv * v_arr
        axpby!(-v2, temp_U2.data, v_arr, tempd1d1.data)
        smat_to_svec!(U_prod, tempd1d1, cone.rt2)
    end

    return prod
end

# TODO reduce allocs
function dder3(cone::MatrixEpiPerSquare, dir::AbstractVector)
    cone.hess_aux_updated || update_hess_aux(cone)
    d1 = cone.d1
    d2 = cone.d2
    dim = cone.dim
    fact_Z = cone.fact_Z

    U_idxs = cone.U_idxs
    v_idx = cone.v_idx
    W_idxs = cone.W_idxs
    U = cone.U
    v = cone.point[cone.v_idx]
    W = cone.W

    @views U_dir = Hermitian(svec_to_smat!(zero(U.data), dir[U_idxs], cone.rt2))
    v_dir = dir[v_idx]
    @views W_dir = vec_copyto!(zero(W), dir[W_idxs])

    dder3 = cone.dder3
    U_dder3 = view(dder3, U_idxs)
    W_dder3 = view(dder3, W_idxs)

    v2 = 2 * v
    vd2 = 2 * v_dir
    Zi = cone.Zi
    ZiW = cone.ZiW # TODO AKA tau
    ZiU = fact_Z \ U # TODO cache
    ZiWd = fact_Z \ W_dir
    ZiUd = fact_Z \ U_dir
    ZiUZi = cone.ZiUZi
    ZiUZiUZi = Hermitian(ZiUZi * ZiU', :U)
    ZiUZi2v = Hermitian(ZiUZi - v2 * ZiUZiUZi, :U)
    WdWZi = W_dir * ZiW'
    WdZiW = W_dir' * ZiW
    UdZiW = U_dir * ZiW
    ZiWdWZi = fact_Z \ WdWZi
    ZiWdWZi2 = Hermitian(ZiWdWZi + ZiWdWZi', :U)
    ZiUdZiW = fact_Z \ UdZiW
    ZiUZiWdWZi = ZiU * ZiWdWZi
    ZiUZiUdZiW = ZiU * ZiUdZiW + ZiUd * cone.ZiUZiW
    ZiWdWZiUZi = ZiWdWZi * ZiU'
    ZiWdWZiUZi2 = ZiWdWZiUZi + ZiWdWZiUZi' + ZiUZiWdWZi'
    ZiUdZi = Hermitian(ZiUd / fact_Z, :U)
    ZiUZiUdZi = ZiU * ZiUdZi
    ZiUZiUdZi2 = Hermitian(ZiUZiUdZi + ZiUZiUdZi')
    ZiUdZiWdWZi = ZiUd * ZiWdWZi + ZiWdWZi * ZiUd'
    WtZiWI = cone.WtZiW + I
    ZiWdWtZiWI = ZiWd * WtZiWI
    vdZiUZiUZiW = vd2 * ZiUZiUZi * W
    WdWtZiWI = W_dir * WtZiWI
    ZiUZiWdWZiWI = ZiUZi * WdWtZiWI + ZiWdWZiUZi2 * W
    vZiUZiUdZi2 = v * ZiUZiUdZi2 - ZiUdZi

    Utemp = vd2 * (-vd2 * ZiUZi2v + ZiWdWZi2 - v2 * (ZiUZiWdWZi + ZiWdWZiUZi2 -
        2 * vZiUZiUdZi2)) + v2 * (ZiWdWZi * WdWZi + WdWZi' * ZiWdWZi2 +
        ZiWdWtZiWI * ZiWd' + v2 * (v2 * ZiUd * ZiUdZi - ZiUdZiWdWZi - ZiUdZiWdWZi'))
    smat_to_svec!(U_dder3, Utemp, cone.rt2)

    v_Wd_dot = -4 * (v * ZiUZiUdZiW + vdZiUZiUZiW) + ZiUZiWdWZiWI + 2 * ZiUdZiW
    dder3[v_idx] = v_dir * (-8 * dot(ZiUZi2v, U_dir) + v_dir *
        (8 * real(dot(ZiUZiUZi, U)) - (d1 - 1) / v / v / v)) +
        4 * v * real(dot(vZiUZiUdZi2, U_dir)) + 2 * real(dot(v_Wd_dot, W_dir))

    Wtemp = 4 * v_dir * (ZiUdZiW - v2 * ZiUZiUdZiW + ZiUZiWdWZiWI - vdZiUZiUZiW) +
        4 * v * (ZiUdZiW * WdZiW + ZiWdWZi * UdZiW + WdWZi' * ZiUdZiW +
        ZiUdZi * WdWtZiWI - v2 * ZiUd * ZiUdZiW) + -2 * (ZiW * WdZiW * WdZiW +
        WdWZi' * ZiWdWtZiWI + ZiWdWtZiWI * WdZiW + ZiWd * WdZiW' * WtZiWI)
    vec_copyto!(W_dder3, Wtemp)

    return dder3
end

# compute an element of the real spectral Kronecker in-place
function spectral_kron_element!(
    skr::AbstractMatrix{T},
    i::Int,
    j::Int,
    a::T,
    b::T,
    ) where {T <: Real}
    @inbounds skr[i, j] = a + b
    return skr
end

# compute an element of the complex spectral Kronecker in-place
function spectral_kron_element!(
    skr::AbstractMatrix{T},
    i::Int,
    j::Int,
    a::Complex{T},
    b::Complex{T},
    ) where {T <: Real}
    apb = a + b
    amb = a - b
    @inbounds begin
        skr[i, j] = real(apb)
        skr[i + 1, j] = -imag(amb)
        skr[i, j + 1] = imag(apb)
        skr[i + 1, j + 1] = real(amb)
    end
    return skr
end


text:
include("run_mms.jl")

let
    N0 = 17
    for sbp_order in (2, 4, 6)
        @show sbp_order
        main(
            sbp_order,
            4,
            N0;
            characteristic_method = true,
            cfl = 1 / 2,
            friction = (V) -> 128 * asinh(V),
            tspan = (0.0, 1.0),
            do_output = true,
        )
    end
end


text:
function clean(phone_number)
    justdigits = filter(isdigit, phone_number)
    oneremoved = justdigits[1] == '1' ? justdigits[2:end] : justdigits
    (occursin(r"[2-9]..[2-9]......",  oneremoved) && length(oneremoved)) == 10 ? oneremoved : nothing
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
using ParallelAccelerator

@acc function mult()
    w = 5#1000
    h = 5#1000
    A = fill(1,h,w)
    B = fill(1,h,w)
    C = fill(0,h,w)
    C += A*B
    return C
end

function main()
    tic()
    mult()
    #D = mult()
    println("SELFTIMED ", toq())

    #println(D)

end

main()

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
using QuantileRegression
using DataFrames
cd(Pkg.dir("QuantileRegression")*"/examples")
df = readtable("engel.csv")

out_ip = qreg(foodexp~income, df, IP()) # or just qreg(foodexp~income, df)
out_irls = qreg(foodexp~income, df, IRLS())

println("...... Interior point ......")
println(out_ip)
println(".......... IRLS ............")
println(out_irls)


text:
module TestParseRelations
using Test

using Catlab.WiringDiagrams.UndirectedWiringDiagrams
using Catlab.Programs.ParseRelations
using Catlab.Programs.ParseRelations: RelationDiagram, TypedRelationDiagram

# Untyped relations

parsed = @relation (x,z) where (x,y,z) begin
  R(x,y)
  S(y,z)
end

d = RelationDiagram(2)
add_box!(d, 2, name=:R); add_box!(d, 2, name=:S)
add_junctions!(d, 3, variable=[:x,:y,:z])
set_junction!(d, [1,2,2,3])
set_junction!(d, [1,3], outer=true)
@test parsed == d

parsed = @relation ((x,z) where (x,y,z)) -> (R(x,y); S(y,z))
@test parsed == d

parsed = @relation function (x,z) where (x,y,z); R(x,y); S(y,z) end
@test parsed == d

# Typed relations

parsed = @relation (x,y,z) where (x::X, y::Y, z::Z, w::W) begin
  R(x,w)
  S(y,w)
  T(z,w)
end
d = TypedRelationDiagram([:X,:Y,:Z])
add_box!(d, [:X,:W], name=:R)
add_box!(d, [:Y,:W], name=:S)
add_box!(d, [:Z,:W], name=:T)
add_junctions!(d, [:X,:Y,:Z,:W], variable=[:x,:y,:z,:w])
set_junction!(d, [1,4,2,4,3,4])
set_junction!(d, [1,2,3], outer=true)
@test parsed == d

end


text:
# numerical range
# This test is based on S&P chap. 43
using Pseudospectra
using Test

@testset "Num. range" begin
    η = 0.015
    A = Pseudospectra.advdiff(100,η)
    opts = Dict{Symbol,Any}(:showfov => true, :fov_npts => 200,
                            :ax => [-60,20,-50,50])

    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)

    @test isapprox(maximum(real(ps_data.ps_dict[:fov])),-η*π^2,atol=1e-6)
end