text:
"A AbstractModel specifies all parameters of a system."
abstract type AbstractModel{T,N,P} <: AbstractParameter{T} end
abstract type AbstractODEModel{T,N,P} <: AbstractModel{T,N,P} end
abstract type AbstractModelwithDelay{T,N,P} <: AbstractModel{T,N,P} end
abstract type AbstractNoisyModel{T,N,P} <: AbstractModel{T,N,P} end
abstract type AbstractNoisySpaceAction{T,N} <: AbstractSpaceAction{T,N} end
abstract type AbstractSimulation{T} <: AbstractParameter{T} end
abstract type AbstractExecution{T,SIM} end
# FIXME shouldn't the second type of full execution be the type of the solution?
abstract type AbstractFullExecution{T,SIM} <: AbstractExecution{T,SIM} end
export AbstractFullExecution, AbstractModel, AbstractSimulation

#function (am::Type{<:AbstractModel})(fallback_args...; fallback_kwargs...)
#    @warn """Model $am undefined!
#    ---------------------
#    $fallback_args
#
#    $fallback_kwargs
#    ---------------------
#    """
#    missing
#end

struct FailedSimulation{T} <: AbstractSimulation{T} end
struct FailedExecution{T,S<:FailedSimulation{T}} <: AbstractExecution{T,S}
    sim::S
end

# export AbstractNoisyModel, NoisyInputModel, NoisyInputAction, AbstractODEModel

# struct NoisyInputModel{T,N,P,M<:AbstractModel{T,N,P}} <: AbstractNoisyModel{T,N,P}
#     model::M
#     noise_process::NoiseProcess
# end
# function Base.getproperty(wnm::NoisyInputModel, sym::Symbol)
#     if sym ∈ [:noise_process, :model]
#         return getfield(wnm, sym)
#     else
#         return getproperty(getfield(wnm, :model), sym)
#     end
# end

# struct NoisyInputAction{T,N,SA<:AbstractSpaceAction{T,N}} <: AbstractNoisySpaceAction{T,N}
#     space_action::SA
#     noise_process::NoiseProcess
# end
# function (act::NoisyInputAction)(du, u, p, t, W)
#     du .= W
#     act.space_action(du, u, p, t)
# end

# function (nim::NoisyInputModel)(args...)
#     inner_fn = nim.model(args...)
#     NoisyInputAction(inner_fn, nim.noise_process)
# end

n_populations(::AbstractModel{T,N,P}) where {T,N,P} = P

function initial_value(::AbstractModel{T,N,P}, space::AbstractSpace{T,N}) where {T,N,P}
    init_val = AxisArray(population_repeat(zeros(space), P), (coordinate_axes(space)..., 1:P))
    return init_val
end

"A Simulation holds an AbstractModel to be solved, the space on which to solve it, the time for which to solve it, the initial value, and various solver options."
struct Simulation{
        T,
        M<:AbstractModel{T},
        S<:AbstractSpace{T},
        IV<:AbstractArray{T},
        ALG,
        DT<:Union{T,Nothing},
        SV_IDX<:Union{AbstractArray,Nothing},
        CB<:Union{Tuple{Function,NamedTuple},Tuple{DECallback,NamedTuple},Nothing,Vector{Tuple{Function,NamedTuple}}},
        GR<:Function} <: AbstractSimulation{T}
    model::M
    space::S
    tspan::Tuple{T,T}
    initial_value::IV
    algorithm::ALG
    dt::DT
    save_idxs::SV_IDX
    callback::CB
    global_reduction::GR
    solver_options
end
function Simulation(model::M; space::S, tspan, initial_value::IV=initial_value(model,space),
                    algorithm::ALG=nothing, dt::DT=nothing, save_idxs::SV_IDX=nothing, callback::CB=nothing, global_reduction::GR=identity,
        opts...) where {T,N,P,M<:AbstractModel{T,N,P}, S<:AbstractSpace{T,N},IV,ALG,DT,SV_IDX,CB,GR}
    save_idxs = parse_save_idxs(space, P, save_idxs)
    return Simulation{T,M,S,IV,ALG,DT,typeof(save_idxs),CB,GR}(model, space, tspan, initial_value, algorithm, dt, save_idxs, callback, global_reduction, opts)
end
function Simulation(model::Missing; kwargs...)
    return FailedSimulation{Missing}()
end

"An Execution holds a Simulation and the solution obtained by running the Simulation."
struct Execution{T,S<:AbstractSimulation{T},D<:DESolution} <: AbstractFullExecution{T,S}
    simulation::S
    solution::D
end
struct ReducedExecution{T,ST,S<:AbstractSimulation{T},SV<:SavedValues{T,ST}} <: AbstractExecution{T,S}
    simulation::S
    saved_values::SV
end
struct AugmentedExecution{T,ST,S<:AbstractSimulation{T},D<:DESolution,SV<:SavedValues{T,ST}} <: AbstractFullExecution{T,S}
    simulation::S
    solution::D
    saved_values::SV
end
make_execution(s::S, sol::DESolution) where {T,S <: Simulation{T}} = Execution(s,sol)
make_execution(s::S, sv::SavedValues) where {T,S <: Simulation{T}} = ReducedExecution(s,sv)
make_execution(s::S, (sol,sv)::Tuple{<:DESolution,<:SavedValues}) where {T,S <: Simulation{T}} = AugmentedExecution(s,sol,sv)
export ReducedExecution, AugmentedExecution,make_execution
function execute(s::Simulation)
    return make_execution(s, solve(s))
end
function execute(s::FailedSimulation)
    return FailedExecution(s)
end

coordinates(sim::Simulation) = coordinates(space(sim))
coordinates(ex::AbstractExecution) = coordinates(space(ex))
timepoints(ex::AbstractFullExecution) = ex.solution.t
reduced_space(sim::Simulation) = subsample(sim.space, sim.save_idxs)
reduced_space(ex::AbstractExecution) = reduced_space(ex.simulation)
frame_xs(exec::AbstractExecution{T,<:Simulation{T,<:M}}) where {T,M<:AbstractModel{T,1}} = [x[1] for x in reduced_space(exec).arr]
frame_xs(sim::Simulation) = [x[1] for x in reduced_space(sim).arr]
origin_idx(sim::Simulation) = origin_idx(sim.space)
origin_idx(ex::AbstractExecution) = origin_idx(ex.simulation)
Base.extrema(ex::AbstractFullExecution) = extrema.(ex.solution.u) |> extrema_arr -> reduce(((max1, min1), (max2, min2)) -> (max(max1, max2), min(min1, min2)), extrema_arr)
Base.maximum(ex::AbstractFullExecution) = maximum(maximum.(ex.solution.u))
Base.minimum(ex::AbstractFullExecution) = minimum(minimum.(ex.solution.u))
stimulus_center(mod::AbstractModel) = center(mod.stimulus)
stimulus_center(sim::Simulation) = stimulus_center(sim.model)
export stimulus_center

"""
    make_system_mutator(simulation)

Construct the differential function to be provided to the ODE solver.
"""
function make_system_mutator(sim::SIM) where {SIM <: Simulation}
    sim.model(sim.space)
end

"""
    generate_problem(simulation)

Return an ODEProblem of the `simulation.model` with time span specified by `simulation.solver`.
"""

function generate_problem(simulation::Simulation{T,<:AbstractODEModel}, p::NamedTuple, callback::DECallback) where {T}
    system_fn! = make_system_mutator(simulation)# simulation.model(simulation.space)
    ode_fn = convert(ODEFunction{true}, system_fn!)
    return ODEProblem(ode_fn, simulation.initial_value, simulation.tspan, p;callback=callback)
end
function generate_problem(simulation::Simulation{T,<:AbstractODEModel}, callback::Nothing, p::NamedTuple) where {T}
    system_fn! = make_system_mutator(simulation)# simulation.model(simulation.space)
    ode_fn = convert(ODEFunction{true}, system_fn!)
    return ODEProblem(ode_fn, simulation.initial_value, simulation.tspan, p)
end

# function generate_problem(simulation::Simulation{T,<:AbstractNoisyModel}; callback=nothing) where {T}
#     system_fn! = make_system_mutator(simulation)# simulation.model(simulation.space)
#     return RODEProblem(system_fn!, simulation.initial_value, simulation.tspan, noise=simulation.model.noise_process, noise_prototype=zeros(size(simulation.initial_value)...), callback=callback)
# end

# TODO: Add history functionality
# function generate_problem(simulation::Simulation{T,MwD}) where {T, MwD<:AbstractModelwithDelay}
#     system_mutator! = make_system_mutator(simulation)
#     return DDEProblem(system_mutator!, simulation.initial_value, history(simulation), simulation.tspan)
# end
parse_save_idxs(::Any, P, save_idx) = save_idx
function parse_save_idxs(space::AbstractSpace, P, subsampler::Union{AbstractSubsampler,AbstractArray{<:AbstractSubsampler}})
	one_pop_coordinates = coordinate_indices(space, subsampler)
	population_coordinates(one_pop_coordinates, P)
end
export handle_callback
function handle_callback(sim::Simulation)
    return handle_callback(sim, sim.callback)
end
function handle_callback(sim::Simulation, cb::CB) where {CB<:Nothing}
    return ((;), nothing)
end
function handle_callback(sim::Simulation, unhandled_cbs::CB) where {CB<:Vector{<:Tuple}}
    cb_nts = map(unhandled_cbs) do unhandled_cb
        nt::NamedTuple, cb::DECallback = handle_callback(sim, unhandled_cb)
    end
    grand_cb = CallbackSet([cb_nt[2] for cb_nt in cb_nts]...)
    grand_nt = merge([cb_nt[1] for cb_nt in cb_nts]...) # FIXME should reduce
    return (grand_nt, grand_cb)
end

# TODO this could probably all be better served by dispatch; but that gets difficult
function _solve(simulation::Simulation, alg; dt=nothing, solver_options...)
    if dt != nothing
        solver_options = (solver_options..., dt=dt)
    end
    p, callback = handle_callback(simulation)
    problem = generate_problem(simulation, p, callback)
    solve(problem, alg; solver_options...)
end
function _solve(simulation::Simulation, alg, ensemble_alg; prob_func::Function, output_func=(sol, i) -> (simulation.global_reduction(sol), false), reduction=(u,data,i) -> (append!(u,data), false), dt=nothing, u_init=[], solver_options...)
    if dt != nothing
        solver_options = (solver_options..., dt=dt)
    end
    p, callback = handle_callback(simulation)
    initial_problem = generate_problem(simulation, p, callback)
    ensemble_problem = EnsembleProblem(initial_problem;
                    output_func=output_func,
                    prob_func=prob_func,
                    reduction=reduction,
                    u_init=u_init)
                    @show solver_options
    solve(ensemble_problem, alg, ensemble_alg; solver_options...)
end

function solve(simulation::Simulation, args...; kwargs...)
    sol = _solve(simulation, simulation.algorithm, args...; dt=simulation.dt, save_idxs=simulation.save_idxs, simulation.solver_options..., kwargs...)
    return sol
end


text:
using WavePropagation
using Test

@testset "WavePropagation.jl" begin
    # Write your tests here.
end

ext: jl


text:
using BrowseTables, Tables, Test
using BrowseTables: escape_string, write_tags, make_cell, write_html

@testset "HTML escapes" begin
    @test escape_string("&<>\"test") == "&<>"test"
end

@testset "tag writing" begin
    io = IOBuffer()
    write_tags(io, "foo"; attributes = (bar = "baz", )) do io
        print(io, "test")
    end
    @test String(take!(io)) == "<foo bar=\"baz\">test</foo>"
end

"For testing `write_html` and representations."
repr_html(args...) = (io = IOBuffer(); write_html(io, args...); String(take!(io)))

const opt = TableOptions()

"For testing cell representations. Uses `<td>` for tags."
repr_cell(x, o = opt) = repr_html(:td, make_cell(o, x))

@testset "cell formatting" begin
    @test repr_cell(1) == "1	"
    if VERSION < v"1.2-"
        @test repr_cell(π) == "π = 3.1415926535897...	"
    else
        @test repr_cell(π) == "π	"
    end
    @test repr_cell(missing) == "missing	"
    @test repr_cell(nothing) == "nothing	"
    @test repr_cell("string<>") == "string<>	"
end

