text:
const error1 = 1e-2 # relative error for adaptive integration
const Nc = 3.0
const Λ = 10^5.0; #UV sacle
const δ = 0.1;
const cmax = 10.0^16
const a = 0.05; # regulator coeff

const kbar = 128.989/0.955;
const pmin = 1.0;
const pmax = Λ

const N_pgrid = 200 #number of points along p axis

const N_pstep = (log(pmax) - log(pmin)) / (N_pgrid - 1)
const pgrid = exp.(collect(log(pmin):N_pstep:log(pmax)))

# -5.1285020065
# -5.1585*10^7

# -5.13800708*10^7

# -5.001256120152*10^7

const inihatza = fill(0.9999999999999999, N_pgrid) #initial condition for hatza, we need transform to za later
const iniza = @. inihatza * (pgrid^2 + deltamT2) / (pgrid^2) #initial condition for za
# iniza[1] = (deltamT2 / (pmin^2)) * inihatza[1] #we modify the first point of za to avoid the non zero pmin effect
const inizc = fill(0.9999999999999999, N_pgrid)#initial condition for zc
const iniλcA = fill(0.9999999999999999, N_pgrid)
const iniλ3A = fill(0.9999999999999999, N_pgrid)

# Plots.plot(pgrid, [iniza.^-1],xaxis=:log)
# Plots.plot(pgrid, [@. round((pgrid^2 * (pgrid^2 + deltamT2)^-1 * iniza)^-1; sigdigits=5)],xaxis=:log)

# Nx = 128 # number of gauss point along ps
# Ny = 8 # number of gauss point along cosθ

# nodesx, weightsx = cu(gausslegendre(Nx)) # gauss legendre point and weights along ps
# nodesy, weightsy = cu(gausslegendre(Ny)) # gauss legendre point and weights along cosθ

# zakLamTfun = extrapolate(interpolate((pgrid,), iniza, Gridded(Linear())), Flat()) # Interpolate initial condition of za
# zckLamTfun = extrapolate(interpolate((pgrid,), inizc, Gridded(Linear())), Flat()) # Interpolate initial condition of zc

zakLamTfun = Spline1D(pgrid, iniza)
zckLamTfun = Spline1D(pgrid, inizc)
λcALamTfun = Spline1D(pgrid, iniλcA)
λ3ALamTfun = Spline1D(pgrid, iniλ3A)

# λ3Afun = extrapolate(
#     interpolate(
#         (readdlm("pgrid.txt")[:, 1], readdlm("kout.txt")[:, 1]),
#         readdlm("lambda3a.dat"),
#         Gridded(Linear()),
#     ),
#     Flat(),
# )
#
# λ3Afun = Spline2D(
#     readdlm("pgrid.txt")[:, 1],
#     readdlm("kout.txt")[:, 1],
#     readdlm("lambda3a.dat"),
# )

# λcAfun = extrapolate(
#     interpolate(
#         (readdlm("pgrid.txt")[:, 1], readdlm("kout.txt")[:, 1]),
#         readdlm("lambdacca.dat"),
#         Gridded(Linear()),
#     ),
#     Flat(),
# )
# λcAfun = Spline2D(
#     readdlm("pgrid.txt")[:, 1],
#     readdlm("kout.txt")[:, 1],
#     readdlm("lambdacca.dat"),
# )

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
"""
    writeloghist(fname,data;opts...)

Write a log-scaled histogram to a the file `fname`.
"""
function writeloghist(fname,data;opts...)
    x,y = loghistogram(data;opts...)
    x1 = truncfloat.(x,1)
    writedlm(fname,zip(x1,y)," ")
end

"""
    loghistogram(data;opts...)

Returns a tuple, the coordinate and the ordinates for a log scale
histogram. `opts...` are passed to `fit(Histogram,...`.
The output `res` can be written to a file with `writedlm("f.txt",zip(res...))`.
This function is in lieu of a proper log-scaled histogram.
"""
function loghistogram(data;opts...)
    hist = fit(Histogram,log.(data); opts..., closed=:left)
    collect(exp.((hist.edges)[1])), hist.weights
end

"""
    append(a::Vector)::Vector

Flatten a vector of vectors, all of the same
eltype.
"""
# Don't restrict the input type, so we can use iterators.
function append(a)::Vector
    tp = eltype(a[1])
    vout = tp[]
    for x in a
        for y in x
            push!(vout,y)
        end
    end
    vout
end

"""
    arrcountmap(a1::Vector, a2::Vector{T}) where T<:Integer

Sum the counts in each element of `a2` for each occurrence of an item
in all of the elements of `a1`. Return a countmap of these sums for each of the
unique elements in the arrays in `a1`.

Example:
```julia
julia> d = Dict( ["a","b"] => 1, ["b", "c"] => 2 );

julia> arrcountmap(keys(d), values(d))
DataStructures.Accumulator{String,Int64} with 3 entries:
  "c" => 2
  "b" => 3
  "a" => 1
```
"""
# This allows iterators for `a1` and `a2` to be used
function arrcountmap(a1, a2)::DataStructures.Accumulator
    acc = counter(a1 |> eltype |> eltype)
    for (teams,nclicks) in zip(a1,a2)
        for team in teams
            push!(acc,team,nclicks)
        end
    end
    acc
end

"""
    arrcountmap(d::AbstractDict)

Pass the keys and values of `d` as the two arrays to `arrcountmap`.
"""
arrcountmap(d::AbstractDict) = arrcountmap(keys(d),values(d))

# Catches errors, but prevents iterators from working
# function arrcountmap(a1::Vector, a2::Vector{T})::DataStructures.Accumulator where T<:Integer
#     acc = counter(a1 |> eltype |> eltype)
#     for (teams,nclicks) in zip(a1,a2)
#         for team in teams
#             push!(acc,team,nclicks)
#         end
#     end
#     acc
# end

"""
    subst!(arr,from,to)::Void

Substitute each occurrence of `from` in the values of collection `arr` with `to`.
"""
function subst!(arr,from,to)::Void
    for i in linearindices(arr)
    @inbounds  if arr[i] == from
            arr[i] = to
        end
    end
end

"""
    subst!(d::AbstractDict,from,to)::Void

Substitute each occurrence of `from` in the values of `d` with `to`.
"""
function subst!(d::AbstractDict,from,to)::Void
    for (k,v) in d
         if d[k] == from
            d[k] = to
        end
    end
end

"""
    keepkeys(dict::AbstractDict, keylist::AbstractArray)

Return a copy of `dict` keeping only keys in `keylist`.
It makes sense to broadcast over the first argument more than
it does the second, but this does not work at the moment.
"""
keepkeys(dict::AbstractDict, keylist::AbstractArray) = keepkeys(dict,typeof(dict),keylist)

"""
    keepkeys(dict::AbstractDict, keylist::AbstractArray, T::DataType)

Specify the the type `T` of the returned dictionary.
"""
keepkeys(dict::AbstractDict, keylist::AbstractArray, T::DataType) = T(k => dict[k] for k in keylist)

"""
    arraytodict(arr::AbstractArray, pkey; dicttype=Dict)

Construct a dictionary of arrays from the array of dictionaries `arr`.
The keys of the output dictionary are the values correpsonding to
the key `pkey` in each dictionary in `arr`. Each value in the
output dictionary is an array of dictionaries
each with the same value for key `pkey`.

This function could use a better name.
"""
function arraytodict(arr::AbstractArray, pkey; dicttype=Dict)
    kt = typeof(arr[1][pkey])
    dout = dicttype{kt,Any}() # could make this more specific
    for obj in arr
        nkey = obj[pkey]
        haskey(dout,nkey) ? nothing : dout[nkey] = Any[]
        push!(dout[nkey],obj)
    end
    return dout
end

"""
    fields(x)

Unpack the fields in struct `x` into a `Tuple`.
"""
fields(x) = ((getfield(x, field) for field in fieldnames(typeof(x)))...,)

ext: jl


text:
# Defines a differential equation describing changes in concentration with respect to time
# due to advection and reaction in a hyporheic zone with a power law residence time distribution
# and first-order reaction kinetics.

  macro advectReact(Δ, bc)
    return :(
        function advectionReaction!(du,u,p,t)
            τ₋, τ₊, α, k, q, v = p
             integral = zeros(length(u))
            if(τ₋<t<τ₊)
                integral = QuadGK.quadgk(τ -> ((τ^(-α) - τ₊^(-α))*(-α+1)/(τ₊^(-α+1)-τ₋^(-α+1)-τ₊^(-α)*(τ₊ - τ₋)*(1-α))) * h(p,t-τ) * exp(-k*τ),τ₋, t, τ₊)[1]
            else
                integral = QuadGK.quadgk(τ -> ((τ^(-α) - τ₊^(-α))*(-α+1)/(τ₊^(-α+1)-τ₋^(-α+1)-τ₊^(-α)*(τ₊- τ₋)*(1-α))) * h(p,t-τ) * exp(-k*τ), τ₋, τ₊)[1]
            end
            du = -v*Δ*bc*u + q*(integral - u)
          end
    )
  end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
