text:
using ModelingToolkit
using OrdinaryDiffEq
using SteadyStateDiffEq
using SciMLBase

"""
Find steady states of the ODE system `sys` by randomly
"""
function find_steady_states(;
    S=ZERO_SIGNAL, proteins=STRESSED, params=Dict(),
    trajectories=10, batch_size=trajectories, ntarget=trajectories,
    solver=SSRootfind(), ensemble_method=EnsembleThreads()
)

    @named sys = RtgMTK(S; proteinlevels=proteins)

    statesmap = Dict(k => i for (i, k) in enumerate(states(sys)))

    # Create a set of initial conditions respecting conservation relationships for the ensemble
    function make_u0()
        u0 = zeros(length(statesmap))

        # Randomly split total into `n` fractions
        function _frac(n::Int, total=1.0)
            total .* diff([zero(total); sort([rand(typeof(total)) for _ in 1:(n-1)]); one(total)])
        end

        # Split ΣMks into 3
        frac = _frac(3, proteins[ΣMks])
        # u0[statesmap[Mks]] = frac[1]
        u0[statesmap[BmhMks]] = frac[2]
        u0[statesmap[Rtg2Mks_c]] = frac[3]

        # Remainder of Bmh
        # u0[statesmap[Bmh]] = proteins[ΣBmh] - u0[statesmap[BmhMks]]

        # Split ΣRtg2 into 2
        frac = _frac(2, proteins[ΣRtg2] - u0[statesmap[Rtg2Mks_c]])
        # u0[statesmap[Rtg2I_c]] = frac[1]
        u0[statesmap[Rtg2A_c]] = frac[2]

        # Split ΣRtg1 into 6
        frac = _frac(6, proteins[ΣRtg1])
        u0[statesmap[Rtg13A_c]] = frac[1]
        u0[statesmap[Rtg13I_c]] = frac[2]
        # u0[statesmap[Rtg1_c]] = frac[3]
        u0[statesmap[Rtg1_n]] = frac[4]
        u0[statesmap[Rtg13A_n]] = frac[5]
        u0[statesmap[Rtg13I_n]] = frac[6]

        # Split ΣRtg3 not in Rtg13 complex into 4
        frac = _frac(4, proteins[ΣRtg3] - u0[statesmap[Rtg13A_c]] - u0[statesmap[Rtg13I_c]] - u0[statesmap[Rtg13A_n]] - u0[statesmap[Rtg13I_n]])

        # u0[statesmap[Rtg3I_c]] = frac[1]
        u0[statesmap[Rtg3A_c]] = frac[2]
        u0[statesmap[Rtg3A_n]] = frac[3]
        u0[statesmap[Rtg3I_n]] = frac[4]

        return u0
    end

    """Rerun when the solution is invalid (NaNs and negatives)"""
    function output_func(sol, i; negtol=-1e-6)
        rerun = any(isnan.(sol.u)) || any(sol.u .< negtol)
        clamp!(sol.u, 0.0, Inf)
        return (sol, rerun)
    end

    function prob_func(prob, i, repeat)
        # Roll a random initial condition
        remake(prob, u0=make_u0())
    end

    """Only save unique solutions"""
    function reduction(us, batch, I; norm_atol=1e-4)
        for res in batch
            isUnique = true
            for existing in us
                if isapprox(existing.u, res.u, atol=norm_atol)
                    isUnique = false
                    break
                end
            end

            if isUnique
                us = push!(us, res)
            end
        end
        return (us, length(us) >= ntarget)
    end

    prob = SteadyStateProblem(sys, resting_u0(sys), params)
    ensprob = EnsembleProblem(prob; output_func, prob_func, reduction)
    return solve(ensprob, solver, ensemble_method; trajectories, batch_size)
end


text:
include("operators.jl")

using GaussQuadrature
using OrdinaryDiffEq
using LinearAlgebra
using UnicodePlots

struct SpectralOP{T}
  x  :: Array{T, 1}
  ω  :: Array{T, 1}
  D  :: Array{T, 2}
  e0 :: Array{T, 1}
  eN :: Array{T, 1}
  function SpectralOP{T}(N) where T
    Nq = N + 1
    x, ω = legendre(T, Nq, both)
    x = (x .+ 1) / 2
    ω /= 2
    D = derivative(x)
    e0 = [1; zeros(T, N)]
    eN = [zeros(T, N); 1]
    new{T}(x, ω, D, e0, eN)
  end
end

