text:
function read_rates(path, opts)

    types = Dict(
        :Field => String,
        :Tank => String,
        :Reservoir => String,
        :Date => Date,
        :Qoil => Float64,
        :Qwat => Float64,
        :Qliq => Float64,
        :Qinj => Float64,
        :Pres => Float64,
        :Source_resp => String,
        :Pbhp_prod => Float64,
        :Pbhp_inj => Float64,
        :Source_bhp => String,
        :Wellwork => String,
        :Wcut => Float64,
        :Total_mobility => Float64,
        :Wres => Float64,
        :Wbhp_prod => Float64,
        :Wbhp_inj => Float64,
        :Pres_min => Float64,
        :Pres_max => Float64
    )
    # TODO: Forced column conversion to 'PooledArray' removed (albeit memory inefficient)
    opts_csv = Dict(
        :dateformat => opts["dateformat"],
        :delim => opts["delim"],
        :types => types,
        :pool => false
    )
    df = CSV.File(opencsv(path); opts_csv...) |> DataFrame

    # Removing weights for missed pressure measurements
    @with df begin
        @. :Wres[ismissing(:Pres)] = missing
        @. :Wbhp_prod[ismissing(:Pbhp_prod)] = missing
        @. :Wbhp_inj[ismissing(:Pbhp_inj)] = missing
    end

    # Sort by dates
    sort!(df, [:Field, :Tank, :Date])

    return df
end

function read_params(path, opts)

    types = Dict(
        :Field => String,
        :Tank => String,
        :Neighb => String,
        :Reservoir => String,
        :Parameter => String,
        :Date => Date,
        :Init_value => Float64,
        :Min_value => Float64,
        :Max_value => Float64,
        :alpha => Float64
    )
    # TODO: Forced column conversion to 'PooledArray' removed (albeit memory inefficient)
    opts_csv = Dict(
        :dateformat => opts["dateformat"],
        :delim => opts["delim"],
        :types => types,
        :pool => false
    )
    df = CSV.File(opencsv(path); opts_csv...) |> DataFrame

    # Fill in the missing values
    crit = ismissing.(df.Neighb)
    df[crit, :Neighb] .= df[crit, :Tank]
    # TODO: To convert 'Union{Missing, String}' to 'String'
    df.Neighb = convert.(String, df.Neighb)

    # Convert 'String' to 'Symbol'
    df.Parameter = Symbol.(df.Parameter)
    replace!(df.Parameter, PSYMS...)

    # Sort by dates
    @eachrow!(df, (:Tank, :Neighb) = sort!([:Tank, :Neighb]))
    sort!(df, [:Parameter, :Field, :Tank, :Neighb, :Date])

    return df
end

function mark_null_params!(df::AbstractDataFrame, name::Symbol, isnull)
    df_view = getparams(df, Val(name))
    df_view[!, :Ignore] .|= @with df_view @byrow begin
        rows = UnitRange(:Jstart::Int, :Jstop::Int)
        all(isnull[rows, :Istart::Int])
    end
    df_view[!, :Const] .|= df_view.Ignore
    return df
end

function set_null_weights!(weight, df::AbstractDataFrame, name::Symbol)
    df_view = getparams(df, Val(name), Val(:const))
    @with df_view @byrow begin
        rows = UnitRange(:Jstart::Int, :Jstop::Int)
        # TODO: Multiplication selected for 'missing propagation',
        # i.e. (missing * number => missing)
        weight[rows, :Istart::Int] .*= zero(eltype(weight))
    end
    return weight
end

function process_params!(df_params::AbstractDataFrame, df_rates::AbstractDataFrame)

    # Numbering tank names and dates
    numberof(data) = @_ data |>
                    unique(__) |>
                    Dict(__ .=> 1:length(__))
    tanks = numberof(df_rates.Tank::Vector{String})
    dates = numberof(df_rates.Date::Vector{Date})
    N = length(dates)

    # Add the tank number and the beginning of the parameter period
    @transform! df_params begin
        :Istart = getindex.(Ref(tanks), :Tank::Vector{String})
        :Jstart = getindex.(Ref(dates), :Date::Vector{Date})
    end

    # For connection parameters, set the connection number
    df = getparams(df_params, Val(:Tconn))
    @with df begin
        connections = numberof(eachrow([:Tank :Neighb]::Matrix{String}))
        :Istart .= getindex.(Ref(connections), eachrow([:Tank :Neighb]))
    end

    # Adding the end of the parameter period
    df = groupby(df_params, [:Parameter, :Field, :Tank, :Neighb])
    @transform!(df, :Jstop = [:Jstart[2:end] .- 1; N])

    # TODO: Additionally, we split (if required) the intervals of fixing 'Jinj',
    # so that within each interval the value of 'λ' would not change.
    # This makes it much easier to differentiate the objective function with respect to 'P'
    split_params!(df_params, :λ, :Jinj)

    # Marking parameters that do not require fitting
    @transform!(df_params, :Const = :Min_value .== :Max_value)
    @transform!(df_params, :Ignore = :Const .& (:Init_value .≠ :Min_value))

    # Some parameters can be immediately excluded from the fitting
    crits = @with df_rates begin (
        λ = (:Qinj .== 0)::BitVector,
        Jp = (@. (ismissing(:Pbhp_prod)
                | ismissing(:Wbhp_prod)
                | (:Wbhp_prod == 0)
                | (:Qliq == 0)))::BitVector,
        Jinj = (@. (ismissing(:Pbhp_inj)
                | ismissing(:Wbhp_inj)
                | (:Wbhp_inj == 0)
                | (:Qinj == 0)))::BitVector,
    ) end
    for (key, values) ∈ pairs(crits)
        mark_null_params!(df_params, key, reshape(values, N, :))
    end

    # Forcing weights of measurements for ignored parameters to zero (Ignore == true)
    crits = @with df_rates begin (
        Jp = :Wbhp_prod,
        Jinj = :Wbhp_inj,
    ) end
    for (key, values) ∈ pairs(crits)
        set_null_weights!(reshape(values, N, :), df_params, key)
    end

    # Instead of well's PI we add well's geometric factor to fit
    crits = @with df_rates begin (
        Jp = :Total_mobility,
    ) end
    for (key, values) ∈ pairs(crits)
        add_geom_factor!(df_params, key, reshape(values, N, :))
    end

    return df_params, df_rates
end

