text:
const HASH_FILENAME = "check.sha256"

"""
Return a decoded JSON file.
"""
function read_json(filename::AbstractString)
    return open(filename) do io
        return JSON.parse(io)
    end
end

"""
Return the SHA 256 hash of a file.
"""
function compute_sha256(filename::AbstractString)
    return open(filename) do io
        return bytes2hex(SHA.sha256(io))
    end
end

"""
Return the key for the given value
"""
function find_key_with_value(d, value)
    for (k, v) in d
        v == value && return k
    end
    error("dict does not have value == $value")
end

function compute_file_hash(path::String, files::Vector{String})
    data = Dict("files" => [])
    for file in files
        file_path = joinpath(path, file)
        file_info = Dict("filename" => file_path, "hash" => compute_sha256(file_path))
        push!(data["files"], file_info)
    end

    open(joinpath(path, HASH_FILENAME), "w") do io
        write(io, JSON.json(data))
    end
end

function read_file_hashes(path)
    data = open(joinpath(path, HASH_FILENAME), "r") do io
        JSON.parse(io)
    end

    return data["files"]
end

function check_kwargs(input_kwargs, valid_set::Array{Symbol}, function_name::String)
    if isempty(input_kwargs)
        return
    else
        for (key, value) in input_kwargs
            if !(key in valid_set)
                throw(ArgumentError("keyword argument $(key) is not a valid input for $(function_name)"))
            end
        end
    end
    return
end

# writing a dictionary of dataframes to files

function write_data(vars_results::Dict, save_path::String; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if :duals in keys(kwargs)
        name = "dual_"
    elseif :params in keys(kwargs)
        name = "parameter_"
    else
        name = ""
    end
    if file_type == Feather || file_type == CSV
        for (k, v) in vars_results
            file_path = joinpath(save_path, "$name$k.$(lowercase("$file_type"))")
            if isempty(vars_results[k])
                @debug "$name$k is empty, not writing $file_path"
            else
                file_type.write(file_path, vars_results[k])
            end
        end
    end
end

# writing a dictionary of dataframes to files and appending the time

function write_data(
    vars_results::Dict,
    time::DataFrames.DataFrame,
    save_path::AbstractString;
    kwargs...,
)
    file_type = get(kwargs, :file_type, Feather)
    for (k, v) in vars_results
        var = DataFrames.DataFrame()
        if file_type == CSV && size(time, 1) == size(v, 1)
            var = hcat(time, v)
        else
            var = v
        end
        file_path = joinpath(save_path, "$(k).$(lowercase("$file_type"))")
        file_type.write(file_path, var)
    end
end

function write_data(
    data::DataFrames.DataFrame,
    save_path::AbstractString,
    file_name::String;
    kwargs...,
)
    if isfile(save_path)
        save_path = dirname(save_path)
    end
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        file_path = joinpath(save_path, "$(file_name).$(lowercase("$file_type"))")
        file_type.write(file_path, data)
    end
    return
end

function write_optimizer_log(optimizer_log::Dict, save_path::AbstractString)
    JSON.write(joinpath(save_path, "optimizer_log.json"), JSON.json(optimizer_log))
end

function write_data(base_power::Float64, save_path::String)
    JSON.write(joinpath(save_path, "base_power.json"), JSON.json(base_power))
end

function _jump_value(input::JuMP.VariableRef)
    return JuMP.value(input)
end

function _jump_value(input::PJ.ParameterRef)
    return PJ.value(input)
end

function _jump_value(input::JuMP.ConstraintRef)
    return JuMP.dual(input)
end

function axis_array_to_dataframe(input_array::JuMP.Containers.DenseAxisArray{Float64})
    if length(axes(input_array)) == 1
        result = Vector{Float64}(undef, length(first(input_array.axes)))

        for t in input_array.axes[1]
            result[t] = input_array[t]
        end

        return DataFrames.DataFrame(var = result)

    elseif length(axes(input_array)) == 2
        result = Array{Float64, length(input_array.axes)}(
            undef,
            length(input_array.axes[2]),
            length(input_array.axes[1]),
        )
        names = Array{Symbol, 1}(undef, length(input_array.axes[1]))

        for t in input_array.axes[2], (ix, name) in enumerate(input_array.axes[1])
            result[t, ix] = input_array[name, t]
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(input_array)) == 3
        extra_dims = sum(length(axes(input_array)[2:(end - 1)]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(input_array)[1]))

        for i in input_array.axes[2]
            third_dim = collect(fill(i, size(input_array)[end]))
            result = Array{Float64, 2}(
                undef,
                length(last(input_array.axes)),
                length(first(input_array.axes)),
            )
            for t in last(input_array.axes),
                (ix, name) in enumerate(first(input_array.axes))

                result[t, ix] = input_array[name, i, t]
            end
            res = DataFrames.DataFrame(hcat(third_dim, result))
            result_df = vcat(result_df, res)
        end

        return DataFrames.names!(result_df, names)

    else
        error("Dimension Number $(length(axes(input_array))) not Supported")
    end
end

function axis_array_to_dataframe(input_array::JuMP.Containers.DenseAxisArray{})
    if length(axes(input_array)) == 1
        result = Vector{Float64}(undef, length(first(input_array.axes)))
        for t in input_array.axes[1]
            result[t] = _jump_value(input_array[t])
        end

        return DataFrames.DataFrame(var = result)
    elseif length(axes(input_array)) == 2
        result = Array{Float64, length(input_array.axes)}(
            undef,
            length(input_array.axes[2]),
            length(input_array.axes[1]),
        )
        names = Array{Symbol, 1}(undef, length(input_array.axes[1]))

        for t in input_array.axes[2], (ix, name) in enumerate(input_array.axes[1])
            result[t, ix] = _jump_value(input_array[name, t])
            names[ix] = Symbol(name)
        end

        return DataFrames.DataFrame(result, names)

    elseif length(axes(input_array)) == 3
        extra_dims = sum(length(axes(input_array)[2:(end - 1)]))
        extra_vars = [Symbol("S$(s)") for s in 1:extra_dims]
        result_df = DataFrames.DataFrame()
        names = vcat(extra_vars, Symbol.(axes(input_array)[1]))

        for i in input_array.axes[2]
            third_dim = collect(fill(i, size(input_array)[end]))
            result = Array{Float64, 2}(
                undef,
                length(last(input_array.axes)),
                length(first(input_array.axes)),
            )
            for t in last(input_array.axes),
                (ix, name) in enumerate(first(input_array.axes))

                result[t, ix] = _jump_value(input_array[name, i, t])
            end
            res = DataFrames.DataFrame(hcat(third_dim, result))
            result_df = vcat(result_df, res)
        end
        return DataFrames.names!(result_df, names)
    else
        @warn("Dimension Number $(length(axes(input_array))) not Supported, returning empty DataFrame")
        return DataFrames.DataFrame()
    end
end

function axis_array_to_dataframe(input_array::JuMP.Containers.SparseAxisArray)
    column_names = unique([(k[1], k[3]) for k in keys(input_array.data)])
    result_df = DataFrames.DataFrame()
    array_values = Vector{Vector{Float64}}(undef, length(column_names))
    for (ix, col) in enumerate(column_names)
        res = values(filter(v -> first(v)[[1, 3]] == col, input_array.data))
        array_values[ix] = PSI._jump_value.(res)
    end
    return DataFrames.DataFrame(array_values, Symbol.(column_names))
end

# this ensures that the time_stamp is not double shortened
function find_var_length(es::Dict, e_list::Array)
    return size(es[Symbol(splitext(e_list[1])[1])], 1)
end

function shorten_time_stamp(time::DataFrames.DataFrame)
    time = time[1:(size(time, 1) - 1), :]
    return time
end

""" Returns the correct container spec for the selected type of JuMP Model"""
function container_spec(m::M, axs...) where {M <: JuMP.AbstractModel}
    return JuMP.Containers.DenseAxisArray{JuMP.variable_type(m)}(undef, axs...)
end

""" Returns the correct container spec for the selected type of JuMP Model"""
function sparse_container_spec(m::M, axs...) where {M <: JuMP.AbstractModel}
    indexes = Base.Iterators.product(axs...)
    contents = Dict{eltype(indexes), Any}(indexes .=> 0)
    return JuMP.Containers.SparseAxisArray(contents)
end

function middle_rename(original::Symbol, split_char::String, addition::String)
    parts = split(String(original), split_char)
    parts[1] = parts[1] * "_" * addition
    return Symbol(join(parts, split_char))
end

"Replaces the string in `char` with the string`replacement`"
function replace_chars(s::String, char::String, replacement::String)
    return replace(s, Regex("[$char]") => replacement)
end

"Removes the string `char` from the original string"
function remove_chars(s::String, char::String)
    return replace_chars(s::String, char::String, "")
end

function get_available_components(::Type{T}, sys::PSY.System) where {T <: PSY.Component}
    return PSY.get_components(T, sys, x -> PSY.get_available(x))
end

function get_available_components(
    ::Type{PSY.RegulationDevice{T}},
    sys::PSY.System,
) where {T <: PSY.Component}
    return PSY.get_components(
        PSY.RegulationDevice{T},
        sys,
        x -> (PSY.get_available(x) && PSY.has_service(x, PSY.AGC)),
    )
end


text:
using ReusePatterns
using Test

supertypes(Int)
supertypes(String)
supertypes(Vector)
supertypes(Vector{String})
supertypes(Array)
supertypes(UnionAll)
supertypes(Union{Type,UnionAll,DataType})
supertypes(Tuple{Int,Int})

@quasiabstract struct Rectangle1
    width::Real
    height::Real
end

@quasiabstract struct Rectangle2{T<:Real}
    width::T
    height::T
end

@quasiabstract struct Rectangle3{T} <: Number
    width::T
    height::T
end

@quasiabstract struct Rectangle4{T<:Real} <: Number
    width::T
    height::T
end

@quasiabstract struct Box1 <: Rectangle1
    depth::Real
end

@quasiabstract struct Box2{T<:Real} <: Rectangle2{T}
    depth::T
end

@quasiabstract struct Box3{T} <: Rectangle3{T}
    depth::T
end

@quasiabstract struct Box4{T<:Real} <: Rectangle4{T}
    depth::T
end

#_____________________________________________________________________
#                            Alice's code
#
using Statistics, ReusePatterns

abstract type AbstractPolygon end

@quasiabstract mutable struct Polygon <: AbstractPolygon
    x::Vector{Float64}
    y::Vector{Float64}
end

# Retrieve the number of vertices, and their X and Y coordinates
vertices(p::Polygon) = length(p.x)
coords_x(p::Polygon) = p.x
coords_y(p::Polygon) = p.y

# Move, scale and rotate a polygon
function move!(p::Polygon, dx::Real, dy::Real)
    p.x .+= dx
    p.y .+= dy
end

function scale!(p::Polygon, scale::Real)
    m = mean(p.x); p.x = (p.x .- m) .* scale .+ m
    m = mean(p.y); p.y = (p.y .- m) .* scale .+ m
end

function rotate!(p::Polygon, angle_deg::Real)
    θ = float(angle_deg) * pi / 180
    R = [cos(θ) -sin(θ); sin(θ) cos(θ)]
    x = p.x .- mean(p.x)
    y = p.y .- mean(p.y)
    (x, y) = R * [x, y]
    p.x = x .+ mean(p.x)
    p.y = y .+ mean(p.y)
end

#_____________________________________________________________________
#                             Bob's code
#

@quasiabstract mutable struct RegularPolygon <: Polygon

    radius::Float64
end

function RegularPolygon(n::Integer, radius::Real)
    @assert n >= 3
    θ = range(0, stop=2pi-(2pi/n), length=n)
    c = radius .* exp.(im .* θ)
    return RegularPolygon(real(c), imag(c), radius)
end

# Compute length of a side and the polygon area
side(p::RegularPolygon) = 2 * p.radius * sin(pi / vertices(p))
area(p::RegularPolygon) = side(p)^2 * vertices(p) / 4 / tan(pi / vertices(p))