function spectralodefun!(dq, q, p, t)
  o = p.o
  f = p.f
  Nq = length(o.x)
  v,  u  = @views  q[1:Nq],  q[Nq+1:2Nq]
  dv, du = @views dq[1:Nq], dq[Nq+1:2Nq]

  D = o.D
  ω = o.ω
  e0 = o.e0
  dv .= (-D' * (ω .* D * u) - e0 * f(v[1])) ./ ω
  du .= v
end

let
  N = 100
  T = Float64

  tspan = (T(0), T(3//10))

  Nq = N + 1
  o = SpectralOP{T}(N)
  q = [zeros(T, Nq);sin.(2o.x * π)]
  f(v) = asinh(v)

  p = (o=o, f=f)
  prob = ODEProblem(spectralodefun!, q, tspan, p)
  sol = solve(prob, Tsit5();
              saveat=range(tspan[1], stop=tspan[2], length = 4),
              atol = 1e-6, rtol = 1e-3,
              internalnorm=(x, _)->norm(x, Inf))

  for n = 1:length(sol.t)
    display(lineplot(o.x, sol.u[n][Nq+1:2Nq], width=100, xlim=(0,1),
                     title="t = $(sol.t[n])"))
  end
end

ext: jl


text:
function hessfunc{T <: Number}(
    f::Function,
    x::T;
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function h(x::T)
            return finitehess(f, x, direction)
        end
    else
        throw(ArgumentError("method must be :finite"))
    end
    return h
end

function hessfunc!{T <: Number}(
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function h!(H::Matrix{T}, x::Vector{T})
            return finitehess!(H, f, x, direction)
        end
    else
        throw(ArgumentError("method must be :finite"))
    end
    return h!
end

function hessfunc{T <: Number}(
    f::Function,
    x::Vector{T};
    method::Symbol = :finite,
    direction::Symbol = :central,
) # -> Function
    if method == :finite
        function h(x::Vector{T})
            n = length(x)
            H = Array(T, n, n)
            finitehess!(H, f, x, direction)
            return H
        end
    else
        throw(ArgumentError("method must be :finite"))
    end
    return h
end


text:
#=
Contains functions/type prototypes for the mesh.
=#

#= Mesh =#
#Subtypes should contain all variables that define a mesh.
#Should make an effort to follow the naming conventions established
#in existing methods in order to give compatibility between methods
abstract type AbstractMesh end
mutable struct Mesh{T} <: AbstractMesh
    G::Matrix{T}
    D::Matrix{T}
    l::Int
    Δᵐ::T
    Δᵖ::T

    # Override constructor for different default meshes for
    # different poll techniques.
    Mesh(N::Int64) = Mesh{Float64}(N)
    function Mesh{T}(N::Int64) where T
        mesh = new()
        mesh.l = 0
        mesh.Δᵐ = min(1, 4.0^(-mesh.l))
        mesh.Δᵖ = 2.0^(-mesh.l)
        mesh.G = Matrix(I,N,N)
        mesh.D = hcat(Matrix(I,N,N),-Matrix(I,N,N))
        return mesh
    end
end

"""
    MeshUpdate!(mesh::Mesh, improvement_found::Bool)

Implements LTMADS update rule from Audet & Dennis 2006 pg. 203 adapted for progressive
barrier constraints with Audet & Dennis 2009 expression 2.4.
"""
function MeshUpdate!(m::Mesh, ::AbstractPoll, result::IterationOutcome)
    if result == Unsuccessful
        m.l += 1
    elseif result == Dominating && m.l > 0
        m.l -= 1
    elseif result == Improving
        m.l == m.l
    end

    m.Δᵐ = min(1, 4.0^(-m.l))
    m.Δᵖ = 2.0^(-m.l)
end


text:
using Etcetera
using Test

@testset "LongestRepeatedStrings.jl" begin

    @test longest_repeated_substrings([9,1,2,3,5,1,2,9]) == (2, Dict([1,2] => [2, 6]))

    @test longest_repeated_substrings("ananas") == (3, Dict("ana" => [1, 3]))
    @test longest_repeated_substrings("ananas"; overlap_allowed=false) == (2, Dict("an" => [1, 3], "na" => [2, 4]))

    @test longest_repeated_substrings("Mississippi") == (4, Dict("issi" => [2, 5]))
    @test longest_repeated_substrings("Mississippi"; overlap_allowed=false) == (3, Dict("iss" => [2, 5], "ssi" => [3, 6]))

     #=
        only repeating strings longer than 1 are reported.

    =#
    @test longest_repeated_substrings("") == (nothing, Dict())
    @test longest_repeated_substrings("a") == (nothing, Dict())
    @test longest_repeated_substrings("aa") == (nothing, Dict())

    @test longest_repeated_substrings("aaa") == (2, Dict( "aa" => [1,2]))
    @test longest_repeated_substrings("aaa"; overlap_allowed=false) == (nothing, Dict())

    @test longest_repeated_substrings("aaaa") == (3, Dict("aaa" =>[1,2]))
    @test longest_repeated_substrings("aaaa",overlap_allowed=false) == (2, Dict("aa" =>[1,3]))

    @test longest_repeated_substrings("aaaaa",overlap_allowed=false) == (2, Dict("aa" =>[1,3]))

    @test longest_repeated_substrings([1, 2, 3, 1, 2, 3]) == (3, Dict([1, 2, 3] => [1, 4]))
    @test longest_repeated_substrings(["one",2,(3,),"one",2,(3,)] ) == (3, Dict(["one", 2, (3,)] => [1, 4]))

end


text:
using LoopVectorization, LinearAlgebra, Base.Threads
import Base.Threads.@spawn

struct InputBlock{T<:Real,
                  M<:AbstractMatrix{T},
                  TF<:Real,
                  TT<:Real,
                  TKD<:Real} <: AbstractMatrix{T}
    data::M
    # Frequency of channel 1 and N
    f_ch1::TF
    f_chn::TF
    # Time step in seconds
    t_samp::TT
    # Dispersion delay across the bands covered by the input block
    Δkdisp::TKD
end

struct OutputBlock{T<:Real,M<:Matrix{T}} <: AbstractMatrix{T}
    data::M
    # Dispersion delay (s) for first DM trial
    y_min::Int
    # Dispersion delay (s) for last DM trial
    y_max::Int
end

# Methods for custom array types
Base.size(ib::InputBlock) = size(ib.data)
Base.size(ob::OutputBlock) = size(ob.data)
Base.getindex(ib::InputBlock, i::Int, j::Int) = ib.data[i, j]
Base.getindex(ob::OutputBlock, i::Int, j::Int) = ob.data[i, j]

# Default constructors
function InputBlock(data, f_ch1, f_chn, t_samp)
    @assert f_ch1 >= f_chn "Channel 1 must be higher than Channel N"
    Δkdisp = KDM * (f_chn^-2 - f_ch1^-2)
    return InputBlock(data, f_ch1, f_chn, t_samp, Δkdisp)
end

function OutputBlock(ib::InputBlock{T}, y_min, y_max) where {T}
    n_trials = y_max - y_min + 1
    n_samp = size(ib)[1]
    data = zeros(T, n_samp, n_trials)
    return OutputBlock(data, y_min, y_max)
end

# Channel spacing
f_step(ib::InputBlock) = (ib.f_ch1 - ib.f_chn) / size(ib)[2]

# Natural DM step of FastDMTransform for this block
dm_step(ib::InputBlock) = ib.t_samp / ib.Δkdisp

function split(ib::InputBlock)
    # The split occurs at the channel where half the dispersion
    # delay has been accounted for
    f = (0.5 * ib.f_chn^-2 + 0.5 * ib.f_ch1^-2)^-0.5
    i = round(Int, (ib.f_ch1 - f) / f_step(ib))

    # New frequency ranges
    step = f_step(ib)
    head_start, head_end = promote(ib.f_ch1, ib.f_ch1 - i * step)
    tail_start, tail_end = promote(ib.f_ch1 - (i + 1) * step, ib.f_chn)

    # New data slices
    ch = size(ib)[2]
    head_slice = view(ib.data, :, 1:i)
    tail_slice = view(ib.data, :, (i + 1):ch)

    # New blocks
    head = InputBlock(head_slice, head_start, head_end, ib.t_samp)
    tail = InputBlock(tail_slice, tail_start, tail_end, ib.t_samp)

    return head, tail
end

circmod(x, y) = mod(x - 1, y) + 1

# We really want go get rid of all these allocations, it's slowing us down a lot
function transform_recursive(block::InputBlock, y_min::Int, y_max::Int)
    out = OutputBlock(block, y_min, y_max)
    # Base Case
    if size(block)[2] == 1
        out.data[:, 1] = block.data[:, 1]
        return out
    end
    # Split
    head, tail = split(block)
    # Transform
    y_min_head = round(Int, y_min * head.Δkdisp / block.Δkdisp + 0.5)
    y_max_head = round(Int, y_max * head.Δkdisp / block.Δkdisp + 0.5)
    transformed_head_task = @spawn transform_recursive(head, y_min_head, y_max_head)

    y_min_tail = round(Int, y_min * tail.Δkdisp / block.Δkdisp + 0.5)
    y_max_tail = round(Int, y_max * tail.Δkdisp / block.Δkdisp + 0.5)
    transformed_tail_task = @spawn transform_recursive(tail, y_min_tail, y_max_tail)

    n_samp = size(block)[1]
    j_range = 1:n_samp

    transformed_head = fetch(transformed_head_task)
    transformed_tail = fetch(transformed_tail_task)

    # Merge
    @tturbo for y in y_min:y_max
        # yh = delay across head band
        yh = round(Int, y * head.Δkdisp / block.Δkdisp + 0.5)
        # yt = delay across tail band
        yt = round(Int, y * tail.Δkdisp / block.Δkdisp + 0.5)
        # yb = delay at interface between head and tail
        yb = y - yh - yt
        ih = yh - transformed_head.y_min + 1
        it = yt - transformed_tail.y_min + 1
        i = y - out.y_min + 1
        # Update
        for j in j_range
            j_shift = circmod(j + yh - yb, n_samp)
            out.data[j, i] = transformed_head.data[j, ih] +
                             transformed_tail.data[j_shift, it]
        end
    end

    return out
end

"""
    fdmt(data,f_ch1,f_chn,t_samp,dm_min,dm_max)

Computes the Fast DM Transform for dynamic spectra `data`. This spectra must have time in
the first axis and cover frequencies from `f_ch1` to `f_chn` in descending order (in MHz).
The sample spacing in time (seconds) is given by `t_samp` and the transform covers DMs `dm_min`
to `dm_max`
"""
function fdmt(data::AbstractMatrix, f_ch1::Real, f_chn::Real, t_samp::Real,
              dm_min::Real, dm_max::Real)
    @assert dm_min >= 0 "Minimum DM must be zero"
    @assert dm_max >= dm_min "Maximum DM must be greater than the minimum"

    # Build input block
    block = InputBlock(data, f_ch1, f_chn, t_samp)
    dm_step = t_samp / block.Δkdisp

    # Convert DMs to delays in sample space
    y_min = trunc(Int, dm_min / dm_step)
    y_max = ceil(Int, dm_max / dm_step)

    # Perform transformation
    output = transform_recursive(block, y_min, y_max)
    return output.data
end

export fdmt


text:
using Glob

if isempty(Glob.glob("CEDMF.so"))

    using Pkg
    using PackageCompiler
    using CalibrateEDMF

    cedmf = pkgdir(CalibrateEDMF)
    pkgs = [:CalibrateEDMF]
    append!(pkgs, [Symbol(v.name) for v in values(Pkg.dependencies()) if v.is_direct_dep])

    create_sysimage(
        pkgs;
        sysimage_path = "CEDMF.so",
        # Caltech Central CPU architecture, `native` leads to issues as well.
        # This one works most of the time, but needs a failsafe mechanism.
        cpu_target = "skylake-avx512",
        precompile_execution_file = joinpath(cedmf, "test", "runtests.jl"),
    )

    # Other cpu_target options
    #cpu_target = "native", # Default
    #cpu_target = PackageCompiler.default_app_cpu_target(),
    #cpu_target = "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)", # Same as Julia Base
    #cpu_target = "generic",

end


text:
using Test
using LogisticRegression: predict

W = ones(1, 1)
b = zeros(1)
y1 = predict(W, b, [-1.0])
y2 = predict(W, b, [ 1.0])
@test y1 == [0]
@test y2 == [1]


text:
# Bit Vectors
# ===========

"""
A bit vector type to represent hash values.
"""
abstract type BitVec end

# Concrete types.
primitive type BitVec64  <: BitVec  64 end
primitive type BitVec128 <: BitVec 128 end
primitive type BitVec256 <: BitVec 256 end
primitive type BitVec512 <: BitVec 512 end

using Core: Intrinsics

const BaseInts = Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64}

# Define bit operations.
for T in [:BitVec64, :BitVec128, :BitVec256, :BitVec512]
    @eval begin
        $T(x::BaseInts) = convert($T, x)
        Base.convert(::Type{$T}, x::BaseInts) = Intrinsics.zext_int($T, x)
        Base.zero(::Type{$T}) = convert($T, 0)
        Base.:&(bv1::$T, bv2::$T) = Intrinsics.and_int(bv1, bv2)
        Base.:|(bv1::$T, bv2::$T) = Intrinsics.or_int(bv1, bv2)
        Base.:⊻(bv1::$T, bv2::$T) = Intrinsics.xor_int(bv1, bv2)
        Base.:<<(bv::$T, n::BaseInts) = Intrinsics.shl_int(bv, n)
        Base.:>>(bv::$T, n::BaseInts) = Intrinsics.lshr_int(bv, n)
        Base.getindex(bv::$T, i::Integer) = Intrinsics.trunc_int(Bool, bv >> (i - 1))
        Base.rem(bv::$T, ::Type{U}) where {U<:BaseInts} = Intrinsics.trunc_int(U, bv)
        Base.count_ones(bv::$T) = Intrinsics.trunc_int(Int, Intrinsics.ctpop_int(bv))
    end
end

function Base.rand(mt::MersenneTwister, ::Type{T}) where {T<:BitVec}
    bv::T = 0
    for _ in 1:div(sizeof(T), 8)
        bv = (bv << 64) | T(rand(mt, UInt64))
    end
    return bv
end


text:
# Convenience type to allow dispatch on solution objects that were created by Trixi
#
# This is a union of a Trixi-specific DiffEqBase.ODESolution and of Trixi's own
# TimeIntegratorSolution.
#
# Note: This is an experimental feature and may be changed in future releases without notice.
const TrixiODESolution = Union{ODESolution{T, N, uType, uType2, DType, tType, rateType, P} where
    {T, N, uType, uType2, DType, tType, rateType, P<:ODEProblem{uType_, tType_, isinplace, P_, F_} where
     {uType_, tType_, isinplace, P_<:AbstractSemidiscretization, F_<:ODEFunction{true, typeof(rhs!)}}}, TimeIntegratorSolution}

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# This file holds plotting types which can be used for both Plots.jl and Makie.jl.

# This abstract type is used to derive PlotData types of different dimensions; but still allows to share some functions for them.
abstract type AbstractPlotData{NDIMS} end

# Define additional methods for convenience.
# These are defined for AbstractPlotData, so they can be used for all types of plot data.
Base.firstindex(pd::AbstractPlotData) = first(pd.variable_names)
Base.lastindex(pd::AbstractPlotData) = last(pd.variable_names)
Base.length(pd::AbstractPlotData) = length(pd.variable_names)
Base.size(pd::AbstractPlotData) = (length(pd),)
Base.keys(pd::AbstractPlotData) = tuple(pd.variable_names...)

function Base.iterate(pd::AbstractPlotData, state=1)
  if state > length(pd)
    return nothing
  else
    return (pd.variable_names[state] => pd[pd.variable_names[state]], state + 1)
  end
end

"""
    Base.getindex(pd::AbstractPlotData, variable_name)

Extract a single variable `variable_name` from `pd` for plotting with `Plots.plot`.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
function Base.getindex(pd::AbstractPlotData, variable_name)
  variable_id = findfirst(isequal(variable_name), pd.variable_names)

  if isnothing(variable_id)
    throw(KeyError(variable_name))
  end

  return PlotDataSeries(pd, variable_id)
end

Base.eltype(pd::AbstractPlotData) = Pair{String, PlotDataSeries{typeof(pd)}}

"""
    PlotData2D

Holds all relevant data for creating 2D plots of multiple solution variables and to visualize the
mesh.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
struct PlotData2DCartesian{Coordinates, Data, VariableNames, Vertices} <: AbstractPlotData{2}
  x::Coordinates
  y::Coordinates
  data::Data
  variable_names::VariableNames
  mesh_vertices_x::Vertices
  mesh_vertices_y::Vertices
  orientation_x::Int
  orientation_y::Int
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pd::PlotData2DCartesian)
  @nospecialize pd # reduce precompilation time

  print(io, "PlotData2DCartesian{",
            typeof(pd.x), ",",
            typeof(pd.data), ",",
            typeof(pd.variable_names), ",",
            typeof(pd.mesh_vertices_x),
            "}(<x>, <y>, <data>, <variable_names>, <mesh_vertices_x>, <mesh_vertices_y>)")
end

# holds plotting information for UnstructuredMesh2D and DGMulti-compatible meshes
struct PlotData2DTriangulated{DataType, NodeType, FaceNodeType, FaceDataType, VariableNames, PlottingTriangulation} <: AbstractPlotData{2}
  x::NodeType # physical nodal coordinates, size (num_plotting_nodes x num_elements)
  y::NodeType
  data::DataType
  t::PlottingTriangulation
  x_face::FaceNodeType
  y_face::FaceNodeType
  face_data::FaceDataType
  variable_names::VariableNames
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pd::PlotData2DTriangulated)
  @nospecialize pd # reduce precompilation time

  print(io, "PlotData2DTriangulated{",
            typeof(pd.x), ", ",
            typeof(pd.data), ", ",
            typeof(pd.x_face), ", ",
            typeof(pd.face_data), ", ",
            typeof(pd.variable_names),
            "}(<x>, <y>, <data>, <plot_triangulation>, <x_face>, <y_face>, <face_data>, <variable_names>)")
end

"""
    PlotData1D

Holds all relevant data for creating 1D plots of multiple solution variables and to visualize the
mesh.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
struct PlotData1D{Coordinates, Data, VariableNames, Vertices} <: AbstractPlotData{1}
  x::Coordinates
  data::Data
  variable_names::VariableNames
  mesh_vertices_x::Vertices
  orientation_x::Integer
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pd::PlotData1D)
  print(io, "PlotData1D{",
            typeof(pd.x), ",",
            typeof(pd.data), ",",
            typeof(pd.variable_names), ",",
            typeof(pd.mesh_vertices_x),
            "}(<x>, <data>, <variable_names>, <mesh_vertices_x>)")
end

# Auxiliary data structure for visualizing a single variable
#
# Note: This is an experimental feature and may be changed in future releases without notice.
struct PlotDataSeries{PD<:AbstractPlotData}
  plot_data::PD
  variable_id::Int
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pds::PlotDataSeries)
  @nospecialize pds # reduce precompilation time

  print(io, "PlotDataSeries{", typeof(pds.plot_data), "}(<plot_data>, ",
        pds.variable_id, ")")
end

# Generic PlotMesh wrapper type.
struct PlotMesh{PD<:AbstractPlotData}
  plot_data::PD
end

# Show only a truncated output for convenience (the full data does not make sense)
function Base.show(io::IO, pm::PlotMesh)
  @nospecialize pm # reduce precompilation time

  print(io, "PlotMesh{", typeof(pm.plot_data), "}(<plot_data>)")
end

"""
    getmesh(pd::AbstractPlotData)

Extract grid lines from `pd` for plotting with `Plots.plot`.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
getmesh(pd::AbstractPlotData) = PlotMesh(pd)

"""
    PlotData2D(u, semi [or mesh, equations, solver, cache];
               solution_variables=nothing,
               grid_lines=true, max_supported_level=11, nvisnodes=nothing,
               slice=:xy, point=(0.0, 0.0, 0.0))

Create a new `PlotData2D` object that can be used for visualizing 2D/3D DGSEM solution data array
`u` with `Plots.jl`. All relevant geometrical information is extracted from the semidiscretization
`semi`. By default, the primitive variables (if existent) or the conservative variables (otherwise)
from the solution are used for plotting. This can be changed by passing an appropriate conversion
function to `solution_variables`.

If `grid_lines` is `true`, also extract grid vertices for visualizing the mesh. The output
resolution is indirectly set via `max_supported_level`: all data is interpolated to
`2^max_supported_level` uniformly distributed points in each spatial direction, also setting the
maximum allowed refinement level in the solution. `nvisnodes` specifies the number of visualization
nodes to be used. If it is `nothing`, twice the number of solution DG nodes are used for
visualization, and if set to `0`, exactly the number of nodes in the DG elements are used.

When visualizing data from a three-dimensional simulation, a 2D slice is extracted for plotting.
`slice` specifies the plane that is being sliced and may be `:xy`, `:xz`, or `:yz`.
The slice position is specified by a `point` that lies on it, which defaults to `(0.0, 0.0, 0.0)`.
Both of these values are ignored when visualizing 2D data.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.

# Examples
```julia
julia> using Trixi, Plots

julia> trixi_include(default_example())
[...]

julia> pd = PlotData2D(sol)
PlotData2D(...)

julia> plot(pd) # To plot all available variables

julia> plot(pd["scalar"]) # To plot only a single variable

julia> plot!(getmesh(pd)) # To add grid lines to the plot
```
"""

PlotData2D(u_ode, semi; kwargs...) = PlotData2D(wrap_array_native(u_ode, semi),
                                                mesh_equations_solver_cache(semi)...;
                                                kwargs...)

# Redirect `PlotDataTriangulated2D` constructor.
PlotData2DTriangulated(u_ode, semi; kwargs...) = PlotData2DTriangulated(wrap_array_native(u_ode, semi),
                                                                        mesh_equations_solver_cache(semi)...;
                                                                        kwargs...)

# Create a PlotData2DCartesian object for TreeMeshes on default.
PlotData2D(u, mesh::TreeMesh, equations, solver, cache; kwargs...) = PlotData2DCartesian(u, mesh::TreeMesh, equations, solver, cache; kwargs...)

# Create a PlotData2DTriangulated object for any type of mesh other than the TreeMesh.
PlotData2D(u, mesh, equations, solver, cache; kwargs...) = PlotData2DTriangulated(u, mesh, equations, solver, cache; kwargs...)

# Create a PlotData2DCartesian for a TreeMesh.
function PlotData2DCartesian(u, mesh::TreeMesh, equations, solver, cache;
                             solution_variables=nothing,
                             grid_lines=true, max_supported_level=11, nvisnodes=nothing,
                             slice=:xy, point=(0.0, 0.0, 0.0))
  @assert ndims(mesh) in (2, 3) "unsupported number of dimensions $ndims (must be 2 or 3)"
  solution_variables_ = digest_solution_variables(equations, solution_variables)

  # Extract mesh info
  center_level_0 = mesh.tree.center_level_0
  length_level_0 = mesh.tree.length_level_0
  leaf_cell_ids = leaf_cells(mesh.tree)
  coordinates = mesh.tree.coordinates[:, leaf_cell_ids]
  levels = mesh.tree.levels[leaf_cell_ids]

  unstructured_data = get_unstructured_data(u, solution_variables_, mesh, equations, solver, cache)
  x, y, data, mesh_vertices_x, mesh_vertices_y = get_data_2d(center_level_0, length_level_0,
                                                             leaf_cell_ids, coordinates, levels,
                                                             ndims(mesh), unstructured_data,
                                                             nnodes(solver), grid_lines,
                                                             max_supported_level, nvisnodes,
                                                             slice, point)
  variable_names = SVector(varnames(solution_variables_, equations))

  orientation_x, orientation_y = _get_orientations(mesh, slice)

  return PlotData2DCartesian(x, y, data, variable_names, mesh_vertices_x, mesh_vertices_y,
                             orientation_x, orientation_y)
end

"""
    PlotData2D(sol; kwargs...)

Create a `PlotData2D` object from a solution object created by either `OrdinaryDiffEq.solve!` (which
returns a `DiffEqBase.ODESolution`) or Trixi's own `solve!` (which returns a
`TimeIntegratorSolution`).

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
PlotData2D(sol::TrixiODESolution; kwargs...) = PlotData2D(sol.u[end], sol.prob.p; kwargs...)

# Also redirect when using PlotData2DTriangulate.
PlotData2DTriangulated(sol::TrixiODESolution; kwargs...) = PlotData2DTriangulated(sol.u[end], sol.prob.p; kwargs...)

# If `u` is an `Array{<:SVectors}` and not a `StructArray`, convert it to a `StructArray` first.
function PlotData2D(u::Array{<:SVector, 2}, mesh, equations, dg::DGMulti, cache;
                    solution_variables=nothing, nvisnodes=2*nnodes(dg))
  nvars = length(first(u))
  u_structarray = StructArray{eltype(u)}(ntuple(_->zeros(eltype(first(u)), size(u)), nvars))
  for (i, u_i) in enumerate(u)
    u_structarray[i] = u_i
  end

  # re-dispatch to PlotData2D with mesh, equations, dg, cache arguments
  return PlotData2D(u_structarray, mesh, equations, dg, cache;
                    solution_variables=solution_variables, nvisnodes=nvisnodes)
end

# constructor which returns an `PlotData2DTriangulated` object.
function PlotData2D(u::StructArray, mesh, equations, dg::DGMulti, cache;
                    solution_variables=nothing, nvisnodes=2*nnodes(dg))

  rd = dg.basis
  md = mesh.md

  # Vp = the interpolation matrix from nodal points to plotting points
  @unpack Vp = rd
  interpolate_to_plotting_points!(out, x) = mul!(out, Vp, x)

  solution_variables_ = digest_solution_variables(equations, solution_variables)
  variable_names = SVector(varnames(solution_variables_, equations))

  num_plotting_points = size(Vp, 1)
  nvars = nvariables(equations)
  uEltype = eltype(first(u))
  u_plot = StructArray{SVector{nvars, uEltype}}(ntuple(_->zeros(uEltype, num_plotting_points, md.num_elements), nvars))

  for e in eachelement(mesh, dg, cache)
    # interpolate solution to plotting nodes element-by-element
    StructArrays.foreachfield(interpolate_to_plotting_points!, view(u_plot, :, e), view(u, :, e))

    # transform nodal values of the solution according to `solution_variables`
    transform_to_solution_variables!(view(u_plot, :, e), solution_variables_, equations)
  end

  # interpolate nodal coordinates to plotting points
  x_plot, y_plot = map(x->Vp * x, md.xyz) # md.xyz is a tuple of arrays containing nodal coordinates

  # construct a triangulation of the reference plotting nodes
  t = reference_plotting_triangulation(rd.rstp) # rd.rstp = reference coordinates of plotting points

  x_face, y_face, face_data = mesh_plotting_wireframe(u, mesh, equations, dg, cache;
                                                      nvisnodes=nvisnodes)

  return PlotData2DTriangulated(x_plot, y_plot, u_plot, t, x_face, y_face, face_data, variable_names)
end

# specializes the PlotData2D constructor to return an PlotData2DTriangulated for any type of mesh.
function PlotData2DTriangulated(u, mesh, equations, dg::DGSEM, cache;
                                solution_variables=nothing, nvisnodes=2*polydeg(dg))

  @assert ndims(mesh) == 2 "Input must be two-dimensional."

  n_nodes_2d = nnodes(dg)^ndims(mesh)
  n_elements = nelements(dg, cache)

  # build nodes on reference element (seems to be the right ordering)
  r, s = reference_node_coordinates_2d(dg)

  # reference plotting nodes
  if nvisnodes == 0 || nvisnodes === nothing
    nvisnodes = polydeg(dg) + 1
  end
  plotting_interp_matrix = plotting_interpolation_matrix(dg; nvisnodes=nvisnodes)

  # create triangulation for plotting nodes
  r_plot, s_plot = (x->plotting_interp_matrix*x).((r, s)) # interpolate dg nodes to plotting nodes

  # construct a triangulation of the plotting nodes
  t = reference_plotting_triangulation((r_plot, s_plot))

  # extract x,y coordinates and solutions on each element
  uEltype = eltype(u)
  nvars = nvariables(equations)
  x = reshape(view(cache.elements.node_coordinates, 1, :, :, :), n_nodes_2d, n_elements)
  y = reshape(view(cache.elements.node_coordinates, 2, :, :, :), n_nodes_2d, n_elements)
  u_extracted = StructArray{SVector{nvars, uEltype}}(ntuple(_->similar(x, (n_nodes_2d, n_elements)), nvars))
  for element in eachelement(dg, cache)
    sk = 1
    for j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, element)
      u_extracted[sk, element] = u_node
      sk += 1
    end
  end

  # interpolate to volume plotting points
  xplot, yplot = plotting_interp_matrix*x, plotting_interp_matrix*y
  uplot = StructArray{SVector{nvars, uEltype}}(map(x->plotting_interp_matrix*x,
                                                   StructArrays.components(u_extracted)))

  xfp, yfp, ufp = mesh_plotting_wireframe(u_extracted, mesh, equations, dg, cache; nvisnodes=nvisnodes)

  # convert variables based on solution_variables mapping
  solution_variables_ = digest_solution_variables(equations, solution_variables)
  variable_names = SVector(varnames(solution_variables_, equations))

  transform_to_solution_variables!(uplot, solution_variables_, equations)
  transform_to_solution_variables!(ufp, solution_variables_, equations)

  return PlotData2DTriangulated(xplot, yplot, uplot, t, xfp, yfp, ufp, variable_names)
end

# Wrapper struct to indicate that an array represents a scalar data field. Used only for dispatch.
struct ScalarData{T}
  data::T
end

ScalarPlotData2D(u, semi::AbstractSemidiscretization; kwargs...) =
  ScalarPlotData2D(u, mesh_equations_solver_cache(semi)...; kwargs...)

# Returns an `PlotData2DTriangulated` which is used to visualize a single scalar field
function ScalarPlotData2D(u, mesh, equations, dg::DGMulti, cache;
                          variable_name=nothing, nvisnodes=2*nnodes(dg))

  rd = dg.basis
  md = mesh.md

  # Vp = the interpolation matrix from nodal points to plotting points
  @unpack Vp = rd

  # interpolate nodal coordinates and solution field to plotting points
  x_plot, y_plot = map(x->Vp * x, md.xyz) # md.xyz is a tuple of arrays containing nodal coordinates
  u_plot = Vp * u

  # construct a triangulation of the reference plotting nodes
  t = reference_plotting_triangulation(rd.rstp) # rd.rstp = reference coordinates of plotting points

  # Ignore face data when plotting `ScalarPlotData2D`, since mesh lines can be plotted using
  # existing functionality based on `PlotData2D(sol)`.
  x_face, y_face, face_data = mesh_plotting_wireframe(ScalarData(u), mesh, equations, dg, cache;
                                                      nvisnodes=2*nnodes(dg))

  # wrap solution in ScalarData struct for recipe dispatch
  return PlotData2DTriangulated(x_plot, y_plot, ScalarData(u_plot), t,
                                x_face, y_face, face_data, variable_name)
end

function ScalarPlotData2D(u, mesh, equations, dg::DGSEM, cache; variable_name=nothing, nvisnodes=2*nnodes(dg))

  n_nodes_2d = nnodes(dg)^ndims(mesh)
  n_elements = nelements(dg, cache)

  # build nodes on reference element (seems to be the right ordering)
  r, s = reference_node_coordinates_2d(dg)

  # reference plotting nodes
  if nvisnodes == 0 || nvisnodes === nothing
    nvisnodes = polydeg(dg) + 1
  end
  plotting_interp_matrix = plotting_interpolation_matrix(dg; nvisnodes=nvisnodes)

  # create triangulation for plotting nodes
  r_plot, s_plot = (x->plotting_interp_matrix*x).((r, s)) # interpolate dg nodes to plotting nodes

  # construct a triangulation of the plotting nodes
  t = reference_plotting_triangulation((r_plot, s_plot))

  # extract x,y coordinates and reshape them into matrices of size (n_nodes_2d, n_elements)
  x = view(cache.elements.node_coordinates, 1, :, :, :)
  y = view(cache.elements.node_coordinates, 2, :, :, :)
  x, y = reshape.((x, y), n_nodes_2d, n_elements)

  # interpolate to volume plotting points by multiplying each column by `plotting_interp_matrix`
  x_plot, y_plot = plotting_interp_matrix * x, plotting_interp_matrix * y
  u_plot = plotting_interp_matrix * reshape(u, size(x))

  # Ignore face data when plotting `ScalarPlotData2D`, since mesh lines can be plotted using
  # existing functionality based on `PlotData2D(sol)`.
  x_face, y_face, face_data = mesh_plotting_wireframe(ScalarData(u), mesh, equations, dg, cache;
                                                      nvisnodes=2*nnodes(dg))

  # wrap solution in ScalarData struct for recipe dispatch
  return PlotData2DTriangulated(x_plot, y_plot, ScalarData(u_plot), t,
                                x_face, y_face, face_data, variable_name)
end

"""
    PlotData1D(u, semi [or mesh, equations, solver, cache];
               solution_variables=nothing, nvisnodes=nothing)

Create a new `PlotData1D` object that can be used for visualizing 1D DGSEM solution data array
`u` with `Plots.jl`. All relevant geometrical information is extracted from the semidiscretization
`semi`. By default, the primitive variables (if existent) or the conservative variables (otherwise)
from the solution are used for plotting. This can be changed by passing an appropriate conversion
function to `solution_variables`.

`nvisnodes` specifies the number of visualization nodes to be used. If it is `nothing`,
twice the number of solution DG nodes are used for visualization, and if set to `0`,
exactly the number of nodes in the DG elements are used.

When visualizing data from a two-dimensional simulation, a 1D slice is extracted for plotting.
`slice` specifies the axis along which the slice is extracted and may be `:x`, or `:y`.
The slice position is specified by a `point` that lies on it, which defaults to `(0.0, 0.0)`.
Both of these values are ignored when visualizing 1D data.
This applies analogously to three-dimensonal simulations, where `slice` may be `xy:`, `:xz`, or `:yz`.

Another way to visualize 2D/3D data is by creating a plot along a given curve.
This is done with the keyword argument `curve`. It can be set to a list of 2D/3D points
which define the curve. When using `curve` any other input from `slice` or `point` will be ignored.

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
PlotData1D(u_ode, semi; kwargs...) = PlotData1D(wrap_array_native(u_ode, semi),
                                                mesh_equations_solver_cache(semi)...;
                                                kwargs...)

function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
                    solution_variables=nothing, nvisnodes=nothing,
                    slice=:x, point=(0.0, 0.0, 0.0), curve=nothing)

  solution_variables_ = digest_solution_variables(equations, solution_variables)
  variable_names = SVector(varnames(solution_variables_, equations))

  original_nodes = cache.elements.node_coordinates
  unstructured_data = get_unstructured_data(u, solution_variables_, mesh, equations, solver, cache)

  orientation_x = 0 # Set 'orientation' to zero on default.

  if ndims(mesh) == 1
    x, data, mesh_vertices_x = get_data_1d(original_nodes, unstructured_data, nvisnodes)
    orientation_x = 1
  elseif ndims(mesh) == 2
    if curve !== nothing
      x, data, mesh_vertices_x = unstructured_2d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, curve, mesh, solver, cache)
    else
      x, data, mesh_vertices_x = unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, point)
    end
  else # ndims(mesh) == 3
    if curve !== nothing
      x, data, mesh_vertices_x = unstructured_3d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, curve, mesh, solver, cache)
    else
      x, data, mesh_vertices_x = unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, point)
    end
  end

  return PlotData1D(x, data, variable_names, mesh_vertices_x,
                    orientation_x)