@testset "rudimentary check for whole file" begin
    filename = tempname() * ".html"
    @info "generating HTML for table" filename
    tb = (a = 1:2, b = [missing, 3])
    @test Tables.schema(tb) ≢ nothing # if this fails, interface changed, rewrite test
    write_html_table(filename, tb)
    html = read(filename, String)
    schema_html = "
#	a	b
"
    @test occursin("
$(schema_html)
", html)
    @test occursin("
$(schema_html)
", html)
    @test occursin("
\n
1	1	" *
                   "missing
\n
" *
                   "2	2	3
\n
", html)
end

@testset "HTMLTable" begin
    tb = (a = 1:2, b = [missing, 3])
    html = sprint(show, "text/html", HTMLTable(tb))
    @test occursin("
", html)
    @test occursin("
", html)
    @test occursin("
", html)
end


text:
using Plots, LaTeXStrings
gr()
fontFamily = "Times"
default(titlefont=(18, fontFamily), guidefont=(18, fontFamily),
      tickfont=(16, fontFamily),legendfont=(13, fontFamily),
      bg_legend = nothing,framestyle = :semi, grid = true)
      # another helpful keyword is:  size= (600, 400)
# ...............................................................#
x = LinRange(0,2π,100)
plot(x, sin.(x), c = :red,label = L"sin(x)", title = "Sine",
    yticks = ([-1,0,1], [L"-1", L"0", L"1"]), ylab = L"f(x)",
    xlab = L"x", xticks = ([0,π/2, π, 3π/2, 2π],
    [L"0", L"\pi/2", L"\pi", L"3\pi/2", L"2\pi"]))
savefig("figSin2.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
module sparsematrix

typealias SparseMatrix SparseMatrixCSC{Complex128}

function gemm!(alpha::Complex128, M::SparseMatrix, B::Matrix{Complex128}, beta::Complex128, result::Matrix{Complex128})
    @inbounds for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for j=1:M.m
        for colindex = 1:M.n
            m2 = alpha*B[colindex, j]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i], j] += M.nzval[i]*m2
            end
        end
    end
    nothing
end

function gemm!(alpha::Complex128, B::Matrix{Complex128}, M::SparseMatrix, beta::Complex128, result::Matrix{Complex128})
    @inbounds for j=1:size(result,2), i=1:size(result,1)
        result[i,j] *= beta
    end
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            mi = M.nzval[i]*alpha
            mrowvali = M.rowval[i]
            for j=1:size(result,1)
                result[j,colindex] += mi*B[j,mrowvali]
            end
        end
    end
    nothing
end

function gemv!(alpha::Complex128, M::SparseMatrix, v::Vector{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n
        vj = alpha*v[colindex]
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[M.rowval[i]] += M.nzval[i]*vj
        end
    end
    nothing
end

function gemv!(alpha::Complex128, v::Vector{Complex128}, M::SparseMatrixCSC{Complex128}, beta::Complex128, result::Vector{Complex128})
    @inbounds for i=1:length(result)
        result[i] *= beta
    end
    @inbounds for colindex = 1:M.n
        for i=M.colptr[colindex]:M.colptr[colindex+1]-1
            result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
        end
    end
    nothing
end

end # module


text:
@testset "Profiles" begin

σ = rand()
g = Gaussian(σ,1)
@test g(σ) ≈ exp(-1)/sqrt(π*σ^2)

@test SpaceTimeFields.radius(g) == σ
@test SpaceTimeFields.strength(g) == 1

σ = 0.2
A = 1
dg = DGaussian(σ,A)
x = 0.1
@test dg(x) == -2*A*x/sqrt(π)/σ^3*exp(-x^2/σ^2)

end

@testset "SpatialFields" begin

σ = 0.2
A = 1

g = EmptySpatialField()
@test g(rand(2)...) == 0.0

g = SpatialGaussian(σ,σ,0.0,0.5,1)
@test g(0,0.5+σ) ≈ g(σ,0.5) ≈ g(-σ,0.5) ≈ g(0,0.5-σ) ≈ exp(-1)/(π*σ^2)

u = 1
v = 0
gc = SpatialGaussian(σ,σ,0.0,0.5,1,u,v)
t = 1
@test gc(u*t,0.5+σ+v*t,t) ≈ gc(σ+u*t,0.5+v*t,t) ≈ gc(-σ+u*t,0.5+v*t,t) ≈ gc(0+u*t,0.5-σ+v*t,t) ≈ exp(-1)/(π*σ^2)

g = EmptySpatialField()
for x in [-0.5,0,0.5], y in [-0.5,0,0.5]
  g += SpatialGaussian(0.2,0.5,x,y,1)
end
@test g(0.5,0) ≈ 5.5357580598289235

g2 = -g
@test g2(0.5,0) ≈ -5.5357580598289235

end

@testset "Derivatives of gaussians" begin

  σx = 0.5
  σy = 0.5
  x0 = 0
  y0 = 0
  A = 1
  dgaussx = SpatialGaussian(σx,σy,x0,y0,A,deriv=1)
  dgaussy = SpatialGaussian(σx,σy,x0,y0,A,deriv=2)
  x = 0.1
  y = 0.2
  @test dgaussx(0.1,0.2) == dgaussy(0.2,0.1) ≈ -2*A*x/π/σx^3/σy*exp(-x^2/σx^2)*exp(-y^2/σy^2)

end


text:
"""
Return a dictionary with the file names as keys and
field values as values.
"""
function features{T<:AbstractMolecularContainer}(x::T)
    Dict([(k, getfield(x, k)) for k in fieldnames(x)])
end


text:
using DataFrames
using SLiDE

"""
    src/parse/check_data.jl [see file for function documentation]
These functions that can be used to compare DataFrames with the same column names.
Unfortunately, this currently only works for DataFrames.
1. compare_summary() - row-by row comparison of keys and values.
2. compare_keys() - are non-value columns consistent?
3. compare_values() - are values consistent within a given tolerance?
"""

<<<<<<< HEAD
println("\nEX 4: COMPARE DATA")
println("  Initializing DataFrames to compare...")

=======
>>>>>>> dev
# Make three different DataFrames to compare.
N = 3
dfa = DataFrame(yr = sort(repeat(2020-(N-1):2020, outer=[2])),
                r = repeat(["co","wi"], outer=[N]),
                value = Float64.(1:N*2))
dfb = copy(dfa)
dfc = copy(dfa)

# Change keys.
dfb[2,:r] = "md"
dfc[3,:r] = "Co"
dfc = edit_with(copy(dfa), Drop(:r, "co", "=="))

# Change values.
dfb[end-1,:value] *= 2.     # Values will not be equal regardless of tolerance.
dfc[end,:value] *= 1.001    # Values could be considered equal depending on tolerance.

# Apply functions.
df_lst = copy.([dfa,dfb,dfc])
inds = [:a,:b,:c];

# Print the three input DataFrames and results of the three comparison functions.
[(println("\ndf", ind, " = "); display(df)) for (df, ind) in zip(df_lst, inds)]

# SUMMARY: Print a summary comparison of the two DataFrames
df6_summary = compare_summary(df_lst, inds)
df2_summary = compare_summary(df_lst, inds; tol = 1E-2)

# EQUAL KEYS? Print a DataFrame summarizing difference between descriptive,
# "key" DataFrame columns.
df_keys = compare_keys(df_lst, inds)

# EQUAL VALUES? Print the rows of the summary DataFrame where values are unequal.
df6_values = compare_values(df_lst, inds)
df2_values = compare_values(df_lst, inds; tol = 1E-2)
println("")

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
mutable struct MetaGraph{T} <: AbstractMetaNexus{T}
    core::Graph{T}
    nprops::Dict{T,Dict{Symbol,Any}}
    aprops::Dict{Arc{T},Dict{Symbol,Any}}
end

MetaGraph{T}() where {T} = MetaGraph(Graph{T}())
function MetaGraph(core::Graph{T}) where {T}
    MetaGraph{T}(core, Dict{T,Dict{Symbol,Any}}(), Dict{Arc{T},Dict{Symbol,Any}}())
end
is_directed(::MetaGraph) = false
Graph(g::MetaGraph{T}) where {T} = g.core

function set_props!(g::MetaGraph{T}, a::Arc{T}, d::Dict) where {T}
    (has_arc(g, a) && has_arc(g, reverse(a))) || return false
    !_hasdict(g, a) ? g.aprops[a] = d : merge!(g.aprops[a], d)
    !_hasdict(g, reverse(a)) ? g.aprops[reverse(a)] = d : merge!(g.aprops[reverse(a)], d)
    return true
end

function rem_prop!(g::MetaGraph{T}, a::Arc{T}, prop::Symbol) where {T}
    delete!(g.aprops[a], prop)
    delete!(g.aprops[reverse(a)], prop)
end

function clear_props!(g::MetaGraph{T}, a::Arc{T}) where {T}
    _hasdict(g, a) && delete!(g.aprops, a)
    _hasdict(g, reverse(a)) && delete!(g.aprops, reverse(a))
end


text:
using Thermo
using Test

@testset "Thermo.jl" begin
    # Write your tests here.
end


text:
using SafeTestsets, Test

@safetestset "SafeTestsets Tests" begin
    @safetestset "Tests" begin
        include("tests.jl")
        @isdefined(a) == true
    end

    @test @isdefined(a) == false

    @safetestset MyTests = "Tests" begin
        include("tests.jl")
        @isdefined(a) == true
    end

    @test @isdefined(a) == false
    @test @isdefined(MyTests) == true
end


text:
module TestExtraSimulateKeywordArguments

using TinyModia
using DifferentialEquations
@usingModiaPlot

# Register extra simulate! keyword arguments
registerExtraSimulateKeywordArguments([:dummy1, :dummy2, :dummy3])

function inputSignal(instantiatedModel, t)
    if isInitial(instantiatedModel)
        dict = get_extraSimulateKeywordArgumentsDict(instantiatedModel)
        dummy1 = get(dict, :dummy1, 1.0)
        dummy2 = get(dict, :dummy2, true)
        dummy3 = get(dict, :dummy3, "this is dummy3")
        println("\nExtra simulate! keyword arguments: dummy1 = $dummy1, dummy2 = $dummy2, dummy3 = \"$dummy3\"\n")
    end

    sin(t)
end

FirstOrder = Model(
    T = 0.2,
    x = Var(init=0.3),
    equations = :[u = inputSignal(instantiatedModel, time/u"s"),
                  T * der(x) + x = u,
                  y = 2*x]
)

firstOrder = @instantiateModel(FirstOrder, logCode=false)

simulate!(firstOrder, Tsit5(), stopTime = 10, log=false, dummy1 = 2.0, dummy3 = "value changed",
          requiredFinalStates = [-0.3617373025974107])

plot(firstOrder, ["u", "y"])

end


text:
module ElsaPlots

using Elsa, Makie, GeometryTypes, LinearAlgebra, SparseArrays
import AbstractPlotting: default_theme, Plot, plot!, plot, to_value
using Base.Iterators: take, cycle
using Statistics: mean
using Elsa: vertices, Hamiltonian, Lattice, Bandstructure, padright, siterange

export plot, scale!

include("tools.jl")
include("plotlattice.jl")
include("plotbands.jl")

end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_m7_opt_rho6_4(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_m7_opt_rho6_4!(A_copy)
end

@inline function exp_mono_m7_opt_rho6_4!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 T2k10
    # Computing Bb2 = x*I+x*A
    coeff1=0.1976456832538621
    coeff2=0.45972372405481887
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.1976456832538621
    coeff2=0.45972372405481887
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.039937522468321524
    coeff2=0.1584087562239318
    coeff3=0.9306183490351906
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.3035616192483739
    coeff2=0.4669306299429674
    coeff3=0.04404867971366269
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.012066871713697488
    coeff2=-0.00802148596442477
    coeff3=0.3230204538240676
    coeff4=0.9638185168573307
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.011627318576143263
    coeff2=0.49862155431380795
    coeff3=0.07051226319340288
    coeff4=0.002052240135641998
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0014318669453225418
    coeff2=-0.0029780792038992215
    coeff3=0.004322688752720727
    coeff4=0.013673005502563394
    coeff5=0.8216219421030205
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.8270926502847625
    coeff2=0.32439972472942114
    coeff3=0.058368756227684025
    coeff4=0.0031897875289976277
    coeff5=3.234759791150232e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.00011888017236476145
    coeff2=-0.00010785272062853496
    coeff3=0.004179296166690874
    coeff4=-0.13039438638797307
    coeff5=0.8475672106088131
    coeff6=0.8145974867536132
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.10691367390390938
    coeff2=0.4936148415929095
    coeff3=0.10652748687630387
    coeff4=0.008329230127484285
    coeff5=0.00021583757875408033
    coeff6=-1.0968604718569746e-7
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-4.500140737568134e-6
    coeff2=1.9301987005076668e-6
    coeff3=4.30074986368367e-5
    coeff4=-0.0061294802470251075
    coeff5=0.0688513777099278
    coeff6=-0.12038470670165373
    coeff7=0.9948034699390058
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-9.34103418432391e-8
    coeff2=0.49812886964499026
    coeff3=0.11904681716211153
    coeff4=0.010576806709956765
    coeff5=0.0003949349245182783
    coeff6=8.856309878298438e-6
    coeff7=7.553245375703866e-9
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9839939777638327
    coeff2=0.40889409510333663
    coeff3=0.41640804341987003
    coeff4=0.15933029209747282
    coeff5=0.03184541272812379
    coeff6=0.007635550247315368
    coeff7=0.0008615033049494069
    coeff8=3.134369854398072e-5
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9839939777638327
    coeff2=0.40889409510333663
    coeff3=0.41640804341987003
    coeff4=0.15933029209747282
    coeff5=0.03184541272812379
    coeff6=0.007635550247315368
    coeff7=0.0008615033049494069
    coeff8=3.134369854398072e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.0002507256194653783
    coeff2=2.0145781949958903e-5
    coeff3=-5.521866470188974e-5
    coeff4=-4.6492324041568414e-5
    coeff5=-1.364785597087004e-5
    coeff6=-4.045094239671855e-6
    coeff7=-5.508925657804225e-7
    coeff8=-2.3628259481592187e-8
    coeff9=0.9995948986117282
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using RDatasets, DataFrames, InteractiveFixedEffectModels, Base.Test
df = dataset("plm", "Cigar")
df[:pState] =  pool(df[:State])
df[:pYear] =  pool(df[:Year])

for method in [:gauss_seidel,  :dogleg, :levenberg_marquardt]
	precision = 1e-1
	println(method)
	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 1), :Sales, df , method =  method, maxiter = 10_000) ;
	@test_approx_eq_eps abs(result.augmentdf[:factors1][1]) 0.18662770198472406 precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings1][1]) 587.2272 precision
	@test_approx_eq_eps result.augmentdf[:residuals][1] -15.6928 precision

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), :Sales,  df, method =  method, maxiter = 10_000) ;
	@test_approx_eq_eps abs(result.augmentdf[:factors1][1])  0.18662770198472406 precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings1][1])  587.227 precision
	@test_approx_eq_eps result.augmentdf[:residuals][1] 2.16611 precision

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 1), Sales ~ 1 |> pState, df, method = method, maxiter = 10_000) ;
	@test_approx_eq_eps abs(result.augmentdf[:factors1][1])   0.17636 precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings1][1])  20.176432452716522 precision
	@test_approx_eq_eps result.augmentdf[:residuals][1]  -10.0181 precision
	@test_approx_eq_eps result.augmentdf[:pState][1]  107.4766 precision

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ 1 |> pState, df, method =  method, maxiter = 10_000) ;
	@test_approx_eq_eps abs(result.augmentdf[:factors2][1]) 0.244 precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings2][1]) 49.7943 precision
	@test_approx_eq_eps result.augmentdf[:residuals][1]  2.165319 precision
	@test_approx_eq_eps result.augmentdf[:pState][1]  107.47666 precision

	result = fit(InteractiveFixedEffectModel(:pState, :pYear, 2), Sales ~ 1 |> pState, df, method =  method, subset = (df[:State] .<= 30), maxiter = 10_000) ;
	@test size(result.augmentdf, 1) == size(df, 1)
	@test_approx_eq_eps abs(result.augmentdf[:factors1][1])  0.20215 precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings1][1]) 29.546 precision
	@test_approx_eq_eps abs(result.augmentdf[:factors2][1])  0.250768  precision
	@test_approx_eq_eps abs(result.augmentdf[:loadings2][1])   43.578 precision
	@test_approx_eq_eps result.augmentdf[:residuals][1]  2.9448  precision
	@test_approx_eq_eps result.augmentdf[:pState][1] 107.4766 precision