using Jute
using BenchmarkTools

include("TestUtils.jl")

include("builtin_fixtures.test.jl")
include("fixture_factory.test.jl")
include("fixtures.test.jl")
include("options.test.jl")
include("performance.test.jl")
include("reporting.test.jl")
include("run_testcase.test.jl")
include("runtests.test.jl")
include("testcases.test.jl")
include("utils.test.jl")

exit(runtests())


text:
module SimpleTypePrint

export show_type, repr_type, config_type_display

"""
# Examples
```jldoctest
julia> using SimpleTypePrint: print_multiple
       print_multiple(show, stdout, 1:10, "{", ",", "}")
{1,2,3,4,5,6,7,8,9,10}
```
"""
function print_multiple(f, io::IO, xs, left="{", sep=",", right="}")
    n = length(xs)::Int
    print(io, left)
    for (i, x) in enumerate(xs)
        f(io,x)
        i < n && print(io, sep)
    end
    print(io, right)
end

"""
    show_type(io, type; kwargs...)
"""
function show_type(io::IO, @nospecialize(ty::Type); max_depth::Int = 3, short_type_name::Bool = true)
    t_var_scope::Set{Symbol} = Set{Symbol}()

    function rec(x::DataType, d)
        # to be consistent with the Julia Compiler.
        (x === Tuple) && return Base.show_type_name(io, x.name)
        if x.name === Tuple.name
            n = length(x.parameters)::Int
            if n > 3 && all(i -> (x.parameters[1] === i), x.parameters)
                # Print homogeneous tuples with more than 3 elements compactly as NTuple
                print(io, "NTuple{", n, ',')
                rec(x.parameters[1], d-1)
                print(io, "}")
                return
            end
        end

        # hanlde function type specially
        is_func = startswith(string(nameof(x)), "#")
        if short_type_name && !is_func
             print(io, nameof(x))
        else
            Base.show_type_name(io, x.name)
        end

        if !isempty(x.parameters)
            if d ≤ 1
                print(io, "{...}")
            else
                print_multiple((_,p) -> rec(p, d-1), io, x.parameters)
            end
        end
    end

    function rec(x::Union, d)
        if x.a isa DataType && Core.Compiler.typename(x.a) === Core.Compiler.typename(DenseArray)
            T, N = x.a.parameters
            if x == StridedArray{T,N}
                print(io, "StridedArray")
                print_multiple(show, io, (T,N))
                return
            elseif x == StridedVecOrMat{T}
                print(io, "StridedVecOrMat")
                print_multiple(show, io, (T,))
                return
            elseif StridedArray{T,N} <: x
                print(io, "Union")
                elems = vcat(StridedArray{T,N},
                    Base.uniontypes(Core.Compiler.typesubtract(x, StridedArray{T,N})))
                print_multiple((_, p) -> rec(p, d-1), io, elems)
                return
            end
        end
        print(io, "Union")
        print_multiple((_, p) -> rec(p, d-1), io, Base.uniontypes(x))
    end

    function show_tv_def(tv::TypeVar, d)
        function show_bound(io::IO, @nospecialize(b))
            parens = isa(b,UnionAll) && !Base.print_without_params(b)
            parens && print(io, "(")
            rec(b, d-1)
            parens && print(io, ")")
        end
        lb, ub = tv.lb, tv.ub
        if lb !== Base.Bottom
            if ub === Any
                Base.show_unquoted(io, tv.name)
                print(io, ">:")
                show_bound(io, lb)
            else
                show_bound(io, lb)
                print(io, "<:")
                Base.show_unquoted(io, tv.name)
            end
        else
            Base.show_unquoted(io, tv.name)
        end
        if ub !== Any
            print(io, "<:")
            show_bound(io, ub)
        end
        nothing
    end

    function rec(x::UnionAll, d, var_group::Vector{TypeVar}=TypeVar[])
        # rename tvar as needed
        var_symb = x.var.name
        if var_symb === :_ || var_symb ∈ t_var_scope
            counter = 1
            while true
                newname = Symbol(var_symb, counter)
                if newname ∉ t_var_scope
                    newtv = TypeVar(newname, x.var.lb, x.var.ub)
                    x = UnionAll(newtv, x{newtv})
                    break
                end
                counter += 1
            end
        end
        var_symb = x.var.name

        push!(var_group, x.var)
        push!(t_var_scope, var_symb)
        if x.body isa UnionAll
            # current var group continues
            rec(x.body, d, var_group)
        else
            # current var group ends
            rec(x.body, d)
            print(io, " where ")
            if length(var_group) == 1
                show_tv_def(var_group[1], d)
            else
                print_multiple((_, v) -> show_tv_def(v, d), io, var_group)
            end
        end
        delete!(t_var_scope, var_symb)
    end

    function rec(tv::TypeVar, d)
        Base.show_unquoted(io, tv.name)
    end

    function rec(ohter, d)
        @nospecialize(other)
        show(io, ohter)
    end

    rec(ty, max_depth)
end

"""
    repr_type(ty::Type; kwargs...)::String
"""
repr_type(ty::Type; kwargs...) = sprint((io,t) -> show_type(io,t; kwargs...), ty)

"""
    config_type_display(;kwargs...)

Replace `Base.show(io::IO, x::Type)` with the simpler type printing funciton
`show_type`.

See also: `show_type`, `repr_type`.
"""
function config_type_display(;kwargs...)
    @eval Base.show(io::IO, x::Type) = show_type(io, x; $(kwargs)...)
end

"""
# Keyword args
- `max_depth=3`: the maximal type AST depth to show. Type arguments deeper than this value
will be printed as `...`.
- `short_type_name=true`: when set to `true`, will print simple type names without their
corresponding module path. e.g. "Name" instead of "ModuleA.ModuleB.Name". Note that the
shorter name will always be used if the type is visible from the current scope.
"""
show_type, repr_type, config_type_display

end  # module SimpleTypePrint

max_forks_repo_name: MrVPlusOne/SimpleTypePrint


text:
@testset "1319.number-of-operations-to-make-network-connected.jl" begin
    @test make_connected(4, [[1,2],[1,3],[2,3]]) == 1
    @test make_connected(6, [[1, 2], [1, 3], [1, 4], [2, 3], [2, 4]]) == 2
    @test make_connected(6, [[1, 2], [1, 3], [1, 4], [2, 3]]) == -1
    @test make_connected(5, [[1, 2], [1, 3], [4, 5], [3, 4]]) == 0
end


text:
module stream

self.list = ((a, x) -> a.append(x), () -> [])
self.sum = (a, x) -> a + x

; Lazy-computes map() and collect() operations on iterables.
class Stream {
  field fns
  new(field iterable) {
    fns = []
  }
  ; Applies [fn] to all members of the stream.
  ; The order of the iterable is unaffected.
  method map(fn) {
    fns.append(fn)
    return self
  }
  ; Applies [fn] to the each member of the stream.
  ;
  method each(fn) {
    for (_, elt) in iterable {
      result = _apply_fns(elt)
      fn(result)
    }
  }
  method collect(fn, fn_acc=None) {
    acc = if fn_acc then fn_acc() else _apply_fns(iterable[0])
    for (i, elt) in iterable {
      if i == 0 and ~fn_acc {
        continue
      }
      acc = fn(acc, _apply_fns(elt))
    }
    return acc
  }
  method _apply_fns(elt) {
    result = elt
    for (_, func) in fns {
      result = func(result)
    }
    return result
  }
}


text:
using Documenter
using Meshing

makedocs(
    sitename = "Meshing",
    format = Documenter.HTML(),
    modules = [Meshing],
    pages = ["Index" => "index.md",
             "API" => "api.md",
             "Examples" => "examples.md",
             "Internals" => "internals.md"]
)

deploydocs(
    repo = "github.com/JuliaGeometry/Meshing.jl.git",
)


text:
mutable struct SimulationInternal
    sim_files_dir::String
    store_dir::String
    logs_dir::String
    models_dir::String
    recorder_dir::String
    results_dir::String
    run_count::OrderedDict{Int, OrderedDict{Int, Int}}
    date_ref::OrderedDict{Int, Dates.DateTime}
    status::RunStatus
    build_status::BuildStatus
    simulation_state::SimulationState
    store::Union{Nothing, SimulationStore}
    recorders::Vector{Symbol}
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    cache_size_mib::Int
    min_cache_flush_size_mib::Int
end