end

function PlotData1D(u, mesh, equations, solver, cache;
                    solution_variables=nothing, nvisnodes=nothing,
                    slice=:x, point=(0.0, 0.0, 0.0), curve=nothing)

  solution_variables_ = digest_solution_variables(equations, solution_variables)
  variable_names = SVector(varnames(solution_variables_, equations))

  original_nodes = cache.elements.node_coordinates
  unstructured_data = get_unstructured_data(u, solution_variables_, mesh, equations, solver, cache)

  orientation_x = 0 # Set 'orientation' to zero on default.

  if ndims(mesh) == 1
    x, data, mesh_vertices_x = get_data_1d(original_nodes, unstructured_data, nvisnodes)
    orientation_x = 1
  else # ndims(mesh) == 2
    pd = PlotData2DTriangulated(u, mesh, equations, solver, cache; solution_variables, nvisnodes)
    x, data, mesh_vertices_x = unstructured_2d_to_1d_curve(pd, curve, slice, point)
  end

  return PlotData1D(x, data, variable_names, mesh_vertices_x,
                    orientation_x)
end

"""
    PlotData1D(sol; kwargs...)

Create a `PlotData1D` object from a solution object created by either `OrdinaryDiffEq.solve!` (which
returns a `DiffEqBase.ODESolution`) or Trixi's own `solve!` (which returns a
`TimeIntegratorSolution`).
!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
PlotData1D(sol::TrixiODESolution; kwargs...) = PlotData1D(sol.u[end], sol.prob.p; kwargs...)

function PlotData1D(time_series_callback::TimeSeriesCallback, point_id::Integer)
  @unpack time, variable_names, point_data = time_series_callback

  n_solution_variables = length(variable_names)
  data = Matrix{Float64}(undef, length(time), n_solution_variables)
  reshaped = reshape(point_data[point_id], n_solution_variables, length(time))
  for v in 1:n_solution_variables
    @views data[:, v] = reshaped[v, :]
  end

  mesh_vertices_x = Vector{Float64}(undef, 0)

  return PlotData1D(time, data, SVector(variable_names), mesh_vertices_x, 0)
end

function PlotData1D(cb::DiscreteCallback{<:Any, <:TimeSeriesCallback}, point_id::Integer)
  return PlotData1D(cb.affect!, point_id)
end

end # @muladd


text:
# conjgrad-tests.jl -
#
# Tests for conjugate gradient algorithm.
#

if ! isdefined(:OptimPackNextGen)
    include("../src/OptimPackNextGen.jl")
end

module ConjGradTests

using Base.Test
using OptimPackNextGen

verb = true
@testset "conjgrad" begin
    for (T, prec) in ((Float16, 5e-2),
                      (Float32, 1e-5),
                      (Float64, 1e-7))
        # FIXME: Due to the randomness, the tests mail fail...
        dims = (70,39)
        H = randn(T, dims)
        H *= T(1/maximum(abs(H)))
        x = randn(T, dims[2])
        y = H*x + T(0.01)*randn(T, dims[1])
        A = H'*H
        A!(dst, src) = (dst[:] = A*src; dst)
        b = H'*y
        x0 = A\b
        err = maximum(abs(conjgrad(A!, b; verb=verb, ftol=1e-16,
                                   maxiter=length(b)) - x0))

        verb && println(T, " --> ", err, "\n")
        @test err ≤ prec
    end
end

end # module


text:
using Test
using JSExpr
using JSExpr: @crawl, deparse, JSTerminal, JSString

@testset "control expressions" begin

    @testset "if expressions" begin
        @test @js(
            if foo
                bar()
            end
        ) == JSString("if (foo) { bar(); }")

        @test @js(
            if foo
                bar()
            else
                spam()
            end
        ) == JSString("if (foo) { bar(); } else { spam(); }")

        @test string(@js(
            if foo
                foo()
            elseif bar
                bar()
            elseif spam
                spam()
            else
                eggs()
            end
        )) == string(
            "if (foo) { foo(); } elseif (bar) { bar(); } ",
            "elseif (spam) { spam(); } else { eggs(); }"
        )
    end

    @testset "while expressions" begin
        @test @js(
            while true
                foo()
            end
        ) == JSString("while (true) { foo(); }")

        @test @js(
            while foo
                if (bar)
                    continue
                end
                break
            end
        ) == JSString("while (foo) { if (bar) { continue; }; break; }");
    end

    @testset "for expressions" begin
        @test @js(
            for i in myiterable
                console.log(i)
            end
        ) == JSString("for (let i of myiterable) { console.log(i); }")

        @test_throws ErrorException JSExpr.crawl(:(
            for i in 1:10, j in i:10
                console.log(i, j)
            end
        ))
    end

end


text:
##############################################################
# constructor for Sem types
##############################################################

function Sem(;
        observed::O = SemObsCommon,
        imply::I = RAM,
        loss::L = SemML,
        diff::D = SemDiffOptim,
        kwargs...) where {O, I, L, D}

    kwargs = Dict{Symbol, Any}(kwargs...)

    set_field_type_kwargs!(kwargs, observed, imply, loss, diff, O, I, D)

    observed, imply, loss, diff = get_fields!(kwargs, observed, imply, loss, diff)

    sem = Sem(observed, imply, loss, diff)

    return sem
end

function SemFiniteDiff(;
        observed::O = SemObsCommon,
        imply::I = RAM,
        loss::L = SemML,
        diff::D = SemDiffOptim,
        has_gradient = false,
        kwargs...) where {O, I, L, D}

    kwargs = Dict{Symbol, Any}(kwargs...)

    set_field_type_kwargs!(kwargs, observed, imply, loss, diff, O, I, D)

    observed, imply, loss, diff = get_fields!(kwargs, observed, imply, loss, diff)

    sem = SemFiniteDiff(observed, imply, loss, diff, has_gradient)

    return sem
end

function SemForwardDiff(;
        observed::O = SemObsCommon,
        imply::I = RAM,
        loss::L = SemML,
        diff::D = SemDiffOptim,
        has_gradient = false,
        kwargs...) where {O, I, L, D}

    kwargs = Dict{Symbol, Any}(kwargs...)

    set_field_type_kwargs!(kwargs, observed, imply, loss, diff, O, I, D)

    observed, imply, loss, diff = get_fields!(kwargs, observed, imply, loss, diff)

    sem = SemForwardDiff(observed, imply, loss, diff, has_gradient)

    return sem
end

##############################################################
# functions
##############################################################

function set_field_type_kwargs!(kwargs, observed, imply, loss, diff, O, I, D)
    kwargs[:observed_type] = O <: Type ? observed : typeof(observed)
    kwargs[:imply_type] = I <: Type ? imply : typeof(imply)
    if loss isa SemLoss
        kwargs[:loss_types] = [lossfun isa SemLossFunction ? typeof(lossfun) : lossfun for lossfun in loss.functions]
    elseif applicable(iterate, loss)
        kwargs[:loss_types] = [lossfun isa SemLossFunction ? typeof(lossfun) : lossfun for lossfun in loss]
    else
        kwargs[:loss_types] = [loss isa SemLossFunction ? typeof(loss) : loss]
    end
    kwargs[:diff_type] = D <: Type ? diff : typeof(diff)
end

# construct Sem fields
function get_fields!(kwargs, observed, imply, loss, diff)
    # observed
    if !isa(observed, SemObs)
        observed = observed(;kwargs...)
    end
    kwargs[:observed] = observed

    # imply
    if !isa(imply, SemImply)
        imply = imply(;kwargs...)
    end

    kwargs[:imply] = imply
    kwargs[:n_par] = n_par(imply)

    # loss
    loss = get_SemLoss(loss; kwargs...)
    kwargs[:loss] = loss

    # diff
    if !isa(diff, SemDiff)
        diff = diff(;kwargs...)
    end

    return observed, imply, loss, diff
end

# construct loss field
function get_SemLoss(loss; kwargs...)
    if loss isa SemLoss
        nothing
    elseif applicable(iterate, loss)
        loss_out = []
        for lossfun in loss
            if isa(lossfun, SemLossFunction)
                push!(loss_out, lossfun)
            else
                lossfun = lossfun(;kwargs...)
                push!(loss_out, lossfun)
            end
        end
        loss = SemLoss(loss_out...; kwargs...)
    else
        if !isa(loss, SemLossFunction)
            loss = SemLoss(loss(;kwargs...); kwargs...)
        end
    end
    return loss
end

##############################################################
# pretty printing
##############################################################

#= function Base.show(io::IO, sem::Sem{O, I, L, D})  where {O, I, L, D}
    lossfuntypes = @. nameof(typeof(sem.loss.functions))
    print(io, "Sem{$(nameof(O)), $(nameof(I)), $lossfuntypes, $(nameof(D))}")
end =#

function Base.show(io::IO, sem::Sem{O, I, L, D})  where {O, I, L, D}
    lossfuntypes = @. string(nameof(typeof(sem.loss.functions)))
    lossfuntypes = "   ".*lossfuntypes.*("\n")
    print(io, "Structural Equation Model \n")
    print(io, "- Loss Functions \n")
    print(io, lossfuntypes...)
    print(io, "- Fields \n")
    print(io, "   observed:  $(nameof(O)) \n")
    print(io, "   imply:     $(nameof(I)) \n")
    print(io, "   diff:      $(nameof(D)) \n")
end

function Base.show(io::IO, sem::SemFiniteDiff{O, I, L, D})  where {O, I, L, D}
    lossfuntypes = @. string(nameof(typeof(sem.loss.functions)))
    lossfuntypes = "   ".*lossfuntypes.*("\n")
    print(io, "Structural Equation Model : Finite Diff Approximation\n")
    print(io, "- Loss Functions \n")
    print(io, lossfuntypes...)
    print(io, "- Fields \n")
    print(io, "   observed:  $(nameof(O)) \n")
    print(io, "   imply:     $(nameof(I)) \n")
    print(io, "   diff:      $(nameof(D)) \n")
end

function Base.show(io::IO, sem::SemForwardDiff{O, I, L, D})  where {O, I, L, D}
    lossfuntypes = @. string(nameof(typeof(sem.loss.functions)))
    lossfuntypes = "   ".*lossfuntypes.*("\n")
    print(io, "Structural Equation Model : Forward Mode Autodiff\n")
    print(io, "- Loss Functions \n")
    print(io, lossfuntypes...)
    print(io, "- Fields \n")
    print(io, "   observed:  $(nameof(O)) \n")
    print(io, "   imply:     $(nameof(I)) \n")
    print(io, "   diff:      $(nameof(D)) \n")
end

function Base.show(io::IO, loss::SemLoss)
    lossfuntypes = @. string(nameof(typeof(loss.functions)))
    lossfuntypes = "   ".*lossfuntypes.*("\n")
    print(io, "SemLoss \n")
    print(io, "- Loss Functions \n")
    print(io, lossfuntypes...)
    print(io, "- Fields \n")
    print(io, "   F:  $(typeof(loss.F))) \n")
    print(io, "   G:  $(typeof(loss.G))) \n")
    print(io, "   H:  $(typeof(loss.H))) \n")
end

function Base.show(io::IO, models::SemEnsemble)

    print(io, "SemEnsemble \n")
    print(io, "- Number of Models: $(models.n) \n")
    print(io, "- Weights: $(round.(models.weights, digits = 2)) \n")
    print(io, "- diff: $(nameof(typeof(models.diff))) \n")

    print(io, "\n", "Models: \n")
    print(io, "=========================", "\n")
    for (model, i) in zip(models.sems, 1:models.n)
        print(io, "----------- ", i, " -----------", "\n")
        print(io, model)
    end
end


text:
#=
Copyright 2021 BlackRock, Inc.

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.
=#

@testset "Zero" begin
    @test zero(PiecewiseQuadratic) ≈ PiecewiseQuadratic(BoundedQuadratic(0.0, 0.0, 0.0))
end

@testset "Simplify" begin

    # Unsimplifiable.
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 3.0, 4.0),
                            BoundedQuadratic(3.0, 4.0, 1.0, 3.0, 4.0)])
    @test simplify(f) ≈ f

    # One non-redundant point.
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 1.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])
    @test simplify(f) ≈ f

    # Redundant point (left end).
    f = PiecewiseQuadratic([BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 5.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 1.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])
    g = PiecewiseQuadratic([BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 1.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])
    @test simplify(f) ≈ g

    # Redundant point (right end).
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0),
                            BoundedQuadratic(4.0, 4.0, 0.0, 0.0, 10.0)])
    g = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])
    @test simplify(f) ≈ g

    # Redundant point (left).
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 50.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 20.0)])
    g = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 20.0)])
    @test simplify(f) ≈ g

    # Redundant point (right).
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 20.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 50.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 4.0)])
    g = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 20.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 4.0)])
    @test simplify(f) ≈ g

    # All points in constructor
    point = BoundedQuadratic(0, 0, 0, 0, 0)
    f = PiecewiseQuadratic(repeat([point], 5))
    @test simplify(f) ≈ PiecewiseQuadratic(point)

    # No redundant point (because of right gap).
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 20.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(4.0, 5.0, 0.0, 0.0, 4.0)])
    @test simplify(f) ≈ f

    # No redundant point (because of left gap).
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 2.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 3.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 20.0)])
    @test simplify(f) ≈ f

    # Complex
    f = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0),
                            BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, 1.0, 0.0),
                            BoundedQuadratic(1.0, 2.0, 0.0, 1.0, 0.0),
                            BoundedQuadratic(3.0, 5.0, 1.0, 1.0, 1)])
    g = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0),
                            BoundedQuadratic(0.0, 2.0, 0.0, 1.0, 0.0),
                            BoundedQuadratic(3.0, 5.0, 1.0, 1.0, 1.0)])
    @test simplify(f) ≈ g

    # Simplify in constructor
    f = PiecewiseQuadratic([BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0),
                            BoundedQuadratic(0.0, 0.0, 0.0, 0.0, 0.0),
                            BoundedQuadratic(0.0, 1.0, 0.0, 1.0, 0.0),
                            BoundedQuadratic(1.0, 2.0, 0.0, 1.0, 0.0),
                            BoundedQuadratic(3.0, 5.0, 1.0, 1.0, 1)]; simplify_result=true)
    @test f ≈ g
end

@testset "Minimize" begin
    f = PiecewiseQuadratic([BoundedQuadratic(1.0, 2.0, 0.0, 0.0, 14.0),
                            BoundedQuadratic(2.0, 3.0, 0.0, 0.0, 24.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 23.0)])
    @test minimize(f) == (1.0, 14.0)
end

@testset "Evaluate" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 3.1, 0.0, 0.0, 1.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])
    @test f(2.0) ≈ 4.0
    @test f.([3.001, 3.05, 3.099]) ≈ ones(3)
    @test f(10.0) == Inf
end

@testset "Sum" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0, 0.0, 4.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 3.0)])

    g = PiecewiseQuadratic([BoundedQuadratic(1.0, 2.0, 0.0, 0.0, 10.0),
                            BoundedQuadratic(2.0, 5.0, 0.0, 0.0, 20.0)])

    h = PiecewiseQuadratic([BoundedQuadratic(1.0, 2.0, 0.0, 0.0, 14.0),
                            BoundedQuadratic(2.0, 3.0, 0.0, 0.0, 24.0),
                            BoundedQuadratic(3.0, 4.0, 0.0, 0.0, 23.0)])

    @test f + g ≈ h
    @test g + f ≈ h
    @test sum([f, g, h]) ≈ 2.0 * h

    bq = BoundedQuadratic(0.0, 10.0, 1.0, 2.0, 4.0)
    exp = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 1.0, 2.0, 8.0),
                              BoundedQuadratic(3.0, 4.0, 1.0, 2.0, 7.0)])
    @test f + bq ≈ exp
end

@testset "Sum with disjoint domain" begin
    f1 = PiecewiseQuadratic([BoundedQuadratic(1.0, 1.0, 0.0)])
    f2 = PiecewiseQuadratic([BoundedQuadratic(-Inf, -1.0, 0.0, 0.0, 0.0)])
    f3 = PiecewiseQuadratic([BoundedQuadratic(1.0, Inf, 0.0, 0.0, 0.0)])
    h = sum([f1, f2, f3])
    h_true = PiecewiseQuadratic(BoundedQuadratic[])
    @test h ≈ h_true
end

@testset "Proximal operator" begin
    # Affine.
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 3.0, 0.0)])
    @test prox(f, 5.0, 1.0) ≈ 2.0

    # Huber.
    μ = 1.0
    f = PiecewiseQuadratic([BoundedQuadratic(-Inf, -μ, 0.0, -μ, -μ^2 / 2),
                            BoundedQuadratic(-μ, μ, 0.5, 0.0, 0.0),
                            BoundedQuadratic(μ, Inf, 0.0, μ, -μ^2 / 2)])

    for u in [-5.0, 1.0, 3.0]
        for ρ in [0.1, 1.0, 10]
            σ = 1 / ρ
            prox_f = u - σ * u / max(abs(u), σ + 1)
            @test prox(f, u, ρ) ≈ prox_f
        end
    end

    # Indicator.
    f = indicator(-1, 1)

    for u in [-5.0, 0, 1.0, 3.0]
        for ρ in [0.1, 1.0, 10]
            prox_f = min(max(u, -1), 1)
            @test prox(f, u, ρ) ≈ prox_f
        end
    end
end

@testset "Shift" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 2, 4)])
    shift!(f, 5.0)
    @test f ≈ PiecewiseQuadratic([BoundedQuadratic(5, 10, 1, -8, 19)])
end

@testset "Scale" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 2, 4)])
    scale!(f, 5.0)
    @test f ≈ PiecewiseQuadratic([BoundedQuadratic(0, 1, 25, 10, 4)])
end

@testset "Tilt" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 2, 4)])
    tilt!(f, 5.0)
    @test f ≈ PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 7, 4)])
end

@testset "Perspective" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 2, 4)])
    perspective!(f, 5.0)
    @test f ≈ PiecewiseQuadratic([BoundedQuadratic(0, 25, 0.2, 2, 20)])
end

@testset "Scalar multiplication" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0.0, 10.0, 1.0, 2.0, 4.0)])
    g = PiecewiseQuadratic([BoundedQuadratic(0.0, 10.0, 3.0, 6.0, 12.0)])
    @test 3.0 * f ≈ g
    @test f * 3.0 ≈ g
    mult_scalar!(f, 3.0)
    @test f ≈ g
end

@testset "Empty" begin
    f = PiecewiseQuadratic()
    g = PiecewiseQuadratic([BoundedQuadratic(0, 5, 1, 2, 4)])
    @test isempty(f)
    @test !isempty(g)
end

@testset "Negative" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 0, 2, 4)])
    g = PiecewiseQuadratic([BoundedQuadratic(0, 5, 0, -2, -4)])
    @test -f ≈ g
end

@testset "Get Plot" begin
    f = PiecewiseQuadratic([BoundedQuadratic(0, 5, 0, 2, 4)])
    N = 10
    x = [-5.0, -3.3333333333333335, -1.6666666666666667, 0.0, 1.6666666666666667,
         3.3333333333333335, 5.0, 6.666666666666667, 8.333333333333334]
    y = [Inf, Inf, Inf, 4.0, 7.333333333333334, 10.666666666666668, 14.0, Inf, Inf]
    @test get_plot(f; N=N) == (x, y)
end

@testset "Is Convex" begin
    #####
    ##### not convex
    #####

    # test single piece that is non convex (never enters loop)
    f = BoundedQuadratic(0.0, 1.0, -1.0, 0.0, 1.0)
    p = PiecewiseQuadratic(f)
    @test(!is_convex(p))

    # test non continuous
    f = BoundedQuadratic(0.0, 1.0, 1.0, 0.0, 1.0)
    g = BoundedQuadratic(5.0, 6.0, 2.0, 0.0, 1.0)
    p = PiecewiseQuadratic([f, g])
    @test(!is_convex(p))

    # test derivatives out of order
    f = BoundedQuadratic(-1.0, 0.0, 0.0, 1.0, 1.0)
    g = BoundedQuadratic(0.0, 1.0, 0.0, -1.0, 1.0)
    p = PiecewiseQuadratic([f, g])
    @test(!is_convex(p))

    # test middle piece not convex
    f = BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0)
    g = BoundedQuadratic(0.0, 1.0, -1.0, 0.0, 0.0)
    h = BoundedQuadratic(1.0, 2.0, 0.0, 1.0, -2.0)
    p = PiecewiseQuadratic([f, g, h])
    @test(!is_convex(p))

    #####
    ##### is convex
    #####

    # test indicator is convex
    f = zero(PiecewiseQuadratic)
    @test(is_convex(f))

    # test single convex piece
    f = BoundedQuadratic(0.0, 1.0, 1.0, 0.0, 1.0)
    p = PiecewiseQuadratic(f)
    @test(is_convex(p))

    # test multi-piece convex
    f1 = BoundedQuadratic(-1.0, 0.0, 0.0, -1.0, 0.0)
    f2 = BoundedQuadratic(0.0, 1.0, 0.0, 0.0, 0.0)
    f3 = BoundedQuadratic(1.0, 2.0, 1.0, -2.0, 1.0)
    f4 = BoundedQuadratic(2.0, 3.0, 0.0, 3.0, -5.0)
    p = PiecewiseQuadratic([f1, f2, f3, f4])
    @test(is_convex(p))
end


text:
###Define functions to be used in testing below

###Test functions for TargetModel
function targetsin!(r::Vector,t::AbstractVector,paras::Vector)
    for i=1:length(t)
        r[i] = sin(paras[1]*t[i])
    end
    r
end

targetsin(t::AbstractVector,paras::Vector) = targetsin!(zeros(eltype(t),length(t)),t,paras)

###Test function for ODEModel
function odesin(t,y,ydot,paras)
    ydot[1] = paras[1]*cos(paras[1]*t)
end

#########################################################################################################

### Test model creation and evaluation functions

timepoints1 = linspace(0.0,10.0,100)
parameters1 = parameters([:a],[1.0],[5.0],[3.0])
measurements1 = data(:function,timepoints1,targetsin,timepoints1,[3.0])
measurements2 = data(:function,timepoints1,targetsin,timepoints1,[2.5])
measurements3 = data(:function,timepoints1,targetsin,timepoints1,[2.0])
noisemodel1 = noise(:gaussian,[0.01])
initial1 = [0.0]

model1 = model(:target,parameters1,measurements1,noisemodel1,targetsin;name="Test")
model2 = model(:target!,parameters1,measurements1,noisemodel1,targetsin!;name="Test!")
model3 = model(:ode,parameters1,measurements1,noisemodel1,odesin,initial1,1,[1];name="Test")

@test_approx_eq_eps evaluate!(model1,[3.0]) measurements(model1) 1e-4
@test_approx_eq_eps evaluate!(model2,[3.0]) measurements(model2) 1e-4
@test_approx_eq_eps evaluate!(model3,[3.0]) measurements(model3) 1e-4

#########################################################################################################

### Test the common initialize! function
s0 = samples(:base,1,1,Float64,Float64)
@test initialize!(trait(:initialize,:default),model1,s0).values == [3.0]
@test (srand(345) ; initialize!(trait(:initialize,:prior),model1,s0).values) == (srand(345) ; [rand(parameters1[1].prior)])

### Test the common interface functions
type TestModel <: AbstractModel
    parameters::Vector
end

@test numparas(TestModel(parameters1)) == 1
@test parameters(TestModel(parameters1)) == parameters1

@test_throws MethodError evaluate!(TestModel(parameters1),[1.0])
@test_throws MethodError dataindex(TestModel(parameters1))
@test_throws MethodError measurements(TestModel(parameters1))
@test_throws MethodError noisemodel(TestModel(parameters1))

########################################################################################################

### Test the geometry calculations
r1 = copy(evaluate!(model1,[3.0]))
r2 = copy(evaluate!(model2,[2.5]))
r3 = copy(evaluate!(model3,[2.0]))

@test_approx_eq loglikelihood(model1,r1) loglikelihood(noisemodel1,datavalues(measurements1),r1)
@test_approx_eq loglikelihood(model2,r2) loglikelihood(noisemodel1,datavalues(measurements1),r2)
@test_approx_eq loglikelihood(model3,r3) loglikelihood(noisemodel1,datavalues(measurements1),r3)

### Test the geometry function for zero-order samples
for m in [model1,model2,model3]
    s = samples(:base,1,3,Float64,Float64)
    copy!(s.values,[3.0,2.0,0.0])
    geometry!(m,s)
    ll = [loglikelihood(m,evaluate!(m,[3.0])) loglikelihood(m,evaluate!(m,[2.0])) -Inf]
    @test_approx_eq s.logprior logprior(parameters1,s.values,Float64)
    @test_approx_eq s.loglikelihood ll
end

###test the show functions
println()
println("====================")
println("Test show() function")
println("====================")
show(model1)
show(model2)
show(model3)
println("====================")
println("End  show() function")
println("====================")
println()


text:
using GeneralizedSDistributions
using Distributions
using Test

@testset "Implementation" begin

    D = Gamma(2, 2)
    G = GSDist(D)

    @testset "Parameter retrieval" begin
        @test_nowarn params(G)
    end

    @testset "Computation of statistics" begin
        @test_nowarn maximum(G)
        @test_nowarn minimum(G)
        @test_nowarn extrema(G)
        @test_nowarn mean(G)
        @test_nowarn var(G)
        @test_nowarn std(G)
        @test_nowarn median(G)
        @test_nowarn modes(G)
        @test_nowarn mode(G)

        @test_throws MethodError skewness(G)

        @test_throws MethodError kurtosis(G)
        @test_throws MethodError kurtosis(G, true)
        @test_throws MethodError kurtosis(G, false)

        @test_throws MethodError isplatykurtic(G)
        @test_throws MethodError isleptokurtic(G)
        @test_throws MethodError ismesokurtic(G)

        @test_throws MethodError entropy(G)

        @test_throws MethodError mgf(G, 3)

        @test_throws MethodError cf(G, 3)
    end

    @testset "Probability evaluation" begin
        x = mode(G)
        y = mode(G) + 1

        @test_nowarn insupport(G, x)

        @test_nowarn pdf(G, x)
        @test_nowarn logpdf(G, x)
        @test_nowarn loglikelihood(G, x)

        @test_nowarn cdf(G, x)
        @test_nowarn logcdf(G, x)
        @test_nowarn logdiffcdf(G, y, x)
        @test_nowarn ccdf(G, x)
        @test_nowarn logccdf(G, x)

        @test_nowarn quantile(G, 0.3)
        @test_nowarn cquantile(G, 0.3)

        @test_nowarn invlogcdf(G, log(0.3))
        @test_nowarn invlogccdf(G, log(0.3))
    end

    @testset "Sampling" begin
        @test_nowarn sampler(G)
        @test_nowarn rand(G, 10)
    end

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    NodeType

An enum describing the possible node types. Each [`Node`](@ref) has a `.index`
field, which should be interpreted as follows:

 * `NODE_CALL_MULTIVARIATE`: the index into `operators.multivariate_operators`
 * `NODE_CALL_UNIVARIATE`: the index into `operators.univariate_operators`
 * `NODE_LOGIC`: the index into `operators.logic_operators`
 * `NODE_COMPARISON`: the index into `operators.comparison_operators`
 * `NODE_MOI_VARIABLE`: the value of `MOI.VariableIndex(index)` in the user's
   space of the model.
 * `NODE_VARIABLE`: the 1-based index of the internal vector
 * `NODE_VALUE`: the index into the `.values` field of `Expression`
 * `NODE_PARAMETER`: the index into `data.parameters`
 * `NODE_SUBEXPRESSION`:  the index into `data.expressions`
"""
@enum(
    NodeType,
    # Index into multivariate operators
    NODE_CALL_MULTIVARIATE,
    # Index into univariate operators
    NODE_CALL_UNIVARIATE,
    # Index into logic operators
    NODE_LOGIC,
    # Index into comparison operators
    NODE_COMPARISON,
    # Index is the value of `MOI.VariableIndex`. This is from the original
    # model, and is not consecutive.
    NODE_MOI_VARIABLE,
    # Index of the internal, consecutive, and ordered `MOI.VariableIndex`.
    NODE_VARIABLE,
    # Index is into the list of constants
    NODE_VALUE,
    # Index is into the list of parameters
    NODE_PARAMETER,
    # Index is into the list of subexpressions
    NODE_SUBEXPRESSION,
)