function save_rates!(df::AbstractDataFrame, prob::NonlinearProblem, path, opts)

    @transform! df begin
        :Pres_calc = vec(prob.params.Pcalc')
        :Pbhp_calc = vec(prob.params.Pbhp')
        :Pinj_calc = vec(prob.params.Pinj')
        :Geom_factor = vec(prob.params.Gw')
        :Prod_Index = vec(prob.params.Jp')
        :Inj_index = vec(prob.params.Jinj')
        :Frac_inj = vec(prob.params.λ')
        :Qliq_calc = vec(prob.params.Qliq')
        :Qinj_calc = vec((prob.params.Qinj .* prob.params.λ)')
    end

    open(path, CSV_ENC, "w") do io
        CSV.write(io, df; delim=opts["delim"], dateformat=opts["dateformat"])
    end

    return df
end

function save_params!(df::AbstractDataFrame, fset::FittingSet, targ::TargetFunction, path, opts)

    @transform!(df, :Calc_value = :Init_value)

    # Fill in tank parameters
    df_view = getparams(df, Val(:tanks), Val(:var))
    copyto!(df_view.Calc_value, fset.cache.ybuf)

    # Fill in the parameters of the wells
    df_view = getparams(df, Val(:Gw), Val(:var))
    copyto!(df_view.Calc_value, targ.wviews.Gw.yviews)
    df_view = getparams(df, Val(:Jinj), Val(:var))
    copyto!(df_view.Calc_value, targ.wviews.Jinj.yviews)
    df_view = getparams(df, Val(:Jp), Val(:var))
    copyto!(df_view.Calc_value, targ.wviews.Jp.yviews)

    # Rename and remove extra columns
    replace!(df.Parameter, reverse.(PSYMS)...)
    select!(df, Not([:Istart, :Jstart, :Jstop, :Link, :Const, :Ignore]))

    open(path, CSV_ENC, "w") do io
        CSV.write(io, df; delim=opts["delim"], dateformat=opts["dateformat"])
    end

    return df
end

function split_params!(df::AbstractDataFrame, src::Symbol, dst::Symbol)

    df_src = getparams(df, Val(src))
    df_dst = getparams(df, Val(dst))

    df_new = map(eachrow(df_dst)) do dfr
        # Interval intersection criterion
        crit = @with df_src begin
            @. (dfr.Istart == :Istart) &
                (
                    ((dfr.Jstart >= :Jstart) &
                        (dfr.Jstart <= :Jstop))
                    |
                    ((dfr.Jstart <= :Jstart) &
                        (dfr.Jstop >= :Jstart))
                )
        end
        # Split based on intersection criterion
        @transform! df_src[crit, :] begin
            # Duplicate information partially
            :Parameter = dfr.Parameter
            :Init_value = dfr.Init_value
            :Min_value = dfr.Min_value
            :Max_value = dfr.Max_value
            :alpha = dfr.alpha
            # Correcting the beginning and end of the fixing period
            :Date = [dfr.Date; :Date[2:end]]
            :Jstart = [dfr.Jstart; :Jstart[2:end]]
            :Jstop = [:Jstop[1:end-1]; dfr.Jstop]
            # Link to source
            :Link = axes(df_src, 1)[crit]
        end
    end |> items -> vcat(items...)

    # Forming a new list of fixing intervals
    deleteat!(df, df.Parameter .=== dst)
    append!(df, df_new; cols=:union)

    return df
end

function add_geom_factor!(df::AbstractDataFrame, name::Symbol, mobt)

    # Copy PI to calculate the geometric factor
    df_new = df[df.Parameter .=== name, :]

    # Calculate well geometric factors
    @eachrow! df_new begin
        # Fluid mobility over time
        rows = UnitRange(:Jstart::Int, :Jstop::Int)
        mobt_ = @view mobt[rows, :Istart::Int]

        # We do mobility equals to one if constant PI is required
        :Const && (mobt_ .= one(eltype(mobt)))

        # Current, min and max mobilities
        Minit = first(mobt_)
        Mmin, Mmax = extrema(mobt_)

        # Mapping PI => geometric factor
        :Init_value /= Minit
        :Min_value /= Mmin
        :Max_value /= Mmax

        # Handle the case when Gmin > Gmax
        if :Max_value < :Min_value
            # Correct min/max mobilities that is Gmin == Gmax
            Mmin *= √(:Min_value / :Max_value)
            Mmax *= √(:Max_value / :Min_value)
            clamp!(mobt_, Mmin, Mmax)
            # Calculate new Gmin, Gmax
            :Min_value = :Max_value = √(:Min_value * :Max_value)
        end

        # Putting values within admissible limits
        :Init_value = clamp(:Init_value, :Min_value, :Max_value)
        :Const |= :Min_value .== :Max_value
    end

    # Removing old information (if presented)
    deleteat!(df, df.Parameter .=== :Gw)
    # Add well geometric factors
    df_new[!, :Parameter] .= :Gw
    append!(df, df_new)

    return df
end


text:
module GLVisualize

using ..GLAbstraction
using Makie: RaymarchAlgorithm, IsoValue, Absorption, MaximumIntensityProjection, AbsorptionRGBA, IndexedAbsorptionRGBA

using ..GLMakie.GLFW
using ModernGL
using StaticArrays
using GeometryBasics
using Colors
using Makie
using FixedPointNumbers
using FileIO
using Markdown
using Observables

import Base: merge, convert, show
using Base.Iterators: Repeated, repeated
using LinearAlgebra

import Makie: to_font, glyph_uv_width!
import ..GLMakie: get_texture!, loadshader

const GLBoundingBox = FRect3D

include("visualize_interface.jl")
export visualize # Visualize an object

include(joinpath("visualize", "lines.jl"))
include(joinpath("visualize", "image_like.jl"))
include(joinpath("visualize", "mesh.jl"))
include(joinpath("visualize", "particles.jl"))
include(joinpath("visualize", "surface.jl"))

export CIRCLE, RECTANGLE, ROUNDED_RECTANGLE, DISTANCEFIELD, TRIANGLE

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end # module

ext: jl


text:
using TightlyBound
using Plots

#setup chosen Plots backend, i recommend pyplot
#pyplot()
#gr()

#make the crystal object
#here we choose fcc Al

types=["Al"];

#positions, crystal units
pos=zeros((1,3));

#lattice vectors, in Bohr units currently
A=[ [3.8 3.8 0]; [3.8 0 3.8 ]; [ 0 3.8 3.8]] * 0.529177;

#makes the crystal
c=makecrys(A, pos, types)

println("The crystal is")
println(c)

# Do scf calculation for the energy
# returns energy, tight-binding-crystal object, and a flag (true is good)
energy, tbc, flag = scf_energy(c)

println()
println("The energy is $energy")
println()
println(tbc)

#plot band structure
kpath=[0.0 0.0 0.0;0.5 0 0; 0.5 0.25 0.75; 0.0 0.0 0.0]
knames=["Γ", "X", "W", "Γ"]

plot_bandstr(tbc, kpath=kpath, names=knames)

savefig("al_fcc.pdf")


text:
@testset "KuhnPokerEnv" begin

    env = KuhnPokerEnv()
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
"""
        BoundaryFunction{B, X1, X2}(func, parameters=nothing)

A wrapper for the user-defined boundary condition function `func`, on the
boundary specified by symbol `B` and at location `(X1, X2)`, and with `parameters`.

Example
=======

julia> using Oceananigans, Oceananigans.BoundaryConditions, Oceananigans.Fields

julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
(::BoundaryFunction{:z,Cell,Cell,var"#7#8",Nothing}) (generic function with 1 method)

julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);

julia> flux_func(x, y, t, p) = cos(p.k * x) * cos(p.ω * t); # function with parameters

julia> parameterized_u_velocity_flux = BoundaryFunction{:z, Face, Cell}(flux_func, (k=4π, ω=3.0))
(::BoundaryFunction{:z,Face,Cell,typeof(flux_func),NamedTuple{(:k, :ω),Tuple{Float64,Float64}}}) (generic function with 1 method)

julia> top_u_bc = BoundaryCondition(Flux, parameterized_u_velocity_flux);
"""
struct BoundaryFunction{B, X1, X2, F, P} <: Function
    func :: F
    parameters :: P

    function BoundaryFunction{B, X1, X2}(func, parameters=nothing) where {B, X1, X2}
        B ∈ (:x, :y, :z) || throw(ArgumentError("The boundary B at which the BoundaryFunction is
                                                to be applied must be either :x, :y, or :z."))
        return new{B, X1, X2, typeof(func), typeof(parameters)}(func, parameters)
    end
end

@inline call_boundary_function(func, ξ, η, t, ::Nothing) = func(ξ, η, t)
@inline call_boundary_function(func, ξ, η, t, parameters) = func(ξ, η, t, parameters)

@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, clock, state) where {Y, Z} =
    call_boundary_function(bc.func, ynode(Y, j, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, clock, state) where {X, Z} =
    call_boundary_function(bc.func, xnode(X, i, grid), znode(Z, k, grid), clock.time, bc.parameters)

@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, clock, state) where {X, Y} =
    call_boundary_function(bc.func, xnode(X, i, grid), ynode(Y, j, grid), clock.time, bc.parameters)

#####
##### Convenience constructors
#####

""" Returns the location of `loc` on the boundary `:x`, `:y`, or `:z`. """
boundary_location(::Val{:x}, loc) = (loc[2], loc[3])
boundary_location(::Val{:y}, loc) = (loc[1], loc[3])
boundary_location(::Val{:z}, loc) = (loc[1], loc[2])

""" Returns a boundary function at on the boundary `B` at the appropriate
    location for a tracer field. """
TracerBoundaryFunction(B, args...) = BoundaryFunction{B, Cell, Cell}(args...)

""" Returns a boundary function at on the boundary `B` at the appropriate
    location for u, the x-velocity field. """
function UVelocityBoundaryFunction(B, args...)
    loc = boundary_location(Val(B), (Face, Cell, Cell))
    return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

""" Returns a boundary function at on the boundary `B` at the appropriate
    location for v, the y-velocity field. """
function VVelocityBoundaryFunction(B, args...)
    loc = boundary_location(Val(B), (Cell, Face, Cell))
    return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

""" Returns a boundary function at on the boundary `B` at the appropriate
    location for w, the z-velocity field. """
function WVelocityBoundaryFunction(B, args...)
    loc = boundary_location(Val(B), (Cell, Cell, Face))
    return BoundaryFunction{B, loc[1], loc[2]}(args...)
end

"""
    TracerBoundaryCondition(bctype, B, args...)

Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to a tracer on the boundary `B`, which is one of `:x, :y, :z`.
The boundary function has the signature

    `func(ξ, η, t)`

where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
TracerBoundaryCondition(bctype, B, args...) =
    BoundaryCondition(bctype, TracerBoundaryFunction(B, args...))

"""
    UVelocityBoundaryCondition(bctype, B, args...)

Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `u`, the `x`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature

    `func(ξ, η, t)`

where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
UVelocityBoundaryCondition(bctype, B, args...) =
    BoundaryCondition(bctype, UVelocityBoundaryFunction(B, args...))

"""
    VVelocityBoundaryCondition(bctype, B, args...)

Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `v`, the `y`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature

    `func(ξ, η, t)`

where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
VVelocityBoundaryCondition(bctype, B, args...) =
    BoundaryCondition(bctype, VVelocityBoundaryFunction(B, args...))

"""
    VVelocityBoundaryCondition(bctype, B, args...)

Returns a `BoundaryCondition` of type `bctype`, that applies the function
`func` to `w`, the `z`-velocity field, on the boundary `B`, which is one
of `:x, :y, :z`. The boundary function has the signature

    `func(ξ, η, t)`

where `t` is time, and `ξ` and `η` are coordinates along the
boundary, eg: `(y, z)` for `B = :x`, `(x, z)` for `B = :y`, or
`(x, y)` for `B = :z`.
"""
WVelocityBoundaryCondition(bctype, B, args...) =
    BoundaryCondition(bctype, WVelocityBoundaryFunction(B, args...))


text:
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")

# 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")

# 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


text:
test_ex(x, y) = @test prettify(x) == prettify(y)
test_ex(x::AbstractArray, y::AbstractArray) = foreach(test_ex, x, y)
test_ex(x::WrapperNode, y) = @test prettify(to_expr(x)) == prettify(y)

test(T, f, name, ex) = test_ex(T(f(name)), ex)
test(T, arg, ex) = test(T, identity, arg, ex)

@testset "Wrapper" begin
    @testset "Wrapper utilities" begin
        include("wrapper/utilities/exprs.jl")
        include("wrapper/utilities/naming_conventions.jl")
    end

    @testset "Wrapper generation" begin
        include("wrapper/config.jl")
    end

    @testset "Code generation" begin
        @testset "Structs" begin
            include("wrapper/codegen/structs/low_level.jl")
            include("wrapper/codegen/structs/high_level.jl")
            include("wrapper/codegen/structs/unions.jl")
        end

        include("wrapper/codegen/constants.jl")
        include("wrapper/codegen/enums/enums.jl")
        include("wrapper/codegen/enums/bitmasks.jl")
        include("wrapper/codegen/handles.jl")

        @testset "Functions" begin
            include("wrapper/codegen/functions/low_level.jl")
            include("wrapper/codegen/functions/high_level.jl")
            include("wrapper/codegen/functions/overloads.jl")
        end

        include("wrapper/codegen/aliases.jl")
        include("wrapper/codegen/docs.jl")
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Documenter, ProjB

makedocs(;
    modules=[ProjB],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/ProjB.jl/blob/{commit}{path}#L{line}",
    sitename="ProjB.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/tkf/ProjB.jl",
)


text:
############################################################################################
## Make model
#=
kernel = mcmckernel_pf[1]
references = references_pf[1]
=#

@testset "ParticleMetropolis" begin
    obj = deepcopy(myobjective)
    for kernel in mcmckernel_pf
        mcmcdefault = MCMCDefault(; stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
        mcmckernel = MCMC(
            _rng,
            kernel,
            myobjective_mcmc,
            mcmcdefault
        )
        for references in references_pf
            pfdefault = ParticleFilterDefault(referencing = references,)
            pfkernel = ParticleFilter(
                _rng,
                myobjective_pf,
                pfdefault
            )
            ## Check constructor
            pmcmc_c = ParticleMetropolis(ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, (:μ, :σ, :p), mcmcdefault))
            @test keys(obj.model.val) == BaytesPMCMC.get_sym(pmcmc_c)

            mcmcsym = :μ
            PMCMC(ParticleMetropolis, ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, mcmcsym, mcmcdefault))
            pmc = ParticleMetropolis(ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, mcmcsym, mcmcdefault))
            BaytesPMCMC.get_sym(pmc)
            ## Initialize pmcmc kernel
            pmcmckernel = PMCMC(_rng, ParticleMetropolis, pfkernel, mcmckernel)
            _vals, _diag = propose!(_rng, pmcmckernel, obj.model, obj.data)
            ## Check if constructor can be used
            pfconstructor = BaytesFilters.ParticleFilterConstructor(
                :latent,
                ParticleFilterDefault(referencing = references,)
            )
            mcmcconstructor = BaytesMCMC.MCMCConstructor(
                kernel,
                keys(myobjective_mcmc.tagged.parameter),
                MCMCDefault(; stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            )
            constructor = PMCMCConstructor(
                ParticleMetropolis,
                pfconstructor,
                mcmcconstructor,
                PMCMCDefault()
            )
            constructor(_rng, obj.model, obj.data, 1.0, SampleDefault())
            results([_diag], pmcmckernel, 2, [.1, .2, .5, .8, .9])
            BaytesPMCMC.result!(pmcmckernel, BaytesPMCMC.get_result(pmcmckernel))
            generate_showvalues(_diag)()

            ## Postprocessing
            @test _diag isa infer(_rng, BaytesPMCMC.AbstractDiagnostics, pmcmckernel, obj.model, obj.data)
            @test _diag.base.prediction isa infer(_rng, pmcmckernel, obj.model, obj.data)
            @test !( predict(_rng, pmcmckernel.kernel.mcmc, obj) isa infer(_rng, pmcmckernel, obj.model, obj.data) )
            @test predict(_rng, pmcmckernel.kernel.pf, obj) isa infer(_rng, pmcmckernel, obj.model, obj.data)
            @test predict(_rng, pmcmckernel.kernel, obj) isa typeof(predict(_rng, pmcmckernel, obj))
            @test _diag.base.ℓobjective == _diag.pf.base.ℓobjective

        end
    end
end

@testset "ParticleGibbs" begin
    obj = deepcopy(myobjective)
    for kernel in mcmckernel_pfa
        mcmcdefault = MCMCDefault(; stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
        mcmckernel = MCMC(
            _rng,
            kernel,
            myobjective_mcmc,
            mcmcdefault
        )
        for references in references_pfa
            pfdefault = ParticleFilterDefault(referencing = references,)
            pfkernel = ParticleFilter(
                _rng,
                myobjective_pf,
                pfdefault
            )
            ## Check constructor
            pmcmc_c = ParticleGibbs(ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, (:μ, :σ, :p), mcmcdefault))
            @test keys(obj.model.val) == BaytesPMCMC.get_sym(pmcmc_c)

            mcmcsym = :μ
            PMCMC(ParticleGibbs, ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, mcmcsym, mcmcdefault))
            pmc = ParticleGibbs(ParticleFilterConstructor(:latent, pfdefault), MCMCConstructor(kernel, mcmcsym, mcmcdefault))
            BaytesPMCMC.get_sym(pmc)

            ## Initialize pmcmc kernel
            pmcmckernel = PMCMC(_rng, ParticleGibbs, pfkernel, mcmckernel)
            _vals, _diag = propose!(_rng, pmcmckernel, obj.model, obj.data)
            ## Check if constructor can be used
            pfconstructor = BaytesFilters.ParticleFilterConstructor(
                :latent,
                ParticleFilterDefault(referencing = references,)
            )
            mcmcconstructor = BaytesMCMC.MCMCConstructor(
                kernel,
                keys(myobjective_mcmc.tagged.parameter),
                MCMCDefault(; stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            )
            constructor = PMCMCConstructor(
                ParticleGibbs,
                pfconstructor,
                mcmcconstructor,
                PMCMCDefault()
            )
            constructor(_rng, obj.model, obj.data, 1., SampleDefault())
            results([_diag], pmcmckernel, 2, [.1, .2, .5, .8, .9])
            BaytesPMCMC.result!(pmcmckernel, BaytesPMCMC.get_result(pmcmckernel))
            generate_showvalues(_diag)()

            ## Postprocessing
            @test _diag isa infer(_rng, BaytesPMCMC.AbstractDiagnostics, pmcmckernel, obj.model, obj.data)
            @test _diag.base.prediction isa infer(_rng, pmcmckernel, obj.model, obj.data)
            @test !( predict(_rng, pmcmckernel.kernel.mcmc, obj) isa infer(_rng, pmcmckernel, obj.model, obj.data) )
            @test predict(_rng, pmcmckernel.kernel.pf, obj) isa infer(_rng, pmcmckernel, obj.model, obj.data)
            @test predict(_rng, pmcmckernel.kernel, obj) isa typeof(predict(_rng, pmcmckernel, obj))
            @test _diag.base.ℓobjective == _diag.pf.base.ℓobjective
        end
    end
end


text:
module IOCapture
using Logging

"""
    IOCapture.capture(f; rethrow=Any, color=false)

Runs the function `f` and captures the `stdout` and `stderr` outputs without printing them
in the terminal. Returns an object with the following fields:

* `.value :: Any`: return value of the function, or the error exception object on error
* `.output :: String`: captured `stdout` and `stderr`
* `.error :: Bool`: set to `true` if `f` threw an error, `false` otherwise
* `.backtrace :: Vector`: array with the backtrace of the error if an error was thrown

The behaviour can be customized with the following keyword arguments:

* `rethrow`:

  When set to `Any` (default), `capture` will rethrow any exceptions thrown
  by evaluating `f`.

  To only throw on a subset of possible exceptions pass the exception type
  instead, such as `InterruptException`. If multiple exception types may need
  to be thrown then pass a `Union{...}` of the types. Setting it to `Union{}`
  will capture all thrown exceptions. Captured exceptions will be returned via
  the `.value` field, and will also set `.error` and `.backtrace` accordingly.

* `color`: if set to `true`, `capture` inherits the `:color` property of `stdout` and
  `stderr`, which specifies whether ANSI color/escape codes are expected. This argument is
  only effective on Julia v1.6 and later.

# Extended help

`capture` works by temporarily redirecting the standard output and error streams
(`stdout` and `stderr`) using `redirect_stdout` and `redirect_stderr` to a temporary
buffer, evaluate the function `f` and then restores the streams. Both the captured text
output and the returned object get captured and returned:

```jldoctest
julia> c = IOCapture.capture() do
           println("test")
       end;

julia> c.output
"test\\n"
```

This approach does have some limitations -- see the README for more information.

**Exceptions.** Normally, if `f` throws an exception, `capture` simply re-throws it with
`rethrow`. However, by setting `rethrow` to `false`, it is also possible to capture
errors, which then get returned via the `.value` field. Additionally, `.error` is set to
`true`, to indicate that the function did not run normally, and the `catch_backtrace` of the
exception is returned via `.backtrace`.

As mentioned above, it is also possible to set `rethrow` to `InterruptException`. This will
make `capture` rethrow only `InterruptException`s. This is useful when you want to capture
all the exceptions, but allow the user to interrupt the running code with `Ctrl+C`.

**Recommended pattern.** The recommended way to refer to `capture` is by fully qualifying
the function name with `IOCapture.capture`. This is also why the package does not export
the function. However, if a shorter name is desired, we recommend renaming the function
when importing:

```julia
using IOcapture: capture as iocapture
```

This avoids the function name being too generic.
"""
function capture(f; rethrow::Type=Any, color::Bool=false)
    # Original implementation from Documenter.jl (MIT license)
    # Save the default output streams.
    default_stdout = stdout
    default_stderr = stderr

    # Redirect both the `stdout` and `stderr` streams to a single `Pipe` object.
    pipe = Pipe()
    Base.link_pipe!(pipe; reader_supports_async = true, writer_supports_async = true)
    @static if VERSION >= v"1.6.0-DEV.481" # https://github.com/JuliaLang/julia/pull/36688
        pe_stdout = IOContext(pipe.in, :color => get(stdout, :color, false) & color)
        pe_stderr = IOContext(pipe.in, :color => get(stderr, :color, false) & color)
    else
        pe_stdout = pipe.in
        pe_stderr = pipe.in
    end
    redirect_stdout(pe_stdout)
    redirect_stderr(pe_stderr)
    # Also redirect logging stream to the same pipe
    logger = ConsoleLogger(pe_stderr)

    # Bytes written to the `pipe` are captured in `output` and eventually converted to a
    # `String`. We need to use an asynchronous task to continously tranfer bytes from the
    # pipe to `output` in order to avoid the buffer filling up and stalling write() calls in
    # user code.
    output = IOBuffer()
    buffer_redirect_task = @async while !eof(pipe)
        write(output, readavailable(pipe))
    end

    # Run the function `f`, capturing all output that it might have generated.
    # Success signals whether the function `f` did or did not throw an exception.
    result, success, backtrace = with_logger(logger) do
        try
            f(), true, Vector{Ptr{Cvoid}}()
        catch err
            err isa rethrow && Base.rethrow(err)
            # If we're capturing the error, we return the error object as the value.
            err, false, catch_backtrace()
        finally
            # Force at least a single write to `pipe`, otherwise `readavailable` blocks.
            println()
            # Restore the original output streams.
            redirect_stdout(default_stdout)
            redirect_stderr(default_stderr)
            close(pipe)
            wait(buffer_redirect_task)
        end
    end
    (
        value = result,
        output = chomp(String(take!(output))),
        error = !success,
        backtrace = backtrace,
    )
end

end


text:
using Test

using CompScienceMeshes
using BEAST

const e0 = point(0.0,0.0,0.0)
const e1 = point(1.0,0.0,0.0)
const e2 = point(0.0,1.0,0.0)
const e3 = point(0.0,0.0,1.0)

p = simplex(
    [
        point(0.0,0.0,0.0),
        point(1.0,0.0,0.0)
    ],
    Val{1}
)

q = simplex(
    [
        point(0.0,0.0,0.0),
        point(0.5,0.0,0.0)
    ],
    Val{1}
)

f = BEAST.LagrangeRefSpace{Float64,1,2,2}()
x = neighborhood(p, [0.0])
v = f(x)

@test v[1].value == 0
@test v[2].value == 1

@test v[1].derivative == -1
@test v[2].derivative == +1

Q = restrict(f, p, q)
@test Q == [
    1.0 0.5
    0.0 0.5]

# Test restriction of RT elements
ni, no = 6, 7;
ui = transpose([triangleGaussA[ni] triangleGaussB[ni] ]);
uo = transpose([triangleGaussA[no] triangleGaussB[no] ]);
wi = triangleGaussW[ni];
wo = triangleGaussW[no];
# universe = Universe(1.0, ui, wi, uo, wo);

p = simplex([2*e0,2*e1,2*e2], Val{2})
x = neighborhood(p,[0.5, 0.5])
ϕ = BEAST.RTRefSpace{Float64}()
v = ϕ(x)

Q = restrict(ϕ, p, p)
@test Q == Matrix(LinearAlgebra.I, 3, 3)

q = simplex([e0+e1, 2*e1, e1+e2], Val{2})
Q = restrict(ϕ, p, q)
@test Q == [
    2 -1 0
    0  1 0
    0 -1 2] // 4


text:
using Oceananigans, Printf

include("utils.jl")

function makeplot(axs, model)
    sca(axs[1])
    plotxzslice(model.tracers.S, cmap="RdBu_r")
    title("Passive tracer")

    sca(axs[2])
    plotxzslice(model.tracers.T, cmap="RdBu_r")
    title("Buoyancy")

    axs[1].axis("off"); axs[1].set_aspect(1)
    axs[2].axis("off"); axs[2].set_aspect(1)
    axs[1].tick_params(left=false, labelleft=false, bottom=false, labelbottom=false)
    axs[2].tick_params(left=false, labelleft=false, bottom=false, labelbottom=false)
    tight_layout()

    return nothing
end

#
# Parameters
#

      α = 2                 # aspect ratio
      n = 1                 # resolution multiple
     Ra = 1e6               # Rayleigh number
Nx = Ny = 16n * α           # horizontal resolution
Lx = Ly = 1.0 * α           # horizontal extent
     Nz = 16n               # vertical resolution
     Lz = 1.0               # vertical extent
     Pr = 0.7               # Prandtl number
      a = 1e-1              # noise amplitude for initial condition
     Δb = 1.0               # buoyancy differential

# Rayleigh and Prandtl determine transport coefficients
ν = sqrt(Δb * Pr * Lz^3 / Ra)
κ = ν / Pr

#
# Model setup
#

arch = CPU()
@hascuda arch = GPU() # use GPU if it's available

model = Model(
     arch = arch,
        N = (Nx, Ny, Nz),
        L = (Lx, Ly, Lz),
        ν = ν,
        κ = κ,
      eos = LinearEquationOfState(βT=1., βS=0.),
constants = PlanetaryConstants(g=1., f=0.)
)

# Constant buoyancy boundary conditions on "temperature"
model.boundary_conditions.T.z.top = BoundaryCondition(Value, 0.0)
model.boundary_conditions.T.z.bottom = BoundaryCondition(Value, Δb)

# Force salinity as a passive tracer (βS=0)
S★(x, z) = exp(4z) * sin(2π/Lx * x)
FS(grid, u, v, w, T, S, i, j, k) = 1/10 * (S★(grid.xC[i], grid.zC[k]) - S[i, j, k])
model.forcing = Forcing(FS=FS)

ArrayType = typeof(model.velocities.u.data)
Δt = 0.01 * min(model.grid.Δx, model.grid.Δy, model.grid.Δz)^2 / ν

#prefix = @sprintf("rayleigh_benard_Ra%d_", Ra)
#nc_writer = NetCDFOutputWriter(dir=".", prefix=prefix, frequency=100)
#push!(model.output_writers, nc_writer)

#
# Initial condition setup for creating regression test data
#

ξ(z) = a * rand() * z * (Lz + z) # noise, damped at the walls
b₀(x, y, z) = (ξ(z) - z) / Lz

x, y, z = model.grid.xC, model.grid.yC, model.grid.zC
x, y, z = reshape(x, Nx, 1, 1), reshape(y, 1, Ny, 1), reshape(z, 1, 1, Nz)

model.tracers.T.data .= ArrayType(b₀.(x, y, z))

fig, axs = subplots(nrows=2)

@printf("""
    Crunching Rayleigh-Benard convection with

        N : %d, %d, %d
       Ra : %.0e

    Let's spin the gears.

""", Nx, Ny, Nz, Ra)

for i = 1:100
    walltime = @elapsed time_step!(model, 1000, Δt)

    makeplot(axs, model)
    gcf()

    wb = buoyancy_flux(model)
    @printf("i: %d, t: %.2e, CFL: %.4f, Nuʷᵇ: %.3f, wall: %s\n", model.clock.iteration,
                model.clock.time, cfl(Δt, model), 1 + mean(wb.data) * Lz^2 / (κ*Δb),
                prettytime(1e9*walltime))
end


text:
try
  @eval using PackageCompiler
catch e
  using Pkg
  Pkg.add("PackageCompiler")
  @eval using PackageCompiler
end
sysim_dir = "$(ENV["HOME"])/.julia/sysimages"
!isdir(sysim_dir) && mkdir(sysim_dir)
using PackageCompiler
create_sysimage(:ITensors,
                sysimage_path="$sysim_dir/sys_itensors.so",
                precompile_execution_file="precompile_itensors.jl")


text:
import HorizonSideRobots: move!
import HorizonSideRobots: isborder
import HorizonSideRobots: putmarker!
import HorizonSideRobots: temperature
import HorizonSideRobots: ismarker

HS = [Nord, Ost, Sud, West]
mvs = []

mutable struct Coord
    x::Int
    y::Int
end

Coord() = Coord(0,0)

struct XYRobot
    robot::Robot
    coord::Coord
    XYRobot(r,(x,y))=new(r,Coord(x,y))
    XYRobot(r)=new(r,Coord())
end

function move!(robot::XYRobot, side::HorizonSide)
    move!(robot.robot, side)
    move!(robot.coord, side)
    println(robot.coord.x, " ", robot.coord.y)
end

isborder(robot::XYRobot,  side::HorizonSide) = isborder(robot.robot, side)
putmarker!(robot::XYRobot) = putmarker!(robot.robot)
ismarker(robot::XYRobot) = ismarker(robot.robot)
temperature(robot::XYRobot) = temperature(robot.robot)
get_coord(robot::XYRobot) = get_coord(robot.coord)

function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side)+2,4))
end

function move!(coord::Coord, side::HorizonSide)
    if side==Nord
        coord.y += 1
    elseif side==Sud
        coord.y -= 1
    elseif side==Ost
        coord.x += 1
    elseif side==West
        coord.x -= 1
    end
end

get_coord(coord::Coord) = (coord.x, coord.y)

function go_to_corner_Sud_West!(r::XYRobot)
    while !isborder(r, Sud) || !isborder(r, West)
        while !isborder(r, Sud)
            move!(r, Sud)
        end
        while !isborder(r, West)
            move!(r, West)
        end
    end
end

function move_if_coord_notzero(r::XYRobot)
    for i in [1, 2, 3, 4]
        while (r.coord.x!=0 && r.coord.y!=0)
            move!(r, HS[i])
        end
        if r.coord.x!=0 || r.coord.y!=0
            putmarker!(r)
        end
        while !isborder(r, HS[i])
            move!(r, HS[i])
        end
    end
end

function go_to_corner!(r::XYRobot, mvs::AbstractArray)
    while isborder(r, Sud) == false || isborder(r, West) == false
        moveIfPossible!(r, mvs, Sud)
        moveIfPossible!(r, mvs, West)
        println(mvs)
    end
end

function moveIfPossible!(r::XYRobot, mvs::AbstractArray, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
        push!(mvs, side)
    end
end

function return_back!(r::XYRobot, mvs::AbstractArray)
    println(mvs)
    while length(mvs) > 0
        last_move = pop!(mvs)
        move!(r, inverse(last_move))
    end
end

function main(r::Robot)
    r=XYRobot(r, (0, 0))
    go_to_corner!(r, mvs)
    move_if_coord_notzero(r)
    return_back!(r, mvs)
end

(1, 2, 3, 4)


text:
struct Bias{T}
    array::Vector{Vector{T}}
    function Bias{T}(bias_vec_of_vec::Vector{Vector{T}}) where {T}
        return new(copy(bias_vec_of_vec))
    end
end

function Bias(bias_vec_of_vec::Vector{Vector{T}}) where {T}
    return Bias{T}(bias_vec_of_vec)
end

function Bias(n_neuronns_per_layers::Vector{Int}, T = Float64)
    # for layers = 2..N_layers
    N_layers = length(n_neuronns_per_layers)
    bias_vec_of_vec = Vector{Vector{T}}(undef, N_layers-1)
    for ll = 1:N_layers-1
        tmp = randn(T, n_neuronns_per_layers[ll+1])
        bias_vec_of_vec[ll] = tmp
    end
    return Bias(bias_vec_of_vec)
end

function subarray(bias::Bias; layer = 1)
    return bias.array[layer]
end

struct WeightMatrix{T}
    # n_row = 1..n_neurons in the next layer
    # n_col = 1..n_neurons in the prev layer
    array::Vector{Matrix{T}}

    function WeightMatrix{T}(weightmat_vec::Vector{Matrix{T}}) where {T}
        return new(weightmat_vec)
    end
end

function WeightMatrix(weightmat_vec::Vector{Matrix{T}}) where {T}
    return WeightMatrix{T}(weightmat_vec)
end

function WeightMatrix(n_neuronns_per_layers::Vector{Int}, T = Float64)
    # for between layers = 1..N_layers
    N_layers   = length(n_neuronns_per_layers)
    weight_mat = Vector{Matrix{T}}(undef, N_layers-1)
    for ll = 1:N_layers-1
        n_prev         = n_neuronns_per_layers[ll]   # n_col
        n_next         = n_neuronns_per_layers[ll+1] # n_row
        tmp            = randn(T, n_next, n_prev)
        weight_mat[ll] = tmp
    end

    return WeightMatrix(weight_mat)
end

function submatrix(weight_mat::WeightMatrix; layer = 1)
    return weight_mat.array[layer]
end

function zeros!(weight_mat::WeightMatrix)
    L = length(weight_mat.array)
    for ii = 1:L
        v = weight_mat.array[ii]
        T = eltype(v)
        for jj in eachindex(v)
            v[jj] = zero(T)
        end
    end
end

struct IntermediateVectorResult{T}
    array::Vector{Vector{T}}
    function IntermediateVectorResult{T}(vec_of_vec::Vector{Vector{T}}) where {T}
        return new(copy(vec_of_vec))
    end
end

function IntermediateVectorResult(vec_of_vec::Vector{Vector{T}}) where {T}
    return IntermediateVectorResult{T}(vec_of_vec)
end

function IntermediateVectorResult(n_neuronns_per_layers::Vector{Int}, T = Float64)
    # for layers = 2..N_layers
    N_layers   = length(n_neuronns_per_layers)
    vec_of_vec = Vector{Vector{T}}(undef, N_layers-1)
    for ll = 1:N_layers-1
        # tmp            = Vector{T}(undef, n_neuronns_per_layers[ll+1])
        tmp            = zeros(T, n_neuronns_per_layers[ll+1])
        vec_of_vec[ll] = tmp
    end
    return IntermediateVectorResult(vec_of_vec)
end

function subarray(res::IntermediateVectorResult; layer = 1)
    return res.array[layer]
end

function zeros!(res::IntermediateVectorResult)
    L = length(res.array)
    for ii = 1:L
        v = res.array[ii]
        T = eltype(v)
        for jj in eachindex(v)
            v[jj] = zero(T)
        end
    end
end

struct PerformanceStat{T}
    perf::Vector{T}
    min::T
    max::T
    mean::T
    sigma::T
    function PerformanceStat{T}(perf::Vector{T}) where {T}
        min   = minimum(perf)
        max   = maximum(perf)
        s     = sum(perf)
        N     = length(perf)
        mean  = s / N
        v     = 1/(N-1) * sum( (perf .- mean).^2 )
        sigma = sqrt(v)
        return new(perf, min, max, mean, sigma)
    end
end

function PerformanceStat(perf::Vector{T}) where {T}
    return PerformanceStat{T}(perf)
end

function Base.show(io::IO, perfstat::PerformanceStat)
    println("Performance stat:")
    println("    - min:   $(perfstat.min)")
    println("    - mean:  $(perfstat.mean)")
    println("    - max:   $(perfstat.max)")
    println("    - sigma: $(perfstat.sigma)")

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module PkgServer

using Pkg
using HTTP
using Base.Threads: Event, @spawn
import Base: fetch
using Random
using LibGit2
using FilesystemDatastructures
using SimpleBufferStream
using JSON3, StructTypes
using Sockets
using Sockets: InetAddr
using Dates
using Tar
using Gzip_jll

include("task_utils.jl")
include("resource.jl")
include("meta.jl")
include("dynamic.jl")
mutable struct RegistryMeta
    # Upstream registry URL (e.g. "https://github.com/JuliaRegistries/General")
    upstream_url::String
    # The hash (for each flavor) for each registry
    hashes::Dict{String,String}

    function RegistryMeta(url::String)
        # Check to ensure this path actually exists
        if !url_exists(url)
            throw(ArgumentError("Invalid unreachable registry '$(url)'"))
        end

        # Auto-detect a repository that doesn't have `.git` at the end but could
        git_url = string(url, ".git")
        if !endswith(url, ".git") && url_exists(git_url)
            url = git_url
        end
        return new(url, Dict{String,String}())
    end
end

struct ServerConfig
    root::String
    listen_addr::InetAddr
    listen_server::Ref{Any}
    cache::SizeConstrainedFileCache
    registries::Dict{String,RegistryMeta}
    storage_servers::Vector{String}
    dotflavors::Vector{String}

    # Default server config constructor
    function ServerConfig(; listen_addr = InetAddr(ip"127.0.0.1", 8000),
                            storage_root = "/tmp/pkgserver",
                            registries = Dict(
                                "23338594-aafe-5451-b93e-139f81909106" =>
                                RegistryMeta("https://github.com/JuliaRegistries/General")
                            ),
                            storage_servers = [
                                "https://us-east.storage.juliahub.com",
                                "https://kr.storage.juliahub.com",
                            ],
                            dotflavors = [
                                ".eager",
                                ".conservative",
                            ],
                            keep_free=3*1024^3)
        # Right now, the only thing we store in `static/` is `/registries`
        mkpath(joinpath(storage_root, "static"))
        # Downloads get stored into `temp`
        mkpath(joinpath(storage_root, "temp"))
        # Files get stored into `cache`
        mkpath(joinpath(storage_root, "cache"))
        return new(
            storage_root,
            listen_addr,
            Ref(nothing),
            SizeConstrainedFileCache(
                joinpath(storage_root, "cache"),
                TargetSizeKeepFree(keep_free),
                DiscardLRU()
            ),
            registries,
            sort!(storage_servers),
            dotflavors,
        )
    end
end

# Initialize a default config that will get overridden by `start()` below
global config = ServerConfig()

function __init__()
    # Set default HTTP useragent
    HTTP.setuseragent!("PkgServer (HTTP.jl)")

    # Record our starting time
    global time_start = now()
end

function start(;kwargs...)
    @info("Loading initial Server Config")
    # Create global config object
    global config = ServerConfig(;kwargs...)
    flavorless_mode = config.dotflavors == [""]

    # Update registries first thing
    @info("Performing initial registry update")
    initial_update_changed = any(update_registries.(config.dotflavors))
    if !flavorless_mode
        # Remove any old `registries` files from when we were a flavorless server
        rm(joinpath(config.root, "static", "registries"); force=true)

        if !initial_update_changed
            @warn("Flavorless storage servers detected, falling back to flavorless mode")
            flavorless_mode = true
            empty!(config.dotflavors)
            push!(config.dotflavors, "")
        end
    else
        # Remove any old `registries.$(flavor)` files from when we were a flavorfull server
        for f in readdir(joinpath(config.root, "static"); join=true)
            if startswith(f, "registries.")
                rm(f; force=true)
            end
        end
    end
    if !initial_update_changed && !any(update_registries.(config.dotflavors))
        error("Unable to get initial registry update!")
    end
    global last_registry_update = now()

    # Experimental.@sync throws if _any_ of the tasks fail
    Base.Experimental.@sync begin
        global registry_update_task = @spawn begin
            while true
                sleep(1)
                @try_printerror begin
                    forget_failures()
                    update_registries.(config.dotflavors)
                    last_registry_update = now()
                end
            end
        end

        # Registry watchdog; if `last_registry_update` doesn't change
        # for more than 20 minutes, we kill ourselves so that we can restart
        max_time_lag = Second(20 * 60)
        global registry_watchdog_task = @spawn begin
            while true
                time_lag = now() - last_registry_update
                if time_lag > max_time_lag
                    task_result = try fetch(registry_update_task); catch err; err end
                    @error "registry update watchdog timer tripped" time_lag max_time_lag task_result
                    exit(1)
                end
                sleep(max_time_lag.value)
            end
        end

        listen_server = Sockets.listen(config.listen_addr)
        config.listen_server[] = listen_server
        @info("server listening", config.listen_addr)
        num_requests = 0
        HTTP.listen(config.listen_addr.host, config.listen_addr.port; server=listen_server, max_connections=10) do http
            num_requests += 1
            resource = http.message.target
            request_id = HTTP.header(http, "X-Request-ID", "")
            # If the user is asking for `/meta`, generate the requisite JSON object and send it back
            if resource == "/meta"
                serve_meta(http)
                return
            end
            if resource == "/meta/stats"
                serve_meta_stats(http)
                return
            end
            if resource == "/meta/siblings"
                serve_siblings(http)
                return
            end
            if resource == "/meta/children"
                # PkgServers don't have children, but to be good neighbors,
                # we return an empty list here, to make it easier to recurse.
                serve_children(http)
                return
            end
            if resource == "/meta/parents"
                serve_parents(http)
                return
            end
            if resource == "/robots.txt"
                serve_robots_txt(http)
                return
            end

            if resource == "/registries" && !flavorless_mode
                # If they're asking for just "/registries", inspect headers to figure
                # out which registry flavor they actually want, and if none is given,
                # give them `conservative` by default, unless they are self-reporting
                # as a CI bot, in which case we'll always point them to `eager`.
                ci = any([v == "t" for (k, v) in filter(!isempty, split(HTTP.header(http, "Julia-CI-Variables", ""), ";"))])
                flavor = HTTP.header(http, "Julia-Registry-Preference", ci ? "eager" : "conservative")
                serve_redirect(http, "/registries.$(flavor)")
                return
            end

            if occursin(r"^/registries(\.[a-z]+$)?", resource)
                filepath = joinpath(config.root, "static", basename(resource))
                if !isfile(filepath)
                    HTTP.setstatus(http, 404)
                    startwrite(http)
                    return
                end
                open(filepath) do io
                    serve_file(http, io, "text/plain")
                end
                return
            end

            # If the user asked for something that is an actual resource, send it directly
            if occursin(resource_re, resource)
                # If the resource already exists locally, yay!  Serve it and quit.
                resource_path = resource_filepath(resource)
                io = try_open(resource_path)
                if io !== nothing
                    hit!(config.cache, resource[2:end])
                    serve_file(http, io, "application/x-gzip")
                    close(io)
                    return
                end

                # If it doesn't exist locally, let's request a fetch on that resource.
                # This will return either `nothing` (e.g. resource does not exist) or
                # a `DownloadState` that represents a partial download.
                dl_state = fetch_resource(resource, request_id)
                if dl_state !== nothing
                    HTTP.setheader(http, "X-Cache-Miss" => "miss")
                    stream_path = temp_resource_filepath(resource)
                    # Wait until `stream_path` is created
                    while !isfile(stream_path) && !istaskdone(dl_state.dl_task)
                        sleep(0.001)
                    end
                    # Try to serve `stream_path` file
                    stream_io = try_open(stream_path)
                    if stream_io !== nothing
                        serve_file(http, stream_io, "application/x-gzip";
                                   content_length=dl_state.content_length,
                                   dl_task=dl_state.dl_task)
                        close(stream_io)
                        return
                    end

                    # If we couldn't open `stream_path`, it may be because the file finished
                    # downloading since we last checked 20 lines ago.  Check again.
                    io = try_open(resource_path)
                    if io !== nothing
                        serve_file(http, io, "application/x-gzip")
                        close(io)
                        return
                    end
                end
            end

            # If a user asks for an `Artifacts.toml` for a particular artifact, we'll generate
            # one for them on the fly.
            m = match(artifact_toml_re, resource)
            if m !== nothing
                artifact_resource = m.captures[1]
                artifact_name = m.captures[2]

                # If we don't actually have the artifact locally, we need to fetch it first.
                artifact_path = resource_filepath(artifact_resource)
                if !isfile(artifact_path)
                    dl_state = fetch_resource(artifact_resource, request_id)
                    if dl_state !== nothing
                        # We actually need to wait for the download to finish so that we can
                        # calculate the SHA256 hash of the tarball.
                        wait(dl_state.dl_task)
                    end
                end

                # We check again, because it's possible the above download failed.
                if isfile(artifact_path)
                    serve_artifact_toml(http, artifact_path, artifact_name)
                    return
                end
            end

            if occursin(r"^/*$", resource)
                path = joinpath(dirname(@__DIR__), "static", "index.html")
                io = try_open(path)
                if io !== nothing
                    serve_file(http, io, "text/html")
                    close(io)
                    return
                end
            end

            global total_misses += 1
            HTTP.setstatus(http, 404)
            startwrite(http)
        end
    end
end

# precompilation
include(joinpath(dirname(@__DIR__), "deps", "precompile.jl"))
if get(ENV, "PKGSERVER_GENERATING_PRECOMPILE", nothing) === nothing
    _precompile_()
end

end # module

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


text:
include("test_1d.jl")
include("test_2d.jl")
include("test_3d.jl")
include("test_4d.jl")
# include("test_5d.jl")

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
## Derivatives
function wst_S1_deriv(image::Array{Float64,2}, filter_hash::Dict, FFTthreads::Int=1)

    # Use nthread threads for FFT -- but for Nx<512 nthread=1 is fastest.  Overhead?
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate output array
    dS1dα  = zeros(Float64, Nx, Nx, Nf)

    ## 1st Order
    im_fd = fft(image)

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]
    zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P = plan_ifft(im_fd)  # P is an operator, P*im is ifft(im)

    # Loop over filters
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i

        zarr[f_i] = (f_v.*f_v) .* im_fd[f_i]
        dS1dα[:,:,f] = 2 .* real.(P*zarr)

        zarr[f_i] .= 0   # reset zarr for next loop
    end
    return dS1dα
end

function wst_S20_deriv(image::Array{Float64,2}, filter_hash::Dict, FFTthreads::Int=1)

    # Use nthread threads for FFT -- but for Nx<512 nthread=1 is fastest.  Overhead?
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate output array
    dS20dα  = zeros(Float64, Nx, Nx, Nf, Nf)

    # allocate image arrays for internal use
    im_rdc = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
    im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
    im_fd  = fft(image)

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]
    zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

    # make a FFTW "plan" a complex array, both forward and inverse transform
    P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
    P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

    # Loop over filters
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i

        zarr[f_i] = f_v .* im_fd[f_i]
        Z_λ = P_ifft*zarr  # complex valued ifft of zarr
        zarr[f_i] .= 0   # reset zarr for next loop
        im_rdc[:,:,f] = Z_λ
        im_rd[:,:,f]  = abs.(Z_λ)
    end

    zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
    for f2 = 1:Nf
        f_i = f_ind[f2]  # CartesianIndex list for filter
        f_v = f_val[f2]  # Values for f_i
        uvec = im_rdc[:,:,f2] ./ im_rd[:,:,f2]
        for f1 = 1:Nf
            temp = P_fft*(im_rd[:,:,f1].*uvec)
            zarr[f_i] = f_v .* temp[f_i]

            Z1dZ2 = real.(P_ifft*zarr)
            #  It is possible to do this with rifft, but it is not much faster...
            #   Z1dZ2 = myrealifft(zarr)
            dS20dα[:,:,f1,f2] += Z1dZ2
            dS20dα[:,:,f2,f1] += Z1dZ2
            zarr[f_i] .= 0   # reset zarr for next loop
        end
    end
    return dS20dα
end

function wst_S20_deriv_sum(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}, FFTthreads::Int=1)
    # Sum over (f1,f2) filter pairs for S20 derivative.  This is much faster
    #   than calling wst_S20_deriv() because the sum can be moved inside the FFT.
    # Use FFTthreads threads for FFT -- but for Nx<512 FFTthreads=1 is fastest.  Overhead?
    # On Cascade Lake box, 4 is good for 2D, 8 or 16 for 3D FFTs
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    (Nf, )    = size(filter_hash["filt_index"])

    # allocate image arrays for internal use
    Uvec   = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
    im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
    im_fd  = fft(image)

    # unpack filter_hash
    f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_hash["filt_value"]
    zarr    = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

    # make a FFTW "plan" a complex array, both forward and inverse transform
    P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
    P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

    # Loop over filters
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i

        zarr[f_i] = f_v .* im_fd[f_i]
        Z_λ = P_ifft*zarr  # complex valued ifft of zarr
        zarr[f_i] .= 0     # reset zarr for next loop
        im_rd[:,:,f] = abs.(Z_λ)
        Uvec[:,:,f]  = Z_λ ./ im_rd[:,:,f]
    end

    zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
    for f2 = 1:Nf
        f_i = f_ind[f2]  # CartesianIndex list for filter
        f_v = f_val[f2]  # Values for f_i

        Wtot = reshape( reshape(im_rd,Nx*Nx,Nf)*wt[:,f2], Nx, Nx)
        temp = P_fft*(Wtot.*Uvec[:,:,f2])
        zarr[f_i] .+= f_v .* temp[f_i]
    end
    ΣdS20dα = real.(P_ifft*zarr)

    return ΣdS20dα
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
include("types.jl")
include("functions.jl")
include("units.jl")
include("gamefuns.jl")
include("newfuns.jl")

function distance(loc1,loc2=(0,0,2))
	oloc=loc1.-loc2
	m=max(abs(oloc[1]),abs(oloc[2]))
	return max(m,abs(oloc[1]+oloc[2]))+oloc[3]
end
function save(game)
	dic=Dict()
	dic[:name]=game.name
	dic[:shells]=game.board.shells
	dic[:initlocs]=game.board.initlocs
	dic[:colind]=game.colind
	dic[:colmax]=game.colmax
	dic[:colock]=game.colock
	dic[:delete]=game.delete
	dic[:sequence]=[]
	for entry in game.sequence
		if isa(entry,Unit)
			unit=entry
		#elseif !isa(entry,Symbol) && length(entry)>1 && isa(entry[2],Unit)
		#	unit=entry[2]
		else
			push!(dic[:sequence],entry)
			continue
		end
		udic=Dict()
		udic[:name]=unit.name
		udic[:color]=unit.color
		udic[:loc]=unit.loc
		push!(dic[:sequence],udic)
	end
	dir=joinpath(homedir(),"weilianqi","saves",game.name)
	if !ispath(dir)
		touch(dir)
	end
	io=open(dir,"a+")
	write(io,string(dic),"\n")
	close(io)
	println("Saved $(length(game.units)) units at $dir")
end
function loadsequence!(game::Game,seq,originoffset=(0,0,0))
	for entry in seq
		if entry==:harvest
			harvest!(game)
		elseif entry[1]==:expand
			expandboard!(game,entry[2]...,false)
		else
			(loc,unit)=entry
			loco=loc.+originoffset
			placeunit!(game,unit)
		end
	end
	return "<3"
end
function loadic!(game,dic::Dict,originoffset=(0,0,0))
	#game.board=newboard(dic[:shells],dic[:initlocs]) #seeeegfault
	if haskey(dic,:colind)
		game.colind=dic[:colind]
	end
	game.color=game.colors[game.colind]
	if haskey(dic,:colmax)
		game.colmax=dic[:colmax]
	end
	if haskey(dic,:colock)
		game.colock=dic[:colock]
	end
	set_gtk_property!(game.gui[:colockcheck],:active,game.colock)
	if haskey(dic,:delete)
		game.delete=dic[:delete]
	end
	set_gtk_property!(game.gui[:deletecheck],:active,game.delete)
	if haskey(dic,:initlocs) #&& game.board.initlocs!=dic[:initlocs]
		pushfirst!(dic[:sequence],(:expand,[dic[:shells],dic[:initlocs]])) #let's hope this works
	end
	for entry in dic[:sequence]
		if isa(entry,Dict)
			if entry[:name]=="spawn";entry[:name]="queen";end
			unit=newunit(entry)
			unit.loc=unit.loc.+originoffset
			placeunit!(game,unit)
		elseif entry==:harvest
			harvest!(game)
		elseif entry[1]==:expand
			for ili in 1:length(entry[2][2])
				entry[2][2][ili]=entry[2][2][ili].+originoffset
			end
			expandboard!(game,entry[2][1],entry[2][2],false) #loadicing
		elseif entry[1]==:delete
			removeunit!(game,entry[2])
		end
	end
	sync!(game)
	return game
end
function loadic(dic::Dict,originoffset=(0,0,0))
	game=newgame(dic[:name],[dic[:shells],dic[:initlocs]],sequence=[])
	loadic!(game,dic,originoffset)
	return game
end
function string2game(lastlineofsavefile)
	dic=eval(parse(lastlineofsavefile))
	return loadic(dic)
end
function loadgame!(game::Game,name::String,originoffset=(0,0,0),backtrack::Integer=0)
	path=joinpath(homedir(),"weilianqi","saves",name)
	lines=readlines(path)
	dic=eval(Meta.parse(lines[end-backtrack]))
	loadic!(game,dic,originoffset)
	return game
end
function loadgame(name::String,backtrack::Integer=0)
	path=joinpath(homedir(),"weilianqi","saves",name)
	lines=readlines(path)
	dic=eval(Meta.parse(lines[end-backtrack]))
	if isa(dic,Dict)
		game=loadic(dic)
	else
		game=newgame(name)
		loadsequence!(game,dic)
	end
	sync!(game)
	return game
end


text:
function insertionSort(A)
    for j in 2:length(A)
        key = A[j]
        i = j-1
        while(i>0 && A[i]<key)
        A[i+1] = A[i]
        i = i-1
        end
        A[i+1] = key

    end
    return A
end


text:
#
# Special case for three tensors
#

function compute_cost(external_dims::Tuple{Int,Int,Int}, internal_dims::Tuple{Int,Int,Int})
  dim11, dim22, dim33 = external_dims
  dim12, dim23, dim31 = internal_dims
  cost12 = dim11 * dim22 * dim12 * dim23 * dim31
  return cost12 + dim11 * dim22 * dim33 * dim31 * dim23
end

function three_tensor_contraction_sequence(which_sequence::Int)::Vector{Any}
  @assert 1 ≤ which_sequence ≤ 3
  return if which_sequence == 1
    Any[3, [1, 2]]
  elseif which_sequence == 2
    Any[1, [2, 3]]
  else
    Any[2, [3, 1]]
  end
end

function optimal_contraction_sequence(is1, is2, is3)
  N1 = length(is1)
  N2 = length(is2)
  N3 = length(is3)
  dim2 = dim(is2)
  dim3 = dim(is3)
  dim11 = 1
  dim12 = 1
  dim31 = 1
  @inbounds for n1 in 1:N1
    i1 = is1[n1]
    n2 = findfirst(==(i1), is2)
    if isnothing(n2)
      n3 = findfirst(==(i1), is3)
      if isnothing(n3)
        dim11 *= dim(i1)
        continue
      end
      dim31 *= dim(i1)
      continue
    end
    dim12 *= dim(i1)
  end
  dim23 = 1
  @inbounds for n2 in 1:length(is2)
    i2 = is2[n2]
    n3 = findfirst(==(i2), is3)
    if !isnothing(n3)
      dim23 *= dim(i2)
    end
  end
  dim22 = dim2 ÷ (dim12 * dim23)
  dim33 = dim3 ÷ (dim23 * dim31)
  external_dims1 = (dim11, dim22, dim33)
  internal_dims1 = (dim12, dim23, dim31)
  external_dims2 = (dim22, dim33, dim11)
  internal_dims2 = (dim23, dim31, dim12)
  external_dims3 = (dim33, dim11, dim22)
  internal_dims3 = (dim31, dim12, dim23)
  cost1 = compute_cost(external_dims1, internal_dims1)
  cost2 = compute_cost(external_dims2, internal_dims2)
  cost3 = compute_cost(external_dims3, internal_dims3)
  mincost, which_sequence = findmin((cost1, cost2, cost3))
  sequence = three_tensor_contraction_sequence(which_sequence)
  return sequence
end


text:
#=
Challenge #252 [Easy] Sailors and monkeys and coconuts, oh my!

Description

A number of sailors (let's call it N) are stranded on an island with a huge pile
of coconuts and a monkey. During the night, each sailor (in turn) does the following
without the others knowing:
    1. He takes one N'th (e.g. if N=5, one fifth) of the coconuts in the pile and
    hides them
    2. The division leaves one coconut left over, which is given to the monkey.
In the morning, they split the remaining coconuts between them. This time the
split is even. There's nothing left over for the monkey.
Your task: Given the number of sailors (N), how many coconuts were in the pile
to begin with (lowest possible number)?

Formal inputs/outputs

Input

The input is a single number: N, the number of sailors. This number is a whole
number that is greater than or equal to 2.

Output

The output is a single number: the number of coconuts in the original pile.
Sample input/output

Input:
5
Output:
3121
Sample solution for 5 sailors: https://jsfiddle.net/722gjnze/8/

Credit

This challenge was originally suggested on /r/dailyprogrammer_ideas by
/u/TinyLebowski (prior to some changes by me). Have a cool challenge idea?
Hop on over to /r/dailyprogrammer_ideas to tell everyone about it!

( source: https://redd.it/43ouxy )
=#

how_many_coconut( N::Integer=2 ) = ( N^N - ( N - 1 ) )

function main( )
    local N::Integer = parse( Int, chomp( readline( STDIN ) ) )

    N >= 2 ? println( how_many_coconut( N ) ) : error( "Please, type a valid number: greater or equal to 2" )
end

main( )

max_stars_repo_name: farmfirst/dailyprogrammer

max_issues_repo_name: farmfirst/dailyprogrammer

max_forks_repo_name: farmfirst/dailyprogrammer


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, Pkg

name = "oneTBB"
version = v"2021.5.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/oneapi-src/oneTBB/archive/refs/tags/v2021.5.0.tar.gz", "e5b57537c741400cf6134b428fc1689a649d7d38d9bb9c1b6d64f092ea28178a"),
    DirectorySource("./bundled"),
]

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

# We can't do Link-Time-Optimization with Clang, disable it.
atomic_patch -p1 ../patches/clang-no-lto.patch

if [[ ${target} == *-linux-musl* ]]; then
    # Adapt patch from https://github.com/oneapi-src/oneTBB/pull/203
    atomic_patch -p1 ../patches/musl.patch

elif [[ ${target} == *-mingw* ]]; then
    #derived from https://github.com/oneapi-src/oneTBB/commit/ce476173772f289c66ba98089618c1ff767ecea4, can hopefully be removed next release
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/lowercase-windows-include.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/NOMINMAX-defines.patch
    # `CreateSemaphoreEx` requires at least Windows Vista/Server 2008:
    # https://docs.microsoft.com/en-us/windows/win32/api/winbase/nf-winbase-createsemaphoreexa
    export CXXFLAGS="-D_WIN32_WINNT=0x0600"
fi

mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DTBB_TEST=OFF \
    -DTBB_EXAMPLES=OFF \
    ..
make -j${nproc}
make install
"""

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

# Disable platforms unlikely to work
filter!(p -> arch(p) ∉ ("armv6l", "armv7l"), platforms)

#i686 mingw fails with errors about _control87
filter!(p -> !Sys.iswindows(p) || arch(p) != "i686", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libtbbmalloc", :libtbbmalloc),
    LibraryProduct("libtbbmalloc_proxy", :libtbbmalloc_proxy),
    LibraryProduct(["libtbb", "libtbb12"], :libtbb),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]
a = ssm.a
c = ssm.c
q = ssm.q
r = ssm.r
x10 = ssm.m1
P10 = ssm.P1

###########################
print("Testing acqr_ksmoother_s: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = acqr_kfilter_s(y,a,c,q,r,x10)
xtN,PtN,Pt1tN = acqr_ksmoother_s(a,xtt,Ptt,xtt1,Ptt1)
if sum(xtN.^2) > 1e-6
    println("OK")
end


text:
# Recurrent inference machines for seismic imaging
# Same notation as Putzky and Welling paper (2017)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using JUDI.TimeModeling, JUDI4Flux, LinearAlgebra, InvertibleNetworks
using PyPlot, Random, JLD, Flux
import Flux.Optimise.update!

# Training data
D = load("/data/pwitte3/models/overthrust_images_train.jld")
ntrain = length(D["m"])

# Crop images and models to 400 x 120 (must be evenly dividable by 4)
for j=1:ntrain
    D["m"][j] = D["m"][j][1:400, 1:120]
    D["m0"][j] = D["m0"][j][1:400, 1:120]
    D["dm"][j] = D["dm"][j][1:400, 1:120]
end

# Use one sample to set up operators
m = D["m"][1]
m0 = D["m0"][1]
dm = vec(D["dm"][1])

# Set up model structure
n = size(m0)
d = (25., 25.)
o = (0., 0.)

# Setup info and model structure
nsrc = 1
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 301
xrec = range(400f0, stop=9600f0, length=nxrec)
yrec = 0f0
zrec = range(250f0, stop=250f0, length=nxrec)

# receiver sampling and recording time
time = 2000f0   # receiver recording time [ms]
dt = 4f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dt, t=time, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
num_simsource = 21
xsrc = range(500f0, stop=9500f0, length=num_simsource)
ysrc = range(0f0, stop=0f0, length=num_simsource)
zsrc = range(20f0, stop=20f0, length=num_simsource)

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dt, t=time)

# setup wavelet
f0 = 0.015f0     # MHz
wavelet = zeros(Float32, srcGeometry.nt[1], num_simsource)
for j=1:num_simsource
    wavelet[:, j] = ricker_wavelet(time, dt, f0) * randn(1)[1]/sqrt(num_simsource)
end
q = judiVector(srcGeometry, wavelet)
q_sim = zeros(Float32, size(wavelet))

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model0)
info = Info(prod(n), nsrc, ntComp)

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

# Return data as julia array
opt = Options(return_array=true)

# Setup operators
F0 = judiModeling(info, model0; options=opt)
Pr = judiProjection(info, recGeometry)
Ps = judiProjection(info, srcGeometry)

# Modeling operator and Jacobian
F0 = Pr*F0*Ps'  # initial model
J = judiJacobian(F0, q)

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

# Dimensions
nx1 = n[1]
nx2 = n[2]
nx_in = 32
nx_hidden = 64
batchsize = 1

#ny1 = recGeometry.nt[1]
#ny2 = nxrec
#ny_in = 1
#ny_hidden = 64

maxiter = 8
Ψ(η) = identity(η)

SLIM = Array{AffineCouplingLayerSLIM}(undef, maxiter)
Params = Array{Parameter}(undef, 0)

# i-SLIM coupling layer
for j=1:maxiter
    SLIM[j] = AffineCouplingLayerSLIM(nx1, nx2, nx_in, nx_hidden, batchsize, Ψ; logdet=false, permute=true)
    #SLIM[j] = ConditionalLayerSLIM(nx1, nx2, nx_in, nx_hidden, ny1, ny2, ny_channel, ny_hidden, batchsize; type="affine")
    global Params = cat(Params, get_params(SLIM[j]); dims=1)
end

function forward(X, Y, J)
    for j=1:maxiter
        X = SLIM[j].forward(X, Y, J)
    end
    return X
end

function backward(dX, X, Y, J)
    for j=maxiter:-1:1
        dX, X = SLIM[j].backward(dX, X, Y, J)[[1,3]]
    end
    return dX
end

# Objective function
function loss(Y, J, x_true)

    # Initiliaze w/ zeros
    X0 = zeros(Float32, nx1, nx2, nx_in, batchsize)

    # Forward pass
    X_ = forward(X0, Y, J)

    # Residual and function value
    dX = zeros(Float32, nx1, nx2, nx_in, batchsize)
    dX[:,:,1:1,:] = X_[:,:,1:1,:] - reshape(x_true, nx1, nx2, 1, batchsize)
    f = .5f0*norm(dX)^2

    # Backward pass (set gradients)
    backward(dX, X_, Y, J)

    return f
end

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

# Get network parameters and overwrite w/ saved values
iter_start = 1
if iter_start > 1
    P_save = load(join(["../data/network_islim_iter_8_params_iteration_", string(iter_start), ".jld"]))["P"]
    for j=1:length(Params)
        Params[j].data = P_save[j].data
    end
end

if iter_start > 1
    opt = load(join(["../data/network_islim_iter_8_params_iteration_", string(iter_start), ".jld"]))["opt"]
else
    opt = Flux.ADAM(1f-3)
end

# Optimization parameters
train_iter = 1000
indices = randperm(ntrain)

# Training loop
for j=iter_start+1:train_iter

    # Draw image + velocity from training data
    i = indices[j]
    print("Source no: ", i, "\n")
    x_true = D["dm"][i]
    m0 = D["m0"][i]

    # Draw random source
    for k=1:num_simsource
        q_sim[:,k] = ricker_wavelet(time, dt, f0) * randn(Float32, 1)[1]/sqrt(1f0*num_simsource)
    end

    # Generate observed data on the fly
    J.model.m = m0
    J.source[1] = q_sim
    Y = J*vec(x_true)

    # Evaluate objective and gradients
    @time f = loss(Y, J, x_true)
    print("Iteration: ", j, "; f(x) = ", f, "\n")

    # Update weights
    for p in Params
        update!(opt, p.data, p.grad)
    end
    for j=1:length(SLIM)
        clear_grad!(SLIM[j])
    end

    # Save intermediate results
    if mod(j, 100) == 0
        save(join(["../data/network_islim_iter_8_params_iteration_", string(j), ".jld"]), "P", Params, "opt", opt)
    end
end

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
##################
# Trie #
##################

struct Trie{K,V} <: ChoiceMap
    leaf_nodes::Dict{K,V}
    internal_nodes::Dict{K,Trie{K,V}}
end

Trie{K,V}() where {K,V} = Trie(Dict{K,V}(), Dict{K,Trie{K,V}}())

# copy constructor for something supported read-only trie interface
function Trie(other)
    trie = Trie{Any,Any}()
    for (key, value) in get_leaf_nodes(other)
        set_leaf_node!(trie, key, value)
    end
    for (key, node) in get_internal_nodes(other)
        sub_trie = Trie(node)
        set_internal_node!(trie, key, sub_trie)
    end
    trie
end

# TODO come up with a better printing method (and something nice for Jupyter
# notebooks)
import JSON
Base.println(trie::Trie) = JSON.print(trie, 4)

# invariant: all internal nodes are nonempty
Base.isempty(trie::Trie) = isempty(trie.leaf_nodes) && isempty(trie.internal_nodes)
get_leaf_nodes(trie::Trie) = trie.leaf_nodes
get_internal_nodes(trie::Trie) = trie.internal_nodes

function Base.values(trie::Trie)
    iterators = convert(Vector{Any}, collect(map(values, values(trie.internal_nodes))))
    push!(iterators, values(trie.leaf_nodes))
    Iterators.flatten(iterators)
end

function has_internal_node(trie::Trie, addr)
    haskey(trie.internal_nodes, addr)
end

function has_internal_node(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        has_internal_node(trie.internal_nodes[first], rest)
    else
        false
    end
end

function get_internal_node(trie::Trie, addr)
    trie.internal_nodes[addr]
end

function get_internal_node(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        get_internal_node(trie.internal_nodes[first], rest)
    else
        throw(KeyError(trie, addr))
    end
end

function set_internal_node!(trie::Trie{K,V}, addr, new_node::Trie{K,V}) where {K,V}
    if !isempty(new_node)
        trie.internal_nodes[addr] = new_node
    end
end

function set_internal_node!(trie::Trie{K,V}, addr::Pair, new_node::Trie{K,V}) where {K,V}
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        node = trie.internal_nodes[first]
    else
        node = Trie{K,V}()
        trie.internal_nodes[first] = node
    end
    set_internal_node!(node, rest, new_node)
end

function delete_internal_node!(trie::Trie, addr)
    delete!(trie.internal_nodes, addr)
    return isempty(trie.leaf_nodes) && isempty(trie.internal_nodes)
end

function delete_internal_node!(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        node = trie.internal_nodes[first]
        if delete_internal_node!(node, rest, new_node)
            delete!(trie.internal_nodes, first)
        end
    end
    return isempty(trie.leaf_nodes) && isempty(trie.internal_nodes)
end

function has_leaf_node(trie::Trie, addr)
    haskey(trie.leaf_nodes, addr)
end

function has_leaf_node(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        has_leaf_node(trie.internal_nodes[first], rest)
    else
        false
    end
end

function get_leaf_node(trie::Trie, addr)
    trie.leaf_nodes[addr]
end

function get_leaf_node(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        get_leaf_node(trie.internal_nodes[first], rest)
    else
        throw(KeyError(trie, addr))
    end
end

function set_leaf_node!(trie::Trie, addr, value)
    trie.leaf_nodes[addr] = value
end

function set_leaf_node!(trie::Trie{K,V}, addr::Pair, value) where {K,V}
    (first, rest) = addr
    if haskey(trie.internal_nodes, first)
        node = trie.internal_nodes[first]
    else
        node = Trie{K,V}()
        trie.internal_nodes[first] = node
    end
    node = trie.internal_nodes[first]
    set_leaf_node!(node, rest, value)
end

function delete_leaf_node!(trie::Trie, addr)
    delete!(trie.leaf_nodes, addr)
    return isempty(trie.leaf_nodes) && isempty(trie.internal_nodes)
end

function delete_leaf_node!(trie::Trie, addr::Pair)
    (first, rest) = addr
    if haskey(trie.internal_nodes[first])
        node = trie.internal_nodes[first]
        if delete_leaf_node!(node, rest)
            delete!(trie.internal_nodes, first)
        end
    end
    return isempty(trie.leaf_nodes) && isempty(trie.internal_nodes)
end

Base.setindex!(trie::Trie, value, addr) = set_leaf_node!(trie, addr, value)

function Base.merge!(a::Trie{K,V}, b::Trie{K,V}) where {K,V}
    merge!(a.leaf_nodes, b.leaf_nodes)
    for (key, a_sub) in a.sub
        if haskey(b.sub, key)
            b_sub = b.sub[key]
            merge!(a_sub, b_sub)
        end
    end
    for (key, b_sub) in b.sub
        if !haskey(a.sub, key)
            a.sub[key] = b_sub
        end
    end
    a
end

get_address_schema(::Trie) = DynamicSchema()

Base.haskey(trie::Trie, key) = has_leaf_node(trie, key)

Base.getindex(trie::Trie, key) = get_leaf_node(trie, key)

export Trie
export set_internal_node!
export delete_internal_node!
export set_leaf_node!
export delete_leaf_node!


text:
module TestStepwise
using Test
using CompressedSensing: gaussian_data, perturb, rmp, foba

@testset "Stepwise Regression Algorithms" begin
    n, m, k = 32, 64, 3
    A, x, b = gaussian_data(n, m, k)
    δ = 1e-2
    y = perturb(b, δ)

    @testset "Relevance Matching Pursuit" begin
        xrmp = rmp(A, b, k)
        # noiseless
        @test xrmp.nzind == x.nzind
        @test xrmp.nzval ≈ x.nzval

        # noisy
        xrmp = rmp(A, y, δ)
        @test xrmp.nzind == x.nzind
        @test isapprox(xrmp.nzval, x.nzval, atol = 2δ)

        # multiple outer loop iterations
        xrmp = rmp(A, y, δ, 3)
        @test xrmp.nzind == x.nzind
        @test isapprox(xrmp.nzval, x.nzval, atol = 2δ)
    end

    @testset "FoBa" begin
        # noiseless
        xfoba = foba(A, b, δ)
        @test xfoba.nzind == x.nzind
        @test xfoba.nzval ≈ x.nzval

        # noisy
        xfoba = foba(A, y, δ)
        @test xfoba.nzind == x.nzind
        @test isapprox(xfoba.nzval, x.nzval, atol = 2δ)
    end
end

end