function SimulationInternal(
    steps::Int,
    models::SimulationModels,
    sim_dir::String,
    name::String,
    output_dir::Union{Nothing, String},
    recorders,
    console_level::Logging.LogLevel,
    file_level::Logging.LogLevel,
    cache_size_mib=1024,
    min_cache_flush_size_mib=MIN_CACHE_FLUSH_SIZE_MiB,
)
    count_dict = Dict{Int, Dict{Int, Int}}()

    for s in 1:steps
        count_dict[s] = Dict{Int, Int}()
        model_count = length(get_decision_models(models))
        for st in 1:model_count
            count_dict[s][st] = 0
        end
        if get_emulation_model(models) !== nothing
            count_dict[s][model_count + 1] = 0
        end
    end

    base_dir = joinpath(sim_dir, name)
    mkpath(base_dir)

    output_dir = _get_output_dir_name(base_dir, output_dir)
    simulation_dir = joinpath(base_dir, output_dir)
    if isdir(simulation_dir)
        error("$simulation_dir already exists. Delete it or pass a different output_dir.")
    end

    sim_files_dir = joinpath(simulation_dir, "simulation_files")
    store_dir = joinpath(simulation_dir, "data_store")
    logs_dir = joinpath(simulation_dir, "logs")
    models_dir = joinpath(simulation_dir, "problems")
    recorder_dir = joinpath(simulation_dir, "recorder")
    results_dir = joinpath(simulation_dir, RESULTS_DIR)

    for path in (
        simulation_dir,
        sim_files_dir,
        logs_dir,
        models_dir,
        recorder_dir,
        results_dir,
        store_dir,
    )
        mkpath(path)
    end

    unique_recorders = Set(REQUIRED_RECORDERS)
    foreach(x -> push!(unique_recorders, x), recorders)

    return SimulationInternal(
        sim_files_dir,
        store_dir,
        logs_dir,
        models_dir,
        recorder_dir,
        results_dir,
        count_dict,
        Dict{Int, Dates.DateTime}(),
        RunStatus.NOT_READY,
        BuildStatus.EMPTY,
        SimulationState(),
        nothing,
        collect(unique_recorders),
        console_level,
        file_level,
        cache_size_mib,
        min_cache_flush_size_mib,
    )
end

function _get_output_dir_name(path, output_dir)
    if output_dir !== nothing
        # The user wants a custom name.
        return output_dir
    end

    # Return the next highest integer.
    output_dir = 1
    for name in readdir(path)
        if occursin(r"^\d+$", name)
            num = parse(Int, name)
            if num >= output_dir
                output_dir = num + 1
            end
        end
    end

    return string(output_dir)
end

function configure_logging(internal::SimulationInternal, file_mode)
    return IS.configure_logging(
        console=true,
        console_stream=stderr,
        console_level=internal.console_level,
        file=true,
        filename=joinpath(internal.logs_dir, SIMULATION_LOG_FILENAME),
        file_level=internal.file_level,
        file_mode=file_mode,
        tracker=nothing,
        set_global=false,
    )
end

function register_recorders!(internal::SimulationInternal, file_mode)
    for name in internal.recorders
        IS.register_recorder!(name; mode=file_mode, directory=internal.recorder_dir)
    end
end

function unregister_recorders!(internal::SimulationInternal)
    for name in internal.recorders
        IS.unregister_recorder!(name)
    end
end


text:
@testset "Model" begin
    input = DataFrame(; a=1, b=2)

    model = Model(df -> df.a + 2 * df.b, :c)
    @test isa(model, Model)

    evaluate!(model, input)
    @test input == DataFrame(; a=1, b=2, c=5)
    @test model(input) == [5]
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
#   Unit "csp.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019, 2020
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements Common Spatial Pattern (CSP) filters,
#   that is, filters based on the generalized eigendecomposition