function scale!(p::RegularPolygon, scale::Real)
    invoke(scale!, Tuple{supertype(RegularPolygon), typeof(scale)}, p, scale) # call "super" method
    p.radius *= scale        # update internal state
end

# Attach a label to a polygon
mutable struct Named{T} <: AbstractPolygon
    polygon::T
    name::String
end
Named{T}(name, args...; kw...) where T = Named{T}(T(args...; kw...), name)
name(p::Named) = p.name

# Forward methods from `Named` to `Polygon`
@forward (Named, :polygon) RegularPolygon

#_____________________________________________________________________
#                           Charlie's code
#

# Here I use `Gnuplot.jl`, but any other package would work...
#=
using Gnuplot
Gnuplot.setverb(false)
@gp "set size ratio -1" "set grid" "set key bottom right" xr=(-1.5, 2.5) :-
=#

# Methods to plot a Polygon
plotlabel(p::Polygon) = "Polygon (" * string(length(p.x)) * " vert.)"
plotlabel(p::RegularPolygon) = "RegularPolygon (" * string(vertices(p)) * " vert., area=" * string(round(area(p) * 100) / 100) * ")"
plotlabel(p::Named) = name(p) * " (" * string(vertices(p)) * " vert., area=" * string(round(area(p) * 100) / 100) * ")"
function plot(p::AbstractPolygon; dt=1, color="black")
    x = coords_x(p); x = [x; x[1]]
    y = coords_y(p); y = [y; y[1]]
    title = plotlabel(p)
    # @gp :- x y "w l tit '$title' dt $dt lw 2 lc rgb '$color'"
end

# Finally, let's have fun with the shapes!
line = Polygon([0., 1.], [0., 1.])
triangle = RegularPolygon(3, 1)
square = Named{RegularPolygon}("Square", 4, 1)

plot(line, color="black")
plot(triangle, color="blue")
plot(square, color="red")

rotate!(triangle, 90)
move!(triangle, 1, 1)
scale!(triangle, 0.32)
scale!(square, sqrt(2))

plot(triangle, color="blue", dt=2)
plot(square, color="red", dt=2)

# Increase number of vertices
p = Named{RegularPolygon}("Heptagon", 7, 1)
plot(p, color="orange", dt=4)
circle = Named{RegularPolygon}("Circle", 1000, 1)
plot(circle, color="dark-green", dt=4)

ext: jl


text:
using TableTraits
using NamedTuples
import IterableTables

@testset "table tests" begin
  let x = K_Table(a=[1, 2])
    @test x == K_Table(Any[[1,2]], [:a])
    @test xt(kpointer(x)) == XT
    @test ncol(x) == 1
    @test nrow(x) == 2
    @test Q.colnames(x) == [:a]
    @test DataFrames.index(x) == DataFrames.Index(Dict(:a=>1), Symbol[:a])
    @test x[1] == [1, 2]
    @test x[2,1] == 2
    @test (x′ = K(r1(kpointer(x))); x == x′)
    @test begin
      t = K_Table(a=[0, 1, 2], b=[0., 10., 20.])
      show_to_string(t) == """
      3×2 Q.K_Table
      │ Row │ a │ b    │
      ├─────┼───┼──────┤
      │ 1   │ 0 │ 0.0  │
      │ 2   │ 1 │ 10.0 │
      │ 3   │ 2 │ 20.0 │"""
    end
    @test (x = DataFrame(a=1:3, b=[10., 20, 30]); x == K(x));
    @test (x = K_Table(a=Int[]); names(names!(x, [:b])) == [:b])
    @test (x = K_Table(a=Int[]); x[:a] == Int[])
  end
  @test (x = K_Table(@NT(a::Int, b::Float64), 3); size(x) == (3, 2))
  @test begin
    t = K_Table(a=[1])
    Q.coldata(t)[1] = push!(t[1], 2)
    Q.coldata(t)[1] = push!(t[1], 3)
    t[1] == [1, 2, 3]
  end
end

@testset "table iterator tests" begin
  let t = K_Table(a=[0, 1, 2], b=[0., 10., 20.]), ti = getiterator(t)
    @test isiterable(t)
    @test isiterabletable(t)
    @test length(ti) == 3
    @test eltype(ti) == typeof(@NT(a=0, b=0.))
    @test start(ti) == 1
    @test collect(ti) == [
      @NT(a=0, b=0.),
      @NT(a=1, b=10.),
      @NT(a=2, b=20.),
    ]
  end
  @test K_Table([@NT(a=1), @NT(a=2)]) == DataFrame(a=[1,2])
  @test K_Table(@NT(a=i) for i in 1:3 if i != 2) == DataFrame(a=[1,3])
end


text:
@doc raw"""
    GibbsKernel(; lengthscale)

Gibbs Kernel with lengthscale function `lengthscale`.

The Gibbs kernel is a non-stationary generalisation of the squared exponential
kernel. The lengthscale parameter ``l`` becomes a function of
position ``l(x)``.

# Definition

For inputs ``x, x'``, the Gibbs kernel with lengthscale function ``l(\cdot)``
is defined as
```math
k(x, x'; l) = \sqrt{\left(\frac{2 l(x) l(x')}{l(x)^2 + l(x')^2}\right)}
\quad \exp{\left(-\frac{(x - x')^2}{l(x)^2 + l(x')^2}\right)}.
```

For a constant function ``l \equiv c``, one recovers the [`SqExponentialKernel`](@ref)
with lengthscale `c`.

# References

Mark N. Gibbs. "Bayesian Gaussian Processes for Regression and Classication." PhD thesis, 1997

Christopher J. Paciorek and Mark J. Schervish. "Nonstationary Covariance Functions
for Gaussian Process Regression". NeurIPS, 2003

Sami Remes, Markus Heinonen, Samuel Kaski. "Non-Stationary Spectral Kernels". arXiV:1705.08736, 2017

Sami Remes, Markus Heinonen, Samuel Kaski. "Neural Non-Stationary Spectral Kernel". arXiv:1811.10978, 2018
"""
struct GibbsKernel{T} <: Kernel
    lengthscale::T
end

GibbsKernel(; lengthscale) = GibbsKernel(lengthscale)

function (k::GibbsKernel)(x, y)
    lengthscale = k.lengthscale
    lx = lengthscale(x)
    ly = lengthscale(y)
    l = invsqrt2 * hypot(lx, ly)
    kernel = (sqrt(lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)
    return kernel(x, y)
end


text:
## quadrature.jl : nystrom integration methods

"""
Abstract type `QuadratureRule`

The following quadrature rules are implemented:

- `Midpoint::QuadratureRule`: the midpoint rule
- `Trapezoidal::QuadratureRule`: the trapezoidal rule
- `Simpson::QuadratureRule`: Simpson's rule
- `GaussLegendre::QuadratureRule`: Gauss-Legendre quadrature rule
- `EOLE::QuadratureRule`: expansion-optimal linear estimation

See also: [`Midpoint`](@ref), [`Trapezoidal`](@ref), [`Simpson`](@ref), [`GaussLegendre`](@ref), [`EOLE`](@ref)
"""
abstract type QuadratureRule end

"""
    GaussLegendre()

Gauss-Legendre quadrature method.

See also: [`Midpoint`](@ref), [`Trapezoidal`](@ref), [`Simpson`](@ref), [`EOLE`](@ref)
"""
struct GaussLegendre <: QuadratureRule end

"""
    EOLE()

Expansion-optimal linear estimation.

See also: [`Midpoint`](@ref), [`Trapezoidal`](@ref), [`Simpson`](@ref), [`GaussLegendre`](@ref)
"""
struct EOLE <: QuadratureRule end

"""
    Simpson()

Simpson's method.

See also: [`Midpoint`](@ref), [`Trapezoidal`](@ref), [`GaussLegendre`](@ref), [`EOLE`](@ref)
"""
struct Simpson <: QuadratureRule end

"""
    Midpoint()

The midpoint rule.

See also: [`Trapezoidal`](@ref), [`Simpson`](@ref), [`GaussLegendre`](@ref), [`EOLE`](@ref)
"""
struct Midpoint <: QuadratureRule end

"""
    Trapezoidal()

The trapezoidal rule.

See also: [`Midpoint`](@ref), [`Simpson`](@ref), [`GaussLegendre`](@ref), [`EOLE`](@ref)
"""
struct Trapezoidal <: QuadratureRule end

# get nodes and weights of Gauss-Legendre quadrature on [a,b]
function get_nodes_and_weights(n::Integer, a, b, q::GaussLegendre)
    nodes, weights = gausslegendre(n)
    weights = (b-a)/2*weights
    nodes = (b-a)/2*nodes .+ (a+b)/2
    return nodes, weights
end

# get nodes and weights of structured grid on [a,b]
function get_nodes_and_weights(n::Integer, a, b, q::EOLE)
    nodes = n == 1 ? [(a+b)/2] : range(a; stop = b,length = n)
    weights = fill((b-a) / n, n)
    return nodes, weights
end

# get nodes and weights of Simpson's rule on [a,b]
function get_nodes_and_weights(n::Integer, a, b, q::Simpson)
    iseven(n) || begin
        @warn "to use Simpson's rule, n must be even (received $(n)). I will continue with n = $(n+1)"
        n += 1
    end
    Δx = (b-a)/n
    nodes = a:Δx:b
    weights = repeat(2:2:4, outer=Int(n/2))
    weights[1] = 1
    push!(weights,1)
    weights *= Δx/3
    return nodes, weights
end

# get nodes and weights of Midpoint rule on [a,b]
function get_nodes_and_weights(n::Integer, a, b, q::Midpoint)
    Δx = (b-a)/n
    nodes = (a+Δx/2):Δx:(b-Δx/2)
    weights = fill(Δx, size(nodes))
    return nodes, weights
end

# get nodes and weights of Trapezoidal rule on [a,b]
function get_nodes_and_weights(n::Integer, a, b, q::Trapezoidal)
    Δx = (b-a)/n
    nodes = a:Δx:b
    weights = fill(Δx, length(nodes))
    weights[1] /= 2
    weights[end] /= 2
    return nodes, weights
end

## EigenSolvers
"""
Abstract type `AbstractEigenSolver`

The following eigensolvers are implemented:
-`EigenSolver<:AbstractEigenSolver`: eigenvalue decomposition using **eigen**
-`EigsSolver<:AbstractEigenSolver`: eigenvalue decomposition using **eigs**

See also: [`EigenSolver`](@ref), [`EigsSolver`](@ref)
"""
abstract type AbstractEigenSolver end

"""
    EigenSolver()

Eigenvalue decomposition using **eigen**.

See also: [`EigsSolver`](@ref)
"""
struct EigenSolver <: AbstractEigenSolver end

compute(A, n, ::EigenSolver) = eigen(A, sortby=λ->-abs(real(λ)))

"""
    EigsSolver()

Eigenvalue decomposition using **eigs**.

See also: [`EigenSolver`](@ref)
"""
struct EigsSolver <: AbstractEigenSolver end

compute(A, n, ::EigsSolver) = eigs(A, nev=n, ritzvec=true, which=:LM, v0=randn(size(A,1)))


text:
using LinearAlgebra, ProgressMeter, Langevin

cd(@__DIR__)
D = read_dict("linesource.txt")
set = set_setup(D)

ps = set_geometry(D)
#ps = UnstructPSpace("linesource.su2")

# quadrature
quadratureorder = 8
points, weights = legendre_quadrature(quadratureorder)
vs = UnstructVSpace(-1.0, 1.0, length(weights), points, weights)
# material
radiation = Radiation(D[:knudsen], D[:sigmaS], D[:sigmaA])
# initial/boundary conditions
nq = size(points, 1)
f0 = ones(nq) .* 1e-4
w0 = discrete_moments(f0, weights)
ib = IB1F([w0], [w0], f0, [w0], [w0], [w0], f0, [w0])

ks = SolverSet(set, ps, vs, radiation, ib, pwd())

uq = UQ1D(3, 6, -0.1, 0.1)

function init_field(x, y)
    s2 = 0.03^2
    flr = 1e-4
    return max(flr, 1.0 / (4.0 * pi * s2) * exp(-(x^2 + y^2) / 4.0 / s2))
end

ctr = Array{KitBase.ControlVolumeUS1F}(undef, size(ps.cellid, 1))
for i in eachindex(ctr)
    n = Vector{Float64}[]
    for j = 1:3
        push!(
            n,
            KitBase.unit_normal(
                ps.points[ps.facePoints[ps.cellFaces[i, j], 1], :],
                ps.points[ps.facePoints[ps.cellFaces[i, j], 2], :],
            ),
        )

        if dot(ps.faceCenter[ps.cellFaces[i, j], :] .- ps.cellCenter[i, :], n[j]) < 0
            n[j] .= -n[j]
        end
    end

    phi = zeros(nq, uq.nq)
    for j = 1:uq.nq
        phi[:, j] .=
            init_field(ps.cellCenter[i, 1], ps.cellCenter[i, 2]) * (1.0 + uq.pceSample[j])
    end

    w = zeros(uq.nq)
    for i = 1:uq.nq
        w[i] = sum(weights .* phi[:, i])
    end
    dx = [
        KitBase.point_distance(
            ps.cellCenter[i, :],
            ps.points[ps.cellid[i, 1], :],
            ps.points[ps.cellid[i, 2], :],
        ),
        KitBase.point_distance(
            ps.cellCenter[i, :],
            ps.points[ps.cellid[i, 2], :],
            ps.points[ps.cellid[i, 3], :],
        ),
        KitBase.point_distance(
            ps.cellCenter[i, :],
            ps.points[ps.cellid[i, 3], :],
            ps.points[ps.cellid[i, 1], :],
        ),
    ]

    ctr[i] = KitBase.ControlVolumeUS1F(n, ps.cellCenter[i, :], dx, w, w, phi)
end

face = Array{KitBase.Interface2D1F}(undef, size(ps.facePoints, 1))
for i in eachindex(face)
    len = norm(ps.points[ps.facePoints[i, 1], :] .- ps.points[ps.facePoints[i, 2], :])
    n = KitBase.unit_normal(
        ps.points[ps.facePoints[i, 1], :],
        ps.points[ps.facePoints[i, 2], :],
    )

    if !(-1 in ps.faceCells[i, :])
        n0 = ps.cellCenter[ps.faceCells[i, 2], :] .- ps.cellCenter[ps.faceCells[i, 1], :]
    else
        n0 = zero(n)
    end
    if dot(n, n0) < 0
        n .= -n
    end

    fw = zeros(uq.nq)
    ff = zeros(nq, uq.nq)

    face[i] = KitBase.Interface2D1F(len, n[1], n[2], fw, ff)
end

dt = 1.2 / 150 * ks.set.cfl
nt = ks.set.maxTime ÷ dt |> Int
@showprogress for iter = 1:nt
    @inbounds Threads.@threads for i in eachindex(face)
        velo = vs.u[:, 1] .* face[i].n[1] + vs.u[:, 2] .* face[i].n[2]
        if !(-1 in ps.faceCells[i, :])
            for j = 1:uq.nq
                ff = @view face[i].ff[:, j]

                KitBase.flux_kfvs!(
                    ff,
                    ctr[ps.faceCells[i, 1]].f[:, j],
                    ctr[ps.faceCells[i, 2]].f[:, j],
                    velo,
                    dt,
                )
            end
        end
    end

    @inbounds Threads.@threads for i in eachindex(ctr)
        if ps.cellType[i] == 0
            for j = 1:3
                dirc = sign(dot(ctr[i].n[j], face[ps.cellFaces[i, j]].n))
                @. ctr[i].f -=
                    dirc * face[ps.cellFaces[i, j]].ff * face[ps.cellFaces[i, j]].len /
                    ps.cellArea[i]
            end

            integral = zeros(uq.nq)
            for k = 1:uq.nq
                integral[k] = discrete_moments(ctr[i].f[:, k], vs.weights)
            end
            #integral = KitBase.discrete_moments(ctr[i].f, vs.weights)
            integral ./= 4.0 * π
            for k = 1:uq.nq
                @. ctr[i].f[:, k] += (integral[k] - ctr[i].f[:, k]) * dt
                ctr[i].w[k] = sum(ctr[i].f[:, k] .* vs.weights)
            end
        end
    end
end

sol = zeros(length(ctr), uq.nq)
for i in axes(sol, 1)
    sol[i, :] = ctr[i].w
end

pce = zeros(length(ctr), uq.nr + 1)
for i in axes(sol, 1)
    pce[i, :] = ran_chaos(sol[i, :], uq)
end

meanstd = zeros(length(ctr), 2)
for i in axes(sol, 1)
    meanstd[i, 1] = mean(pce[i, :], uq.op)
    meanstd[i, 2] = std(pce[i, :], uq.op)
end

write_vtk(ks.ps.points, ks.ps.cellid, meanstd)


text:
n = 10 # rand(1:10)
A = matrixdepot("minij", n)
@test issymmetric(A)
@test isposdef(A)
println("'minij' passed test...")


text:
using Luna

a = 13e-6
gas = :Ar
pres = 7.5

τfwhm = 30e-15
λ0 = 800e-9
energy = 1e-6

L = 15e-2

grid = Grid.EnvGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

coren, densityfun = Capillary.gradient(gas, L, pres, 0);
m = Capillary.MarcatiliMode(a, coren, loss=false);
aeff = let m = m
    aeff(z) = Modes.Aeff(m, z=z)
end

energyfun, energyfunω = Fields.energyfuncs(grid)

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_env(PhysData.γ3_gas(gas)),)