"""
    struct Node
        type::NodeType
        index::Int
        parent::Int
    end

A single node in a nonlinear expression tree. Used by
[`Expression`](@ref).

See the MathOptInterface documentation for information on how the nodes and
values form an expression tree.
"""
struct Node
    type::NodeType
    index::Int
    parent::Int
end

"""
    struct Expression
        nodes::Vector{Node}
        values::Vector{Float64}
    end

The core type that represents a nonlinear expression. See the MathOptInterface
documentation for information on how the nodes and values form an expression
tree.
"""
struct Expression
    nodes::Vector{Node}
    values::Vector{Float64}
    Expression() = new(Node[], Float64[])
end

function Base.:(==)(x::Expression, y::Expression)
    return x.nodes == y.nodes && x.values == y.values
end

"""
    struct Constraint
        expression::Expression
        set::Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
            MOI.Interval{Float64},
        }
    end

A type to hold information relating to the nonlinear constraint `f(x) in S`,
where `f(x)` is defined by `.expression`, and `S` is `.set`.
"""
struct Constraint
    expression::Expression
    set::Union{
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
        MOI.Interval{Float64},
    }
end

"""
    ParameterIndex

An index to a nonlinear parameter that is returned by [`add_parameter`](@ref).
Given `data::Model` and `p::ParameterIndex`, use `data[p]` to retrieve
the current value of the parameter and `data[p] = value` to set a new value.
"""
struct ParameterIndex
    value::Int