"""
```
(1)
function csp(C₁ :: SorH, C₂ :: SorH;
             eVar     :: TeVaro = ○,
             eVarC    :: TeVaro = ○,
             eVarMeth :: Function = searchsortedfirst,
             selMeth  :: Symbol = :extremal,
             simple   :: Bool = false)

(2)
function csp(X₁ :: Mat, X₂ :: Mat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX₁   :: Tmean = 0,
             meanX₂   :: Tmean = 0,
             wX₁      :: Tw = ○,
             wX₂      :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarC    :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          selMeth  :: Symbol = :extremal,
          simple   :: Bool = false)

(3)
function csp(𝐗₁::VecMat, 𝐗₂::VecMat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX₁   :: Into = 0,
             meanX₂   :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarC    :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          selMeth  :: Symbol = :extremal,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       w₁       :: Vector = [],
       w₂       :: Vector = [],
       ✓w       :: Bool = true,
       init₁    :: SorHo = nothing,
       init₂    :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)
```

Return a [LinearFilter](@ref) object:

**(1) Common spatial pattern**
with covariance matrices `C_1` and `C_2` of dimension
``n⋅n`` as input. The subscript of the covariance matrices refers to the `dims`
used to compute it (see above).

`eVar`, `eVarC` and `eVarMeth` are keyword optional arguments
for defining the [subspace dimension](@ref) ``p``. Particularly:
-  By default, the two-step procedure described above is used to find the
   solution. In this case `eVarC` is used for defining the subspace dimension of
   the whitening step. If `eVarC=0.0` is passed (not to be confused with
   `eVarC=0` ), the solution will be find by the generalized
   eigenvalue-eigenvector procedure.
- `eVar` is the keyword optional argument for defining the
   [subspace dimension](@ref) ``p`` using the `.arev` vector
   given by [csp.5].
- `eVarMeth` applies to both `eVarC` and `eVar`. The default value is
   `evarMeth=searchsortedfirst`.

if `selMeth=:extremal` (default) use case **a) Separating two classes**
described above is considered. Any other symbol for `selMeth` will instruct
to consider instead the use case **b) Enhance the signal-to-noise ratio**.

If `simple` is set to `true`, ``p`` is set equal to ``n``
and only the fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

**(2) Common spatial pattern**
with data matrices `X₁` and `X₂` as input.

`X₁` and `X₂` are real or complex data matrices.

`covEst`, `dims`, `meanX₁`, `meanX₂`,  `wX₁` and `wX₂` are optional
keyword arguments to regulate the estimation of the
covariance matrices ``(C_1, C_2)`` of (`X₁`, `X₂`).
Particularly (See [covariance matrix estimations](@ref)),
- `meanX₁` is the `mean` argument for data matrix `X₁`.
- `meanX₂` is the `mean` argument for data matrix `X₂`.
- `wX₁` is the `w` argument for estimating a weighted covariance matrix for `X₁`.
- `wX₂` is the `w` argument for estimating a weighted covariance matrix for `X₂`.
- `covEst` applies to the estimations of both covariance matrices.

Once the two covariance matrices ``C_1`` and ``C_2`` estimated,
method (1) is invoked with optional keyword arguments
`eVar`, `eVarC`, `eVarMeth`, `selMeth` and `simple`.
See method (1) for details.

**(3) Common spatial pattern**
with two vectors of data matrices
`𝐗₁` and `𝐗₂` as input.

`𝐗₁` and `𝐗₂` do not need to hold the same number
of matrices and the number of samples in the matrices they contain
is arbitrary.

`covEst`, `dims`, `meanX₁` and `meanX₂` are optional
keyword arguments to regulate the estimation of the
covariance matrices for all matrices in `𝐗₁` and `𝐗₂`.
See method (2) and [covariance matrix estimations](@ref).

A mean covariance matrix is computed separatedly from the covariance matrices
computed from the data matrices in `𝐗₁` and `𝐗₂`,
using optional keywords arguments
`metric`, `w₁`, `w₂`, `✓w`,
`init₁`, `init₂`, `tol` and `verbose`. Particularly
(see [mean covariance matrix estimations](@ref)),
- `w₁` are the weights for the covariance matrices computed from `𝐗₁`,
- `w₂` are the weights for the covariance matrices computed from `𝐗₂`,
- `init₁` is the initialization for the mean of the covariance matrices computed from `𝐗₁`,
- `init₂` is the initialization for the mean of the covariance matrices computed from `𝐗₂`.
By default, the arithmetic mean is computed.

**See also:** [CSTP](@ref), [PCA](@ref), [AJD](@ref), [mAJD](@ref).

**Examples:**

```
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
t, n=50, 10
X1=genDataMatrix(n, t)
X2=genDataMatrix(n, t)
Cx1=Symmetric((X1*X1')/t)
Cx2=Symmetric((X2*X2')/t)
C=Cx1+Cx2
cC=csp(Cx1, Cx2; simple=true)
Dx1=cC.F'*Cx1*cC.F
@test norm(Dx1-Diagonal(Dx1))+1≈1.
Dx2=cC.F'*Cx2*cC.F
@test norm(Dx2-Diagonal(Dx2))+1≈1.
@test cC.F'*C*cC.F≈I
@test norm(Dx1-(I-Dx2))+1≈1.

# Method (1) complex
t, n=50, 10
X1c=genDataMatrix(ComplexF64, n, t)
X2c=genDataMatrix(ComplexF64, n, t)
Cx1c=Hermitian((X1c*X1c')/t)
Cx2c=Hermitian((X2c*X2c')/t)
Cc=Cx1c+Cx2c
cCc=csp(Cx1c, Cx2c; simple=true)
Dx1c=cCc.F'*Cx1c*cCc.F
@test norm(Dx1c-Diagonal(Dx1c))+1. ≈ 1.
Dx2c=cCc.F'*Cx2c*cCc.F
@test norm(Dx2c-Diagonal(Dx2c))+1. ≈ 1.
@test cCc.F'*Cc*cCc.F≈I
@test norm(Dx1c-(I-Dx2c))+1. ≈ 1.

# Method (2) real
c12=csp(X1, X2, simple=true)
Dx1=c12.F'*Cx1*c12.F
@test norm(Dx1-Diagonal(Dx1))+1≈1.
Dx2=c12.F'*Cx2*c12.F
@test norm(Dx2-Diagonal(Dx2))+1≈1.
@test c12.F'*C*c12.F≈I
@test norm(Dx1-(I-Dx2))+1≈1.
@test cC==c12

# Method (2) complex
c12c=csp(X1c, X2c, simple=true)
Dx1c=c12c.F'*Cx1c*c12c.F
@test norm(Dx1c-Diagonal(Dx1c))+1. ≈ 1.
Dx2c=c12c.F'*Cx2c*c12c.F
@test norm(Dx2c-Diagonal(Dx2c))+1. ≈ 1.
@test c12c.F'*Cc*c12c.F≈I
@test norm(Dx1c-(I-Dx2c))+1. ≈ 1.
@test cCc==c12c

# Method (3) real
# CSP of the average covariance matrices
k=10
Xset=[genDataMatrix(n, t) for i=1:k]
Yset=[genDataMatrix(n, t) for i=1:k]

c=csp(Xset, Yset)

# ... selecting subspace dimension allowing an explained variance = 0.9
c=csp(Xset, Yset; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
c=csp(Xset, Yset; meanX₁=nothing, meanX₂=nothing, eVar=0.9)

# csp on the average of the covariance and cross-covariance matrices
# computed along dims 1
c=csp(Xset, Yset; dims=1, eVar=0.9)

# name of the filter
c.name

using Plots
# plot regularized accumulated eigenvalues
plot(c.arev)

# plot the original covariance matrices and the transformed counterpart
# example when argument `selMeth` is `extremal` (default): 2-class separation
 cC=csp(Cx1, Cx2)
 Cx1Max=maximum(abs.(Cx1));
 h1 = heatmap(Cx1, clim=(-Cx1Max, Cx1Max), title="Cx1", yflip=true, c=:bluesreds);
 h2 = heatmap(cC.F'*Cx1*cC.F, clim=(0, 1), title="F'*Cx1*F", yflip=true, c=:amp);
 Cx2Max=maximum(abs.(Cx2));
 h3 = heatmap(Cx2, clim=(-Cx2Max, Cx2Max), title="Cx2", yflip=true, c=:bluesreds);
 h4 = heatmap(cC.F'*Cx2*cC.F, clim=(0, 1), title="F'*Cx2*F", yflip=true, c=:amp);
 CMax=maximum(abs.(C));
 h5 = heatmap(C, clim=(-CMax, CMax), title="Cx1+Cx2", yflip=true, c=:bluesreds);
 h6 = heatmap(cC.F'*C*cC.F, clim=(0, 1), title="F'*(Cx1+Cx2)*F", yflip=true, c=:amp);
 📈=plot(h1, h3, h5, h2, h4, h6, size=(800,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCSP1.png")
```

 ![Figure CSP1](assets/FigCSP1.png)

```
# example when argument `selMeth` is different from `extremal`: enhance snr
 cC=csp(Cx1, Cx2; selMeth=:enhaceSNR)
 Cx1Max=maximum(abs.(Cx1));
 h1 = heatmap(Cx1, clim=(-Cx1Max, Cx1Max), title="Cx1", yflip=true, c=:bluesreds);
 h2 = heatmap(cC.F'*Cx1*cC.F, clim=(0, 1), title="F'*Cx1*F", yflip=true, c=:amp);
 Cx2Max=maximum(abs.(Cx2));
 h3 = heatmap(Cx2, clim=(-Cx2Max, Cx2Max), title="Cx2", yflip=true, c=:bluesreds);
 h4 = heatmap(cC.F'*Cx2*cC.F, clim=(0, 1), title="F'*Cx2*F", yflip=true, c=:amp);
 CMax=maximum(abs.(C));
 h5 = heatmap(C, clim=(-CMax, CMax), title="Cx1+Cx2", yflip=true, c=:bluesreds);
 h6 = heatmap(cC.F'*C*cC.F, clim=(0, 1), title="F'*(Cx1+Cx2)*F", yflip=true, c=:amp);
 📉=plot(h1, h3, h5, h2, h4, h6, size=(800,400))
# savefig(📉, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCSP2.png")

```

 ![Figure CSP2](assets/FigCSP2.png)

```
# Method (3) complex
# CSP of the average covariance matrices
k=10
Xsetc=[genDataMatrix(ComplexF64, n, t) for i=1:k]
Ysetc=[genDataMatrix(ComplexF64, n, t) for i=1:k]

cc=csp(Xsetc, Ysetc)

# ... selecting subspace dimension allowing an explained variance = 0.9
cc=csp(Xsetc, Ysetc; eVar=0.9)

# ... subtracting the mean from the matrices in Xset and Yset
cc=csp(Xsetc, Ysetc; meanX₁=nothing, meanX₂=nothing, eVar=0.9)

# csp on the average of the covariance and cross-covariance matrices
# computed along dims 1
cc=csp(Xsetc, Ysetc; dims=1, eVar=0.9)

# name of the filter
cc.name
```
"""
function csp(C₁ :: SorH, C₂ :: SorH;
             eVar     :: TeVaro = ○,
             eVarC    :: TeVaro = ○,
             eVarMeth :: Function = searchsortedfirst,
             selMeth  :: Symbol = :extremal,
             simple   :: Bool = false)

  # checks
  size(C₁, 1)==size(C₁, 2) || throw(ArgumentError(📌*", csp function: Matrix `C₁` must be square"))
  size(C₂, 1)==size(C₂, 2) || throw(ArgumentError(📌*", csp function: Matrix `C₂` must be square"))
  size(C₁)==size(C₂) || throw(ArgumentError(📌*", csp function: Matrices `C₁` and `C₂` must have the same size"))
  eVar isa Int && eVarC isa Int && eVar>eVarC && throw(ArgumentError(📌*", csp function: `eVar` cannot be larger than `eVarC`"))

  args=("Common Spatial Pattern", false)

  if eVarC≠○ && eVarC≈0. # use gevd, which also actually whitens C1+C2
     λ, U = eig(C₁, C₁+C₂)
     λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

     simple ? LF(U, inv(U), Diagonal(λ), ○, ○, ○, args...) :
     begin
        eVar, D, U, p, arev=_ssdcsp!(eVar, λ, U, _minDim(C₁, C₂), eVarMeth, selMeth) # subspace dimension
        LF(U, pinv(U), D, eVar, λ, arev, args...)
     end
  else
     w=whitening(C₁+C₂; eVar=eVarC, eVarMeth=eVarMeth)

     # alert user if eVar passed as an integer exceeds the whitening dim
     whiteDim=size(w.F, 2)
     if eVar isa Int && eVar>whiteDim
        @warn(📌*", csp function: the whitening step reduced the rank to $(whiteDim); `eVar` has been lowered to this value.")
        eVar=whiteDim
     end

     λ, U = eig(Hermitian(w.F'*C₁*w.F)) # get evd of whitened C1
     # Hermitian is necessary for complex data
     λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

     simple ? LF(w.F*U, U'*w.iF, Diagonal(λ), ○, ○, ○, args...) :
     begin
        eVar, D, U, p, arev=_ssdcsp!(eVar, λ, U, _minDim(C₁, C₂), eVarMeth, selMeth) # subspace dimension
        LF(w.F*U, U'*w.iF, D, eVar, λ, arev, args...)
     end
  end
end

function csp(X₁ :: Mat, X₂ :: Mat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX₁   :: Tmean = 0,
             meanX₂   :: Tmean = 0,
             wX₁      :: Tw = ○,
             wX₂      :: Tw = ○,
          eVar     :: TeVaro = ○,
          eVarC    :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          selMeth  :: Symbol = :extremal,
          simple   :: Bool = false)

   dims===○ && (dims=_set_dims(X₁, X₂))
   _check_data(X₁, dims, covEst, meanX₁, wX₁)===○ && return
   _check_data(X₂, dims, covEst, meanX₂, wX₂)===○ && return

   C₁=_cov(X₁, covEst, dims, meanX₁, wX₁)
   C₂=_cov(X₂, covEst, dims, meanX₂, wX₂)

   csp(C₁, C₂; eVar=eVar, eVarC=eVarC, eVarMeth=eVarMeth,
       selMeth=selMeth, simple=simple)
end

function csp(𝐗₁::VecMat, 𝐗₂::VecMat;
             covEst   :: StatsBase.CovarianceEstimator = SCM,
             dims     :: Into = ○,
             meanX₁   :: Into = 0,
             meanX₂   :: Into = 0,
          eVar     :: TeVaro = ○,
          eVarC    :: TeVaro = ○,
          eVarMeth :: Function = searchsortedfirst,
          selMeth  :: Symbol = :extremal,
          simple   :: Bool = false,
       metric   :: Metric = Euclidean,
       w₁       :: Vector = [],
       w₂       :: Vector = [],
       ✓w       :: Bool = true,
       init₁    :: SorHo = nothing,
       init₂    :: SorHo = nothing,
       tol      :: Real = 0.,
       verbose  :: Bool = false)

   Metric==VonNeumann && throw(ArgumentError(📌*", csp function: A solution for the mean is not available for the Von Neumann metric. Use another metric as `metric` argument"))
   dims===○ && (dims=_set_dims(𝐗₁, 𝐗₂))
   _check_data(𝐗₁, dims, covEst, meanX₁, ○)===○ && return
   _check_data(𝐗₂, dims, covEst, meanX₂, ○)===○ && return

   𝐂₁= _cov(𝐗₁; covEst=covEst, dims=dims, meanX=meanX₁)
   𝐂₂= _cov(𝐗₂; covEst=covEst, dims=dims, meanX=meanX₂)

   csp(mean(metric, 𝐂₁;
            w = w₁, ✓w = ✓w,
            init = init₁===○ ? ○ : Hermitian(init₁), #just init₁ here when you upfate PodDefManifold
            tol = tol, verbose = verbose),
       mean(metric, 𝐂₂;
            w = w₂, ✓w = ✓w,
            init = init₂===○ ? ○ : Hermitian(init₂), #just init₂ here when you upfate PodDefManifold
            tol = tol, verbose = verbose);
       eVar=eVar, eVarC=eVarC, eVarMeth=eVarMeth,
       selMeth=selMeth, simple=simple)
end

"""
```
(1)
function cstp( X :: Mat, C₍₁₎ :: SorH, C₍₂₎ :: SorH;
               eVar     :: TeVaro = ○,
               eVarC    :: TeVaro = ○,
               eVarMeth :: Function = searchsortedfirst,
               simple   :: Bool = false)

(2)
function cstp( 𝐗::VecMat;
               covEst   :: StatsBase.CovarianceEstimator = SCM,
               meanXd₁  :: Into = 0,
               meanXd₂  :: Into = 0,
            eVar     :: TeVaro = ○,
            eVarC    :: TeVaro = ○,
            eVarMeth :: Function = searchsortedfirst,
            simple   :: Bool = false,
         metric   :: Metric = Euclidean,
         w        :: Vector = [],
         ✓w       :: Bool = true,
         init1    :: SorHo = nothing,
         init2    :: SorHo = nothing,
         tol      :: Real = 0.,
         verbose  :: Bool = false)
```

Return a [LinearFilter](@ref) object:

(1)
**Common spatio-temporal pattern**
with ``n⋅m`` mean data matrix `X`,
``m⋅m`` covariance matrices `C₍₁₎` and ``n⋅n`` covariance matrix `C₍₂₎` as input.

`eVar`, `eVarC` and `eVarMeth` are keyword optional arguments
for defining the [subspace dimension](@ref) ``p``. Particularly:
-  `eVarC` is used for defining the subspace dimension of
   the whitening step. The default is 0.999.
- `eVar` is the keyword optional argument for defining the
   [subspace dimension](@ref) ``p`` using the `.arev` vector
   given by [cstp.5]. The default is given in [cstp.6] here above.
- `eVarMeth` applies to both `eVarC` and `eVar`. The default value is
   `evarMeth=searchsortedfirst`.

If `simple` is set to `true`, ``p`` is set equal to ``n``
and only the fields `.F` and `.iF` are written in the constructed object.
This option is provided for low-level work when you don't need to define
a subspace dimension or you want to define it by your own methods.

(2)
**Common spatio-temporal pattern**
with a set of ``k`` data matrices `𝐗` as input.

The ``k`` matrices in `𝐗` are real or complex data matrices.
They must all have the same size.

`covEst`, `meanXd₁` and `meanXd₂` are optional
keyword arguments to regulate the estimation of the
covariance matrices of the data matrices in `𝐗`,
to be used to compute the mean covariance matrices in [cstp.2] here above.
See [covariance matrix estimations](@ref).
`meanXd₁` and `meanXd₂` are the means along dimension 1 and 2, respectively,
of the data matrices in `𝐗`.

The mean covariance matrices ``C_{(1)}`` and ``C_{(1)}`` in [cstp.2]
are computed using optional keywords arguments
`metric`, `w`, `✓w`, `init1`, `init2`, `tol` and `verbose`,
which allow to compute non-Euclidean means.
Particularly (see [mean covariance matrix estimations](@ref)),
- `init1` is the initialization for ``C_{(1)}``,
- `init2` is the initialization for ``C_{(2)}``.
By default, the arithmetic means [cstp.2] are computed.

Once the two covariance matrices ``C_{(1)}`` and ``C_{(2)}`` estimated,
method (1) is invoked with optional keyword arguments
`eVar`, `eVarC`, `eVarMeth` and `simple`.
See method (1) for details.

**Examples:**

```
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

# Method (1) real
t, n, k=10, 20, 10
Xset = [genDataMatrix(t, n) for i = 1:k]
Xfixed=randn(t, n)./1
for i=1:length(Xset) Xset[i]+=Xfixed end
C1=Hermitian( mean((X'*X)/t for X∈Xset) )
C2=Hermitian( mean((X*X')/n for X∈Xset) )
Xbar=mean(Xset)
c=cstp(Xbar, C1, C2; simple=true)
@test c.F[1]'*C2*c.F[1]≈I
@test c.F[2]'*C1*c.F[2]≈I
Z=c.F[1]'*Xbar*c.F[2]
n=minimum(size(Z))
@test norm(Z[1:n, 1:n]-Diagonal(Z[1:n, 1:n]))+1. ≈ 1.
cX=cstp(Xset; simple=true)
@test c==cX

# Method (1) complex
t, n, k=10, 20, 10
Xcset = [genDataMatrix(ComplexF64, t, n) for i = 1:k]
Xcfixed=randn(ComplexF64, t, n)./1
for i=1:length(Xcset) Xcset[i]+=Xcfixed end
C1c=Hermitian( mean((Xc'*Xc)/t for Xc∈Xcset) )
C2c=Hermitian( mean((Xc*Xc')/n for Xc∈Xcset) )
Xcbar=mean(Xcset)
cc=cstp(Xcbar, C1c, C2c; simple=true)
@test cc.F[1]'*C2c*cc.F[1]≈I
@test cc.F[2]'*C1c*cc.F[2]≈I
Zc=cc.F[1]'*Xcbar*cc.F[2]
n=minimum(size(Zc))
@test norm(Zc[1:n, 1:n]-Diagonal(Zc[1:n, 1:n]))+1. ≈ 1.
cXc=cstp(Xcset; simple=true)
@test cc==cXc

# Method (2) real
c=cstp(Xset)

# ... selecting subspace dimension allowing an explained variance = 0.9
c=cstp(Xset; eVar=0.9)

# ... giving weights `w` to the covariance matrices
c=cstp(Xset; w=abs2.(randn(k)), eVar=0.9)

# ... subtracting the means
c=cstp(Xset; meanXd₁=nothing, meanXd₂=nothing, w=abs2.(randn(k)), eVar=0.9)

# explained variance
c.eVar

# name of the filter
c.name

using Plots
# plot the original covariance matrices and the transformed counterpart
c=cstp(Xset)

C1Max=maximum(abs.(C1));
 h1 = heatmap(C1, clim=(-C1Max, C1Max), title="C1", yflip=true, c=:bluesreds);
 D1=c.F[1]'*C2*c.F[1];
 D1Max=maximum(abs.(D1));
 h2 = heatmap(D1, clim=(0, D1Max), title="F[1]'*C2*F[1]", yflip=true, c=:amp);
 C2Max=maximum(abs.(C2));
 h3 = heatmap(C2, clim=(-C2Max, C2Max), title="C2", yflip=true, c=:bluesreds);
 D2=c.F[2]'*C1*c.F[2];
 D2Max=maximum(abs.(D2));
 h4 = heatmap(D2, clim=(0, D2Max), title="F[2]'*C1*F[2]", yflip=true, c=:amp);

XbarMax=maximum(abs.(Xbar));
 h5 = heatmap(Xbar, clim=(-XbarMax, XbarMax), title="Xbar", yflip=true, c=:bluesreds);
 DX=c.F[1]'*Xbar*c.F[2];
 DXMax=maximum(abs.(DX));
 h6 = heatmap(DX, clim=(0, DXMax), title="F[1]'*Xbar*F[2]", yflip=true, c=:amp);
 📈=plot(h1, h3, h5, h2, h4, h6, size=(800,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigCSTP.png")

```

 ![Figure CSTP](assets/FigCSTP.png)

```

# Method (2) complex
cc=cstp(Xcset)

# ... selecting subspace dimension allowing an explained variance = 0.9
cc=cstp(Xcset; eVar=0.9)

# ... giving weights `w` to the covariance matrices
cc=cstp(Xcset; w=abs2.(randn(k)), eVar=0.9)

# ... subtracting the mean
cc=cstp(Xcset; meanXd₁=nothing, meanXd₂=nothing,
        w=abs2.(randn(k)), eVar=0.9)

# explained variance
c.eVar

# name of the filter
c.name

```

"""
function cstp( X :: Mat, C₍₁₎ :: SorH, C₍₂₎ :: SorH;
               eVar     :: TeVaro = ○,
               eVarC    :: TeVaro = ○,
               eVarMeth :: Function = searchsortedfirst,
               simple   :: Bool = false)

   d₍₂₎, d₍₁₎, dₓ=size(C₍₂₎, 1), size(C₍₁₎, 1), size(X)
   (d₍₁₎==dₓ[2] && d₍₂₎==dₓ[1]) || throw(ArgumentError(📌*", cstp function: For n⋅m matrix X, matrix C₍₁₎ must be m⋅m and matrix C₍₂₎ n⋅n"))
   eVar isa Int && eVarC isa Int && eVar>eVarC && throw(ArgumentError(📌*", cstp function: `eVar` cannot be larger than `eVarC`"))
   args=("Common Spatio-Temporal Pattern", false)

   kwargs=(eVar=eVarC, eVarMeth=eVarMeth, simple=false)
   t=whitening(C₍₁₎; kwargs...)
   s=whitening(C₍₂₎; kwargs...)

   # alert user if eVar passed as an integer exceeds the minimum whitening dim
   whiteDim=min(size(t.F, 2), size(t.F, 2))
   if eVar isa Int && eVar>whiteDim
      @warn(📌*", cstp function: the whitening step reduced the rank to $(whiteDim); `eVar` has been lowered to this value.")
      eVar=whiteDim
   end

   U, λ, V = svd(s.F'*X*t.F; full=true)
   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF([s.F*U, t.F*V], [U'*s.iF, V'*t.iF], Diagonal(λ), ○, ○, ○, args...) :
   begin
     eVar, D, U, V, p, arev=_ssdcstp!(eVar, λ, U, Matrix(V), _minDim(X), eVarMeth) # subspace dimension
     LF([s.F*U, t.F*V], [U'*s.iF, V'*t.iF], D, eVar, λ, arev, args...)
   end
end

# no dims argument since it is the same for cstp
function cstp( 𝐗::VecMat;
               covEst   :: StatsBase.CovarianceEstimator = SCM,
               meanXd₁  :: Into = 0,
               meanXd₂  :: Into = 0,
            eVar     :: TeVaro = ○,
            eVarC    :: TeVaro = ○,
            eVarMeth :: Function = searchsortedfirst,
            simple   :: Bool = false,
         metric   :: Metric = Euclidean,
         w        :: Vector = [],
         ✓w       :: Bool = true,
         init1     :: SorHo = nothing,
         init2     :: SorHo = nothing,
         tol      :: Real = 0.,
         verbose  :: Bool = false)

   Metric==VonNeumann && throw(ArgumentError(📌*", cstp function: A solution for the mean is not available for the Von Neumann metric. Use another metric as `metric` argument"))
   covEst==SCM && metric ∉ (Euclidean, Wasserstein) && throw(ArgumentError(📌*", cstp function: Only the Euclidean and Wasserstein `metric` can be used if the covariance estimator is `SCM`"))
   _check_data(𝐗, 1, covEst, meanXd₁, ○)===○ && return
   _check_data(𝐗, 2, covEst, meanXd₂, ○)===○ && return

   𝐂₍₁₎=_cov(𝐗; covEst=covEst, dims = 1, meanX = meanXd₁)
   𝐂₍₂₎=_cov(𝐗; covEst=covEst, dims = 2, meanX = meanXd₂)

   cstp(PosDefManifold.fVec(mean, Vector{Matrix}(𝐗)), # multi-threaded Euclidean mean
        mean(metric, 𝐂₍₁₎;
             w = w, ✓w = ✓w,
             init = init1===○ ? ○ : Hermitian(init1), #just init here when you update PodDefManifold
             tol = tol, verbose = verbose),
        mean(metric, 𝐂₍₂₎;
             w = w, ✓w = ✓w,
             init = init2===○ ? ○ : Hermitian(init2), #just init here when you update PodDefManifold.jl
             tol = tol, verbose = verbose),
        eVar     = eVar,
        eVarC    = eVarC,
        eVarMeth = eVarMeth,
        simple   = simple)
end


text:
using ADCME
using PyPlot
using ForwardDiff
using NiLang, NiLang.AD

# ADCME
function myloss_adcme(b; n=101)
    h = 1/(n-1)
    x = LinRange(0,1,n)[2:end-1]
    f = @. (4*(2 + x - x^2))

    u = trisolve(-b/h^2 * ones(n-3), 2b/h^2+1 * ones(n-2), -b/h^2 * ones(n-3), f)
    ue = u[div(n+1,2)] # extract values at x=0.5

    return (ue-1.0)^2
end

function benchmark_adcme(n)
    b = Variable(1.0)
    loss = myloss_adcme(b; n = n)
    g = gradients(loss, b)
    sess = Session(); init(sess)
    run(sess, g)
    ts_fwd = 0.0
    ts_bwd = 0.0
    for i = 1:11
        init(sess)
        d1 = @timed run(sess, loss)
        (i>1) && (ts_fwd += d1[2])
        init(sess)
        d2 = @timed run(sess, g)
        (i>1) && (ts_bwd += d2[2])
        @info "adcme", n, d1[2], d2[2]
    end
    ts_fwd /= 10
    ts_bwd /= 10
    return ts_fwd, ts_bwd
end

# ForwardDiff
function trisolve!(a::T, b::T, c::T, D::AbstractVector, X::AbstractVector) where T
    N = length(X)
    D = copy(D)
    B = zeros(T, N)
    @inbounds B[1] = b
    @inbounds for i = 2:N
        w = a / B[i-1]
        B[i] = b - w * c
        D[i] = D[i] - w * D[i - 1]
    end
    @inbounds X[N] = D[N] / B[N]
    @inbounds for i = N-1:-1:1
        X[i] = (D[i] - c * X[i + 1]) / B[i]
    end
    return X
end

function myloss(b::T; n=101) where T
    h = 1/(n-1)
    x = LinRange(0,1,n)[2:end-1]
    f = @. T(4*(2 + x - x^2))

    u = trisolve!(-b/h^2, 2b/h^2+1, -b/h^2, f, zeros(T, n-2))
    ue = u[div(n+1,2)] # extract values at x=0.5

    return (ue-1.0)^2
end

function benchmark_forwarddiff(n)
    ts_fwd = 0.0
    ts_bwd = 0.0
    for i = 1:11
        d1 = @timed myloss(10.0, n = n)
        (i>1) && (ts_fwd += d1[2])
        d2 = @timed myloss(ForwardDiff.Dual(10.0, 1.0),  n = n)
        (i>1) && (ts_bwd += d2[2])
        @info "forwarddiff", n, d1[2], d2[2]
    end
    ts_fwd /= 10
    ts_bwd /= 10
    ts_fwd, ts_bwd
end

# NiLang
@i function i_trisolve!(a::T, b::T, c::T, D!::AbstractVector, X!::AbstractVector, B!::AbstractVector) where T
    @invcheckoff @inbounds begin
    B![1] += b
    for i = 2:length(X!)
        @routine begin
            w ← zero(T)
            w += a / B![i-1]
        end
        B![i] += b
        B![i] -= w * c
        D![i] -= w * D![i - 1]
        ~@routine
    end
    X![end] += D![end] / B![end]
    for i = length(X!)-1:-1:1
        @routine begin
            anc ← zero(T)
            anc += D![i]
            anc -= c * X![i + 1]
        end
        X![i] += anc / B![i]
        ~@routine
    end
    end
end

@i function i_myloss!(loss::T, f!::AbstractVector{T}, u!::AbstractVector{T},
              b_cache!::AbstractVector{T}, b::T) where T
    @invcheckoff @inbounds begin
    n ← length(f!) + 2
    h ← zero(T)
    h += 1 / (n-1)
    for i=1:length(f!)
        @routine begin
            @zeros T xi anc
            xi += i * h
            anc += 2 + xi
            anc -= xi ^ 2
        end
        f![i] += 4 * anc
        ~@routine
    end

    @routine begin
        @zeros T h2 factor_a factor_b factor_c
        h2 += h^2
        factor_a -= b / h2
        factor_c += factor_a
        factor_b -= 2 * factor_a
        factor_b += 1
    end
    i_trisolve!(factor_a, factor_b, factor_c, f!, u!, b_cache!)
    ~@routine
    @routine begin
        ue ← zero(T)
        ue += u![div(n+1,2)] # extract values at x=0.5
        ue -= 1
    end
    loss += ue^2
    ~@routine
    h -= 1 / (n-1)
    end
end

function benchmark_nilang(n)
    ts_fwd = 0.0
    ts_bwd = 0.0
    for i = 1:11
        d1 = @timed i_myloss!(0.0, zeros(n-2), zeros(n-2), zeros(n-2), 10.0)
        (i>1) && (ts_fwd += d1[2])
        d2 = @timed Grad(i_myloss!)(Val(1), 0.0, zeros(n-2), zeros(n-2), zeros(n-2), 10.0)
        (i>1) && (ts_bwd += d2[2])
        @info "nilang", n, d1[2], d2[2]
    end
    ts_fwd /= 10
    ts_bwd /= 10
    ts_fwd, ts_bwd
end

tsf, tsb = zeros(3,6), zeros(3,6)
for (k,n) in enumerate([101, 1001,10001, 100001, 1000001, 10000001])
    tsf[1,k], tsb[1,k] = benchmark_adcme(n)
    tsf[2,k], tsb[2,k] = benchmark_forwarddiff(n)
    tsf[3,k], tsb[3,k] = benchmark_nilang(n)
end

figure(figsize=(10,4))
subplot(121)
plot([101, 1001,10001, 100001, 1000001, 10000001], tsf[1,:], "o-", label="ADCME")
plot([101, 1001,10001, 100001, 1000001, 10000001], tsf[2,:], "o-", label="Julia")
plot([101, 1001,10001, 100001, 1000001, 10000001], tsf[3,:], "o-", label="NiLang")
title("Forward Computation")
legend()
xlabel("\$n\$")
ylabel("Time (seconds)")
grid("on")
subplot(122)
plot([101, 1001,10001, 100001, 1000001, 10000001], tsb[1,:], "o-", label="ADCME")
plot([101, 1001,10001, 100001, 1000001, 10000001], tsb[2,:], "o-", label="ForwardDiff")
plot([101, 1001,10001, 100001, 1000001, 10000001], tsb[3,:], "o-", label="NiLang")
legend()
title("Gradient Back-Propagation")
xlabel("\$n\$")
ylabel("Time (seconds)")
grid("on")
savefig("benchmark.png")


text:
import StaticArrays
function similar_type(SAT::Type{<:StaticArrays.SVector},
                      size::StaticArrays.Size, T::Type)
    StaticArrays.similar_type(SAT, T, size)
end
# Called when mixing polyhedra with one using StaticArrays and the other not.
similar_type(::Type{<:StaticArrays.SVector}, ::Int, T::Type) = Vector{T}
similar_type(::Type{<:Vector}, ::StaticArrays.Size, T::Type) = Vector{T}

# The dimension is Int for Vector, SparseVector and
# StaticArrays.Size for StaticArrays.SVector
# This allows similar_type to be compiled as only one method for Vector
# and SparseVector and be type stable for StaticArrays.SVector
"""
    FullDim(p)::FullDim

Similar to [`fulldim`](@ref) but used for type stability with the vector type.
If the vector type is `StaticArrays.SVector` then it returns a
`StaticArrays.Size`.
"""
const FullDim = Union{Int, StaticArrays.Size}

FullDim(::Type{<:AbstractVector}) = -1 # Shouldn't hurt as it will not be used
FullDim(v::AbstractVector) = length(v)
FullDim(v::Union{StaticArrays.SVector{N}, Type{<:StaticArrays.SVector{N}}}) where N = StaticArrays.Size(v)

function FullDim_convert(::Type{StaticArrays.Size{N}}, d::Int) where N
    @assert N[1] == d
    return StaticArrays.Size{N}()
end
function FullDim_convert(::Type{StaticArrays.Size{N}},
                         d::StaticArrays.Size{N}) where N
    return d
end
FullDim_convert(::Type{Int}, d::Int) = d
function FullDim_convert(::Type{Int}, d::StaticArrays.Size{N}) where N
    return N[1]
end

"""
    fulldim(rep::Rep)::Int

Returns the dimension of the space in which polyhedron, representation,
element or vector is defined. That is, a straight line in a 3D space has
`fulldim` 3 even if its dimension is 1.
"""
fulldim(p) = fulldim(FullDim(p))

function fulldim end

fulldim(N::Integer) = convert(Int, N)
fulldim(::StaticArrays.Size{N}) where N = N[1]

neg_fulldim(N::Int) = -N
# FIXME May need generated function to make it type stable
neg_fulldim(::StaticArrays.Size{N}) where N = StaticArrays.Size{(-N[1],)}()

sum_fulldim(N1::Int, N2::Int) = N1 + N2
# FIXME May need generated function to make it type stable
sum_fulldim(::StaticArrays.Size{N1}, ::StaticArrays.Size{N2}) where {N1, N2} = StaticArrays.Size{(N1[1] + N2[1],)}()
# Called in cartesian product of two polyhedra, only the second using static arrays.
sum_fulldim(N1::Int, ::StaticArrays.Size{N2}) where N2 = N1 + N2[1]


text:
#=

A call to Resynthesizer using medium data:
large image and medium mask.

For testing mallocs.

Usage:

>julia
julia> include("medium.jl")
julia> @time test(image, mask)

OR

>julia --track-allocation=user
and the usual process for analyzing mallocs...
julia> include("medium.jl")
julia> test(image, mask)
julia>using Profile
julia>Profile.clear_malloc_data()
julia> test(image, mask)
=#

# For "load" local  files
using FileIO
# The tested item
include("../../src/apps/inpaint.jl")

# mock test data
# global so it persists and out of the test
image = load("/home/bootch/git/resynthesizerJulia/test/data/in/lighthouse.png")
# medium mask
mask = falses(size(image))
# 10x20 mask
mask[50:60,300:320] .= true

function test(image, mask)
    inpaint(image, mask)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using Unitful, UnitfulAstro

# Given
sellingPrice = 167 # Dollars
materialsAndParts = 25 # Dollars per unit
directLabor = 2 * 20
fixedCost = 1400000

# Find
# a.) overhead expenses are charged at
overheadExpenses = 0.8 * directLabor
totalManufacturingCost = materialsAndParts + directLabor + overheadExpenses

# b.) Breakeven
m = sellingPrice - totalManufacturingCost
b = -fixedCost
x = -b / m

# c.) Profit/unit
units = 30000
y = m * units + b
profitPerUnit = y / units

# d.)
x = 15000 # when y = 0, what should
P = -b / x + totalManufacturingCost


text:
families_registry = Dict{DataType, Vector{Family}}()

available_families(element::BIMShape) = available_families(element.family)
available_families(family::Family) = available_families(typeof(family))
available_families(family_type::DataType) = families_registry[family_type]

register_family(family::Family) = push!(get!(families_registry, typeof(family), Family[]), family)
export available_families, register_family

#=
We define standard family types
=#

for thickness in [11, 15, 30]
  register_family(
    wall_family_element(
      default_wall_family(),
      name="Wall_$(thickness)",
      thickness=thickness))
end

for thickness in [20, 22, 25, 30]
  register_family(
    slab_family_element(
      default_slab_family(),
      name="Slab_$(thickness)",
      thickness=thickness))
end


text:
using Pkgx
using Documenter

DocMeta.setdocmeta!(Pkgx, :DocTestSetup, :(using Pkgx); recursive=true)

makedocs(;
    modules=[Pkgx],
    authors="Thomas A. Christensen II <25492070+MillironX@users.noreply.github.com> and contributors",
    repo="https://github.com/MillironX/Pkgx.jl/blob/{commit}{path}#{line}",
    sitename="Pkgx.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MillironX.github.io/Pkgx.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MillironX/Pkgx.jl",
    devbranch="master",
)


text:
include("writeMat.jl")

nxs = 12
SP, X = PolynomialRing(Nemo.ZZ,["X$i" for i in 1:nxs])

############# M1.11
#############

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.0111.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[12] X[2]+X[4];
       X[1]+X[5] X[11];
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev11.1111.mpl")

############# M1.01
#############

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.0111.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[12] X[2]+X[4];
       0 0;
       X[1]+X[6] X[12]];

toMaple(net, X, "M1M2nonrev01.1111.mpl")

############# M1.10
#############

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.0111.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[11] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       X[1]+X[5] X[11];
       0 0];

toMaple(net, X, "M1M2nonrev10.1111.mpl")

############# M1.00
#############

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[1]+X[4] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       0 0;
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.0111.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1000.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1001.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1010.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[2]+X[4] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       0 0;
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1011.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1100.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[1]+X[6] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       0 0;
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1101.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[2]+X[6] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       0 0;
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1110.mpl")

net = [X[2] X[1];
       X[1] X[2];
       X[3] X[5];
       X[5] X[4];
       X[4] X[6];
       X[4] X[3];
       X[6] X[5];
       X[7] X[1]+X[3];
       X[8] X[2]+X[3];
       X[9] X[1]+X[5];
       X[10] X[2]+X[5];
       X[1]+X[4] X[7];
       X[2]+X[4] X[8];
       X[1]+X[6] X[9];
       X[2]+X[6] X[10];
       X[1]+X[5] X[2]+X[3];
       X[1]+X[6] X[2]+X[4];
       0 0;
       0 0];

toMaple(net, X, "M1M2nonrev00.1111.mpl")

max_stars_repo_name: LauraBMo/CRNTscripts

max_issues_repo_name: LauraBMo/CRNTscripts

max_forks_repo_name: LauraBMo/CRNTscripts


text:
module PicoConsole

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Pico Console" => Ahorn.EntityPlacement(
        Maple.PicoConsole
    )
)

sprite = "objects/pico8Console"

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

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.PicoConsole, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)