end


text:
function main()
	N = 10_001
	i = 1
	p = 2
	while i <= N
		if isprime(p)
			println(p)
			i += 1
		end
		p += 1
	end
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
import FractionalTimeDG
using PyPlot

rmax = 5
M = 201
τ = range(-1, 1, length=M)
Ψ = zeros(rmax, M)
FractionalTimeDG.legendre_polys!(Ψ, τ)

figure(1)
fmt = [":", "--", "-.", "-"]
for r = 1:rmax-1
    plot(τ, Ψ[r+1,:]-Ψ[r,:], fmt[r])
end
grid(true)
xlabel(L"$\tau$")
s = LaTeXString[latexstring("\$r=", r, "\$") for r = 1:rmax-1]
legend(s)
savefig("fig1.eps")


text:
using CSDP
using MathOptInterface
const MOI = MathOptInterface
using SemidefiniteOptInterface
const SOI = SemidefiniteOptInterface
m = CSDP.CSDPSolverInstance()
SOI.initinstance!(m, [-2], 0)
SOI.setobjectivecoefficient!(m, -1., 1, 1, 1)
SOI.setobjectivecoefficient!(m, -1., 1, 2, 2)
MOI.optimize!(m)
@show m.status


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

name = "GtkSourceView"
version = v"4.4.0"