end

"""
    ExpressionIndex

An index to a nonlinear expression that is returned by [`add_expression`](@ref).

Given `data::Model` and `ex::ExpressionIndex`, use `data[ex]` to
retrieve the corresponding [`Expression`](@ref).
"""
struct ExpressionIndex
    value::Int
end

"""
    ConstraintIndex

An index to a nonlinear constraint that is returned by [`add_constraint`](@ref).

Given `data::Model` and `c::ConstraintIndex`, use `data[c]` to
retrieve the corresponding [`Constraint`](@ref).
"""
struct ConstraintIndex
    value::Int
end

"""
    Model()

The core datastructure for representing a nonlinear optimization problem.

It has the following fields:
 * `objective::Union{Nothing,Expression}` : holds the nonlinear objective
   function, if one exists, otherwise `nothing`.
 * `expressions::Vector{Expression}` : a vector of expressions in the model.
 * `constraints::OrderedDict{ConstraintIndex,Constraint}` : a map from
   [`ConstraintIndex`](@ref) to the corresponding [`Constraint`](@ref). An
   `OrderedDict` is used instead of a `Vector` to support constraint deletion.
 * `parameters::Vector{Float64}` : holds the current values of the parameters.
 * `operators::OperatorRegistry` : stores the operators used in the model.
"""
mutable struct Model
    objective::Union{Nothing,Expression}
    expressions::Vector{Expression}
    constraints::OrderedDict{ConstraintIndex,Constraint}
    parameters::Vector{Float64}
    operators::OperatorRegistry
    # This is a private field, used only to increment the ConstraintIndex.
    last_constraint_index::Int64
    function Model()
        return new(
            nothing,
            Expression[],
            OrderedDict{ConstraintIndex,Constraint}(),
            Float64[],
            OperatorRegistry(),
            0,
        )
    end