end


text:
using Pkg3
using Pkg3.Types
if Base.isdeprecated(Main, :Test)
    using Test
else
    using Base.Test
end

function temp_pkg_dir(fn::Function)
    local project_path
    try
        project_path = joinpath(tempdir(), randstring())
        withenv("JULIA_ENV" => project_path) do
            fn()
        end
    finally
        rm(project_path, recursive=true, force=true)
    end
end

# Tests for Example.jl fail on master,
# so let's use another small package
# in the meantime
const TEST_PKG = "Crayons"

temp_pkg_dir() do
    Pkg3.add(TEST_PKG; preview = true)
    @test_warn "not in project" Pkg3.API.rm("Example")
    Pkg3.add(TEST_PKG)
    @eval import $(Symbol(TEST_PKG))
    Pkg3.up()
    Pkg3.rm(TEST_PKG; preview = true)

    # TODO: Check coverage kwargs
    # TODO: Check that preview = true doesn't actually execute the test
    # by creating a package with a test file that fails.
    Pkg3.test(TEST_PKG)
    Pkg3.test(TEST_PKG; preview = true)

    Pkg3.rm(TEST_PKG)

    try
        Pkg3.add([PackageSpec(TEST_PKG, VersionSpec(v"55"))])
    catch e
        @test contains(sprint(showerror, e), TEST_PKG)
    end

    nonexisting_pkg = randstring(14)
    @test_throws CommandError Pkg3.add(nonexisting_pkg)
    @test_throws CommandError Pkg3.up(nonexisting_pkg)
    @test_warn "not in project" Pkg3.rm(nonexisting_pkg)