linop!, βfun! = LinearOps.make_linop(grid, m, λ0);

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)
statsfun = Stats.default(grid, Eω, m, linop!, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop!, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)


text:
#######################################################
### ------------------ DEFINITIONS -------------------
#######################################################

abstract type Dimension end
abstract type TwoD <: Dimension end
abstract type ThreeD <: Dimension end

struct Square <: TwoD end
struct Circle <: TwoD end
struct Cube <: ThreeD end
struct Sphere <: ThreeD end
struct Cylinder <: ThreeD end

struct Atom{T<:Dimension}
    shape::T
    r::AbstractMatrix
    N::Int64
    sizes::Any
end

get_dimension(atom::Atom{T}) where {T<:TwoD} = 2
get_dimension(atom::Atom{T}) where {T<:ThreeD} = 3

abstract type Pump end
abstract type PlaneWave <: Pump end
abstract type Gaussian <: Pump end

struct PlaneWave2D <: PlaneWave
    direction::Vector
end
struct PlaneWave3D <: PlaneWave
    direction::Vector
end
struct Gaussian2D <: Gaussian
    w₀::Float64
end
struct Gaussian3D <: Gaussian
    w₀::Float64
    function Gaussian3D(w₀)
        λ = 2π/k₀
        if w₀ < 2λ
            @warn "Waist is Too Small, it should be larger than 2λ (~ $(round(2λ, digits=3))).
            \n If the input value is correct, at least make sure that `w₀` is larger than `system size`.
            \n This advice garantee that Transmission curves will have values below 1." maxlog=50
        end
        return new(w₀)
    end
end

struct Laser{T<:Pump}
    pump::T
    s::Float64
    Δ::Float64
end
struct Lasers{T<:Pump}
    pump::Vector{T}
    s::Vector{Float64}
    Δ::Vector{Float64}
end

abstract type Physics end
abstract type Linear <: Physics end
abstract type NonLinear <: Physics end

struct Scalar <: Linear end
struct Vectorial <: Linear end
struct MeanField <: NonLinear end
struct BBGKY <: NonLinear end

struct LinearOptics{T<:Linear}
    physic::T
    atoms::Atom
    laser::Laser
    kernelFunction::Function
    spectrum::Dict
    data::Dict
end

struct NonLinearOptics{T<:NonLinear}
    physic::T
    atoms::Atom
    laser::Laser
    excitations::Dict
    data::Dict
end


text:
module FileFind
	include("find.jl")
end


text:
module FastInverseSqrt

export rsqrt

# See https://en.wikipedia.org/wiki/Fast_inverse_square_root#Overview_of_the_code & https://www.youtube.com/watch?v=p8u_k2LIZyo
function rsqrt(x::Float32)
    xₛ = x
    int32 = reinterpret(UInt32, xₛ)
    int32 = 0x5f3759df - int32 >> 1
    xₛ = reinterpret(Float32, int32)
    xₛ *= 1.5f0 - x * 0.5f0 * xₛ^2
    return xₛ
end
function rsqrt(x::Float64)
    xₛ = x
    int64 = reinterpret(UInt64, xₛ)
    int64 = 0x5fe6eb50c7b537a9 - int64 >> 1  # See https://stackoverflow.com/a/11644533/3260253
    xₛ = reinterpret(Float64, int64)
    xₛ *= 1.5 - x * 0.5 * xₛ^2
    return xₛ
end
rsqrt(x::Real) = rsqrt(float(x))

end


text:
export BFFPSV18

# ===============================================================
# Bogomolov, Forets, Frehse, Podelski, Schilling, Viry. HSCC 2018
# ===============================================================

# dummy functions for option :lazy_inputs_interval
lazy_inputs_interval_always = (k -> true)
lazy_inputs_interval_never = (k -> false)

function ispartition(partition::AbstractVector{<:AbstractVector{Int}})
    current = 1
    for block in partition
        for i in block
            if i != current
                return false
            end
            current += 1
        end
    end
    return true
end

function options_BFFPSV18()
    return OptionSpec[
        # general options
        OptionSpec(:discretization, "forward", domain=String,
            info="model for bloating/continuous time analysis"),
        OptionSpec(:algorithm, "explicit", domain=String, domain_check=(
            v  ->  v in ["explicit", "wrap"]), info="algorithm backend"),
        OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
            domain_check=(v  ->  v > 0.), info="time step"),
        OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=(
            v  ->  length(v) > 0 && all(e -> e > 0, v)),
            info="variables of interest; default: all variables"),
        OptionSpec(:partition, [Int[]],
            domain=AbstractVector{<:AbstractVector{Int}}, domain_check=
            ispartition,
            info="block partition; a block is represented by a vector " *
                 "containing its indices"),

        # discretization options
        OptionSpec(:sih_method, "concrete", domain=String,
            info="method to compute the symmetric interval hull in discretization"),
        OptionSpec(:exp_method, "base", domain=String,
            info="method to compute the matrix exponential"),
        OptionSpec(:assume_sparse, false, domain=Bool,
            info="use an analysis for sparse discretized matrices?"),

        # reachability options
        OptionSpec(:lazy_inputs_interval, lazy_inputs_interval_always,
            domain=Union{Int, Function},
            domain_check=(v  ->  !(v isa Int) || v >= -1),
            info="length of interval in which the inputs are handled as a " *
                 "lazy set (``-1`` for 'never'); may generally also be a " *
                 "predicate over indices; the default corresponds to ``-1``"),

        # approximation options
        OptionSpec(:block_options, nothing, domain=Any,
            info="short hand to set ':block_options_init' and " *
                 "':block_options_iter'"),
        OptionSpec(:block_options_init, nothing, domain=Any,
            info="option for the approximation of the initial states " *
                 "(during decomposition)"),
        OptionSpec(:block_options_iter, nothing, domain=Any,
            info="option for the approximation of the states ``X_k``, ``k>0``"),

        # convenience options
        OptionSpec(:assume_homogeneous, false, domain=Bool,
            info="ignore dynamic inputs during the analysis?"),
        OptionSpec(:eager_checking, true, domain=Bool,
            info="terminate as soon as property violation was detected?"),
    ]
end

function normalization_BFFPSV18!(𝑂::TwoLayerOptions)
    # :lazy_inputs_interval option: convert integers to functions
    if haskey_specified(𝑂, :lazy_inputs_interval)
        v = 𝑂[:lazy_inputs_interval]
        if v isa Int
            if v == -1
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_never
            elseif v == 0
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_always
            else
                𝑂.specified[:lazy_inputs_interval] = (k -> k % v == 0)
            end
        end
    end

    # :block_options options: use convenience option for '_init' and '_iter'
    if !haskey_specified(𝑂, :block_options_init) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_init] = 𝑂[:block_options]
    end
    if !haskey_specified(𝑂, :block_options_iter) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_iter] = 𝑂[:block_options]
    end

    nothing