end

"""
    AbstractAutomaticDifferentiation

An abstract type for extending [`Evaluator`](@ref).
"""
abstract type AbstractAutomaticDifferentiation end

"""
    Evaluator(
        model::Model,
        backend::AbstractAutomaticDifferentiation,
        ordered_variables::Vector{MOI.VariableIndex},
    )

Create `Evaluator`, a subtype of `MOI.AbstractNLPEvaluator`, from `Model`.
"""
mutable struct Evaluator{B} <: MOI.AbstractNLPEvaluator
    # The internal datastructure.
    model::Model
    # The abstract-differentiation backend
    backend::B
    # A vector to store expressions in Julia's Expr form. The eltype is Any
    # because expressions may be constants, and `:(1)` is just `1`.
    julia_expressions::Vector{Any}
    # ordered_constraints is needed because `OrderedDict` doesn't support
    # looking up a key by the linear index.
    ordered_constraints::Vector{ConstraintIndex}
    # Storage for the NLPBlockDual, so that we can query the dual of individual
    # constraints without needing to query the full vector each time.
    constraint_dual::Vector{Float64}
    # Timers
    eval_objective_timer::Float64
    eval_constraint_timer::Float64
    eval_objective_gradient_timer::Float64
    eval_constraint_jacobian_timer::Float64
    eval_hessian_lagrangian_timer::Float64

    function Evaluator(
        model::Model,
        backend::B = nothing,
    ) where {B<:Union{Nothing,MOI.AbstractNLPEvaluator}}
        return new{B}(
            model,
            backend,
            Expr[],
            ConstraintIndex[],
            Float64[],
            0.0,
            0.0,
            0.0,
            0.0,
            0.0,
        )
    end
end

_bound(s::MOI.LessThan) = MOI.NLPBoundsPair(-Inf, s.upper)
_bound(s::MOI.GreaterThan) = MOI.NLPBoundsPair(s.lower, Inf)
_bound(s::MOI.EqualTo) = MOI.NLPBoundsPair(s.value, s.value)
_bound(s::MOI.Interval) = MOI.NLPBoundsPair(s.lower, s.upper)

"""
    MOI.NLPBlockData(evaluator::Evaluator)

Create an [`MOI.NLPBlockData`](@ref) object from an [`Evaluator`](@ref)
object.
"""
function MOI.NLPBlockData(evaluator::Evaluator)
    return MOI.NLPBlockData(
        [_bound(c.set) for (_, c) in evaluator.model.constraints],
        evaluator,
        evaluator.model.objective !== nothing,
    )
end

"""
    ExprGraphOnly() <: AbstractAutomaticDifferentiation

The default implementation of `AbstractAutomaticDifferentiation`. The only
supported feature is `:ExprGraph`.
"""
struct ExprGraphOnly <: AbstractAutomaticDifferentiation end

function Evaluator(model::Model, ::ExprGraphOnly, ::Vector{MOI.VariableIndex})
    return Evaluator(model)
end

"""
    SparseReverseMode() <: AbstractAutomaticDifferentiation

An implementation of `AbstractAutomaticDifferentiation` that uses sparse
reverse-mode automatic differentiation to compute derivatives. Supports all
features in the MOI nonlinear interface.
"""
struct SparseReverseMode <: AbstractAutomaticDifferentiation end

function Evaluator(
    model::Model,
    ::SparseReverseMode,
    ordered_variables::Vector{MOI.VariableIndex},
)
    return Evaluator(model, ReverseAD.NLPEvaluator(model, ordered_variables))
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test,
    MatrixFactorizations, ArrayLayouts, LinearAlgebra, Random, LazyBandedMatrices, StaticArrays
import InfiniteLinearAlgebra: qltail, toeptail, tailiterate, tailiterate!, tail_de, ql_X!,
    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix, InfBandCartesianIndices,
    rightasymptotics, QLHessenberg, ConstRows, PertConstRows, chop, chop!,
    BandedToeplitzLayout, PertToeplitzLayout, TridiagonalToeplitzLayout, BidiagonalToeplitzLayout
import Base: BroadcastStyle
import BlockArrays: _BlockArray
import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
import MatrixFactorizations: QLPackedQ
import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
import LazyArrays: colsupport, MemoryLayout, ApplyLayout, LazyArrayStyle, arguments
import InfiniteArrays: OneToInf, oneto, RealInfinity
import LazyBandedMatrices: BroadcastBandedBlockBandedLayout, BroadcastBandedLayout, LazyBandedLayout

@testset "chop" begin
    a = randn(5)
    b = [a; zeros(5)]
    chop!(b, eps())
    @test b == a

    @test isempty(chop!([0]))

    A = randn(5, 5)
    @test chop([A zeros(5, 2); zeros(2, 5) zeros(2, 2)], eps()) == A

    c = PseudoBlockArray([randn(5); zeros(10)], (blockedrange(1:5),))
    d = chop!(c, 0)
    @test length(d) == 6
end

include("test_infconv.jl")
include("test_infbanded.jl")