end


text:
using Test
using Distributed
using Sockets
using Suppressor
using JSON
using Syslogs
using Memento
using Memento.TestUtils
using TimeZones
using Dates
using Serialization

using Base.CoreLogging: global_logger, min_enabled_level

files = [
    "records.jl",
    "formatters.jl",
    "handlers.jl",
    "loggers.jl",
    "io.jl",
    "concurrency.jl",
    "test.jl",
    "stdlib.jl",
    "ext/json.jl",
]

Sys.isunix() ? push!(files, "ext/syslogs.jl") : nothing
haskey(ENV, "MEMENTO_BENCHMARK") ? files = ["benchmarks.jl"] : nothing

# for loggers.jl
struct TestError <: Exception
    msg
end

# for records.jl
struct SimpleRecord{T} <: Record
    level::String
    msg::String
    date::T
end

SimpleRecord(level, msg) = SimpleRecord(level, msg, Dates.now(tz"UTC"))

struct ConstRecord <: AttributeRecord
end

function Base.getproperty(::ConstRecord, attr::Symbol)
    if attr === :level
        return "error"
    elseif attr === :msg
        return "It's a ConstRecord"
    else
        throw(KeyError(attr))
    end
end

struct SimplestHandler <: Handler{Union{}}
    buf::IOBuffer