end

function validation_BFFPSV18(𝑂)
    # block_options
    for b_options in [:block_options, :block_options_init, :block_options_iter]
        if haskey_specified(𝑂, b_options)
            bo = 𝑂[b_options]
            if bo isa Type{<:LazySet} || bo isa Type{<:AbstractDirections}
                # uniform options
            elseif bo isa Symbol
                # template directions
                option = get(Utils.template_direction_symbols, bo, nothing)
                if option == nothing
                    throw(DomainError(bo, "if the `$b_options` option " *
                        "is a Symbol, it must be one of " *
                        "$(keys(Utils.template_direction_symbols))"))
                end
                𝑂.specified[b_options] = option
            elseif bo isa AbstractVector || bo isa Dict{Int, Any}
                # mapping
            elseif bo isa Real || bo isa Pair{<:UnionAll, <:Real}
                ε = bo isa Real ? bo : bo[2]
                if ε <= 0
                    throw(DomainError(ε, "the `$b_options` option must be " *
                                         "positive"))
                end
            elseif b_options == :block_options_iter && bo == nothing
                # no overapproximation
            else
                throw(DomainError(bo == nothing ? "nothing" : bo,
                    "the `$b_options` option does not accept the given input"))
            end
        end
    end

    nothing
end

"""
    BFFPSV18 <: ContinuousPost

Implementation of the reachability algorithm for purely continuous linear
time-invariant systems using block decompositons by S. Bogomolov, M. Forets,
G. Frehse, A. Podelski, C. Schilling and F. Viry [1].

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Notes

The following options are available:

```julia
$(print_option_spec(options_BFFPSV18()))
```

### Algorithm

We refer to [1] for technical details.

[1] [Reach Set Approximation through Decomposition with Low-dimensional Sets
and High-dimensional Matrices](https://dl.acm.org/citation.cfm?id=3178128).
S. Bogomolov, M. Forets, G. Frehse, A. Podelski, C. Schilling, F. Viry.
HSCC '18 Proceedings of the 21st International Conference on Hybrid Systems:
Computation and Control (part of CPS Week).
"""
struct BFFPSV18 <: ContinuousPost
    options::TwoLayerOptions

    function BFFPSV18(𝑂::Options)
        normalized_𝑂 = validate_and_wrap_options(𝑂, options_BFFPSV18();
            validation=validation_BFFPSV18,
            normalization=normalization_BFFPSV18!)
        return new(normalized_𝑂)
    end
end

# convenience constructor from pairs of symbols
BFFPSV18(𝑂::Pair{Symbol,<:Any}...) = BFFPSV18(Options(Dict{Symbol,Any}(𝑂)))

# default options
BFFPSV18() = BFFPSV18(Options())

init(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂))

function init!(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)
    # state dimension for (purely continuous or purely discrete systems)
    𝑂copy = copy(𝑂)
    𝑂copy[:n] = statedim(𝑆)

    # solver-specific options (adds default values for unspecified options)
    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)

    # :vars option; default: all variables
    if haskey_specified(𝒫.options, :vars)
        𝑂validated[:vars] = 𝒫.options[:vars]
    else
        𝑂validated[:vars] = 1:𝑂validated[:n]
    end

    # :partition option: use 1D blocks
    if haskey_specified(𝒫.options, :partition)
        𝑂validated[:partition] = 𝒫.options[:partition]
    else
        𝑂validated[:partition] = [[i] for i in 1:𝑂validated[:n]]
    end

    # :blocks option (internal only)
    # list of all interesting block indices in the partition
    𝑂validated[:blocks] = compute_blocks(𝑂validated[:vars], 𝑂validated[:partition])

    # :block_options_init & :block_options_iter options:
    # set default according to :partition
    if !haskey_specified(𝒫.options, :block_options_init)
        𝑂validated[:block_options_init] =
            compute_default_block_options(𝑂validated[:partition])
    end
    if !haskey_specified(𝒫.options, :block_options_iter)
        𝑂validated[:block_options_iter] =
            compute_default_block_options(𝑂validated[:partition])
    end

    # Input -> Output variable mapping
    𝑂validated[:inout_map] = inout_map_reach(𝑂validated[:partition], 𝑂validated[:blocks], 𝑂validated[:n])

    if 𝑂validated[:project_reachset]
        𝑂validated[:output_function] = nothing
    else
        𝑂validated[:output_function] = 𝑂validated[:projection_matrix]
    end

    return 𝑂validated
end

"""
    post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂::Options)

Calculate the reachable states of the given initial value problem using `BFFPSV18`.

### Input

- `𝒫` -- post operator of type `BFFPSV18`
- `𝑆` -- sytem, initial value problem for a continuous ODE
- `𝑂` -- algorithm-specific options
"""
function post(𝒫::BFFPSV18, 𝑆::AbstractSystem, 𝑂_input::Options)
    𝑂 = TwoLayerOptions(merge(𝑂_input, 𝒫.options.specified), 𝒫.options.defaults)

    if 𝑂[:mode] == "reach"
        info("Reachable States Computation...")
        @timing begin
            Rsets = reach(𝑆, 𝑂)
            info("- Total")
        end

        # Projection
        if 𝑂[:project_reachset]
            info("Projection...")
            RsetsProj = @timing project(Rsets, 𝑂)
        else
            RsetsProj = Rsets
        end

        return ReachSolution(RsetsProj, 𝑂_input)

    elseif 𝑂[:mode] == "check"
        info("invariants are currently not supported in 'check' mode")

        # Input -> Output variable mapping in property
        property = inout_map_property(𝑂[:property], 𝑂[:partition], 𝑂[:blocks], 𝑂[:n])

        # =================
        # Property checking
        # =================
        info("Property Checking...")
        @timing begin
            answer = check_property(𝑆, property, 𝑂)
            info("- Total")
        end

        if answer == 0
            info("The property is satisfied!")
            return CheckSolution(true, -1, 𝑂_input)
        else
            info("The property may be violated at index $answer," *
                " (time point $(answer * 𝑂[:δ]))!")
            return CheckSolution(false, answer, 𝑂_input)
        end
    else
        error("unsupported mode $(𝑂[:mode])")
    end # mode
end

function compute_blocks(vars, partition)
    blocks = Vector{Int}()
    sizehint!(blocks, length(vars))
    next = 0
    var_idx = 1
    for (i, block) in enumerate(partition)
        next += length(block)
        if vars[var_idx] <= next
            push!(blocks, i)
            var_idx += 1
            while var_idx <= length(vars) && vars[var_idx] <= next
                var_idx += 1
            end
            if var_idx > length(vars)
                break
            end
        end
    end
    @assert var_idx == length(vars) + 1
    sizehint!(blocks, length(blocks))
    return blocks
end

function compute_default_block_options(partition)
    # use Interval for 1D blocks and Hyperrectangle otherwise
    block_options = Vector{Type{<:LazySet}}(undef, length(partition))
    for (i, block) in enumerate(partition)
        block_options[i] = length(block) == 1 ? Interval : Hyperrectangle
    end
    return block_options
end


text:
module NumBorders
    export get_num_borders!

    include("../10/horizonside.jl")
    include("../10/functional_robot.jl")
    include("interface_fold_borders_counter.jl")

    function get_num_borders!(robot)
        robot = interface_save_robot(robot)

        horizont_borders_counter = interface_fold_borders_counter(robot, Nord)
        rectangular_borders = interface_rectangular_borders(robot)
        trajectories = interface_trajectories(robot)

        trajectories.snake!(Ost, Nord) do side
            if robot.isborder(Nord)
                return false
            end
            while horizont_borders_counter.move!(side) || rectangular_borders.move!(side)
            end
        end

       return horizont_borders_counter.get_num()
    end
end # module NumBorders

#-------Исполняемая часть файла:
using .NumBorders
using HorizonSideRobots

robot = Robot("11/horizontal_borders.sit", animate=true)

get_num_borders!(robot) |> println


text:
import JSON

timed_function_call(f, input) = @timed f(input...)

function json_array_dim(a)
    if length(size(a)) > 0
        return 1 + json_array_dim(a[1])
    else
        return 0
    end
end

function json_array_eltype(a)
    if eltype(a) == Any
        return json_array_eltype(a[1])
    else
        return eltype(a)
    end
end

juliafy_json(t) = t
juliafy_json(a::Array) = convert(Array{{json_array_eltype(a), json_array_dim(a)}}, hcat(a...))

tupleit(t) = tuple(t)
tupleit(t::Tuple) = t

input_tuples = JSON.Parser.parsefile("{:s}")

for (i, input_tuple) in enumerate(input_tuples)

    input_tuple = [juliafy_json(elem) for elem in input_tuple]

    output_tuple = $FUNCTION_NAME(input_tuple...) |> tupleit

    open("{:s}.output$i.json", "w") do f
       JSON.print(f, output_tuple)
    end
end

max_stars_repo_name: BinaryFissionGames/lovelace-engine

max_issues_repo_name: BinaryFissionGames/lovelace-engine

max_forks_repo_name: BinaryFissionGames/lovelace-engine


text:
using StrainRecon
using Base.Test

@testset "findM" begin
    n = 3
    m = 5
    nsamp = 5
    w = randomw(n,nsamp)
    Mtrue = randomM(m,n)
    R = 1/(1e-8)*StrainRecon.sdiag(ones(m))
    for i = 1:nsamp
        d = Mtrue*w[:,i]
        param = getStrainReconParam(d,n,R=R)
        M_nsol,pm = findM(param, w[:,i], nsol = 2^n, probM = true)
        for j = 1:2^n
            M = M_nsol[:,:,j]*1
            #@test M*w[:,i]
        end

        @test all(diff(pm,2) .<= 0) # check that subsequent rows of M in nth + 1 matrix are not better than previous matrix

        @test norm(sum(pm,2) .- 1) < 1e-8 # sum of p(M) across matrices at the same row = 1

    end
end


text:
using PowerModels
using .PP2PM

import Ipopt
import Juniper
import JuMP