@testset "∞-block arrays" begin
    @testset "fixed block size" begin
        k = Base.OneTo.(oneto(∞))
        n = Fill.(oneto(∞), oneto(∞))
        @test broadcast(length, k) ≡ map(length, k) ≡ OneToInf()
        @test broadcast(length, n) ≡ map(length, n) ≡ OneToInf()

        b = mortar(Fill([1, 2], ∞))
        @test blockaxes(b, 1) ≡ Block.(OneToInf())
        @test b[Block(5)] == [1, 2]
        @test b[Block.(2:∞)][Block.(2:10)] == b[Block.(3:11)]
        @test exp.(b)[Block.(2:∞)][Block.(2:10)] == exp.(b[Block.(3:11)])

        c = PseudoBlockArray(1:∞, Vcat(2, Fill(3, ∞)))
        @test c[Block.(2:∞)][Block.(2:10)] == c[Block.(3:11)]

        @test length(axes(b, 1)) ≡ ℵ₀
        @test last(axes(b, 1)) ≡ RealInfinity()
        @test Base.BroadcastStyle(typeof(b)) isa LazyArrayStyle{1}

        @test unitblocks(oneto(∞)) ≡ blockedrange(Ones{Int}(∞))
        @test unitblocks(2:∞) == 2:∞

        @test unitblocks(oneto(∞))[Block.(2:∞)] == 2:∞
    end

    @testset "1:∞ blocks" begin
        a = blockedrange(oneto(∞))
        @test axes(a, 1) == a
        o = Ones((a,))
        @test Base.BroadcastStyle(typeof(a)) isa LazyArrayStyle{1}
        b = exp.(a)
        @test axes(b, 1) == a
        @test o .* b isa typeof(b)
        @test b .* o isa typeof(b)
    end

    @testset "concat" begin
        a = unitblocks(1:∞)
        b = exp.(a)
        c = BlockBroadcastArray(vcat, a, b)
        @test length(c) == ∞
        @test blocksize(c) == (∞,)
        @test c[Block(5)] == [a[5], b[5]]

        A = unitblocks(BandedMatrix(0 => 1:∞, 1 => Fill(2.0, ∞), -1 => Fill(3.0, ∞)))
        B = BlockBroadcastArray(hvcat, 2, A, Zeros(axes(A)), Zeros(axes(A)), A)
        @test B[Block(3, 3)] == [A[3, 3] 0; 0 A[3, 3]]
        @test B[Block(3, 4)] == [A[3, 4] 0; 0 A[3, 4]]
        @test B[Block(3, 5)] == [A[3, 5] 0; 0 A[3, 5]]
        @test blockbandwidths(B) == (1, 1)
        @test subblockbandwidths(B) == (0, 0)
        @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix

        C = BlockBroadcastArray(hvcat, 2, A, A, A, A)
        @test C[Block(3, 3)] == fill(A[3, 3], 2, 2)
        @test C[Block(3, 4)] == fill(A[3, 4], 2, 2)
        @test C[Block(3, 5)] == fill(A[3, 5], 2, 2)
        @test blockbandwidths(C) == (1, 1)
        @test subblockbandwidths(C) == (1, 1)
        @test B[Block.(1:10), Block.(1:10)] isa BlockSkylineMatrix
    end

    @testset "KronTrav" begin
        Δ = BandedMatrix(1 => Ones(∞), -1 => Ones(∞)) / 2
        A = KronTrav(Δ - 2I, Eye(∞))
        @test axes(A, 1) isa InfiniteLinearAlgebra.OneToInfBlocks
        V = view(A, Block.(Base.OneTo(3)), Block.(Base.OneTo(3)))
        @test MemoryLayout(V) isa BlockBandedMatrices.BandedBlockBandedLayout

        u = A * [1; zeros(∞)]
        @test u[1:3] == A[1:3, 1]
        @test bandwidths(view(A, Block(1, 1))) == (1, 1)
    end

    @testset "triangle recurrences" begin
        @testset "n and k" begin
            n = mortar(Fill.(oneto(∞), oneto(∞)))
            k = mortar(Base.OneTo.(oneto(∞)))

            @test n[Block(5)] ≡ layout_getindex(n, Block(5)) ≡ view(n, Block(5)) ≡ Fill(5, 5)
            @test k[Block(5)] ≡ layout_getindex(k, Block(5)) ≡ view(k, Block(5)) ≡ Base.OneTo(5)
            @test Base.BroadcastStyle(typeof(n)) isa LazyArrays.LazyArrayStyle{1}
            @test Base.BroadcastStyle(typeof(k)) isa LazyArrays.LazyArrayStyle{1}

            N = 1000
            v = view(n, Block.(Base.OneTo(N)))
            @test view(v, Block(2)) ≡ Fill(2, 2)
            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
            @test @allocated(axes(v)) ≤ 40

            dest = PseudoBlockArray{Float64}(undef, axes(v))
            @test copyto!(dest, v) == v
            @test @allocated(copyto!(dest, v)) ≤ 40

            v = view(k, Block.(Base.OneTo(N)))
            @test view(v, Block(2)) ≡ Base.OneTo(2)
            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
            @test @allocated(axes(v)) ≤ 40
            @test copyto!(dest, v) == v

            @testset "stack overflow" begin
                i = Base.to_indices(k, (Block.(2:∞),))[1].indices
                last(i)
            end

            v = view(k, Block.(2:∞))
            @test Base.BroadcastStyle(typeof(v)) isa LazyArrayStyle{1}
            @test v[Block(1)] == 1:2
            @test v[Block(1)] ≡ k[Block(2)] ≡ Base.OneTo(2)

            @test axes(n, 1) isa BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,OneToInf{Int64}}}
        end

        @testset "BlockHcat copyto!" begin
            n = mortar(Fill.(oneto(∞), oneto(∞)))
            k = mortar(Base.OneTo.(oneto(∞)))

            a = b = c = 0.0
            dat = BlockHcat(
                BroadcastArray((n, k, b, bc1) -> (k + b - 1) * (n + k + bc1) / (2k + bc1), n, k, b, b + c - 1),
                BroadcastArray((n, k, abc, bc, bc1) -> (n + k + abc) * (k + bc) / (2k + bc1), n, k, a + b + c, b + c, b + c - 1)
            )
            N = 1000
            KR = Block.(Base.OneTo(N))
            V = view(dat, Block.(Base.OneTo(N)), :)
            @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(hcat)}
            @test PseudoBlockArray(V)[Block.(1:5), :] == dat[Block.(1:5), :]
            V = view(dat', :, Block.(Base.OneTo(N)))
            @test MemoryLayout(V) isa LazyArrays.ApplyLayout{typeof(vcat)}
            a = dat.arrays[1]'
            N = 100
            KR = Block.(Base.OneTo(N))
            v = view(a, :, KR)
            @time r = PseudoBlockArray(v)
            @test v == r
        end

        @testset "BlockBanded" begin
            a = b = c = 0.0
            n = mortar(Fill.(oneto(∞), oneto(∞)))
            k = mortar(Base.OneTo.(oneto(∞)))
            Dy = BlockBandedMatrices._BandedBlockBandedMatrix((k .+ (b + c))', axes(k, 1), (-1, 1), (-1, 1))
            N = 100
            @test Dy[Block.(1:N), Block.(1:N)] == BlockBandedMatrices._BandedBlockBandedMatrix((k.+(b+c))[Block.(1:N)]', axes(k, 1)[Block.(1:N)], (-1, 1), (-1, 1))
        end

        @testset "Symmetric" begin
            k = mortar(Base.OneTo.(oneto(∞)))
            n = mortar(Fill.(oneto(∞), oneto(∞)))

            dat = BlockHcat(
                BlockBroadcastArray(hcat, float.(k), Zeros((axes(n, 1),)), float.(n)),
                Zeros((axes(n, 1), Base.OneTo(3))),
                Zeros((axes(n, 1), Base.OneTo(3))))
            M = BlockBandedMatrices._BandedBlockBandedMatrix(dat', axes(k, 1), (1, 1), (1, 1))
            Ms = Symmetric(M)
            @test blockbandwidths(M) == (1, 1)
            @test blockbandwidths(Ms) == (1, 1)
            @test Ms[Block.(1:5), Block.(1:5)] == Symmetric(M[Block.(1:5), Block.(1:5)])
            @test Ms[Block.(1:5), Block.(1:5)] isa BandedBlockBandedMatrix

            b = [ones(10); zeros(∞)]
            @test (Ms * b)[Block.(1:6)] == Ms[Block.(1:6), Block.(1:4)]*ones(10)
            @test ((Ms * Ms) * b)[Block.(1:6)] == (Ms * (Ms * b))[Block.(1:6)]
            @test ((Ms + Ms) * b)[Block.(1:6)] == (2*(Ms * b))[Block.(1:6)]

            dat = BlockBroadcastArray(hcat, float.(k), Zeros((axes(n, 1),)), float.(n))
            M = BlockBandedMatrices._BandedBlockBandedMatrix(dat', axes(k, 1), (-1, 1), (1, 1))
            Ms = Symmetric(M)
            @test Symmetric((M+M)[Block.(1:10), Block.(1:10)]) == (Ms+Ms)[Block.(1:10), Block.(1:10)]
        end
    end

    @testset "blockdiag" begin
        D = Diagonal(mortar(Fill.((-(0:∞) - (0:∞) .^ 2), 1:2:∞)))
        x = [randn(5); zeros(∞)]
        x̃ = PseudoBlockArray(x, (axes(D, 1),))
        @test (D*x)[1:10] == (D*x̃)[1:10]
    end

    @testset "sortedunion" begin
        a = cumsum(1:2:∞)
        @test BlockArrays.sortedunion(a, a) ≡ a
        @test BlockArrays.sortedunion([∞], a) ≡ BlockArrays.sortedunion(a, [∞]) ≡ a
        @test BlockArrays.sortedunion([∞], [∞]) == [∞]

        b = Vcat([1, 2], 3:∞)
        c = Vcat(1, 3:∞)
        @test BlockArrays.sortedunion(b, b) ≡ b
        @test BlockArrays.sortedunion(c, c) ≡ c
    end
end

@testset "Algebra" begin
    @testset "BandedMatrix" begin
        A = BandedMatrix(-3 => Fill(7 / 10, ∞), -2 => 1:∞, 1 => Fill(2im, ∞))
        @test A isa BandedMatrix{ComplexF64}
        @test A[1:10, 1:10] == diagm(-3 => Fill(7 / 10, 7), -2 => 1:8, 1 => Fill(2im, 9))

        A = BandedMatrix(0 => Vcat([1, 2, 3], Zeros(∞)), 1 => Vcat(1, Zeros(∞)))
        @test A[1, 2] == 1

        A = BandedMatrix(-3 => Fill(7 / 10, ∞), -2 => Fill(1, ∞), 1 => Fill(2im, ∞))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]'
        @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10])

        A = BandedMatrix(-1 => Vcat(Float64[], Fill(1 / 4, ∞)), 0 => Vcat([1.0 + im], Fill(0, ∞)), 1 => Vcat(Float64[], Fill(1, ∞)))
        @test MemoryLayout(typeof(view(A.data, :, 1:10))) == ApplyLayout{typeof(hcat)}()
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]'
        @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10])

        A = BandedMatrix(-2 => Vcat(Float64[], Fill(1 / 4, ∞)), 0 => Vcat([1.0 + im, 2, 3], Fill(0, ∞)), 1 => Vcat(Float64[], Fill(1, ∞)))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10, 1:10] ≈ (A')[1:10, 1:10] ≈ A[1:10, 1:10]'
        @test At[1:10, 1:10] ≈ transpose(A)[1:10, 1:10] ≈ transpose(A[1:10, 1:10])

        A = _BandedMatrix(Fill(1, 4, ∞), ℵ₀, 1, 2)
        @test A^2 isa BandedMatrix
        @test (A^2)[1:10, 1:10] == (A*A)[1:10, 1:10] == (A[1:100, 1:100]^2)[1:10, 1:10]
        @test A^3 isa ApplyMatrix{<:Any,typeof(*)}
        @test (A^3)[1:10, 1:10] == (A*A*A)[1:10, 1:10] == ((A*A)*A)[1:10, 1:10] == (A*(A*A))[1:10, 1:10] == (A[1:100, 1:100]^3)[1:10, 1:10]

        @testset "∞ x finite" begin
            A = BandedMatrix(1 => 1:∞) + BandedMatrix(-1 => Fill(2, ∞))
            B = _BandedMatrix(randn(3, 5), ℵ₀, 1, 1)

            @test lmul!(2.0, copy(B)')[:, 1:10] == (2B')[:, 1:10]

            @test_throws ArgumentError BandedMatrix(A)
            @test A * B isa MulMatrix
            @test B'A isa MulMatrix

            @test all(diag(A[1:6, 1:6]) .=== zeros(6))

            @test (A*B)[1:7, 1:5] ≈ A[1:7, 1:6] * B[1:6, 1:5]
            @test (B'A)[1:5, 1:7] ≈ (B')[1:5, 1:6] * A[1:6, 1:7]
        end
    end

    @testset "BlockTridiagonal" begin
        A = BlockTridiagonal(Vcat([fill(1.0, 2, 1), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)),
            Vcat([zeros(1, 1)], Fill(zeros(2, 2), ∞)),
            Vcat([fill(1.0, 1, 2), Matrix(1.0I, 2, 2)], Fill(Matrix(1.0I, 2, 2), ∞)))

        @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz
        @test isblockbanded(A)

        @test A[Block.(1:2), Block(1)] == A[1:3, 1:1] == reshape([0.0, 1.0, 1.0], 3, 1)

        @test BlockBandedMatrix(A)[1:100, 1:100] == BlockBandedMatrix(A, (2, 1))[1:100, 1:100] == BlockBandedMatrix(A, (1, 1))[1:100, 1:100] == A[1:100, 1:100]

        @test (A-I)[1:100, 1:100] == A[1:100, 1:100] - I
        @test (A+I)[1:100, 1:100] == A[1:100, 1:100] + I
        @test (I+A)[1:100, 1:100] == I + A[1:100, 1:100]
        @test (I-A)[1:100, 1:100] == I - A[1:100, 1:100]

        @test (A-im*I)[1:100, 1:100] == A[1:100, 1:100] - im * I
        @test (A+im*I)[1:100, 1:100] == A[1:100, 1:100] + im * I
        @test (im*I+A)[1:100, 1:100] == im * I + A[1:100, 1:100]
        @test (im*I-A)[1:100, 1:100] == im * I - A[1:100, 1:100]

        T = mortar(LazyBandedMatrices.Tridiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞)))
        #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
        @test T[Block(2, 2)] == [1 2; 3 4]
        @test_broken T[Block(1, 3)] == Zeros(2, 2)
    end

    @testset "BlockBidiagonal" begin
        B = mortar(LazyBandedMatrices.Bidiagonal(Fill([1 2; 3 4], ∞), Fill([1 2; 3 4], ∞), :U))
        #TODO: copy BlockBidiagonal code from BlockBandedMatrices to LazyBandedMatrices
        @test B[Block(2, 3)] == [1 2; 3 4]
        @test_broken B[Block(1, 3)] == Zeros(2, 2)
    end

    @testset "Fill" begin
        A = _BandedMatrix(Ones(1, ∞), ℵ₀, -1, 1)
        @test 1.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test Zeros(∞) .* A ≡ Zeros(∞, ∞) .* A ≡ A .* Zeros(1, ∞) ≡ A .* Zeros(∞, ∞) ≡ Zeros(∞, ∞)
        @test Ones(∞) .* A isa BandedMatrix{Float64,<:Ones}
        @test A .* Ones(1, ∞) isa BandedMatrix{Float64,<:Ones}
        @test 2.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Fill}
        @test Eye(∞) * A isa BandedMatrix{Float64,<:Ones}
        @test A * Eye(∞) isa BandedMatrix{Float64,<:Ones}

        @test A * A isa BandedMatrix
        @test (A*A)[1:10, 1:10] == BandedMatrix(2 => Ones(8))

        Ã = _BandedMatrix(Fill(1, 1, ∞), ℵ₀, -1, 1)
        @test A * Ã isa BandedMatrix
        @test Ã * A isa BandedMatrix
        @test Ã * Ã isa BandedMatrix

        B = _BandedMatrix(Ones(1, 10), ℵ₀, -1, 1)
        C = _BandedMatrix(Ones(1, 10), 10, -1, 1)
        D = _BandedMatrix(Ones(1, ∞), 10, -1, 1)

        @test (A*B)[1:10, 1:10] == (B*C)[1:10, 1:10] == (D*A)[1:10, 1:10] == D * B == (C*D)[1:10, 1:10] == BandedMatrix(2 => Ones(8))
    end

    @testset "Banded Broadcast" begin
        A = _BandedMatrix((1:∞)', ℵ₀, -1, 1)
        @test 2.0 .* A isa BandedMatrix{Float64,<:Adjoint}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Adjoint}
        @test Eye(∞) * A isa BandedMatrix{Float64,<:Adjoint}
        @test A * Eye(∞) isa BandedMatrix{Float64,<:Adjoint}
        A = _BandedMatrix(Vcat((1:∞)', Ones(1, ∞)), ℵ₀, 0, 1)
        @test 2.0 .* A isa BandedMatrix
        @test A .* 2.0 isa BandedMatrix
        @test Eye(∞) * A isa BandedMatrix
        @test A * Eye(∞) isa BandedMatrix
        b = 1:∞
        @test BroadcastStyle(typeof(b)) isa LazyArrayStyle{1}
        @test BroadcastStyle(typeof(A)) isa BandedStyle
        @test BroadcastStyle(LazyArrayStyle{1}(), BandedStyle()) isa LazyArrayStyle{2}
        @test BroadcastStyle(LazyArrayStyle{2}(), BandedStyle()) isa LazyArrayStyle{2}
        @test bandwidths(b .* A) == (0, 1)

        @test colsupport(b .* A, 1) == 1:1
        @test Base.replace_in_print_matrix(b .* A, 2, 1, "0.0") == " ⋅ "
        @test bandwidths(A .* b) == (0, 1)
        @test A .* b' isa BroadcastArray
        @test bandwidths(A .* b') == bandwidths(A .* b')
        @test colsupport(A .* b', 3) == 2:3

        A = _BandedMatrix(Ones{Int}(1, ∞), ℵ₀, 0, 0)'
        B = _BandedMatrix((-2:-2:-∞)', ℵ₀, -1, 1)
        C = Diagonal(2 ./ (1:2:∞))
        @test bandwidths(A * (B * C)) == (-1, 1)
        @test bandwidths((A * B) * C) == (-1, 1)

        A = _BandedMatrix(Ones{Int}(1, ∞), ℵ₀, 0, 0)'
        B = _BandedMatrix((-2:-2:-∞)', ℵ₀, -1, 1)
        @test MemoryLayout(A + B) isa BroadcastBandedLayout{typeof(+)}
        @test MemoryLayout(2 * (A + B)) isa BroadcastBandedLayout{typeof(*)}
        @test bandwidths(A + B) == (0, 1)
        @test bandwidths(2 * (A + B)) == (0, 1)
    end

    @testset "Triangle OP recurrences" begin
        k = mortar(Base.OneTo.(1:∞))
        n = mortar(Fill.(1:∞, 1:∞))
        @test k[Block.(2:3)] isa BlockArray
        @test n[Block.(2:3)] isa BlockArray
        @test k[Block.(2:3)] == [1, 2, 1, 2, 3]
        @test n[Block.(2:3)] == [2, 2, 3, 3, 3]
        @test blocksize(BroadcastVector(exp, k)) == (ℵ₀,)
        @test BroadcastVector(exp, k)[Block.(2:3)] == exp.([1, 2, 1, 2, 3])
        # BroadcastVector(+,k,n)
    end
    # Multivariate OPs Corollary (3)
    # n = 5
    # BlockTridiagonal(Zeros.(1:∞,2:∞),
    #         (n -> Diagonal(((n+2).+(0:n)))/ (2n + 2)).(0:∞),
    #         Zeros.(2:∞,1:∞))

    @testset "KronTrav" begin
        Δ = BandedMatrix(1 => Ones(∞) / 2, -1 => Ones(∞))
        A = KronTrav(Δ, Eye(∞))
        @test A[Block(100, 101)] isa BandedMatrix
        @test A[Block(100, 100)] isa BandedMatrix
        @test A[Block.(1:5), Block.(1:5)] isa BandedBlockBandedMatrix
        B = KronTrav(Eye(∞), Δ)
        @test B[Block(100, 101)] isa BandedMatrix
        @test B[Block(100, 100)] isa BandedMatrix
        V = view(A + B, Block.(1:5), Block.(1:5))
        @test MemoryLayout(typeof(V)) isa BroadcastBandedBlockBandedLayout{typeof(+)}
        @test arguments(V) == (A[Block.(1:5), Block.(1:5)], B[Block.(1:5), Block.(1:5)])
        @test (A+B)[Block.(1:5), Block.(1:5)] == A[Block.(1:5), Block.(1:5)] + B[Block.(1:5), Block.(1:5)]

        @test blockbandwidths(A + B) == (1, 1)
        @test blockbandwidths(2A) == (1, 1)
        @test blockbandwidths(2 * (A + B)) == (1, 1)

        @test subblockbandwidths(A + B) == (1, 1)
        @test subblockbandwidths(2A) == (1, 1)
        @test subblockbandwidths(2 * (A + B)) == (1, 1)
    end
end

include("test_hessenbergq.jl")
include("test_infql.jl")
include("test_infqr.jl")
include("test_inful.jl")
include("test_infcholesky.jl")


text:
#
println("Da) Test of the PhotoEmission module with ASF from an internally generated initial- and final-state multiplet.")
#
# setDefaults("unit: rate", "a.u.")
setDefaults("unit: rate", "1/s")
setDefaults("print summary: open", "zzz-radiative.sum")
#== wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                        initialConfigs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                        finalConfigs  =[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                        process = JAC.Radiative,
                        processSettings=PhotoEmission.Settings([E1, M1, E2, M2], [JAC.UseCoulomb, JAC.UseBabushkin], true, true,
                           true, Tuple{Int64,Int64}[(5,0), (7,0), (10,0), (11,0), (12,0), (13,0), (14,0), (15,0), (16,0)], 0., 0., 10000. ) )  ==#

wa = Atomic.Computation("Energies and Einstein coefficients for the spectrum Fe X",  Nuclear.Model(26.);
                        initialConfigs = [Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],
                        finalConfigs   = [Configuration("[Ne] 3s^2 3p^5")],
                        process = Radiative,
                        processSettings = PhotoEmission.Settings([E1, M1, E2, M2], [JAC.UseCoulomb, JAC.UseBabushkin],
                                          false, true, false, Tuple{Int64,Int64}[], 0., 0., 10000. ) );
wb = @time( perform(wa) )
setDefaults("print summary: close", "")

## , Configuration("[Ne] 3s^2 3p^4 3d")  Configuration("[Ne] 3s 3p^6")


text:
using SDEModels
using StaticArrays
using Test
import Random: seed!
import LinearAlgebra: norm

ex1 = :(dX = Y*a*dt + b*dW1 + dW2)

@testset "codegen" begin
  @test collect(SDEModels.symbols(ex1)) == [:dX, :Y, :a, :dt, :b, :dW1, :dW2]
  SDEModels.replace_symbols!(ex1, Dict(:dX => :(sin(0)), :dt => :t))
  @test collect(SDEModels.symbols(ex1)) == [:Y, :a, :t, :b, :dW1, :dW2]
  @test SDEModels.cat_expressions([:(a+b), :(c-d)]) == :(StaticArrays.SVector{2}(a+b,c-d))
  @test SDEModels.cat_expressions([:(a+b) :(c-d); :e 1]) == :(StaticArrays.SMatrix{2,2}(a+b, e, c-d, 1))
end

@sde_model TestModel begin
    dX = Y*a*dt + b*dW1 + d*dW2
    dY = Z*c*dt +   dW2 +   dW3
    dZ = X*e*dt + f*dW3 -   dW4
end

@sde_model Wiener dX = dW
@sde_model Deterministic dX = X*dt
@sde_model TimeDependent dX = t*a*X*dt + b*X*dW

@sde_model SpecialCaseModel1 begin
  dx = dw1 + dw2
end

@sde_model SpecialCaseModel2 begin
  dx = dw1
  dy = dw1
end

@sde_model SpecialCaseModel3 begin
  dx = dw1
  dy = dw1 + dw2
  dz = dw2
end

m1 = BlackScholes(0.01, 0.3)
m2 = Heston(0.01, 10.0, 0.3, 0.1, 0.4)
m3 = TestModel(1, 2, 3, 4, 5, 6)
wiener = Wiener()
deterministic = Deterministic()
timedependend = TimeDependent(1, 2)

@testset "model creation" begin
  @test fieldnames(BlackScholes) == (:r, :σ)
  @test fieldnames(Heston) == (:r, :κ, :θ, :σ, :ρ)
  @test fieldnames(TestModel) == (:a, :c, :e, :b, :d, :f)
  @test dim(m1) == (1,1)
  @test dim(m2) == (2,2)
  @test dim(m3) == (3,4)
  @test dim(wiener) == (1,1)
  @test dim(deterministic) == (1,0)
end

@testset "drift" begin
  t0 = 0.0
  t1 = 1.0
  s2 = [100.0, 0.4]
  @test drift(m1, t0, 100.0) ≈ 100.0 * 0.01
  @test drift(m2, t0, s2) ≈ [0.01*100.0, 10.0*(0.3-0.4)]
  @test drift(m3, t0, [0.1, 0.2, 0.3]) ≈ [0.2*1, 0.3*2, 0.1*3]
  @test drift(wiener, t0, 100.0) ≈ 0.0
  @test drift(deterministic, t0, 100.0) ≈ 100.0
  @test drift(timedependend, t0, 3.0) ≠ drift(timedependend, t1, 3.0)
  @test drift(timedependend, 4.0, 3.0) ≈ 12
  @test corrected_drift(m2, t0, s2, 0.65) ≈ drift(m2, t0, s2) - [s2[2]*s2[1]+m2.σ*m2.ρ*s2[1]/2; m2.σ^2/2] * 0.65
end

@testset "diffusion" begin
  t0 = 0.0
  @test diffusion(m1, t0, 100.0) ≈ 100.0 * 0.3
  @test diffusion(m2, t0, [100.0, 0.4]) ≈ [100.0*sqrt(0.4) 0.0; 0.4*0.1*sqrt(0.4) sqrt(1-0.4^2)*0.1*sqrt(0.4)]
  @test diffusion(m3, t0, [0.1, 0.2, 0.3]) ≈ [4 5 0 0; 0 1 1 0; 0 0 6 -1]
  @test diffusion(wiener, t0, 100.0) ≈ 1.0
  @test diffusion(deterministic, t0, 100.0) ≈ 0.0
end

@testset "simulation" begin
  t0 = 0.0
  @test size(simulate(m1, EulerMaruyama(0.01), t0, 100.0, 1000)[1]) == (1001,)
  @test size(simulate(m1, Milstein(0.01), t0, 100.0, 1000)[1]) == (1001,)
  @test size(simulate(m2, EulerMaruyama(0.01), t0, [100.0, 0.4], 500)[1]) == (2, 501)
  @test sample(Deterministic(), EulerMaruyama(1), t0, 1.0, 1) ≈ 2
  # @test sample(Deterministic(), Milstein(1), t0, 1.0, 1) ≈ 2
  @test simulate(deterministic, EulerMaruyama(1.0), t0, 1.0, 5)[1] ≈ [1, 2, 4, 8, 16, 32]
  # @test simulate(deterministic, Milstein(1.0), t0, 1.0, 5)[1] ≈ [1, 2, 4, 8, 16, 32]

  @test size(simulate(SpecialCaseModel1(), EulerMaruyama(0.01), t0, [100.0], 1000)[1]) == (1,1001,)
  @test size(simulate(SpecialCaseModel2(), EulerMaruyama(0.01), t0, [100.0, 100.0], 1000)[1]) == (2,1001)
  @test size(simulate(SpecialCaseModel3(), EulerMaruyama(0.01), t0, [100.0, 100.0, 100.0], 1000)[1]) == (3,1001)
end

@testset "multilevel" begin
  substeps = 4
  coarse = EulerMaruyama(0.1)
  fine = subdivide(coarse, substeps)
  x0 = [100.0, 0.4]
  seed!(0)
  xc, xf, tc, tf = SDEModels._multilevel_simulate(m2, coarse, fine, substeps, 0.0, SVector(x0...), 10, 100)
  seed!(0)
  x, t = simulate(m2, fine, 0.0, SVector(x0...), 40, 100)
  @test t == tf
  @test isapprox(sum(norm.(xf - x).^2), 0, atol=1e-16)
  seed!(0)
  yc, yf = SDEModels._multilevel_sample(m2, coarse, fine, substeps, 0.0, SVector(x0...), 10, 100)
  @test isapprox(sum(norm.(xf[end,:] - yf).^2), 0, atol=1e-16)
  @test isapprox(sum(norm.(xc[end,:] - yc).^2), 0, atol=1e-16)

  ml_scheme = MultilevelScheme(coarse, 4, 0:4)
  seed!(0)
  ml_x1, ml_t = simulate(m2, ml_scheme, 0.0, x0, 1, 1)
  @test length(unique(extrema.(ml_t))) == 1
  seed!(0)
  ml_x2 = sample(m2, ml_scheme, 0.0, x0, 1, 1)
  @test isapprox(sum(norm.(last.(ml_x1) - last.(ml_x2))), 0, atol=1e-13)
end