end

struct FilterHandler <: Handler{Union{}}
    buf::IOBuffer
    filters::Vector{Memento.Filter}
end

FilterHandler(buf::IOBuffer) = FilterHandler(buf, Memento.Filter[])
Base.push!(handler::FilterHandler, filter::Memento.Filter) = push!(handler.filters, filter)
Memento.getfilters(handler::FilterHandler) = handler.filters

function Memento.emit(handler::Union{SimplestHandler, FilterHandler}, record::Record)
    println(handler.buf, record.msg)
end

struct AsyncHandler <: Handler{Union{}}
    buf::IOBuffer
    channel::Channel{Record}

    function AsyncHandler(buf::IOBuffer)
        channel = Channel{Record}(Inf)
        handler = new(buf, channel)

        task = @async process_record!(handler)
        bind(channel, task)

        return handler
    end
end

function process_record!(handler::AsyncHandler)
    while isopen(handler.channel)
        record = take!(handler.channel)
        println(handler.buf, record.msg)
    end
end

function Memento.emit(handler::AsyncHandler, record::Record)
    put!(handler.channel, record)
end

_props(cr::ConstRecord) = (:level => cr.level, :msg => cr.msg)

Base.iterate(cr::ConstRecord, state...) = iterate(_props(cr), state...)