function run_powermodels(json_path)
    # function to run optimal transmission switching (OTS) optimization from powermodels.jl
    pm = PP2PM.load_pm_from_json(json_path)

    juniper_solver = JuMP.with_optimizer(Juniper.Optimizer, nl_solver =
                                JuMP.with_optimizer(Ipopt.Optimizer, tol = 1e-4, print_level = 0))

    result = run_ots(pm, ACPPowerModel, juniper_solver,
                        setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

max_stars_repo_name: Zamwell/pandapower

max_issues_repo_name: Zamwell/pandapower

max_forks_repo_name: Zamwell/pandapower


text:
using FeedbackParticleFilters, Distributions, LinearAlgebra, Random, PDMats, FillArrays

println("Testing BrownianMotionTorus.jl:")

@testset "BrownianMotionTorus.jl" begin
    @testset "$(d)d" for d in 1:3
        print("  inner constructor")
        model = BrownianMotionTorus(d)

        print("  method state_dim")
        print("-")
        @test state_dim(model) == d
        println("DONE")

        print("  method noise_dim")
        print("-")
        @test noise_dim(model) == d
        println("DONE")

        print("  method initialize")
        x0 = initialize(model)
        print("-")
        @test all(0 .<= x0 .<= 2π)
        println("DONE")

        print("  calling instance of BrownianMotionTorus")
        x0 = ones(d)
        x1 = model(x0, 0.01)
        print("-")
        @test x1 != x0
        print("-")
        @test all(0 .<= x0 .<= 2π)
        xx0 = hcat(ones(d), zeros(d))
        xx1 = model(xx0, 0.01)
        print("-")
        @test xx1 != xx0
        print("-")
        @test all(0 .<= x0 .<= 2π)
        println("DONE")

        print("  method propagate!")
        x0 = ones(d)
        x1 = copy(x0)
        propagate!(x1, model, 0.01)
        print("-")
        @test x1 != x0
        xx0 = hcat(ones(d), zeros(d))
        xx1 = copy(xx0)
        propagate!(xx1, model, 0.01)
        print("-")
        @test xx1 != xx0
        println("DONE")
    end
end; #BrownianMotionTorus.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using LazyBandedMatrices, FillArrays, BandedMatrices, BlockBandedMatrices, BlockArrays, ArrayLayouts, LazyArrays, Test
import BlockBandedMatrices: isbandedblockbanded, BandedBlockBandedStyle, BandedLayout
import LazyBandedMatrices: KronTravBandedBlockBandedLayout, BroadcastBandedLayout, arguments, FillLayout, OnesLayout, call
import BandedMatrices: BandedColumns

@testset "Banded kron" begin
    A = brand(5,5,2,2)
    B = brand(2,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (5,4)
    @test Matrix(K) == kron(Matrix(A), Matrix(B))

    A = brand(3,4,1,1)
    B = brand(3,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (7,2)
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B))
    K = kron(B,A)
    @test Matrix(K) ≈ kron(Matrix(B), Matrix(A))

    K = kron(A, B')
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B'))
    K = kron(A', B)
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A'), Matrix(B))
    K = kron(A', B')
    K isa BandedMatrix
    @test Matrix(K) ≈ kron(Matrix(A'), Matrix(B'))

    A = brand(5,6,2,2)
    B = brand(3,2,1,0)
    K = kron(A,B)
    @test K isa BandedMatrix
    @test bandwidths(K) == (12,4)
    @test Matrix(K) ≈ kron(Matrix(A), Matrix(B))

    n = 10; h = 1/n
    D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))
    D_xx = kron(D², Eye(n))
    @test D_xx isa BandedMatrix
    @test bandwidths(D_xx) == (10,10)
    D_yy = kron(Eye(n), D²)
    @test D_yy isa BandedMatrix
    @test bandwidths(D_yy) == (1,1)
    Δ = D_xx + D_yy
    @test Δ isa BandedMatrix
    @test bandwidths(Δ) == (10,10)

    @testset "#87" begin
        @test kron(Diagonal([1,2,3]), Eye(3)) isa Diagonal{Float64,Vector{Float64}}
    end
end

@testset "DiagTrav" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test DiagTrav(A) == [1, 4, 2, 7, 5, 3]
    A = [1 2 3; 4 5 6]
    @test DiagTrav(A) == [1, 4, 2, 5, 3]
    A = [1 2; 3 4; 5 6]
    @test DiagTrav(A) == [1, 3, 2, 5, 4]

    @test resize!(DiagTrav(A), Block(2)) == [1, 3,2]

    A = DiagTrav(randn(3,3,3))
    @test A[Block(1)] == A[1:1,1,1]
    @test A[Block(2)] == [A.array[2,1,1], A.array[1,2,1], A.array[1,1,2]]
    @test A[Block(3)] == [A.array[3,1,1], A.array[2,2,1], A.array[1,3,1],
                          A.array[2,1,2], A.array[1,2,2], A.array[1,1,3]]
    @test A == [A[Block(1)]; A[Block(2)]; A[Block(3)]]
end

@testset "BlockKron" begin
    n = 4
    Δ = BandedMatrix(1 => Ones(n-1), 0 => Fill(-2,n), -1 => Ones(n-1))
    A = BlockKron(Δ, Eye(n))
    @test isblockbanded(A)
    @test isbandedblockbanded(A)
    @test blockbandwidths(A) == (1,1)
    @test BandedBlockBandedMatrix(A) == A
end

@testset "KronTrav" begin
    @testset "vector" begin
        a = [1,2,3]
        b = [4,5,6]
        c = [7,8]
        @test KronTrav(a,b) == DiagTrav(b*a')
        @test KronTrav(a,c) == DiagTrav(c*a')
        @test KronTrav(c,a) == DiagTrav(a*c')
    end

    @testset "matrix" begin
        A = [1 2; 3 4]
        B = [5 6; 7 8]
        K = KronTrav(A,B)

        X = [9 10; 11 0]
        @test K*DiagTrav(X) == DiagTrav(B*X*A')

    end

    @testset "banded" begin
        n = 4
        Δ = BandedMatrix(1 => Ones(n-1), 0 => Fill(-2,n), -1 => Ones(n-1))
        A = KronTrav(Δ, Eye(n))
        B = KronTrav(Eye(n), Δ)

        X = triu!(randn(n,n))[:,end:-1:1]
        @test A * DiagTrav(X) == DiagTrav(X * Δ')
        @test B * DiagTrav(X) == DiagTrav(Δ * X)

        @test blockbandwidths(A) == blockbandwidths(B) == (1,1)
        @test subblockbandwidths(A) == (1,1)
        @test subblockbandwidths(B) == (0,0)
        @test isblockbanded(A)
        @test isbandedblockbanded(A)
        @test BandedBlockBandedMatrix(A) == A
        @test MemoryLayout(A) isa KronTravBandedBlockBandedLayout

        @test Base.BroadcastStyle(typeof(A)) isa BandedBlockBandedStyle
        @test A+B isa BandedBlockBandedMatrix
        @test A+B == Matrix(A) + Matrix(B)
    end

    @testset "Views" begin
        n = 4
        Δ = BandedMatrix(1 => Ones(n-1), 0 => Fill(-2,n), -1 => Ones(n-1))
        A = KronTrav(Δ, Eye(n))

        V = view(A, Block(2,2))
        @test MemoryLayout(V) isa BroadcastBandedLayout{typeof(*)}
        @test call(V) == *
        @test MemoryLayout.(arguments(V)) isa Tuple{BandedColumns{DenseColumnMajor},BandedColumns{OnesLayout}}
        @test BandedMatrix(V) == V == A[Block(2,2)]

        V = view(A, Block.(2:4), Block.(1:4))
        @test blockbandwidths(V) == (0,2)
        @test subblockbandwidths(V) == (1,1)
        @test BandedBlockBandedMatrix(V) == A[Block.(2:4), Block.(1:4)]
        @test A[Block.(2:4), Block.(1:4)] isa BandedBlockBandedMatrix
    end
end


text:
function setupForcing!(
    mb    :: ModelBlock;
    w_max :: Float64,
)
    ev = mb.ev
    fi = mb.fi
    co = mb.co
    cfg = ev.config
    wksp = co.wksp
    gd = ev.gd
    gd_slab = ev.gd_slab
    amo = co.amo
    amo_slab = co.amo_slab

    if cfg["transform_vector_field"]
        PolelikeCoordinate.project!(
            gd,
            fi.TAUX_east,
            fi.TAUY_north,
            fi.TAUX,
            fi.TAUY;
            direction=:Forward,
            grid=:T,
        )
    else
        fi.TAUX .= fi.TAUX_east
        fi.TAUY .= fi.TAUY_north
    end

    # Setup Ekman Flow
    if cfg["advection_scheme"] == "static"

        fi._u .= 0.0
        fi._v .= 0.0
        fi._w .= 0.0

    elseif cfg["advection_scheme"] in [ "ekman_AGA2020", "ekman_AGA2020_allowU", "ekman_KSC2018", "ekman_CO2012" ]

        Mx_east = nothing
        My_north = nothing

        if cfg["advection_scheme"] in [ "ekman_CO2012", "ekman_KSC2018" ]

            f_sT = co.mtx[:f_sT]
            ϵ_sT = co.mtx[:ϵ_sT]
            invD_sT = co.mtx[:invD_sT]

            if cfg["advection_scheme"] == "ekman_CO2012"
                switch = 1.0
            elseif cfg["advection_scheme"] == "ekman_KSC2018"
                switch = 0.0
            else
                throw(ErrorException("Unexpected scenario. Please check."))
            end

            Mx_east  = (   ϵ_sT .* fi.TAUX_east + f_sT .* fi.TAUY_north  ) .* invD_sT / ρ_sw
            My_north = ( - f_sT .* fi.TAUX_east + ϵ_sT .* fi.TAUY_north * switch  ) .* invD_sT / ρ_sw

        elseif cfg["advection_scheme"] in [ "ekman_AGA2020", "ekman_AGA2020_allowU" ]

            if ! haskey(co.mtx, :ϵ2invβ_sT)
                co.mtx[:ϵ2invβ_sT] = co.mtx[:ϵ_sT].^2 * (gd.R / 2.0 / gd.Ω)
            end

            f_sT = co.mtx[:f_sT]
            ϵ_sT = co.mtx[:ϵ_sT]
            invD_sT = co.mtx[:invD_sT]
            ϵ2invβ_sT = co.mtx[:ϵ2invβ_sT]
            γ = cfg["γ"]

            # First, I need to get the curl. I choose to
            # do the curl using line integral in ocean model
            # space
            curlτ_sT = reshape(
                amo_slab.T_CURLx_T * view(fi.TAUY, :) + amo_slab.T_CURLy_T * view(fi.TAUX, :),
                1, gd.Nx, gd.Ny,
            )

            if cfg["advection_scheme"] == "ekman_AGA2020_allowU"
                switch = 1.0
            elseif cfg["advection_scheme"] == "ekman_AGA2020"
                switch = 0.0
            else
                throw(ErrorException("Unexpected scenario. Please check."))
            end

            Mx_east  = (   f_sT .* fi.TAUY_north  ) .* invD_sT / ρ_sw
            My_north = ( - f_sT .* fi.TAUX_east + γ * ϵ2invβ_sT .* curlτ_sT ) .* co.mtx[:invD_sT] / ρ_sw

        else
            throw(ErrorException("Unexpected scenario. Please check."))
        end

        if Mx_east == nothing || My_north == nothing
            throw(ErrorException("Error: Either or both Mx_east or My_north are not assigned."))
        end

        Mx_sT = getSpace!(wksp, :sT; o=0.0)
        My_sT = getSpace!(wksp, :sT; o=0.0)

        PolelikeCoordinate.project!(
            gd,
            Mx_east,
            My_north,
            Mx_sT,
            My_sT;
            direction=:Forward,
            grid=:T,
        )

        # notice that we drop the z dimension for simplicity in the for loop
        Mx_u = reshape( co.amo_slab.U_interp_T * view(Mx_sT, :), gd.Nx, gd.Ny)
        My_v = reshape( co.amo_slab.V_interp_T * view(My_sT, :), gd.Nx, gd.Ny+1)

        H_Ek = sum(gd.Δz_T[1:cfg["Ekman_layers"], 1, 1])
        H_Rf = sum(gd.Δz_T[(cfg["Ekman_layers"]+1):(cfg["Ekman_layers"]+cfg["Returnflow_layers"]), 1, 1])

        for k = 1:cfg["Ekman_layers"]
            fi.sv[:UVEL][k, :, :] .= Mx_u ./ H_Ek
            fi.sv[:VVEL][k, :, :] .= My_v ./ H_Ek
        end

        for k = (cfg["Ekman_layers"]+1):(cfg["Ekman_layers"]+cfg["Returnflow_layers"])
            fi.sv[:UVEL][k, :, :] .= - Mx_u ./ H_Rf
            fi.sv[:VVEL][k, :, :] .= - My_v ./ H_Rf
        end

    else
        throw(ErrorException("Unknown scheme: " * string(cfg["advection_scheme"])))

    end

    _u = getSpace!(wksp, :U, true)
    _v = getSpace!(wksp, :V, true)
    _w = getSpace!(wksp, :W, true)

    _u[:] = fi._u
    _v[:] = fi._v
    _w[:] = fi._w

    mul!(fi._u, co.amo.U_flowmask_U, _u)
    mul!(fi._v, co.amo.V_flowmask_V, _v)
    mul!(fi._w, co.amo.W_flowmask_W, _w)

    tmp_T  = getSpace!(wksp, :T, true)
    tmp2_T = getSpace!(wksp, :T, true)

    mul!(tmp_T, co.amo.T_DIVx_U, fi._u)
    mul!(tmp_T, co.amo.T_DIVy_V, fi._v, 1.0, 1.0)

    # compute w
    DIVvol_T = reshape( tmp_T, amo.bmo.T_dim...)

    fi.sv[:WVEL][1, :, :] .= 0.0
    for k=1:gd.Nz
        fi.sv[:WVEL][k+1, :, :] .= fi.sv[:WVEL][k, :, :] + DIVvol_T[k, :, :] * gd.Δz_T[k, 1, 1]
    end

    CFL_break = 0
    violate_w_max = 0.0
    violate_w_min = 0.0
    for (i, w) in enumerate(fi._w)
        if w > w_max
            CFL_break += 1
            fi._w[i] = w_max
            violate_w_max = max(violate_w_max, w)
        elseif w < - w_max
            CFL_break += 1
            fi._w[i] = - w_max
            violate_w_min = min(violate_w_min, w)
        end
    end
    if CFL_break != 0
        println("CFL condition breaks in $(CFL_break) grid points. Arbitrarily cap w = ±$(w_max). Violated w (min, max) = ($(violate_w_min), $(violate_w_max))")
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Test
using SymbolicRegression

include("test_params.jl")

## Test Base.print
options = Options(;
    default_params..., binary_operators=(+, *, /, -), unary_operators=(cos, sin)
)

f = (x1, x2, x3) -> (sin(cos(sin(cos(x1) * x3) * 3.0) * -0.5) + 2.0) * 5.0

tree = f(Node("x1"), Node("x2"), Node("x3"))

s = repr(tree)
true_s = "((sin(cos(sin(cos(x1) * x3) * 3.0) * -0.5) + 2.0) * 5.0)"

@test s == true_s

# TODO: Next, we test that custom varMaps work:

EquationSearch(
    randn(Float32, 3, 10),
    randn(Float32, 10);
    options=options,
    varMap=["v1", "v2", "v3"],
    niterations=0,
    multithreading=true,
)

s = repr(tree)
true_s = "((sin(cos(sin(cos(v1) * v3) * 3.0) * -0.5) + 2.0) * 5.0)"
@test s == true_s


text:
@testset "unwind function" begin
    iterators = ([1, 2], ["a","b"], ["x", "y", "z"])
    @test MLJTuning.unwind(iterators...) ==
        [1  "a"  "x";
         2  "a"  "x";
         1  "b"  "x";
         2  "b"  "x";
         1  "a"  "y";
         2  "a"  "y";
         1  "b"  "y";
         2  "b"  "y";
         1  "a"  "z";
         2  "a"  "z";
         1  "b"  "z";
         2  "b"  "z"]
end

@test MLJTuning.delete((x=1, y=2, z=3), :x, :z) == (y=2,)

@testset "signature of measure" begin
    measures = [accuracy, confmat, misclassification_rate]
    @test MLJTuning.signature.(measures) == [-1, 0, 1]
end

true


text:
struct Match
    rule::AbstractRule
    # the rhs pattern to instantiate
    pat_to_inst::Union{Nothing,Pattern}
    # the substitution
    sub::Sub
    # the id the matched the lhs
    id::EClassId
end

const MatchesBuf = Vector{Match}

function cached_ids(g::EGraph, p::Pattern)# ::Vector{Int64}
    collect(keys(g.classes))
end

# FIXME
function cached_ids(g::EGraph, p::PatTerm)
    # println("pattern $p, $(p.head)")
    # println("all ids")
    # keys(g.classes) |> println
    # println("cached symbols")
    # cached = get(g.symcache, p.head, Set{Int64}())
    # println("symbols where $(p.head) appears")
    # appears = Set{Int64}()
    # for (id, class) ∈ g.classes
    #     for n ∈ class
    #         if n.head == p.head
    #             push!(appears, id)
    #         end
    #     end
    # end
    # # println(appears)
    # if !(cached == appears)
    #     @show cached
    #     @show appears
    # end

    collect(keys(g.classes))
    # cached
    # get(g.symcache, p.head, [])
end

# function cached_ids(g::EGraph, p::PatLiteral)
#     get(g.symcache, p.val, [])
# end

function (r::SymbolicRule)(g::EGraph, id::EClassId)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, r.right, sub, id)
end

function (r::DynamicRule)(g::EGraph, id::EClassId)
    ematch(g, r.ematch_program, id) .|> sub -> Match(r, nothing, sub, id)
end

function (r::BidirRule)(g::EGraph, id::EClassId)
    vcat(ematch(g, r.ematch_program_l, id) .|> sub -> Match(r, r.right, sub, id),
        ematch(g, r.ematch_program_r, id) .|> sub -> Match(r, r.left, sub, id))
end

macro maybethreaded(x, body)
    esc(quote
        if $x
            Threads.@threads $body
        else
            $body
        end
    end)
end

"""
Returns an iterator of `Match`es.
"""
function eqsat_search!(egraph::EGraph, theory::Vector{<:AbstractRule},
    scheduler::AbstractScheduler, report; threaded=false)
    match_groups = Vector{Match}[]
    function pmap(f, xs)
        # const propagation should be able to optimze one of the branch away
        if threaded
            # # try to divide the work evenly between threads without adding much overhead
            # chunks = Threads.nthreads() * 10
            # basesize = max(length(xs) ÷ chunks, 1)
            # ThreadsX.mapi(f, xs; basesize=basesize)
            ThreadsX.map(f, xs)
        else
            map(f, xs)
        end
    end

    inequalities = filter(Base.Fix2(isa, UnequalRule), theory)
    # never skip contradiction checks
    append_time = TimerOutput()
    for rule ∈ inequalities
        @timeit report.to repr(rule) begin
            ids = cached_ids(egraph, rule.left)
            rule_matches = pmap(i -> rule(egraph, i), ids)
            @timeit append_time "appending matches" begin
                append!(match_groups, rule_matches)
            end
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    for rule ∈ other_rules
        @timeit report.to repr(rule) begin
            # don't apply banned rules
            if !cansearch(scheduler, rule)
                # println("skipping banned rule $rule")
                continue
            end

            ids = cached_ids(egraph, rule.left)
            rule_matches = pmap(i -> rule(egraph, i), ids)

            n_matches = sum(length, rule_matches)
            can_yield = inform!(scheduler, rule, n_matches)
            if can_yield
                @timeit append_time "appending matches" begin
                    append!(match_groups, rule_matches)
                end
            end
        end
    end

    # @timeit append_time "appending matches" begin
    #     result = reduce(vcat, match_groups) # this should be more efficient than multiple appends
    # end
    merge!(report.to, append_time, tree_point=["Search"])

    return Iterators.flatten(match_groups)
    # return result
end


text:
# This file establishes the default initial conditions, boundary conditions and sources
# for the heldsuarez_problem experiment, following:
#
# - Held, I. M. and Suarez, M. J.: A proposal for the intercomparison
# of the dynamical cores of atmospheric general circulation models,
# B. Am. Meteorol. Soc., 75, 1825–1830, 1994.
#
# - Thatcher, D. R. and Jablonowski, C.: A moist aquaplanet variant of the
# Held–Suarez test for atmospheric model dynamical cores, Geosci. Model Dev.,
# 9, 1263–1292, 2016.

# Override default CLIMAParameters for consistency with literature on this case
using ClimateMachine.Thermodynamics

using CLIMAParameters.Planet

struct HeldSuarezProblem{BC, ISP, ISA, WP, BS, MP} <: AbstractAtmosProblem
    boundarycondition::BC
    init_state_prognostic::ISP
    init_state_auxiliary::ISA
    perturbation::WP
    base_state::BS
    moisture_profile::MP
end
function HeldSuarezProblem(;
    boundarycondition = (AtmosBC(), AtmosBC()),
    perturbation = nothing,
    base_state = nothing,
    moisture_profile = nothing,
)
    # Set up defaults
    if isnothing(perturbation)
        perturbation = DeterministicPerturbation()
    end
    if isnothing(base_state)
        base_state = HeldSuarezBaseState()
    end
    if isnothing(moisture_profile)
        moisture_profile = MoistLowTropicsMoistureProfile()
    end

    problem = (
        boundarycondition,
        init_gcm_experiment!,
        (_...) -> nothing,
        perturbation,
        base_state,
        moisture_profile,
    )
    return HeldSuarezProblem{typeof.(problem)...}(problem...)
end

problem_name(::HeldSuarezProblem) = "HeldSuarez"

setup_source(::HeldSuarezProblem) =
    (Gravity(), Coriolis(), held_suarez_forcing!)

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
"""
# Choose

  Displays a number menu in the matlab console.

## Syntax

   ichoose=gnChoose(header,itemlist)

## Description

 An OM_JPL Function.

 Displays a numbered menu on the metlab command window that allows the user
  to make a selection of one or more items.   The item number/s is returned.

 A user may cancel the selection, upon which [] is returned.

## INPUTS

   header          a title displayed at the top of the menu list
   itemlist        a cell array of strings to be displayed a numbered list

## OUTPUTS

   ichoose        the index of the selected menu item

## EXAMPLES

   o={ ...
       'Start MTB', ...
       'Start XWaves', ...
       '2013 Jobs', ...
       '2014 Jobs', ...
      };
    so=gnChoose('',o);

## Example
 To demonstrate the functionality, an example is shown below:

   t = 'Job type';
   o={ ...
       'Start MTB', ...
       'Start XWaves', ...
       '2013 Jobs', ...
       '2014 Jobs', ...
      };
    so=gnChoose(t,o);

## Revisions

 Written by: Jason McConochie
 Rev 1, 21-Dec-2016  Jason McConochie
"""
function gnChoose(header::AbstractString,itemlist::Tuple)

   write(STDOUT,"\n\n$(header)\n\nChoose:\n")
   n=length(itemlist)
   for i=1:n
       write(STDOUT,"$(i). $(itemlist[i])\n")
   end
   write(STDOUT,"$(n+1). Cancel\n")
   write(STDOUT,'\n')
   function input(prompt::AbstractString="")
        print(prompt)
        return chomp(readline())
   end

   ichoose=n+2
   while (ichoose > n+1) | (ichoose < 1)
       ichoose=parse(Int,input("Select Option: "))
   end
   return ichoose

end

"""
  maskReplace!

  Replaces values in array original with values in array replacements if replaceMask is true

  Usage: maskReplace!(original,replacements,replaceMask)

  Example:

```
  julia> original=[1,2,3,4,5]
  julia> replacements=[0,0,4,0,1]
  julia> replaceMask=[false,true,true,false,true]
  julia> maskReplace!(original,replacements,replaceMask)
  julia> original
  5-element Array{Int64,1}:
    1
    0
    4
    4
    1
  ```
"""
function maskReplace!{T<:Number}(original::Array{T,1},replacements::Array{T,1},replaceMask::Array{Bool,1})
  for i in eachindex(replaceMask)
  original[i]=!replaceMask[i]*original[i] + replaceMask[i]*replacements[i]
  end
end

"""
Array based Dates.value function definition
"""
function Dates.value(in::Array{DateTime,1})::Array{Int64,1}
   out=zeros(Int64,size(in,1))
   for i=1:length(in)
   out[i]=Dates.value(in[i])
   end
   return out
end

function gmAngDiff(this_angle,relative_to)
  ang1=this_angle
  ang2=relative_to
  #Make all angles to be between 0 and 360 degrees
  ang1=mod(ang1,360)
  ang2=mod(ang2,360)
  #Calculate angular difference as angle1-angle2
  angdiff=ang1-ang2
  #Identify angles > 180. These need to be sutracted from 360
  if angdiff>180
    angdiff=angdiff-360
  end
  #Identify angles <-180. These need to be sutracted from 360
  ismall=find(angdiff<=-180)
  if angdiff<=-180
    angdiff=angdiff+360
  end
  return angdiff
end

"""
function [minl,maxl,ticks]=ascale(minv,maxv,sType);
 Function will calculate nice axis limits based on the minimum and maximum values
  in the vectors minv and maxv

 Several options for how the limits are set is available based on the value of type
 's' - allows some extra space at either end of the minv and maxv

 Revision 1, 27Dec2010 Jason McConochie (original code by Luciano Mason)
"""
function ascale(minv::Float64,maxv::Float64,sType::AbstractString)
   #minv=real(minv)
   #maxv=real(maxv)
   #minv[isnan.(minv)]=[]
   #maxv[isnan.(maxv)]=[]
   #if isempty(minv) | isempty(maxv)
   #  minl=-1
   #  maxl=1
   #  ticks=[-1,0,1]
   #  return
   #end
   #if isnan(minv) | isnan(maxv)
   #  minl=-1
   #  maxl=1
   #  ticks=[-1,0,1]
   #  return
   #end

   # Vector of Number of tick marks wanted
   di=[5,6,7,8,9]
   m=length(di)
   # Vector of step sizes ie. 1,10,100 2,4,6,8 5,0.5,50 etc
   i=[1,2,5]
   n=length(i)
   if(abs(minv - maxv)<1e-8)
     minv=minv-0.1*abs(minv)
     maxv=maxv+0.1*abs(maxv)
   end
   if( (minv==0) & (maxv==0) )
     minv=-1
     maxv=1
   end

   # Make vector of intervals based on di
   range=abs(maxv-minv)
   # ai=range./di;
   # Work in log space
   li=log10(i)
   try
     lai=log10(range)-log10(di)
   catch
     minl=minv
     maxl=maxv
     ticks=[minl,maxl]
   end
   li=repmat(li',m,1)
   lai=repmat(lai,1,n)
   # Find best combination of Number of ticks/Interval
   r=abs(rem(lai-li,1))
   #r=abs(lai-li)
   minr=minimum(r[:])
   #r=rem(r,1);minr=rem(minr,1);

   imin,jmin=findn(r.==minr)
   #println(['No. of tick gaps=' num2str(di(imin)) ' Step size=' num2str(i(jmin))]);
   i=imin[end]
   j=jmin[end]

   di=lai[i,j]-li[i,j]
   fdi=floor(abs(di))*sign(di)
   interval=10.^(fdi+li[i,j])

   if sType=="e"
     llim=floor(minv/interval)*interval
     ulim=ceil(maxv/interval)*interval
     ticks=[x for x=llim:interval:ulim]
     minl=ticks[1]
     maxl=ticks[end]
   end
   if sType=="i"
     llim=ceil(minv/interval)*interval
     ulim=floor(maxv/interval)*interval
     ticks=[x for x=llim:interval:ulim]
     minl=minv
     maxl=maxv
   end
   if sType=="s"
     llim=ceil((minv-0.05*range)/interval)*interval
     ulim=floor((maxv+0.05*range)/interval)*interval
     ticks=[x for x=llim:interval:ulim]
     minl=minv-0.05*range
     maxl=maxv+0.05*range
   end

   return minl,maxl,ticks
end

max_stars_repo_name: jadamc/MetoceanTools

max_issues_repo_name: jadamc/MetoceanTools

max_forks_repo_name: jadamc/MetoceanTools


text:
using BinaryBuilder

cuda_version = v"11.0.2"   # NOTE: could be less specific

script = raw"""
mkdir -p ${libdir} ${prefix}/include

cd ${WORKSPACE}/srcdir
if [[ ${target} == *-linux-gnu ]]; then
    cd libcutensor
    find .

    install_license license.pdf

    mv lib/11/libcutensor.so* ${libdir}
    mv include/* ${prefix}/include
elif [[ ${target} == x86_64-w64-mingw32 ]]; then
    apk add p7zip
    7z x libcutensor*.exe
    cd libcutensor

    cd libcutensor
    find .

    install_license license.pdf

    mv lib/11/cutensor.dll ${libdir}
    mv include/* ${prefix}/include
fi
"""

include("../common.jl")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
using Test
include("../src/thermo.jl")

@testset "Thermodynamic Functions" begin
    @test isapprox(temperature(290, 1004), 290)
    @test isapprox(temperature(290, 1004, celsius=true), 16.85)

    @test isapprox(pressure(1004), 1000)

    @test isapprox(wsat([270,290], [1000, 1004]), [0.0028296776778632955, 0.012165289848931712])

    @test isapprox(rh([0.002,0.002], [270,290], [1000, 1004]), [70.6794281075225, 16.4402165902823])

    @test isapprox(mslp(290, 950, 1000), 1009.07759)
end


text:
@testset "$C" for C in (Lab, Luv)
    @testset "$C over $C: $T" for T in (Float64, Float32)
        @test BlendNormal(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(60, 50, -40)
        @test BlendNormal(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(40, -50, 60)
        @test BlendMultiply(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(54, 98.75, 51.875)
        @test BlendMultiply(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(4, -62.1875, 15.9375)
        @test BlendScreen(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(96, 98.75, 51.875)
        @test BlendScreen(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(46, -62.1875, 15.9375)
        @test BlendOverlay(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(92, 98.75, 51.875)
        @test BlendOverlay(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(8, -62.1875, 15.9375)
        @test BlendHardLight(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(92, 98.75, 8.125)
        @test BlendHardLight(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(8, -62.1875, 44.0625)
        @test BlendSoftLight(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(90.97366596101027, 88.88878188659973, 55.9814453125)
        @test BlendSoftLight(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(8.2, -44.3896484375, -1.6912879754125099)
        @test BlendHue(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(90, 80, -64)
        @test BlendHue(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(10, -25, 30)
        @test BlendSaturation(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(90, 50, 43.75)
        @test BlendSaturation(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(10, -40, -60)
        @test BlendColor(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(90, 50, -40)
        @test BlendColor(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(10, -50, 60)
        @test BlendLuminosity(C{T}(90,  80,  70), C{T}(60,  50, -40)) ≈ C{T}(60, 80, 70)
        @test BlendLuminosity(C{T}(10, -20, -30), C{T}(40, -50,  60)) ≈ C{T}(40, -20, -30)
    end

    A = coloralpha(C)
    @testset "$A over $C: $T" for T in (Float64, Float32)
        @test BlendNormal(C{T}(90, 80, 70), A{T}(60, 50, -40, 0.6)) ≈ C{T}(72, 62, 4)
    end

    @testset "$C over $A: $T" for T in (Float64, Float32)
        @test BlendNormal(A{T}(90, 80, 70, 0.6), C{T}(60, 50, -40)) ≈ A{T}(60, 50, -40, 1)
    end

    @testset "$A over $A: $T" for T in (Float64, Float32)
        @test BlendNormal(A{T}(90, 80, 70, 0.6), A{T}(60, 50, -40, 0.6)) ≈ A{T}(68.57142857142857, 58.57142857142857, -8.57142857142857, 0.84)
    end
end

@testset "$C" for C in (LCHab, LCHuv)
    @testset "$C over gray $C: $T" for T in (Float64, Float32)
        @test blend(C{T}(40, 0, 100), C{T}(60, 100, 200), opacity=0.5) ≈ C{T}(50, 50, 200)
    end
end


text:
@testset "data and model node io" begin
    an = ArrayNode(ones(2,5))
    @test repr(an) == "ArrayNode(2×5 Array with Float64 elements)"
    @test repr(an; context=:compact => true) == "ArrayNode"
    @test repr(MIME("text/plain"), an) ==
        """
        2×5 ArrayNode{Matrix{Float64}, Nothing}:
         1.0  1.0  1.0  1.0  1.0
         1.0  1.0  1.0  1.0  1.0"""

    anm = reflectinmodel(an)
    @test repr(anm) == "ArrayModel(Dense(2, 10))"
    @test repr(anm; context=:compact => true) == "ArrayModel"
    @test repr(MIME("text/plain"), anm) == repr(anm)

    bn = BagNode(an, bags([1:2, 3:5]))
    @test repr(bn) == "BagNode with 2 obs"
    @test repr(bn; context=:compact => true) == "BagNode"
    @test repr(MIME("text/plain"), bn) ==
        """
        BagNode with 2 obs
          └── ArrayNode(2×5 Array with Float64 elements)"""

    bnm = reflectinmodel(bn, d -> Dense(d, 10), SegmentedMean)
    @test repr(bnm) == "BagModel … ↦ SegmentedMean(10) ↦ ArrayModel(Dense(10, 10))"
    @test repr(bnm; context=:compact => true) == "BagModel"
    @test repr(MIME("text/plain"), bnm) ==
        """
        BagModel … ↦ SegmentedMean(10) ↦ ArrayModel(Dense(10, 10))
          └── ArrayModel(Dense(2, 10))"""

    wbn = WeightedBagNode(an, bags([1:2, 3:5]), rand(5))
    @test repr(wbn) == "WeightedBagNode with 2 obs"
    @test repr(wbn; context=:compact => true) == "WeightedBagNode"
    @test repr(MIME("text/plain"), wbn) ==
        """
        WeightedBagNode with 2 obs
          └── ArrayNode(2×5 Array with Float64 elements)"""

    wbnm = reflectinmodel(wbn)
    @test repr(wbnm) == "BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))"
    @test repr(wbnm; context=:compact => true) == "BagModel"
    @test repr(MIME("text/plain"), wbnm) ==
        """
        BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))
          └── ArrayModel(Dense(2, 10))"""

    pn = ProductNode((; bn, wbn))
    @test repr(pn) == "ProductNode with 2 obs"
    @test repr(pn; context=:compact => true) == "ProductNode"
    @test repr(MIME("text/plain"), pn) ==
        """
        ProductNode with 2 obs
          ├─── bn: BagNode with 2 obs
          │          └── ArrayNode(2×5 Array with Float64 elements)
          └── wbn: WeightedBagNode with 2 obs
                     └── ArrayNode(2×5 Array with Float64 elements)"""

    pnm = reflectinmodel(pn, d -> Dense(d, 10), BagCount ∘ SegmentedMean)
    @test repr(pnm) == "ProductModel … ↦ ArrayModel(Dense(20, 10))"
    @test repr(pnm; context=:compact => true) == "ProductModel"
    @test repr(MIME("text/plain"), pnm) ==
        """
        ProductModel … ↦ ArrayModel(Dense(20, 10))
          ├─── bn: BagModel … ↦ BagCount(SegmentedMean(10)) ↦ ArrayModel(Dense(11, 10))
          │          └── ArrayModel(Dense(2, 10))
          └── wbn: BagModel … ↦ BagCount(SegmentedMean(10)) ↦ ArrayModel(Dense(11, 10))
                     └── ArrayModel(Dense(2, 10))"""

    ln = LazyNode{:Sentence}(["a", "b", "c", "d"])
    @test repr(ln) == "LazyNode{Sentence} with 4 obs"
    @test repr(ln; context=:compact => true) == "LazyNode"
    @test repr(MIME("text/plain"), ln) == repr(ln)

    lnm = reflectinmodel(ln)
    @test repr(lnm) == "LazyModel{Sentence}"
    @test repr(lnm; context=:compact => true) == "LazyModel"
    @test repr(MIME("text/plain"), lnm) ==
        """
        LazyModel{Sentence}
          └── BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))
                └── ArrayModel(Dense(2053, 10))"""
end

@testset "aggregation io" begin
    a = SegmentedMax(3)
    @test repr(a) == "SegmentedMax(3)"
    @test repr(a; context=:compact => true) == "SegmentedMax"
    @test repr(MIME("text/plain"), a) == "SegmentedMax(ψ = Float32[0.0, 0.0, 0.0])"

    a = AggregationStack(
                    SegmentedMean(2),
                    SegmentedMax(2),
                    AggregationStack(SegmentedPNorm(zeros(2) |> f32, ones(2) |> f32, -ones(2) |> f32)),
                    AggregationStack(SegmentedLSE(zeros(2) |> f32, ones(2) |> f32))
                   )
    @test repr(a) == "[SegmentedMean(2); SegmentedMax(2); SegmentedPNorm(2); SegmentedLSE(2)]"
    @test repr(a; context=:compact => true) == "AggregationStack"
    @test repr(MIME("text/plain"), a) ==
        """
        AggregationStack:
         SegmentedMean(ψ = Float32[0.0, 0.0])
         SegmentedMax(ψ = Float32[0.0, 0.0])
         SegmentedPNorm(ψ = Float32[0.0, 0.0], ρ = Float32[1.0, 1.0], c = Float32[-1.0, -1.0])
         SegmentedLSE(ψ = Float32[0.0, 0.0], ρ = Float32[1.0, 1.0])"""

    bc = BagCount(SegmentedMean(3))
    @test repr(bc) == "BagCount(SegmentedMean(3))"
    @test repr(bc; context=:compact => true) == "BagCount"
    @test repr(MIME("text/plain"), bc) == "BagCount(SegmentedMean(ψ = Float32[0.0, 0.0, 0.0]))"

    bc = BagCount(SegmentedMeanMax(3))
    @test repr(bc) == "BagCount([SegmentedMean(3); SegmentedMax(3)])"
    @test repr(bc; context=:compact => true) == "BagCount"
    @test repr(MIME("text/plain"), bc) ==
        """
        BagCount(AggregationStack:
         SegmentedMean(ψ = Float32[0.0, 0.0, 0.0])
         SegmentedMax(ψ = Float32[0.0, 0.0, 0.0])
        )"""
end

@testset "maybe hot vector io" begin
    x = MaybeHotVector(1, 3)
    @test repr(x) == "MaybeHotVector(i = 1, l = 3)"
    @test repr(x; context=:compact => true) == "3-element MaybeHotVector"
    @test repr(MIME("text/plain"), x) ==
        """
        3-element MaybeHotVector{Int64, Int64, Bool}:
         1
         0
         0"""

    x = MaybeHotVector(missing, 4)
    @test repr(x) == "MaybeHotVector(i = missing, l = 4)"
    @test repr(x; context=:compact => true) == "4-element MaybeHotVector"
    @test repr(MIME("text/plain"), x) ==
        """
        4-element MaybeHotVector{Missing, Int64, Missing}:
         missing
         missing
         missing
         missing"""
end

@testset "maybe hot matrix io" begin
    X = MaybeHotMatrix([1, 3], 3)
    @test repr(X) == "MaybeHotMatrix(I = [1, 3], l = 3)"
    @test repr(X; context=:compact => true) == "3×2 MaybeHotMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        3×2 MaybeHotMatrix{Int64, Int64, Bool}:
         1  0
         0  0
         0  1"""

    X = MaybeHotMatrix([missing, 2], 2)
    @test repr(X) == "MaybeHotMatrix(I = Union{Missing, Int64}[missing, 2], l = 2)"
    @test repr(X; context=:compact => true) == "2×2 MaybeHotMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2×2 MaybeHotMatrix{Union{Missing, Int64}, Int64, Union{Missing, Bool}}:
         missing  false
         missing   true"""

    X = MaybeHotMatrix([missing, missing, missing], 4)
    @test repr(X) == "MaybeHotMatrix(I = [missing, missing, missing], l = 4)"
    @test repr(X; context=:compact => true) == "4×3 MaybeHotMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        4×3 MaybeHotMatrix{Missing, Int64, Missing}:
         missing  missing  missing
         missing  missing  missing
         missing  missing  missing
         missing  missing  missing"""
end

@testset "pre imputing matrix io" begin
    X = PreImputingMatrix(reshape(1:8, 4, 2) |> Matrix)
    @test repr(X) == "PreImputingMatrix(W = [1 5; 2 6; 3 7; 4 8], ψ = [0, 0])"
    @test repr(X; context=:compact => true) == "4×2 PreImputingMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        4×2 PreImputingMatrix{Int64, Matrix{Int64}, Vector{Int64}}:
        W:
         1  5
         2  6
         3  7
         4  8

        ψ:
         0  0"""
end

@testset "post imputing matrix io" begin
    X = PostImputingMatrix(reshape(1:6, 2, 3) |> Matrix)
    @test repr(X) == "PostImputingMatrix(W = [1 3 5; 2 4 6], ψ = [0, 0])"
    @test repr(X; context=:compact => true) == "2×3 PostImputingMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2×3 PostImputingMatrix{Int64, Matrix{Int64}, Vector{Int64}}:
        W:
         1  3  5
         2  4  6

        ψ:
         0
         0"""
end

@testset "imputing dense io" begin
    X = PostImputingMatrix(reshape(1:6, 2, 3) |> Matrix)
    Y = PreImputingMatrix(reshape(1:8, 4, 2) |> Matrix)

    m = preimputing_dense(1, 1)
    @test repr(m) == "[pre_imputing]Dense(1, 1)"
    @test repr(m) == repr(m; context=:compact => true) == repr(MIME("text/plain"), m)

    m = postimputing_dense(1, 1)
    @test repr(m) == "[post_imputing]Dense(1, 1)"
    @test repr(m) == repr(m; context=:compact => true) == repr(MIME("text/plain"), m)
end

@testset "ngram matrix io" begin
    X = NGramMatrix(["a", "b", "c"])
    @test repr(X) == "NGramMatrix(S = [\"a\", \"b\", \"c\"], n = 3, b = 256, m = 2053)"
    @test repr(X; context=:compact => true) == "2053×3 NGramMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2053×3 NGramMatrix{String, Vector{String}, Int64}:
         \"a\"
         \"b\"
         \"c\""""

    X = NGramMatrix([missing, "b", missing] |> PooledArray)
    @test repr(X) == "NGramMatrix(S = Union{Missing, String}[missing, \"b\", missing], n = 3, b = 256, m = 2053)"
    @test repr(X; context=:compact => true) == "2053×3 NGramMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2053×3 NGramMatrix{Union{Missing, String}, PooledVector{Union{Missing, String}, UInt32, Vector{UInt32}}, Union{Missing, Int64}}:
         missing
         \"b\"
         missing"""

    X = NGramMatrix([codeunits("hello"), codeunits("world")])
    @test repr(X) == "NGramMatrix(S = CodeUnits{UInt8, String}" *
        "[[0x68, 0x65, 0x6c, 0x6c, 0x6f], [0x77, 0x6f, 0x72, 0x6c, 0x64]], n = 3, b = 256, m = 2053)"
    @test repr(X; context=:compact => true) == "2053×2 NGramMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2053×2 NGramMatrix{CodeUnits{UInt8, String}, Vector{CodeUnits{UInt8, String}}, Int64}:
         UInt8[0x68, 0x65, 0x6c, 0x6c, 0x6f]
         UInt8[0x77, 0x6f, 0x72, 0x6c, 0x64]"""

    X = NGramMatrix([[1,2,3], missing, missing])
    @test repr(X) == "NGramMatrix(S = Union{Missing, Vector{Int64}}[[1, 2, 3], missing, missing], n = 3, b = 256, m = 2053)"
    @test repr(X; context=:compact => true) == "2053×3 NGramMatrix"
    @test repr(MIME("text/plain"), X) ==
        """
        2053×3 NGramMatrix{Union{Missing, Vector{Int64}}, Vector{Union{Missing, Vector{Int64}}}, Union{Missing, Int64}}:
         [1, 2, 3]
         missing
         missing"""
end


text:
"""
    CausalAnalysis(test, data, result)

The result of a causal analysis.

## Fields

- **`test`**: The parameters of the test.
- **`data`**: The data on which the test was applied.
- **`result`**: The result of applying the relevant causality statistic using
the parameters in `test` to `data`.

## Implements

- [`summarise(f::Function, analysis::CausalAnalysis, args...; kwargs...)`](@ref). Summarise
a causal analysis using some summary statistic.
"""
struct CausalAnalysis{CT <: Union{CausalityTest, MetaCausalityTest}, DT, RT} <: AbstractCausalAnalysis{CT, DT, RT}
    test::CT
    data::DT
    result::RT
end

function Base.show(io::IO, x::AbstractCausalAnalysis)
    T = typeof(x).name

    CT = typeof(x.test)
    RT = typeof(x.result)
    DT = typeof(x.data)
    CT_str = "\n  test: $CT"
    DT_str = "\n  data: $DT"
    RT_str = "\n  result: $(RT)"
    summary = join([CT_str, DT_str, RT_str], "")
    print(io, "$(T)$summary\n")
end

"""
    causality(test::Union{CausalityTest, MetaCausalityTest}, x, y)
    causality(test::Union{CausalityTest, MetaCausalityTest}, x, y, z)

Apply the causality `test` to the provided data series. Providing the
test as the first argument returns a [`CausalAnalysis`](@ref)
instance which stores the test, the data used for the test, and the
result in a single struct.
"""
function causality(test::Union{CausalityTest, MetaCausalityTest}) end

function causality(test::Union{CausalityTest, MetaCausalityTest}, x, y)
    CausalAnalysis(test, (x, y), causality(x, y, test))
end

function causality(test::Union{CausalityTest, MetaCausalityTest}, x, y, z)
    CausalAnalysis(test, (x, y, z), causality(x, y, z, test))
end


text:
### ----------------------------------------------------------------------------
### ---- Basic testing for the Shifted Characteristic Polynomial Evaluation:
### ----------------------------------------------------------------------------

include("dynamic_symtridiagonal.jl")

using LinearAlgebra, Logging, Plots

"""
    Basic test on the characteristic polynomial for the tridiagonal Matrix.
"""
function Test1(n=50)
    mainDiag = rand(n)
    subDiag = rand(n - 1)
    dynamicT = DynamicSymTridiagonal(mainDiag[1])
    T = GetT(dynamicT)
    for Idx in 1:n - 1
        dynamicT(mainDiag[Idx + 1], subDiag[Idx])
    end
    @info "Characteristic poly without shift evaluated to be: "
    println(CharPolyShifted(dynamicT, 0))
    println(det(T))
    @info "Plot out the characteric polynomial over the min, max eigenvalue. "
    minEigen = eigmin(T); maxEigen = eigmax(T)
    grid = LinRange(minEigen, maxEigen, 1024)
    p = grid .|> (x) -> CharPolyShifted(dynamicT, x)
    p2 = grid .|> (x) -> det(T - x*I)
    fig = plot(grid, sign.(p).*log10.(abs.(p) .+ 1))
    display(plot!(fig, grid, sign.(p2).*log10.(abs.(p2) .+ 1)))

return dynamicT end

# Test1();

function Test2(n = 200)
    mainDiag = rand(n)
    subDiag = rand(n - 1)
    dynamicT = DynamicSymTridiagonal(mainDiag[1])
    T = GetT(dynamicT)
    for Idx in 1:n - 1
        dynamicT(mainDiag[Idx + 1], subDiag[Idx])
    end
    # try locating one of the inner eigenvalue of matrix T.
    @info "Trying out the EigenValueLocate function. "
    eigenVals = eigvals(T)
    sort!(eigenVals)
    P(x) = CharPolyShifted(dynamicT, x)
    locatedEig = EigenValueLocate(
        dynamicT,
        (eigenVals[end - 1] + eigenVals[end])/2,
        Inf)
    println("locatedEig: $locatedEig")
    println("Expected: $(eigenVals[end])")
    error = abs(locatedEig - eigenVals[end])
    println("error: $(error)")
    println("This is the lsit of all eigenvalues: ")
    display(eigenVals)
    println("The corresonding eigenvector is: ")
    LocatedEigenVec = eigvecs(T, [locatedEig])
    if error > 1e-10
        @warn "Error is bigger than expected. "
        # might over probe and cause big problem.
    end

return T, locatedEig, LocatedEigenVec end

# T, locatedEig, locatedEigenVec = Test2()
# display(plot(locatedEigenVec))

function Test3(n=2048)
    mainDiag = -2*ones(n)
    subDiag = ones(n - 1)
    dynamicT = DynamicSymTridiagonal(mainDiag[1])
    ChangeVelocityTol!(dynamicT, 1e-4)
    referenceT = SymTridiagonal(mainDiag, subDiag)
    for Idx in 1:n÷16
        dynamicT(mainDiag[Idx + 1], subDiag[Idx])
        EigenvaluesUpdate(dynamicT)
        if dynamicT.converged_this_step |> length >= 1
            println("step $(Idx) has: $(dynamicT.converged_this_step) converged. ")
        end
    end
    @info "Ritz System: "
    display(dynamicT.thetas)
    display(dynamicT.converged)
    println("Converged $(CountConverged(dynamicT))")
    # Plot the converged that are identified and the correct eigenvalues,
    # curious how good cauchy interlace is.

return end

Test3()