# Collection of sources required to build GtkSourceView
sources = [
    "https://download.gnome.org/sources/gtksourceview/4.4/gtksourceview-$(version).tar.xz" =>
    "9ddb914aef70a29a66acd93b4f762d5681202e44094d2d6370e51c9e389e689a",
    "./bundled",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/gtksourceview-*/

# We need to run native `xmllint` and `glib-compile-resources`
apk add libxml2-utils glib-dev

# Don't build broken tests
atomic_patch -p1 ../patches/meson_build_no_tests.patch

mkdir build && cd build
meson .. \
    -Dgir=false \
    --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libgtksourceview-4", :libgtksourceview),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Glib_jll",
    "GTK3_jll",
    "FriBidi_jll",
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using DashBootstrapComponents, DashCoreComponents, DashHtmlComponents

pagination = html_div([
    html_div("Select a page", id = "pagination-contents"),
    dbc_pagination(id = "pagination", max_value = 10),
    html_div("Or set the page dynamically using the slider below"),
    dcc_slider(
        id = "page-change",
        min = 1,
        max = 10,
        step = 1,
        value = 1,
        marks = Dict(
            1 => "1",
            2 => "2",
            3 => "3",
            4 => "4",
            5 => "5",
            6 => "6",
            7 => "7",
            8 => "8",
            9 => "9",
            10 => "10",
        ),
    ),
])

callback!(
    app,
    Output("pagination-contents", "children"),
    Input("pagination", "active_page"),
) do page
    return page > 0 ? "Page selected: $page" : "Select a page"
end;

callback!(app, Output("pagination", "active_page"), Input("page-change", "value")) do value
    return value
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
# Provide implementation of default layers

Provided layers
- DenseLayer
- DenseNoBiasLayer
- VectorFunctionLayer
"""

#using Random, Zygote
#import ..Utils
#import Base.size

# ------------------------------------------------------------------------------
# DenseLayer layer

"""
   DenseLayer

Representation of a layer in the network

# Fields:
* `w`:  Weigths matrix with respect to the input from previous layer or data (n x n pr. layer)
* `wb`: Biases (n)
* `f`:  Activation function
* `df`: Derivative of the activation function
"""
mutable struct DenseLayer <: AbstractLayer
     w::Array{Float64,2}
     wb::Array{Float64,1}
     f::Function
     df::Union{Function,Nothing}
     """
        DenseLayer(f,n,nₗ;w,wb,df)

     Instantiate a new DenseLayer

     # Positional arguments:
     * `nₗ`: Number of nodes of the previous layer
     * `n`:  Number of nodes
     # Keyword arguments:
     * `w`:   Initial weigths with respect to input [default: Xavier initialisation, dims = (n,nₗ)]
     * `wb`:  Initial weigths with respect to bias [default: Xavier initialisation, dims = (n)]
     * `f`:   Activation function [def: `identity`]
     * `df`:  Derivative of the activation function [default: `nothing` (i.e. use AD)]
     * `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

     # Notes:
     - Xavier initialization = `rand(Uniform(-sqrt(6)/sqrt(nₗ+n),sqrt(6)/sqrt(nₗ+n))`

     """
     function DenseLayer(nₗ,n;rng = Random.GLOBAL_RNG,w=rand(rng, Uniform(-sqrt(6)/sqrt(nₗ+n),sqrt(6)/sqrt(nₗ+n)),n,nₗ),wb=rand(rng, Uniform(-sqrt(6)/sqrt(nₗ+n),sqrt(6)/sqrt(nₗ+n)),n),f=identity,df=nothing)
         # To be sure w is a matrix and wb a column vector..
         w  = reshape(w,n,nₗ)
         wb = reshape(wb,n)
         return new(w,wb,f,df)
     end
end

function _zComp(layer::DenseLayer,x)
    w  = layer.w
    wb = layer.wb
    z  = zeros(eltype(x),size(w,1))
    @inbounds for n in axes(w,1)
        zn = zero(eltype(x))
        @simd for nl in axes(x,1)
            zn += w[n,nl] * x[nl]
        end
        zn   += wb[n]
        z[n]  = zn
    end
    return z
end

function forward(layer::DenseLayer,x)
  z =  _zComp(layer,x) #@avx layer.w * x + layer.wb #_zComp(layer,x) #layer.w * x + layer.wb # _zComp(layer,x) #   layer.w * x + layer.wb # testd @avx
  return layer.f.(z)
end

function backward(layer::DenseLayer,x,nextGradient)
   z = _zComp(layer,x) #@avx layer.w * x + layer.wb #_zComp(layer,x) # layer.w * x + layer.wb # _zComp(layer,x) # @avx layer.w * x + layer.wb               # tested @avx
   if layer.df != nothing
       dϵ_dz =  layer.df.(z) .* nextGradient # tested @avx
    else
       dϵ_dz = layer.f'.(z) .* nextGradient # using AD
    end
   dϵ_dI =  layer.w' * dϵ_dz # @avx
end

function getParams(layer::DenseLayer)
  return Learnable((layer.w,layer.wb))
end

function getGradient(layer::DenseLayer,x,nextGradient)
   z      =  _zComp(layer,x) #@avx layer.w * x + layer.wb #  _zComp(layer,x) #layer.w * x + layer.wb # @avx
   if layer.df != nothing
       dϵ_dz = layer.df.(z) .* nextGradient
    else
       dϵ_dz = layer.f'.(z) .* nextGradient # using AD
    end
   dϵ_dw  = dϵ_dz * x' # @avx
   dϵ_dwb = dϵ_dz
   return Learnable((dϵ_dw,dϵ_dwb))
end

function setParams!(layer::DenseLayer,w)
   layer.w  = w.data[1]
   layer.wb = w.data[2]
end
function size(layer::DenseLayer)
    return size(layer.w')
end

# ------------------------------------------------------------------------------
# DenseNoBiasLayer layer

"""
   DenseNoBiasLayer

Representation of a layer without bias in the network

# Fields:
* `w`:  Weigths matrix with respect to the input from previous layer or data (n x n pr. layer)
* `f`:  Activation function
* `df`: Derivative of the activation function
"""
mutable struct DenseNoBiasLayer <: AbstractLayer
     w::Array{Float64,2}
     f::Function
     df::Union{Function,Nothing}
     """
        DenseNoBiasLayer(f,nₗ,n;w,df)

     Instantiate a new DenseNoBiasLayer

        DenseNoBiasLayer(f,nₗ,n;w,df)

     Instantiate a new DenseNoBiasLayer

     # Positional arguments:
     * `nₗ`:  Number of nodes of the previous layer
     * `n`:   Number of nodes
     # Keyword arguments:
     * `w`:   Initial weigths with respect to input [default: Xavier initialisation, dims = (nₗ,n)]
     * `f`:   Activation function [def: `identity`]
     * `df`:  Derivative of the activation function [def: `nothing` (i.e. use AD)]
     * `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]
     # Notes:
     - Xavier initialization = `rand(Uniform(-sqrt(6)/sqrt(nₗ+n),sqrt(6)/sqrt(nₗ,n))`
     """
     function DenseNoBiasLayer(nₗ,n;rng = Random.GLOBAL_RNG,w=rand(rng,Uniform(-sqrt(6)/sqrt(nₗ+n),sqrt(6)/sqrt(nₗ+n)),n,nₗ),f=identity,df=nothing)
         # To be sure w is a matrix and wb a column vector..
         w  = reshape(w,n,nₗ)
         return new(w,f,df)
     end
end

function _zComp(layer::DenseNoBiasLayer,x)
    w  = layer.w
    z  = zeros(eltype(x),size(w,1))
    @inbounds for n in axes(w,1)
        zn = zero(eltype(x))
        @simd for nl in axes(x,1)
            zn += w[n,nl] * x[nl]
        end
        z[n]  = zn
    end
    return z
end

function forward(layer::DenseNoBiasLayer,x)
  z = _zComp(layer,x)
  return layer.f.(z)
end

function backward(layer::DenseNoBiasLayer,x,nextGradient)
   z = _zComp(layer,x)
   if layer.df != nothing
       dϵ_dz = layer.df.(z) .* nextGradient
    else
       dϵ_dz = layer.f'.(z) .* nextGradient # using AD
    end
   dϵ_dI = layer.w' * dϵ_dz
end

function getParams(layer::DenseNoBiasLayer)
  return Learnable((layer.w,))
end

function getGradient(layer::DenseNoBiasLayer,x,nextGradient)
   z      = _zComp(layer,x)
   if layer.df != nothing
       dϵ_dz = layer.df.(z) .* nextGradient
    else
       dϵ_dz = layer.f'.(z) .* nextGradient # using AD
    end
   dϵ_dw  =  dϵ_dz * x'
   return Learnable((dϵ_dw,))
end

function setParams!(layer::DenseNoBiasLayer,w)
   layer.w = w.data[1]
end

function size(layer::DenseNoBiasLayer)
    return size(layer.w')
end

# ------------------------------------------------------------------------------
# VectorFunction layer

"""
   VectorFunctionLayer

Representation of a (weightless) VectorFunction layer in the network. Vector
function layer expects a vector activation function, i.e. a function taking the
whole output of the previous layer an input rather than working on a single node
as "normal" activation functions would do.
Useful for example with the SoftMax function in classification or with the `pool1D`
function to implement a "pool" layer in 1 dimensions.
As it is weightless, it doesn't apply any transformation to the output coming
from the previous layer. It means that the number of nodes must be set to the same
as in the previous layer (and if you are using this for classification, to the
number of classes, i.e. the _previous_ layer must be set equal to the number of
classes in the predictions).

# Fields:
* `nₗ`: Number of nodes in input (i.e. length of previous layer)
* `n`:  Number of nodes in output (automatically inferred in the constructor)
* `f`:  Activation function (vector)
* `df`: Derivative of the (vector) activation function

# Notes:
* The output `size` of this layer is given by the size of the output function,
that not necessarily is the same as the previous layers.
"""
mutable struct VectorFunctionLayer <: AbstractLayer
     nₗ::Int64
     n::Int64
     f::Function
     df::Union{Function,Nothing}
     """
        VectorFunctionLayer(f,n;df)

     Instantiate a new VectorFunctionLayer

     # Positional arguments:
     * `nₗ`: Number of nodes (must be same as in the previous layer)
     # Keyword arguments:
     * `f`:  Activation function [def: `softMax`]
     * `df`: Derivative of the activation function [default: `nothing` (i.e. use AD)]
     * `dummyDataToTestOutputSize`: Dummy data to test the output size [def: `ones(nₗ)`]

     # Notes:
     - If the derivative is provided, it should return the gradient as a (n,n) matrix (i.e. the Jacobian)
     - To avoid recomputing the activation function just to determine its output size,
       we compute the output size once here in the layer constructor by calling the
       activation function with `dummyDataToTestOutputSize`. Feel free to change
       it if it doesn't match with the activation function you are setting
     """
     function VectorFunctionLayer(nₗ;f=softMax,df=nothing,dummyDataToTestOutputSize=ones(nₗ))
         n = length(f(dummyDataToTestOutputSize))
         return new(nₗ,n,f,df)
     end
end

function forward(layer::VectorFunctionLayer,x)
  return layer.f(x)
end

function backward(layer::VectorFunctionLayer,x,nextGradient)
   if layer.df != nothing
       dϵ_dI = layer.df(x)' * nextGradient
    else  # using AD
      j = autoJacobian(layer.f,x; nY=layer.n)
      dϵ_dI = j' * nextGradient
    end
   return dϵ_dI
end

function getParams(layer::VectorFunctionLayer)
  return Learnable(())
end

function getGradient(layer::VectorFunctionLayer,x,nextGradient)
   return Learnable(()) # parameterless layer
end

function setParams!(layer::VectorFunctionLayer,w)
   return nothing
end
function size(layer::VectorFunctionLayer)
    # Output size for the VectorFunctionLayer is given by its activation function
    # We test its length with dummy values
    return (layer.nₗ,layer.n)
end


text:
# gkl.jl
"""
    mutable struct GKLFactorization{TU,TV,S<:Real}

Structure to store a Golub-Kahan-Lanczos (GKL) bidiagonal factorization of a linear map `A`
of the form

```julia
A * V = U * B + r * b'
A' * U = V * B'
```

For a given GKL factorization `fact` of length `k = length(fact)`, the two bases `U` and `V`
are obtained via [`basis(fact, :U)`](@ref basis) and `basis(fact, :V)`. Here, `U` and `V`
are instances of [`OrthonormalBasis{T}`](@ref Basis), with also
`length(U) == length(V) == k` and where `T` denotes the type of vector like objects used in
the problem. The Rayleigh quotient `B` is obtained as [`rayleighquotient(fact)`](@ref) and
is of type `Bidiagonal{S<:Number}` with `size(B) == (k,k)`. The residual `r` is
obtained as [`residual(fact)`](@ref) and is of type `T`. One can also query
[`normres(fact)`](@ref) to obtain `norm(r)`, the norm of the residual. The vector `b` has no
dedicated name but can be obtained via [`rayleighextension(fact)`](@ref). It takes the
default value ``e_k``, i.e. the unit vector of all zeros and a one in the last entry, which
is represented using [`SimpleBasisVector`](@ref).

A GKL factorization `fact` can be destructured as `U, V, B, r, nr, b = fact` with
`nr = norm(r)`.

`GKLFactorization` is mutable because it can [`expand!`](@ref) or [`shrink!`](@ref).
See also [`GKLIterator`](@ref) for an iterator that constructs a progressively expanding
GKL factorizations of a given linear map and a starting vector `u₀`.
"""
mutable struct GKLFactorization{TU,TV,S<:Real}
    k::Int # current Krylov dimension
    U::OrthonormalBasis{TU} # basis of length k
    V::OrthonormalBasis{TV} # basis of length k
    αs::Vector{S}
    βs::Vector{S}
    r::TU
end

Base.length(F::GKLFactorization) = F.k
Base.sizehint!(F::GKLFactorization, n) = begin
    sizehint!(F.U, n)
    sizehint!(F.V, n)
    sizehint!(F.αs, n)
    sizehint!(F.βs, n)
    return F
end
Base.eltype(F::GKLFactorization) = eltype(typeof(F))
Base.eltype(::Type{<:GKLFactorization{<:Any,<:Any,S}}) where {S} = S

# iteration for destructuring into components
Base.iterate(F::GKLFactorization) = (basis(F, :U), Val(:V))
Base.iterate(F::GKLFactorization, ::Val{:V}) = (basis(F, :V), Val(:rayleighquotient))
Base.iterate(F::GKLFactorization, ::Val{:rayleighquotient}) =
    (rayleighquotient(F), Val(:residual))
Base.iterate(F::GKLFactorization, ::Val{:residual}) = (residual(F), Val(:normres))
Base.iterate(F::GKLFactorization, ::Val{:normres}) =
    (normres(F), Val(:rayleighextension))
Base.iterate(F::GKLFactorization, ::Val{:rayleighextension}) =
    (rayleighextension(F), Val(:done))
Base.iterate(F::GKLFactorization, ::Val{:done}) = nothing

"""
        basis(fact::GKLFactorization, which::Symbol)

Return the list of basis vectors of a [`GKLFactorization`](@ref), where `which` should take
the value `:U` or `:V` and indicates which set of basis vectors (in the domain or in the
codomain of the corresponding linear map) should be returned. The return type is an
`OrthonormalBasis{T}`, where `T` represents the type of the vectors used by the problem.
"""
function basis(F::GKLFactorization, which::Symbol)
    length(F.U) == F.k || error("Not keeping vectors during GKL bidiagonalization")
    which == :U || which == :V || error("invalid flag for specifying basis")
    return which == :U ? F.U : F.V
end
rayleighquotient(F::GKLFactorization) =
    Bidiagonal(view(F.αs, 1:F.k), view(F.βs, 1:(F.k-1)), :L)
residual(F::GKLFactorization) = F.r
@inbounds normres(F::GKLFactorization) = F.βs[F.k]
rayleighextension(F::GKLFactorization) = SimpleBasisVector(F.k, F.k)

# GKL iteration for constructing the orthonormal basis of a Krylov subspace.
"""
    struct GKLIterator{F,TU,O<:Orthogonalizer}
    GKLIterator(f, u₀, [orth::Orthogonalizer = KrylovDefaults.orth, keepvecs::Bool = true])

Iterator that takes a general linear map `f::F` and an initial vector `u₀::TU` and generates
an expanding `GKLFactorization` thereof. In particular, `GKLIterator` implements the
[Golub-Kahan-Lanczos bidiagonalization procedure](http://www.netlib.org/utk/people/JackDongarra/etemplates/node198.html).
Note, however, that this implementation starts from a vector `u₀` in the codomain of the
linear map `f`, which will end up (after normalisation) as the first column of `U`.

The argument `f` can be a matrix, a tuple of two functions where the first represents the
normal action and the second the adjoint action, or a function accepting two arguments,
where the first argument is the vector to which the linear map needs to be applied, and the
second argument is either `Val(false)` for the normal action and `Val(true)` for the adjoint
action. Note that the flag is thus a `Val` type to allow for type stability in cases where
the vectors in the domain and the codomain of the linear map have a different type.

The optional argument `orth` specifies which [`Orthogonalizer`](@ref) to be used. The
default value in [`KrylovDefaults`](@ref) is to use [`ModifiedGramSchmidtIR`](@ref), which
possibly uses reorthogonalization steps.

When iterating over an instance of `GKLIterator`, the values being generated are
instances `fact` of [`GKLFactorization`](@ref), which can be immediately destructured into a
[`basis(fact, :U)`](@ref), [`basis(fact, :V)`](@ref), [`rayleighquotient`](@ref),
[`residual`](@ref), [`normres`](@ref) and [`rayleighextension`](@ref), for example as

```julia
for (U, V, B, r, nr, b) in GKLIterator(f, u₀)
    # do something
    nr < tol && break # a typical stopping criterion
end
```

Since the iterator does not know the dimension of the underlying vector space of
objects of type `T`, it keeps expanding the Krylov subspace until the residual norm `nr`
falls below machine precision `eps(typeof(nr))`.

The internal state of `GKLIterator` is the same as the return value, i.e. the corresponding
`GKLFactorization`. However, as Julia's Base iteration interface (using `Base.iterate`)
requires that the state is not mutated, a `deepcopy` is produced upon every next iteration
step.

Instead, you can also mutate the `GKLFactorization` in place, using the following
interface, e.g. for the same example above

```julia
iterator = GKLIterator(f, u₀)
factorization = initialize(iterator)
while normres(factorization) > tol
    expand!(iterator, factorization)
    U, V, B, r, nr, b = factorization
    # do something
end
```

Here, [`initialize(::GKLIterator)`](@ref) produces the first GKL factorization of length 1,
and `expand!(::GKLIterator, ::GKLFactorization)`(@ref) expands the factorization in place.
See also [`initialize!(::GKLIterator, ::GKLFactorization)`](@ref) to initialize in an
already existing factorization (most information will be discarded) and
[`shrink!(::GKLIterator, k)`](@ref) to shrink an existing factorization down to length `k`.
"""
struct GKLIterator{F,TU,O<:Orthogonalizer}
    operator::F
    u₀::TU
    orth::O
    keepvecs::Bool
    function GKLIterator{F,TU,O}(
        operator::F,
        u₀::TU,
        orth::O,
        keepvecs::Bool
    ) where {F,TU,O<:Orthogonalizer}
        if !keepvecs && isa(orth, Reorthogonalizer)
            error("Cannot use reorthogonalization without keeping all Krylov vectors")
        end
        return new{F,TU,O}(operator, u₀, orth, keepvecs)
    end
end
GKLIterator(
    operator::F,
    u₀::TU,
    orth::O = KrylovDefaults.orth,
    keepvecs::Bool = true
) where {F,TU,O<:Orthogonalizer} = GKLIterator{F,TU,O}(operator, u₀, orth, keepvecs)

Base.IteratorSize(::Type{<:GKLIterator}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:GKLIterator}) = Base.EltypeUnknown()

function Base.iterate(iter::GKLIterator)
    state = initialize(iter)
    return state, state
end
function Base.iterate(iter::GKLIterator, state::GKLFactorization)
    nr = normres(state)
    if nr < eps(typeof(nr))
        return nothing
    else
        state = expand!(iter, deepcopy(state))
        return state, state
    end
end

function initialize(iter::GKLIterator; verbosity::Int = 0)
    # initialize without using eltype
    u₀ = iter.u₀
    β₀ = norm(u₀)
    iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero"))
    v₀ = apply_adjoint(iter.operator, u₀)
    α = norm(v₀) / β₀
    Av₀ = apply_normal(iter.operator, v₀) # apply operator
    α² = dot(u₀, Av₀) / β₀^2
    α² ≈ α * α || throw(ArgumentError("operator and its adjoint are not compatible"))
    T = typeof(α²)
    # these lines determines the type that we will henceforth use
    u = mul!(zero(T)*Av₀, u₀, 1 / β₀) # (one(T) / β₀) * u₀
    v = (one(T) / (α * β₀)) * v₀
    if typeof(Av₀) == typeof(u)
        r = rmul!(Av₀, 1 / (α * β₀))
    else
        r = mul!(similar(u), Av₀, 1 / (α * β₀))
    end
    r = axpy!(-α, u, r)
    β = norm(r)

    U = OrthonormalBasis([u])
    V = OrthonormalBasis([v])
    S = real(T)
    αs = S[α]
    βs = S[β]
    if verbosity > 0
        @info "GKL iteration step 1: normres = $β"
    end

    return GKLFactorization(1, U, V, αs, βs, r)
end
function initialize!(iter::GKLIterator, state::GKLFactorization; verbosity::Int = 0)
    U = state.U
    while length(U) > 1
        pop!(U)
    end
    V = empty!(state.V)
    αs = empty!(state.αs)
    βs = empty!(state.βs)

    u = mul!(U[1], iter.u₀, 1 / norm(iter.u₀))
    v = apply_adjoint(iter.operator, u)
    α = norm(v)
    rmul!(v, 1 / α)
    r = apply_normal(iter.operator, v)
    r = axpy!(-α, u, r)
    β = norm(r)

    state.k = 1
    push!(V, v)
    push!(αs, α)
    push!(βs, β)
    state.r = r
    if verbosity > 0
        @info "GKL iteration step 1: normres = $β"
    end

    return state
end
function expand!(iter::GKLIterator, state::GKLFactorization; verbosity::Int = 0)
    βold = normres(state)
    U = state.U
    V = state.V
    r = state.r
    U = push!(U, rmul!(r, 1 / βold))
    v, r, α, β = gklrecurrence(iter.operator, U, V, βold, iter.orth)

    push!(V, v)
    push!(state.αs, α)
    push!(state.βs, β)

    #!iter.keepvecs && popfirst!(state.V) # remove oldest V if not keepvecs

    state.k += 1
    state.r = r
    if verbosity > 0
        @info "GKL iteration step $(state.k): normres = $β"
    end

    return state
end
function shrink!(state::GKLFactorization, k)
    length(state) == length(state.V) ||
        error("we cannot shrink GKLFactorization without keeping vectors")
    length(state) <= k && return state
    U = state.U
    V = state.V
    while length(V) > k + 1
        pop!(U)
        pop!(V)
    end
    pop!(V)
    r = pop!(U)
    resize!(state.αs, k)
    resize!(state.βs, k)
    state.k = k
    state.r = rmul!(r, normres(state))
    return state
end

# Golub-Kahan-Lanczos recurrence relation
function gklrecurrence(
    operator,
    U::OrthonormalBasis,
    V::OrthonormalBasis,
    β,
    orth::Union{ClassicalGramSchmidt,ModifiedGramSchmidt}
)
    u = U[end]
    v = apply_adjoint(operator, u)
    v = axpy!(-β, V[end], v)
    α = norm(v)
    rmul!(v, inv(α))

    r = apply_normal(operator, v)
    r = axpy!(-α, u, r)
    β = norm(r)
    return v, r, α, β
end
function gklrecurrence(
    operator,
    U::OrthonormalBasis,
    V::OrthonormalBasis,
    β,
    orth::ClassicalGramSchmidt2
)
    u = U[end]
    v = apply_adjoint(operator, u)
    v = axpy!(-β, V[end], v) # not necessary if we definitely reorthogonalize next step and previous step
    # v, = orthogonalize!(v, V, ClassicalGramSchmidt())
    α = norm(v)
    rmul!(v, inv(α))

    r = apply_normal(operator, v)
    r = axpy!(-α, u, r)
    r, = orthogonalize!(r, U, ClassicalGramSchmidt())
    β = norm(r)
    return v, r, α, β
end
function gklrecurrence(
    operator,
    U::OrthonormalBasis,
    V::OrthonormalBasis,
    β,
    orth::ModifiedGramSchmidt2
)
    u = U[end]
    v = apply_adjoint(operator, u)
    v = axpy!(-β, V[end], v)
    # for q in V # not necessary if we definitely reorthogonalize next step and previous step
    #     v, = orthogonalize!(v, q, ModifiedGramSchmidt())
    # end
    α = norm(v)
    rmul!(v, inv(α))

    r = apply_normal(operator, v)
    r = axpy!(-α, u, r)
    for q in U
        r, = orthogonalize!(r, q, ModifiedGramSchmidt())
    end
    β = norm(r)
    return v, r, α, β
end
function gklrecurrence(
    operator,
    U::OrthonormalBasis,
    V::OrthonormalBasis,
    β,
    orth::ClassicalGramSchmidtIR
)
    u = U[end]
    v = apply_adjoint(operator, u)
    v = axpy!(-β, V[end], v)
    α = norm(v)
    nold = sqrt(abs2(α) + abs2(β))
    while α < orth.η * nold
        nold = α
        v, = orthogonalize!(v, V, ClassicalGramSchmidt())
        α = norm(v)
    end
    rmul!(v, inv(α))

    r = apply_normal(operator, v)
    r = axpy!(-α, u, r)
    β = norm(r)
    nold = sqrt(abs2(α) + abs2(β))
    while eps(one(β)) < β < orth.η * nold
        nold = β
        r, = orthogonalize!(r, U, ClassicalGramSchmidt())
        β = norm(r)
    end

    return v, r, α, β
end
function gklrecurrence(
    operator,
    U::OrthonormalBasis,
    V::OrthonormalBasis,
    β,
    orth::ModifiedGramSchmidtIR
)
    u = U[end]
    v = apply_adjoint(operator, u)
    v = axpy!(-β, V[end], v)
    α = norm(v)
    nold = sqrt(abs2(α) + abs2(β))
    while eps(one(α)) < α < orth.η * nold
        nold = α
        for q in V
            v, = orthogonalize!(v, q, ModifiedGramSchmidt())
        end
        α = norm(v)
    end
    rmul!(v, inv(α))

    r = apply_normal(operator, v)
    r = axpy!(-α, u, r)
    β = norm(r)
    nold = sqrt(abs2(α) + abs2(β))
    while eps(one(β)) < β < orth.η * nold
        nold = β
        for q in U
            r, = orthogonalize!(r, q, ModifiedGramSchmidt())
        end
        β = norm(r)
    end

    return v, r, α, β
end


text:
function main()
	;
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# Simple ridge regression based upon the sindy-mpc
# repository, see https://arxiv.org/abs/1711.05501
# and https://github.com/eurika-kaiser/SINDY-MPC/blob/master/LICENSE

"""
$(TYPEDEF)
`STLQS` is taken from the [original paper on SINDY](https://www.pnas.org/content/113/15/3932) and implements a
sequentially thresholded least squares iteration. `λ` is the threshold of the iteration.
It is based upon [this matlab implementation](https://github.com/eurika-kaiser/SINDY-MPC/utils/sparsifyDynamics.m).
It solves the following problem
```math
\\argmin_{x} \\frac{1}{2} \\| Ax-b\\|_2 + \\lambda \\|x\\|_2
```

# Fields
$(FIELDS)

# Example
```julia
opt = STLQS()
opt = STLQS(1e-1)
opt = STLQS(Float32[1e-2; 1e-1])
```
## Note
This was formally `STRRidge` and has been renamed.
"""
mutable struct STLSQ{T} <: AbstractOptimizer{T}
    """Sparsity threshold"""
    λ::T

    function STLSQ(threshold::T = 1e-1) where T
        @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite"

        return new{typeof(threshold)}(threshold)
    end

end

Base.summary(::STLSQ) = "STLSQ"

function (opt::STLSQ{T})(X, A, Y, λ::U = first(opt.λ);
    maxiter = maximum(size(A)), abstol::U = eps(eltype(T)),
    progress = nothing, kwargs...) where {T,U}

    smallinds = abs.(X) .<= λ
    biginds = @. ! smallinds[:, 1]

    x_i = similar(X)
    x_i .= X

    xzero = zero(eltype(X))
    obj = xzero
    sparsity = xzero
    conv_measure = xzero

    iters = 0
    converged = false

    _progress = isa(progress, Progress)
    initial_prog = _progress ? progress.counter : 0

    @views while (iters < maxiter) && !converged
        iters += 1

        smallinds .= abs.(X) .<= λ
        X[smallinds] .= xzero

        for j in 1:size(Y, 2)
            @. biginds = ! smallinds[:, j]
            X[biginds, j] .= A[:, biginds] \ Y[:,j]
        end

        conv_measure = norm(x_i .- X, 2)

        if _progress
            obj = norm(Y - A*X, 2)
            sparsity = norm(X, 0, λ)

            ProgressMeter.next!(
            progress;
            showvalues = [
                (:Threshold, λ), (:Objective, obj), (:Sparsity, sparsity),
                (:Convergence, conv_measure)
            ]
            )
        end

        if conv_measure < abstol
            converged = true

            _progress ? (progress.counter = initial_prog + maxiter) : nothing

            #    ProgressMeter.update!(
            #    progress,
            #    initial_prog + maxiter
            #    )
            #end

        else
            x_i .= X
        end
    end

    clip_by_threshold!(X, λ; kwargs...)
    return
end


text:
# Themelis, Stella, Patrinos, "Douglas-Rachford splitting and ADMM
# for nonconvex optimization: Accelerated and Newton-type algorithms",
# arXiv:2005.10230, 2020.
#
# https://arxiv.org/abs/2005.10230

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

function drls_default_gamma(f, muf, Lf, alpha, lambda)
    if muf !== nothing && muf > 0
        return 1 / (alpha * muf)
    end
    if ProximalOperators.is_convex(f)
        return alpha / Lf
    else
        return alpha * (2 - lambda) / (2 * Lf)
    end
end

function drls_C(f, muf, Lf, gamma, lambda, beta)
    a = muf === nothing || muf <= 0 ? gamma * Lf : 1 / (gamma * muf)
    m = ProximalOperators.is_convex(f) ? max(a - lambda / 2, 0) : 1
    return (lambda / ((1 + a)^2) * ((2 - lambda) / 2 - a * m))
end

Base.@kwdef struct DRLSIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,Tg,TH}
    f::Tf = Zero()
    g::Tg = Zero()
    x0::Tx
    alpha::R = real(eltype(x0))(0.95)
    beta::R = real(eltype(x0))(0.5)
    lambda::R = real(eltype(x0))(1)
    muf::Maybe{R} = nothing
    Lf::Maybe{R} = nothing
    gamma::R = drls_default_gamma(f, muf, Lf, alpha, lambda)
    c::R = beta * drls_C(f, muf, Lf, gamma, lambda, beta)
    dre_sign::Int = muf === nothing || muf <= 0 ? 1 : -1
    max_backtracks::Int = 20
    H::TH = LBFGS(x0, 5)
end

Base.IteratorSize(::Type{<:DRLSIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct DRLSState{R,Tx,TH}
    x::Tx
    u::Tx
    v::Tx
    w::Tx
    res::Tx
    res_prev::Tx = zero(x)
    xbar::Tx
    xbar_prev::Tx = zero(x)
    d::Tx = zero(x)
    x_d::Tx = zero(x)
    gamma::R
    f_u::R
    g_v::R
    H::TH
    tau::Maybe{R} = nothing
end

DRE(f_u::Number, g_v::Number, x, u, res, gamma) = f_u + g_v - real(dot(x - u, res)) / gamma + 1 / (2 * gamma) * norm(res)^2

DRE(state::DRLSState) = DRE(state.f_u, state.g_v, state.x, state.u, state.res, state.gamma)

function Base.iterate(iter::DRLSIteration)
    x = copy(iter.x0)
    u, f_u = prox(iter.f, x, iter.gamma)
    w = 2 * u - x
    v, g_v = prox(iter.g, w, iter.gamma)
    res = u - v
    xbar = x - iter.lambda * res
    state = DRLSState(
        x=x, u=u, v=v, w=w, res=res, xbar=xbar, gamma=iter.gamma, f_u=f_u,
        g_v=g_v, H=iter.H,
    )
    return state, state
end

function Base.iterate(iter::DRLSIteration{R}, state::DRLSState) where {R}
    DRE_curr = DRE(state)

    mul!(state.d, iter.H, -state.res)
    state.x_d .= state.x .+ state.d
    state.xbar_prev, state.xbar = state.xbar, state.xbar_prev
    state.res_prev, state.res = state.res, state.res_prev
    state.tau = R(1)
    state.x .= state.x_d

    for k = 1:iter.max_backtracks
        state.x .= state.tau .* state.x_d .+ (1 - state.tau) .* state.xbar_prev

        state.f_u = prox!(state.u, iter.f, state.x, iter.gamma)
        state.w .= 2 .* state.u .- state.x
        state.g_v = prox!(state.v, iter.g, state.w, iter.gamma)
        state.res .= state.u .- state.v

        if k == 1
            update!(iter.H, state.d, state.res - state.res_prev)
        end

        state.xbar .= state.x .- iter.lambda * state.res
        DRE_candidate = DRE(state)

        if iter.dre_sign * DRE_candidate <= iter.dre_sign * DRE_curr - iter.c / iter.gamma * norm(state.res_prev)^2
            break
        end

        state.tau = k < iter.max_backtracks - 1 ? state.tau / 2 : R(0)
    end

    return state, state
end

# Solver

struct DRLS{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::DRLS)(x0; kwargs...)
    stop(state::DRLSState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it,
        state.gamma / state.gamma,
        norm(state.res, Inf),
        (state.tau === nothing ? 0.0 : state.tau)
    )
    iter = DRLSIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.u, state_final.v, num_iters
end

DRLS(; maxit=1_000, tol=1e-8, verbose=false, freq=10, kwargs...) =
    DRLS(maxit, tol, verbose, freq, kwargs)


text:
using ResearchBooks
using Test

@testset "ResearchBooks.jl" begin
    # Write your tests here.
end


text:
"""

Sample the parent label of each event for given parameters.

"""
function sample_parents(events::Array{<:Number, 1}, bg, kappa::Number, kern::Function)
    eventdifferencelist = event_difference_list(events)
    sample_parents(events, bg, kappa, kern, eventdifferencelist)
end

function sample_parents(events::Array{<:Number, 1}, bg, kappa::Number, kern::Function, eventdifferencelist)
    n = length(eventdifferencelist)
    numBackgroundEvents = 0
    childEventsVector = zeros(Int64, n)
    childShiftTime = zeros(Float64, n)
    parents = Array{Int64, 1}(undef, n)

    for (i, eventdiff) in enumerate(eventdifferencelist)
        parent = sample_parent(events[i], bg, kappa, kern, eventdiff)
        parents[i] = parent

        if parent == 0
            numBackgroundEvents += 1
            continue
         end
        childEventsVector[parent] += 1
        childShiftTime[i] = events[i] - events[parent]

    end
    (parents, numBackgroundEvents, childEventsVector,  childShiftTime[childShiftTime .!= 0])
end

function sample_parent(event::Number, bg, kappa::Number, kern::Function, eventdifferences::Array{<:Number, 1})
    probs = event_parent_prob(event, bg, kappa, kern, eventdifferences)
    StatsBase.sample(StatsBase.Weights(probs)) - 1
end


text:
using SafeTestsets

@time begin
    @time @safetestset "Initial" begin include("test.jl") end
end


text:
module ManifoldOps

using DifferentialForms: Forms, bitsign, hodge, invhodge, ⋆
using LinearAlgebra

using ..Funs
using ..Manifolds
using ..Ops

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

# Topological operators

# Boundary

export boundary
function boundary(::Val{Pr}, ::Val{R}, manifold::Manifold{D}) where {R,D}
    @assert 0 < R ≤ D
    # println("[boundary(D=$D, P=Pr, R=$R)]")
    return Op{D,Pr,R - 1,Pr,R}(manifold, get_boundaries(manifold, R).op)
end
function boundary(::Val{Dl}, ::Val{R}, manifold::Manifold{D}) where {R,D}
    @assert 0 ≤ R < D
    # println("[boundary(D=$D, P=Dl, R=$R)]")
    return Op{D,Dl,R + 1,Dl,R}(manifold, get_boundaries(manifold, R + 1).op')
end

boundary(f::Fun{D,P,R}) where {D,P,R} = boundary(Val(P), Val(R), f.manifold) * f

export isboundary
function isboundary(::Val{Pr}, ::Val{R}, manifold::Manifold{D}) where {R,D}
    @assert 0 ≤ R < D
    # println("[isboundary(D=$D, P=Pr, R=$R)]")
    return Op{D,Pr,R,Pr,R}(manifold, map(Bool, get_isboundary(manifold, R)).op)
end

# Derivative (coboundary)

export deriv
function deriv(::Val{P}, ::Val{R}, manifold::Manifold{D}) where {P,R,D}
    P::PrimalDual
    s = P == Pr ? +1 : -1
    @assert 0 ≤ R ≤ D
    @assert 0 ≤ R + s ≤ D
    # println("[deriv(D=$D, P=$P, R=$R)]")
    return adjoint(boundary(Val(P), Val(R + s), manifold))::Op{D,P,R + s,P,R}
end

deriv(f::Fun{D,P,R}) where {D,P,R} = deriv(Val(P), Val(R), f.manifold) * f

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

# Geometric operators

# Hodge dual (circumcentric (diagonal) hodge operator)

export hodge, ⋆
function Forms.hodge(::Val{Pr}, ::Val{R},
                     manifold::Manifold{D,C,S}) where {R,D,C,S}
    @assert 0 ≤ R ≤ D
    # println("[hodge(D=$D, P=Pr, R=$R)]")
    vol = get_volumes(manifold, R)
    dualvol = get_dualvolumes(manifold, R)
    if isempty(vol)
        return zero(Op{D,Dl,R,Pr,R,S}, manifold)
    end
    return Op{D,Dl,R,Pr,R,S}(manifold, Diagonal(dualvol ./ vol))
end

export invhodge
function Forms.invhodge(::Val{Dl}, ::Val{R},
                        manifold::Manifold{D,C,S}) where {P,R,D,C,S}
    R::Int
    @assert 0 ≤ R ≤ D
    # println("[invhodge(D=$D, P=Dl, R=$R)]")
    vol = get_volumes(manifold, R)
    dualvol = get_dualvolumes(manifold, R)
    if isempty(vol)
        return zero(Op{D,Pr,R,Dl,R,S}, manifold)
    end
    return Op{D,Pr,R,Dl,R,S}(manifold, Diagonal(vol ./ dualvol))
end

function Forms.hodge(::Val{Dl}, ::Val{R},
                     manifold::Manifold{D,C,S}) where {R,D,C,S}
    @assert 0 ≤ R ≤ D
    # TODO: sign of metric
    s = bitsign(R * (D - R))
    # s = bitsign(iseven(D) ? 0 : iodd(R))
    return s * invhodge(Val(Dl), Val(R), manifold)
end

function Forms.invhodge(::Val{Pr}, ::Val{R},
                        manifold::Manifold{D,C,S}) where {R,D,C,S}
    @assert 0 ≤ R ≤ D
    # TODO: sign of metric
    s = bitsign(R * (D - R))
    # s = bitsign(iseven(D) ? 0 : iodd(R))
    return s * hodge(Val(Pr), Val(R), manifold)
end

Forms.hodge(f::Fun{D,P,R}) where {D,P,R} = hodge(Val(P), Val(R), f.manifold) * f
function Forms.invhodge(f::Fun{D,P,R}) where {D,P,R}
    return invhodge(Val(P), Val(R), f.manifold) * f
end

# Derivatives

export coderiv
function coderiv(::Val{P}, ::Val{R},
                 manifold::Manifold{D,C,S}) where {P,R,D,C,S}
    P::PrimalDual
    dR = P == Pr ? +1 : -1
    @assert 0 ≤ R ≤ D
    @assert 0 ≤ R - dR ≤ D
    # println("[coderiv(D=$D, P=$P, R=$R)]")
    return (bitsign(R) *
            invhodge(Val(!P), Val(R - dR), manifold) *
            deriv(Val(!P), Val(R), manifold) *
            hodge(Val(P), Val(R), manifold))::Op{D,P,R - dR,P,R,S}
end

coderiv(f::Fun{D,P,R}) where {D,P,R} = coderiv(Val(P), Val(R), f.manifold) * f

# Δ
export laplace
# TODO: memoise laplace operator (and hodge, coderiv, etc.)
function laplace(::Val{P}, ::Val{R},
                 manifold::Manifold{D,C,S}) where {P,R,D,C,S}
    P::PrimalDual
    @assert 0 ≤ R ≤ D
    dR = P == Pr ? +1 : -1
    # println("[laplace(D=$D, P=$P, R=$R)]")
    op = zero(Op{D,P,R,P,R,S}, manifold)
    if 0 ≤ R - dR ≤ D
        op += deriv(Val(P), Val(R - dR), manifold) *
              coderiv(Val(P), Val(R), manifold)
    end
    if 0 ≤ R + dR ≤ D
        op += coderiv(Val(P), Val(R + dR), manifold) *
              deriv(Val(P), Val(R), manifold)
    end
    return op::Op{D,P,R,P,R,S}
end

laplace(f::Fun{D,P,R}) where {D,P,R} = laplace(Val(P), Val(R), f.manifold) * f

# ∫
export integral
integral(f::Fun{D,Pr,D}) where {D} = sum(f.values)
integral(f::Fun{D,Dl,0}) where {D} = sum(f.values)

# ⋅
export dot
function LinearAlgebra.dot(f::Fun{D,Pr,D}, g::Fun{D,Pr,D}) where {D}
    @assert f.manifold ≡ g.manifold
    vol = get_volumes(f.manifold, D)
    return dot(f.values, Diagonal(1 ./ vol), g.values)
end
function LinearAlgebra.dot(f::Fun{D,Dl,0}, g::Fun{D,Dl,0}) where {D}
    @assert f.manifold ≡ g.manifold
    dualvol = get_dualvolumes(f.manifold, 0)
    return dot(f.values, Diagonal(1 ./ dualvol), g.values)
end

# LinearAlgebra.norm(f::Fun) = norm(f.values)
# LinearAlgebra.norm(f::Fun, p::Real) = norm(f.values, p)

end


text:
# Finite Group Theory Routines
module FiniteGroups

using LinearAlgebra, SparseArrays, Base.Threads
import Base: ==, *
import LinearAlgebra: dot
export name, order, class, inclass, mult

"""
    AbstractFiniteGroup

Abstract type of finite groups.
Group elements are labeled by integers.
The identity element are 1 by default.
"""
abstract type AbstractFiniteGroup end

function Base.display(g::AbstractFiniteGroup)
    println("Finite group : $(name(g))")
    println("Group order  : $(order(g))")
    println("Classes      : $(length(class(g)))")
end

include("FiniteGroup.jl")
include("Character.jl")
include("Representation.jl")
include("ProjReps.jl")
include("Decomposition.jl")
include("PointGroups.jl")
include("Permutation.jl")
include("Precompile.jl")

end


text:
function avgfilter1d_setup(n, w = 3)
    xs = rand(SVector{w,UInt8}, n)
    ys = zero(xs)
    return (; ys, xs)
end

function avgfilter1d_seq!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = length(xs) - 2
    for i0 in 0:n-1
        i = firstindex(xs) + i0 + 1
        @inbounds ys[i] = (xs[i-1] .+ xs[i] .+ xs[i+1]) .÷ 3
    end
    return ys
end

function avgfilter1d_threads!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = length(xs) - 2
    Threads.@threads for i0 in 0:n-1
        i = firstindex(xs) + i0 + 1
        @inbounds ys[i] = (xs[i-1] .+ xs[i] .+ xs[i+1]) .÷ 3
    end
    return ys
end

function avgfilter1d_tapir_dac!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = length(xs) - 2
    Tapir.@par for i0 in 0:n-1
        i = firstindex(xs) + i0 + 1
        @inbounds ys[i] = (xs[i-1] .+ xs[i] .+ xs[i+1]) .÷ 3
        @grainsize 131072
    end
    return ys
end

function avgfilter1d_tapir_seq!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = length(xs) - 2
    Tapir.@par seq for i0 in 0:n-1
        i = firstindex(xs) + i0 + 1
        @inbounds ys[i] = (xs[i-1] .+ xs[i] .+ xs[i+1]) .÷ 3
        @grainsize 131072
    end
    return ys
end


text:
target = "gasp/libgasp.$(Libdl.dlext)"

if !isfile(target)
    @static if is_linux() || is_apple()
        repo = LibGit2.clone("https://github.com/kpamnany/gasp", "gasp")
        LibGit2.branch!(repo, "v0.5")
        println("Compiling libgasp...")
        run(`make -C gasp`)
    else
	error("gasp is currently Linux/OS-X only")
    end
end


text:
module FourRoomsUtil

using ..Resampling
import ..Resampling.FourRoomsParams
import ..Resampling: AbstractPolicy
using ..Reproduce
import StatsBase
import StatsBase: Weights
using Random
import ..JuliaRL
import ..build_algorithm_dict

include("FourRooms/policies.jl")

export FourRoomsAgent

include("FourRooms/agent.jl")

export
    env_settings!,
    get_experience,
    policy_settings,
    max_is_dict,
    get_policy

const constructed_policy_settings = Dict(
    "uniform"=>UniformRandomPolicy(4),
    "corners"=>FourCornersVariant(),
    "favor_down"=>FavoredRandomPolicy(4, FourRoomsParams.DOWN, 0.1),
    "random_state_variant"=>RandomStateVariant()
)

const GVFS = Dict(
    "down"=>GVF(FunctionalCumulant((state_t, action_t, state_tp1, action_tp1, preds_tp1)-> (state_tp1 == state_t ? 1.0 : 0.0)),
                StateTerminationDiscount(0.9, (state_t, action_t, state_tp1)->(state_tp1 == state_t)),
                PersistentPolicy(3)),
    "favored_down"=>GVF(FunctionalCumulant((state_t, action_t, state_tp1, action_tp1, preds_tp1)-> (state_tp1 == state_t ? 1.0 : 0.0)),
                        StateTerminationDiscount(0.9, (state_t, action_t, state_tp1)->(state_tp1 == state_t)),
                        FavoredRandomPolicy(4, 3, 0.9))
)

const max_is_dict = Dict(
    "easy"=>1.8,
    "hard"=>9.0,
    "hardest"=>99.0)

function env_settings!(s::ArgParseSettings)
    @add_arg_table s begin
        "--policy"
        help = "The policy settings: $(keys(constructed_policy_settings))"
        range_tester=(p)->(p ∈ keys(constructed_policy_settings))
        required=true
        "--gvf"
        help = "The gvf to learn: $(keys(GVFS))"
        range_tester=(gvf)->(gvf ∈ keys(GVFS))
        required=true
        # "--policy_input"
        # help="Inputs to the policy"
        # nargs='+'
    end
end

# function get_policy(policy)
#     return policy_settings[policy][1], policy_settings[policy][2]
# end

function get_policy(parsed::Dict)
    return constructed_policy_settings[parsed["policy"]]
end

end


text:
struct HyperRectangle{N, T}
    origin::SVector{N, T}
    widths::SVector{N, T}
end

@generated function vertices(rect::HyperRectangle{N}) where N
    verts = Expr[]
    for I in CartesianIndices(tuple([2 for i in 1:N]...))
        push!(verts,
        Expr(:call,
        :SVector,
        [I[i] == 1 ? :(rect.origin[$i]) : :(rect.origin[$i] + rect.widths[$i]) for i in 1:N]...))
    end
    Expr(:call, :TwosArray, verts...)
end

center(rect::HyperRectangle) = rect.origin + 0.5 * rect.widths

convert(::Type{HyperRectangle{N, T2}}, r::HyperRectangle{N, T1}) where {N, T1, T2} =
    HyperRectangle{N, T2}(convert(SVector{N, T2}, r.origin),
                          convert(SVector{N, T2}, r.widths))

@generated function faces(rect::HyperRectangle{N}) where N
    quote
        verts = vertices(rect)
        SVector($(Expr(:tuple, [:(selectdim(verts, $d, $i)) for i in 1:2 for d in 1:N]...)))
    end
end

struct BodyAndFaceCenters{N, T}
    rect::HyperRectangle{N, T}
end

function Base.iterate(b::BodyAndFaceCenters{N, T}, i=0) where {N, T}
    if i > 2 * N
        return nothing
    elseif i == 0
        return center(b.rect), i + 1
    else
        @inbounds f = faces(b.rect)[i]
        return sum(f) ./ length(f), i + 1
    end
end

Base.length(iter::BodyAndFaceCenters{N}) where {N} = 2 * N + 1
Base.eltype(iter::BodyAndFaceCenters{N, T}) where {N, T} = SVector{N, T}

body_and_face_centers(rect::HyperRectangle) = BodyAndFaceCenters(rect)


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using Distributions # to define priors
using Random        # to seed the random number generator
using DataFrames
using CSV

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "read_and_write_data.jl"))
include(joinpath(SRC_DIR, DIR, "transforms.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

# decide if first passage time observations or partially observed diffusion
fptObsFlag = true

# pick dataset
filename = "path_fpt_simpleConjug.csv"

# fetch the data
(df, x0, obs, obs_time, fpt,
      fptOrPartObs) = readData(Val(fptObsFlag), joinpath(OUT_DIR, filename))
param = :complexConjug
θ_init = [10.0, -8.0, 25.0, 0.0, 3.0]
P˟ = FitzhughDiffusion(param, θ_init...)
P̃ = std_aux_laws(FitzhughDiffusionAux, param, θ_init, obs, obs_time, 1)
L = @SMatrix [1. 0.]
Σ = @SMatrix [10^(-10)]

model_setup = DiffusionSetup(P˟, P̃, fptOrPartObs)
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time, fpt)
set_imputation_grid!(model_setup, 1/5000)
set_x0_prior!(model_setup, KnownStartingPt(x0))
initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100))
set_auxiliary!(model_setup; skip_for_save=10^1, adaptive_prop=NoAdaptation())

blocks = create_blocks( ChequeredBlocking(), model_setup.P,
                        (knots=collect(1:length(obs)-2)[1:1:end],
                         ϵ=10^(-10),
                         change_pt=SimpleChangePt(100)) )
mcmc_setup = MCMCSetup(
      Imputation(blocks[1], 0.99, Vern7()),
      Imputation(blocks[2], 0.99, Vern7()),
      ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 5))
                  ),
      ParamUpdate(ConjugateUpdt(), [1,2,3], θ_init, nothing,
                  MvNormal(fill(0.0, 3), diagm(0=>fill(1000.0, 3))),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2,3]))
                  ))

schedule = MCMCSchedule(1*10^4, [[1,3,4],[2,3,4]],
                        (save=3*10^2, verbose=10^2, warm_up=100,
                         readjust=(x->x%100==0), fuse=(x->false)))

Random.seed!(4)
out = mcmc(mcmc_setup, schedule, model_setup)
out, elapsed = @timeit mcmc(mcmc_setup, schedule, model_setup)
display(out.accpt_tracker)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out; truth=[10.0, -8.0, 25.0, 0.0, 3.0])
plot_paths(out; transf=[(x,θ)->x, (x,θ)->conjugToRegular(x, θ[1], 0)],
           obs=(times=obs_time[2:end], vals=obs[2:end], indices=1))


text:
"""
Get flatrate policy that optimizes the objective across a set of training histories
"""

"""
Returns lists of possible revenues and utilizations indexed by actions.

Used to select the best flatrate.
"""
function flatrate_analysis(mdp::PMDP, h::AbstractSimHistory)
    # extract request trace from history
    trace = collect(eachstep(h, "s, info"))
    requests = [rec for rec in trace if rec.s.iₚ !=PMDPs.empty_product_id(mdp)]

    # get data from trace
    request_edges = [[PMDPs.product(mdp, rec.s)...] for rec in requests]
    request_budgets = [rec.info.b for rec in requests]

    # Chech revenue of each flatrate
    r_as = PMDPs.Action[] # array containing total revenue for each possible flatrate
    u_as = Int64[] # array of final capacity for each flatrate
    for flatrate in POMDPs.actions(mdp)
        c_init = PMDPs.pp(mdp).c₀
        c = copy(c_init)
        r_a = 0.
        for i in 1:length(requests)
            if ~PMDPs.sale_impossible(mdp, requests[i].s) && PMDPs.user_buy(flatrate, request_budgets[i])
                c = PMDPs.reduce_capacities(c, product(mdp, requests[i].s))
                r_a +=flatrate
            end
        end
        push!(r_as, r_a)
        push!(u_as, sum(c_init - c))
    end

    return (r_a = r_as, u_a = u_as)
end

function optimize_flatrate_policy(mdp::PMDP, training_histories::AbstractArray{<:AbstractSimHistory})::Tuple

    # row index labels are mdp actions, column index labels are training histories
    R = []
    U = []
    for h in training_histories
        r_as, u_as = flatrate_analysis(mdp, h)
        push!(R, r_as)
        push!(U, u_as)
    end

    return hcat(R...), hcat(U...)
end

"""
Get flatrate policy from a training set of histories

To get policy action, use

    action(flatrate_policy, POMDPs.initialstate(mg))

"""
function get_flatrate_policy(mdp::PMDP, training_histories::Array{<:AbstractSimHistory}; objective=:revenue)::Policy
    R, U = optimize_flatrate_policy(mdp, training_histories)

    if objective==:revenue
        M = R
    elseif objective==:utilization
        M = U
    end

    best_results = sum(M; dims=2)
    best_flatrate = POMDPs.actions(mdp)[findmax(best_results)[2]]

    FunctionPolicy(s->best_flatrate)
end


text:
@testset "Test functionality of System" begin
    sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys"; add_forecasts = false)
    summary(devnull, sys)
    @test get_frequency(sys) == PSY.DEFAULT_SYSTEM_FREQUENCY

    generators = collect(get_components(ThermalStandard, sys))
    generator = get_component(ThermalStandard, sys, get_name(generators[1]))
    @test IS.get_uuid(generator) == IS.get_uuid(generators[1])
    @test_throws(IS.ArgumentError, add_component!(sys, generator))

    generators2 = get_components_by_name(ThermalGen, sys, get_name(generators[1]))
    @test length(generators2) == 1
    @test IS.get_uuid(generators2[1]) == IS.get_uuid(generators[1])
    @test !has_time_series(generators2[1])

    @test isnothing(get_component(ThermalStandard, sys, "not-a-name"))
    @test isempty(get_components_by_name(ThermalGen, sys, "not-a-name"))
    @test_throws(
        IS.ArgumentError,
        get_components_by_name(ThermalStandard, sys, "not-a-name")
    )
    @test isempty(get_components(ThermalStandard, sys, x -> (!get_available(x))))
    @test !isempty(get_available_components(ThermalStandard, sys))
    # Test get_bus* functionality.
    bus_numbers = Vector{Int}()
    for bus in get_components(Bus, sys)
        push!(bus_numbers, bus.number)
        if length(bus_numbers) >= 2
            break
        end
    end

    bus = PowerSystems.get_bus(sys, bus_numbers[1])
    @test bus.number == bus_numbers[1]

    buses = PowerSystems.get_buses(sys, Set(bus_numbers))
    sort!(bus_numbers)
    sort!(buses, by = x -> x.number)
    @test length(bus_numbers) == length(buses)
    for (bus_number, bus) in zip(bus_numbers, buses)
        @test bus_number == bus.number
    end

    @test get_forecast_initial_times(sys) == []
    @test get_time_series_resolution(sys) == Dates.Hour(1)

    # Get time_series with a name and without.
    components = collect(get_components(HydroEnergyReservoir, sys))
    @test !isempty(components)
    component = components[1]
    ts = get_time_series(SingleTimeSeries, component, "max_active_power")
    @test ts isa SingleTimeSeries

    returned_it, returned_len = check_time_series_consistency(sys, SingleTimeSeries)
    @test returned_it == first(TimeSeries.timestamp(get_data(ts)))
    @test returned_len == length(get_data(ts))

    # Test all versions of get_time_series_[array|timestamps|values]
    values1 = get_time_series_array(component, ts)
    values2 = get_time_series_array(SingleTimeSeries, component, "max_active_power")
    @test values1 == values2
    values3 = get_time_series_array(SingleTimeSeries, component, "max_active_power")
    @test values1 == values3

    val = get_time_series_array(SingleTimeSeries, component, "max_active_power")
    @test val isa TimeSeries.TimeArray
    val = get_time_series_timestamps(SingleTimeSeries, component, "max_active_power")
    @test val isa Array
    @test val[1] isa Dates.DateTime
    val = get_time_series_values(SingleTimeSeries, component, "max_active_power")
    @test val isa Array
    @test val[1] isa AbstractFloat

    val = get_time_series_array(component, ts)
    @test val isa TimeSeries.TimeArray
    val = get_time_series_timestamps(component, ts)
    @test val isa Array
    @test val[1] isa Dates.DateTime
    val = get_time_series_values(component, ts)
    @test val isa Array
    @test val[1] isa AbstractFloat

    clear_time_series!(sys)
    @test length(collect(get_time_series_multiple(sys))) == 0
    @test IS.get_internal(sys) isa IS.InfrastructureSystemsInternal
end

@testset "Test handling of bus_numbers" begin
    sys = create_rts_system()
    # TODO: sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")

    @test length(sys.bus_numbers) > 0
    buses = get_components(Bus, sys)
    bus_numbers = sort!([get_number(bus) for bus in buses])
    @test bus_numbers == get_bus_numbers(sys)

    # Remove some components
    remove_components!(x -> get_number(x) ∈ [101, 201], sys, Bus)
    @test length(sys.bus_numbers) == length(bus_numbers) - 2

    # Remove entire type
    remove_components!(sys, Bus)
    @test length(sys.bus_numbers) == 0

    # Remove individually.
    for bus in buses
        add_component!(sys, bus)
    end
    @test length(sys.bus_numbers) > 0
    for bus in buses
        remove_component!(sys, bus)
    end
    @test length(sys.bus_numbers) == 0

    # Remove by name.
    for bus in buses
        add_component!(sys, bus)
    end
    @test length(sys.bus_numbers) > 0
    for bus in buses
        remove_component!(Bus, sys, get_name(bus))
    end
    @test length(sys.bus_numbers) == 0
end

@testset "Test System iterators" begin
    sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")

    i = 0
    for component in iterate_components(sys)
        i += 1
    end

    components = get_components(Component, sys)
    @test i == length(components)

    # Test debugging functions.
    component = first(components)
    uuid = IS.get_uuid(component)
    @test get_name(get_component(sys, uuid)) == get_name(component)
    @test get_name(get_component(sys, string(uuid))) == get_name(component)
end

@testset "Test remove_component" begin
    sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")
    generators = get_components(ThermalStandard, sys)
    initial_length = length(generators)
    @assert initial_length > 0
    gen = collect(generators)[1]

    remove_component!(sys, gen)

    @test isnothing(get_component(typeof(gen), sys, get_name(gen)))
    generators = get_components(typeof(gen), sys)
    @test length(generators) == initial_length - 1

    @test_throws(IS.ArgumentError, remove_component!(sys, gen))

    add_component!(sys, gen)
    remove_component!(typeof(gen), sys, get_name(gen))
    @test isnothing(get_component(typeof(gen), sys, get_name(gen)))

    @assert length(get_components(typeof(gen), sys)) > 0
    remove_components!(sys, typeof(gen))
    @test_throws(IS.ArgumentError, remove_components!(sys, typeof(gen)))

    remove_components!(sys, Area)
    @test isempty(get_components(Area, sys))
    @test isnothing(get_area(collect(get_components(Bus, sys))[1]))
end

@testset "Test missing Arc bus" begin
    sys = System(100)
    line = Line(nothing)
    @test_throws(IS.ArgumentError, add_component!(sys, line))
end

@testset "Test frequency set" begin
    sys = System(100; frequency = 50.0)
    @test get_frequency(sys) == 50.0
end

@testset "Test exported names" begin
    @test IS.validate_exported_names(PowerSystems)
end

@testset "Test system ext" begin
    sys = System(100)
    ext = get_ext(sys)
    ext["data"] = 2
    @test get_ext(sys)["data"] == 2
    clear_ext!(sys)
    @test isempty(get_ext(sys))
end

@testset "Test system checks" begin
    sys = System(100)
    @test_logs (:warn, r"There are no .* Components in the System") match_mode = :any check(
        sys,
    )
end

@testset "Test system units" begin
    sys = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys"; add_forecasts = false)
    @test get_units_base(sys) == "DEVICE_BASE"
    set_units_base_system!(sys, "SYSTEM_BASE")
    @test get_units_base(sys) == "SYSTEM_BASE"
end

@testset "Test add_time_series multiple components" begin
    sys = System(100.0)
    bus = Bus(nothing)
    bus.bustype = BusTypes.REF
    add_component!(sys, bus)
    components = []
    len = 2
    for i in 1:len
        gen = ThermalStandard(nothing)
        gen.name = string(i)
        gen.bus = bus
        add_component!(sys, gen)
        push!(components, gen)
    end

    initial_time = Dates.DateTime("2020-01-01T00:00:00")
    end_time = Dates.DateTime("2020-01-01T23:00:00")
    dates = collect(initial_time:Dates.Hour(1):end_time)
    data = collect(1:24)
    ta = TimeSeries.TimeArray(dates, data, ["1"])
    name = "max_active_power"
    ts = SingleTimeSeries(name = name, data = ta)
    add_time_series!(sys, components, ts)

    for i in 1:len
        component = get_component(ThermalStandard, sys, string(i))
        ts = get_time_series(SingleTimeSeries, component, name)
        @test ts isa SingleTimeSeries
    end
end

@testset "Test set_name! of system component" begin
    sys = System(100.0)
    bus = Bus(nothing)
    bus.bustype = BusTypes.REF
    add_component!(sys, bus)
    set_name!(sys, bus, "new_name")
    @test get_component(Bus, sys, "new_name") === bus
end

@testset "Test forecast parameters" begin
    sys = System(100.0)
    bus = Bus(nothing)
    bus.bustype = BusTypes.REF
    add_component!(sys, bus)
    gen = ThermalStandard(nothing)
    gen.name = "gen"
    gen.bus = bus
    add_component!(sys, gen)

    resolution = Dates.Hour(1)
    initial_time = Dates.DateTime("2020-09-01")
    second_time = initial_time + resolution
    name = "test"
    horizon = 24
    data = SortedDict(initial_time => ones(horizon), second_time => ones(horizon))

    forecast = Deterministic(data = data, name = name, resolution = resolution)
    add_time_series!(sys, gen, forecast)

    @test get_time_series_resolution(sys) == resolution
    @test get_forecast_horizon(sys) == horizon
    @test get_forecast_initial_timestamp(sys) == initial_time
    @test get_forecast_interval(sys) == Dates.Millisecond(second_time - initial_time)
    @test get_forecast_window_count(sys) == 2
    @test Dates.Hour(get_forecast_total_period(sys)) ==
          Dates.Hour(second_time - initial_time) + Dates.Hour(resolution * horizon)
    @test get_forecast_initial_times(sys) == [initial_time, second_time]

    remove_time_series!(sys, typeof(forecast), gen, get_name(forecast))
    @test_throws ArgumentError get_time_series(typeof(forecast), gen, get_name(forecast))
end

@testset "Invalid constructor" begin
    @test_throws IS.DataFormatError System("data.invalid")
end

@testset "Test deepcopy with runchecks" begin
    sys = System(100.0)
    @test get_runchecks(sys)
    @test get_runchecks(deepcopy(sys))
end

@testset "Test deepcopy with runchecks disabled" begin
    sys = System(100.0, runchecks = false)
    @test !get_runchecks(sys)
    sys2 = deepcopy(sys)
    @test sys2 isa System
    @test !get_runchecks(sys)
end

@testset "Test time series counts" begin
    c_sys5 = PSB.build_system(
        PSITestSystems,
        "c_sys5_uc";
        add_forecasts = true,
        force_build = true,
    )
    _, ts_count, forecast_count = get_time_series_counts(c_sys5)
    @test ts_count == 0
    @test forecast_count == 3
end

@testset "Test deepcopy with time series options" begin
    sys = PSB.build_system(
        PSITestSystems,
        "test_RTS_GMLC_sys";
        time_series_in_memory = true,
        force_build = true,
    )
    @test sys.data.time_series_storage isa IS.InMemoryTimeSeriesStorage
    sys2 = deepcopy(sys)
    @test sys2.data.time_series_storage isa IS.InMemoryTimeSeriesStorage
    @test IS.compare_values(sys, sys2)
    # Ensure that the storage references got updated correctly.
    for component in get_components(Component, sys2, x -> has_time_series(x))
        @test component.time_series_container.time_series_storage ===
              sys2.data.time_series_storage
    end

    sys = PSB.build_system(
        PSITestSystems,
        "test_RTS_GMLC_sys";
        time_series_in_memory = false,
        force_build = true,
    )
    @test sys.data.time_series_storage isa IS.Hdf5TimeSeriesStorage
    sys2 = deepcopy(sys)
    @test sys2.data.time_series_storage isa IS.Hdf5TimeSeriesStorage
    @test sys.data.time_series_storage.file_path != sys2.data.time_series_storage.file_path
    @test IS.compare_values(sys, sys2)
    for component in get_components(Component, sys2, x -> has_time_series(x))
        @test component.time_series_container.time_series_storage ===
              sys2.data.time_series_storage
    end
end

@testset "Test with compression enabled" begin
    @test get_compression_settings(System(100.0)) == CompressionSettings(enabled = false)

    settings = CompressionSettings(enabled = true, type = CompressionTypes.BLOSC)
    @test get_compression_settings(System(100.0, compression = settings)) == settings
    @test get_compression_settings(System(100.0, enable_compression = true)) ==
          CompressionSettings(enabled = true)
end

@testset "Test compare_values" begin
    sys1 = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")
    sys2 = PSB.build_system(PSITestSystems, "test_RTS_GMLC_sys")
    gen1 = first(get_components(ThermalStandard, sys1))
    gen2 = first(get_components(ThermalStandard, sys2))
    @test IS.compare_values(gen1, gen2)
    @test IS.compare_values(sys1, sys2)

    set_active_power!(gen1, get_active_power(gen1) + 0.1)
    @test(
        @test_logs(
            (:error, r"not match"),
            match_mode = :any,
            !IS.compare_values(gen1, gen2),
        )
    )
    @test(
        @test_logs(
            (:error, r"not match"),
            match_mode = :any,
            !IS.compare_values(sys1, sys2)
        )
    )
end