@testset "Memento" begin

@testset "Logging" begin
    @testset "Sample Usage" begin
        Memento.config!("info"; fmt="[{date} | {level} | {name}]: {msg}", colorized=false)
        logger1 = getlogger(@__MODULE__)
        debug(logger1, "Something that won't get logged.")
        Memento.info(logger1, "Something you might want to know.")
        Memento.warn(logger1, "This might cause an error.")
        Memento.warn(logger1, ErrorException("A caught exception that we want to log as a warning."))
        @test_throws ErrorException error(logger1, "Something that should throw an error.")
        @test_throws ErrorException error(logger1, ErrorException("A caught exception that we should log and rethrow"))
        logger2 = getlogger("Pkg.Foo.Bar")
    end

    @testset "Logger Hierarchy" begin
        Memento.reset!()
        foo = getlogger("Foo")
        bar = getlogger("Foo.Bar")
        baz = getlogger("Foo.Bar.Baz")
        car = getlogger("Foo.Car")

        for l in (foo, bar, baz, car)
            @test isset(l)
            @test getlevel(l) == "info"
            @test length(gethandlers(l)) == 0
        end

        root_io = IOBuffer()
        baz_io = IOBuffer()

        push!(getlogger(), DefaultHandler(root_io, DefaultFormatter("{name} - {level}: {msg}")))

        msg = "This should propagate and log because all loggers are set to :warn by default."
        expected = "Foo.Car - warn: $msg"
        Memento.warn(car, msg)
        result = String(take!(root_io))
        @test occursin(expected, result)

        msg = "This should not log because debug messages don't have high enough importance."
        debug(bar, msg)
        result = String(take!(root_io))
        @test result == ""

        setlevel!(getlogger(), "debug")
        msg = "This should log to because we've set the root logger to info"
        debug(getlogger(), msg)
        result = String(take!(root_io))
        expected = "root - debug: $msg"
        @test occursin(expected, result)

        msg = "This should not log because `bar` is set too low"
        debug(bar, msg)
        result = String(take!(root_io))
        @test result == ""

        setlevel!(bar, "debug")
        msg = "This should still not log because `foo` is set too low"
        debug(bar, msg)
        result = String(take!(root_io))
        @test result == ""

        # Now if we have the chain of root, foo and bar all set to info then the log record
        # will propagate up from the bar logger to be written to the buffer by the root loggers handler.
        setlevel!(foo, "debug")
        msg = "Now this should log because `bar`, `foo` and the root logger are all set appropriately."
        debug(bar, msg)
        expected = "Foo.Bar - debug: $msg"
        result = String(take!(root_io))
        @test occursin(expected, result)

        # Now we'll test that a child logger can log messages to its handler even if it doesn't
        # propagate to the root logger.
        setlevel!(baz, "debug")
        push!(baz, DefaultHandler(baz_io, DefaultFormatter("{name} - {level}: {msg}")))

        # Reset the parent loggers
        setlevel!.((foo, bar), "info")

        msg = "This should log to `baz_io`, but not `root_io`"
        debug(baz, msg)
        expected = "Foo.Bar.Baz - debug: $msg"
        result = String(take!(baz_io))
        @test occursin(expected, result)

        result = String(take!(root_io))
        @test result == ""

        msg = "This should be logged to both `root_io` and `baz_io`"
        Memento.info(baz, msg)
        expected = "Foo.Bar.Baz - info: $msg"
        result = String(take!(baz_io))
        @test occursin(expected, result)

        result = String(take!(root_io))
        @test occursin(expected, result)
    end

    @testset "Default Levels" begin
        @test allunique(values(Memento._log_levels))
    end
end

# Test files should assume the const _loggers has been reset
for file in files
    Memento.reset!()
    include(file)
end

end


text:
"""
Trapezoidal integration of x Nx1 and Δt 1xM to obtain y NxM.
"""
function trapz(Δt,x)
    y = 1/2 * (x[:,2:end] .+ x[:,1:end-1]) .* Δt
    y = sum(y, dims=2)
end
"""
Trapezoidal cumulative integration of x Nx1 and Δt 1xM to obtain y NxM.
"""
function cumtrapz(Δt,x)
    y = 1/2 * (x[:,2:end] .+ x[:,1:end-1]) .* Δt
    y = cumsum(y, dims=2)
end