text:
# # Tracking transient deformation of a cantilever beam: centered difference

# Source code: [`bending_wave_Ray_expl_cd_tut.jl`](bending_wave_Ray_expl_cd_tut.jl)

# ## Description

# A cantilever beam is given an initial velocity and then at time 0.0 it is
# suddenly stopped by fixing one of its ends. This sends a wave down the beam.

# The beam is modeled as a solid. Trapezoidal rule is used to integrate the
# equations of motion in time. No damping is present.

# ## Goals

# - Show how to create the discrete model for explicit dynamics.
# - Demonstrate  centered difference time stepping.

##
# ## Definitions

# Basic imports.
using LinearAlgebra
using Arpack

# This is the finite element toolkit itself.
using FinEtools

# The linear stress analysis application is implemented in this package.
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule

# Input parameters
E = 205000*phun("MPa");# Young's modulus
nu = 0.3;# Poisson ratio
rho = 7850*phun("KG*M^-3");# mass density
loss_tangent = 0.0001;
frequency = 1/0.0058;
Rayleigh_mass = 2*loss_tangent*(2*pi*frequency);
L = 200*phun("mm");
W = 4*phun("mm");
H = 8*phun("mm");
tolerance = W/500;
vmag = 0.1*phun("m")/phun("SEC");
tend = 0.013*phun("SEC");

##
# ## Create the discrete model

MR = DeforModelRed3D
fens,fes  = H8block(L,W,H, 50,1,4)

geom = NodalField(fens.xyz)
u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

nl = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
setebc!(u, nl, true, 1)
setebc!(u, nl, true, 2)
setebc!(u, nl, true, 3)
applyebc!(u)
numberdofs!(u)

corner = selectnode(fens, nearestto=[0 0 0])
cornerzdof = u.dofnums[corner[1], 3]

material = MatDeforElastIso(MR, rho, E, nu, 0.0)

femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material)
femm = associategeometry!(femm, geom)
K = stiffness(femm, geom, u)
# Assemble the mass matrix as diagonal. The HRZ lumping technique is
# applied through the assembler of the sparse matrix.
hrzass = SysmatAssemblerSparseHRZLumpingSymm(0.0)
femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
M = mass(femm, hrzass, geom, u)

# Form the damping matrix.
C = Rayleigh_mass * M

# Figure out the highest frequency in the model, and use a time step that is
# considerably larger than the period of that highest frequency.
evals, evecs = eigs(K, M; nev=1, which=:LM);
@show dt = 0.99 * 2/sqrt(evals[1]);

# The time stepping loop is protected by `let end` to avoid unpleasant surprises
# with variables getting clobbered by globals.
ts, corneruzs = let dt = dt
    # Initial displacement, velocity, and acceleration.
    U0 = gathersysvec(u)
    v = deepcopy(u)
    v.values[:, 3] .= vmag
    V0 = gathersysvec(v)
    F1 = fill(0.0, length(V0))
    U1 = fill(0.0, length(V0))
    V1 = fill(0.0, length(V0))
    A0 = fill(0.0, length(V0))
    A1 = fill(0.0, length(V0))

    # The times and displacements of the corner will be collected into two vectors
    ts = Float64[]
    corneruzs = Float64[]
    # Let us begin the time integration loop:
    t = 0.0;
    step = 0;
    while t < tend
        push!(ts, t)
        push!(corneruzs, U0[cornerzdof])
        t = t+dt;
        step = step + 1;
        (mod(step,1000)==0) && println("Step $(step): $(t)")
        # Zero out the load
        fill!(F1, 0.0);
        # Initial acceleration
        if step == 1
            A0 = M \ (F1)
        end
        # Update displacement.
        @. U1 = U0 + dt*V0 + (dt^2/2)*A0;
        # Compute updated acceleration.
        A1 .= M \ (-K*U1 + F1)
        # Update the velocities.
        @. V1 = V0 + (dt/2)*(A0 + A1)
        # Switch the temporary vectors for the next step.
        U0, U1 = U1, U0;
        V0, V1 = V1, V0;
        A0, A1 = A1, A0;
        if (t == tend) # Are we done yet?
            break;
        end
        if (t+dt > tend) # Adjust the last time step so that we exactly reach tend
            dt = tend-t;
        end
    end
    ts, corneruzs # return the collected results
end

##
# ## Plot the results

using Gnuplot

@gp "set terminal windows 4 " :-
@gp  :- ts corneruzs./phun("mm") "lw 2 lc rgb 'red' with lines title 'Displacement of the corner' "
@gp  :- "set xlabel 'Time [s]'"
@gp  :- "set ylabel 'Displacement [mm]'"

# The end.
true


text:
## First process examples
const EXAMPLES_OUT = joinpath(@__DIR__, "src", "examples")
ispath(EXAMPLES_OUT) && rm(EXAMPLES_OUT; recursive=true)
mkpath(EXAMPLES_OUT)

# TODO use the general approach of ApproximateGPs.jl
examples = filter!(isdir, readdir(joinpath(@__DIR__, "..", "examples"); join=true))
let script = "using Pkg; Pkg.activate(ARGS[1]); Pkg.instantiate()"
    for example in examples
        if !success(`$(Base.julia_cmd()) -e $script $example`)
            error(
                "project environment of example ",
                basename(example),
                " could not be instantiated",
            )
        end
    end
end
# Run examples asynchronously
processes = let literatejl = joinpath(@__DIR__, "literate.jl")
    map(examples) do example
        return run(
            pipeline(
                `$(Base.julia_cmd()) $literatejl $(basename(example)) $EXAMPLES_OUT`;
                stdin=devnull,
                stdout=devnull,
                stderr=stderr,
            );
            wait=false,
        )::Base.Process
    end
end
# Check that all examples were run successfully
isempty(processes) || success(processes) || error("some examples were not run successfully")

## Build the docs
using AugmentedGPLikelihoods
using Documenter
using DocumenterCitations
using Literate
using Pkg

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

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))

likelihoods = filter!(filename -> endswith(filename, ".md"), readdir(EXAMPLES_OUT))

@info "Available likelihoods: $(likelihoods)"

makedocs(
    bib;
    modules=[AugmentedGPLikelihoods],
    authors="Théo Galy-Fajou <theo.galyfajou@gmail.com> and contributors",
    repo="https://github.com/JuliaGaussianProcesses/AugmentedGPLikelihoods.jl/blob/{commit}{path}#{line}",
    sitename="AugmentedGPLikelihoods.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaGaussianProcesses.github.io/AugmentedGPLikelihoods.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Likelihoods" => joinpath.(Ref("likelihoods"), likelihoods),
        "Examples" => joinpath.(Ref("examples"), likelihoods),
        "Misc" => "misc.md",
        "References" => "references.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/AugmentedGPLikelihoods.jl",
    devbranch="main",
    push_preview=true,
)

ext: jl


text:
include("plotting.jl")

# animate_words("hello", "heIIo"; D = KL(1.0), save_path=abspath(joinpath(@__DIR__, "hello_heIIo.gif")))
# animate_words("Julia", "Strings"; D = Balanced(), normalize_density=true, save_path=abspath(joinpath(@__DIR__, "julia_strings.gif")), duration=4)
# animate_words("Julia", "distances"; D = Balanced(), normalize_density=true, save_path=abspath(joinpath(@__DIR__, "Julia_distances.gif")))

# animate_words("DifferentialEquations", "DifferentiaIEquations"; D = KL(1.0), save_path=abspath(joinpath(@__DIR__, "DifferentiaIEquations.gif")), duration=4)

# animate_words("FIux", "Flux"; D = KL(5.0), save_path=abspath(joinpath(@__DIR__, "FIux.gif")), duration=4)

# animate_words("Julia", "Visual"; D = KL(1.0), save_path=abspath(joinpath(@__DIR__, "FIux.gif")), duration=4)

animate_words("Julia", "visual"; D = Balanced(), normalize_density=true, save_path=abspath(joinpath(@__DIR__, "julia_visual.gif")), duration=4,
total_frames = 50*4,
pause_frames = 50*4 ÷ 3)

# animate_words("Julia", "visual"; D = KL(5), save_path=abspath(joinpath(@__DIR__, "julia_visual_unbalanced.gif")), duration=4)

# julia> calculate_contributions(KL(5.0), word_measure("Flux"), word_measure("FIux"), 0.1)
# (transport_cost = 4.074841771685898, regularization = 496.6252591146831, marginal_1_penalty = 0.1502096716498963, marginal_2_penalty = 1.2417285590962812)


text:
colors_tum = Dict{Symbol, RGB{Float64}}(
    :tum_main_blue     => RGB{Float64}(  0/255, 101/255, 189/255),
    :tum_sec_blue1     => RGB{Float64}(  0/255,  82/255, 147/255),
    :tum_sec_blue2     => RGB{Float64}(  0/255,  51/255,  89/255),
    :tum_sec_gray1     => RGB{Float64}( 51/255,  51/255,  51/255),
    :tum_sec_gray2     => RGB{Float64}(128/255, 128/255, 128/255),
    :tum_sec_gray3     => RGB{Float64}(204/255, 204/255, 204/255),
    :tum_accent_gray   => RGB{Float64}(218/255, 215/255, 203/255),
    :tum_accent_orange => RGB{Float64}(227/255, 114/255,  34/255),
    :tum_accent_green  => RGB{Float64}(162/255, 173/255,   0/255),
    :tum_accent_blue1  => RGB{Float64}(152/255, 198/255, 234/255),
    :tum_accent_blue2  => RGB{Float64}(100/255, 160/255, 200/255),
)


text:
## SciMLBase Trait Definitions

SciMLBase.isautodifferentiable(alg::AbstractMethodOfStepsAlgorithm) = SciMLBase.isautodifferentiable(alg.alg)
SciMLBase.allows_arbitrary_number_types(alg::AbstractMethodOfStepsAlgorithm) = SciMLBase.allows_arbitrary_number_types(alg.alg)
SciMLBase.allowscomplex(alg::AbstractMethodOfStepsAlgorithm) = SciMLBase.allowscomplex(alg.alg)
SciMLBase.isdiscrete(alg::AbstractMethodOfStepsAlgorithm) = SciMLBase.isdiscrete(alg.alg)
SciMLBase.isadaptive(alg::AbstractMethodOfStepsAlgorithm) = SciMLBase.isadaptive(alg.alg)

## DelayDiffEq Internal Traits

isconstrained(alg::AbstractMethodOfStepsAlgorithm{constrained}) where constrained =
    constrained
OrdinaryDiffEq.uses_uprev(alg::AbstractMethodOfStepsAlgorithm, adaptive) = true

OrdinaryDiffEq.isimplicit(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.isimplicit(alg.alg)
OrdinaryDiffEq.isdtchangeable(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.isdtchangeable(alg.alg)
OrdinaryDiffEq.ismultistep(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.ismultistep(alg.alg)
OrdinaryDiffEq.isautoswitch(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.isautoswitch(alg.alg)
OrdinaryDiffEq.get_chunksize(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.get_chunksize(alg.alg)
OrdinaryDiffEq.get_chunksize_int(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.get_chunksize_int(alg.alg)
OrdinaryDiffEq.alg_autodiff(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_autodiff(alg.alg)
OrdinaryDiffEq.alg_difftype(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_difftype(alg.alg)
OrdinaryDiffEq.standardtag(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.standardtag(alg.alg)
OrdinaryDiffEq.concrete_jac(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.concrete_jac(alg.alg)
OrdinaryDiffEq.alg_extrapolates(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_extrapolates(alg.alg)
OrdinaryDiffEq.alg_order(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_order(alg.alg)
OrdinaryDiffEq.alg_maximum_order(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_maximum_order(alg.alg)
OrdinaryDiffEq.alg_adaptive_order(alg::AbstractMethodOfStepsAlgorithm) = OrdinaryDiffEq.alg_adaptive_order(alg.alg)

"""
    iscomposite(alg)

Return if algorithm `alg` is a composite algorithm.
"""
iscomposite(alg) = false
iscomposite(::OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm) = true
iscomposite(alg::AbstractMethodOfStepsAlgorithm) = iscomposite(alg.alg)

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using ThreeBodyProblem
using RecipesBase
using Test

r = rand(Float64)
x,y = circle(r)
@test x[1]^2 + y[1]^2 == r^2

x,y,z = sphere(r)
@test x[1]^2 + y[1]^2 + z[1]^2 == r^2

sys = earth_moon()
# plot(sys)


text:
using LocalBinaryPatterns: average_mode

@testset "local_binary_pattern" begin
    # reference result comes from [1]
    # - [1] T. Ojala, M. Pietikäinen, and D. Harwood, “A comparative study of texture measures with classification based on featured distributions,” _Pattern Recognition_, vol. 29, no. 1, pp. 51–59, Jan. 1996, doi: 10.1016/0031-3203(95)00067-4.
    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [192 64 0; 104 169 27; 40 107 0]

    out = local_binary_pattern(X)
    @test eltype(out) == UInt8
    @test size(out) == (3, 3)
    @test out == ref_out

    # ensure default parameter values are not changed accidently
    @test out == local_binary_pattern(X; rotation=false, uniform_degree=nothing)

    # intensity-invariant
    @test local_binary_pattern(X./255) == out
    # Gray inputs
    @test local_binary_pattern(Gray.(X./255)) == out

    # LBP is only defined for scalar values
    @test_throws MethodError local_binary_pattern(RGB.(Gray.(X./255)))

    # not yet ready for N-dimensional array (although it's doable)
    @test_throws MethodError local_binary_pattern(rand(3, 3, 3))

    @testset "OffsetArrays" begin
        Xo = OffsetArray(X, -1, -1)
        out = local_binary_pattern(Xo)
        @test axes(out) == axes(Xo)
        @test OffsetArrays.no_offset_view(out) == ref_out
    end

    @testset "Rotation Invariant" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [3 1 0; 13 53 27; 5 91 0]

        out = local_binary_pattern(X; rotation=true)
        @test eltype(out) == UInt8
        @test size(out) == (3, 3)
        @test out == ref_out
    end

    @testset "Uniform encoding" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [192 64 0; 9 9 9; 9 9 0]

        out = local_binary_pattern(X; rotation=false, uniform_degree=2)
        @test eltype(out) == UInt8
        @test size(out) == (3, 3)
        @test out == ref_out
    end

    @testset "Rotation Invariant, Uniform encoding" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [3 1 0; 9 9 9; 9 9 0]

        out = local_binary_pattern(X; rotation=true, uniform_degree=2)
        @test eltype(out) == UInt8
        @test size(out) == (3, 3)
        @test out == ref_out
    end
end

@testset "local_binary_pattern, Interpolation-based" begin
    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [1 1 0; 3 2 14; 2 7 0]
    out = local_binary_pattern(X, 4, 1)
    @test eltype(out) == UInt32
    @test size(out) == (3, 3)
    @test out == ref_out

    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [129 1 0; 7 14 124; 6 31 0]
    out = local_binary_pattern(X, 8, 1)
    @test eltype(out) == UInt32
    @test size(out) == (3, 3)
    @test out == ref_out

    @testset "OffsetArrays" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [129 1 0; 7 14 124; 6 31 0]
        Xo = OffsetArray(X, -1, -1)
        out = local_binary_pattern(Xo, 8, 1)
        @test axes(out) == axes(Xo)
        @test OffsetArrays.no_offset_view(out) == ref_out
    end

    # ensure default parameter values are not changed accidently
    out = local_binary_pattern(X, 8, 1)
    @test local_binary_pattern(X, 8, 1, Linear(); rotation=false, uniform_degree=nothing) == out
    @test local_binary_pattern(X, 8, 1, Linear()) == local_binary_pattern(X, 8, 1, BSpline(Linear()))

    # The pattern encoding is different because the offset orders are not the same
    @test local_binary_pattern(X, 8, 1, Constant()) != local_binary_pattern(X)

    @testset "Rotation Invariant" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [129 1 0; 7 7 31; 3 31 0]
        out = local_binary_pattern(X, 8, 1; rotation=true)
        @test eltype(out) == UInt32
        @test size(out) == (3, 3)
        @test out == ref_out
    end

    @testset "Uniform encoding" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [33 1 0; 33 33 33; 33 33 0]
        out = local_binary_pattern(X, 8, 1; uniform_degree=2)
        @test eltype(out) == UInt32
        @test size(out) == (3, 3)
        @test out == ref_out
    end

    @testset "Rotation Invariant, Uniform encoding" begin
        X = [6 7 9; 5 6 3; 2 1 7]
        ref_out = [33 1 0; 33 33 33; 33 33 0]
        out = local_binary_pattern(X, 8, 1; rotation=true, uniform_degree=2)
        @test eltype(out) == UInt32
        @test size(out) == (3, 3)
        @test out == ref_out
    end

    @testset "radius" begin
        X = [3  1  1  2  2
             1  2  1  5  5
             1  1  2  4  5
             3  3  5  1  2
             3  3  5  2  2]
        ref_out = [0   9   13  8   8
                   9   9   13  0   0
                   11  11  9   0   0
                   1   0   0   6   6
                   1   0   0   6   6]
        out = local_binary_pattern(X, 4, 2)
        @test eltype(out) == UInt32
        @test size(out) == (5, 5)
        @test out == ref_out

        @test_nowarn local_binary_pattern(X, 4, 1.5)
        @test_throws ArgumentError local_binary_pattern(X, 4, 0.5)
    end
end

@testset "average mode" begin
    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [192 64 0; 104 169 27; 40 105 1]
    out = local_binary_pattern(LBP(average_mode, (1, 1)), X)
    @test eltype(out) == UInt8
    @test size(out) == (3, 3)
    @test out == ref_out

    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [3 1 0; 13 53 27; 5 45 1]
    out = local_binary_pattern(LBP(average_mode, (1, 1)), X, rotation=true)
    @test eltype(out) == UInt8
    @test size(out) == (3, 3)
    @test out == ref_out

    X = [6 7 9; 5 6 3; 2 1 7]
    ref_out = [1 1 0; 3 6 28; 2 15 0]
    out = local_binary_pattern(LBP(average_mode, (1, 1)), X, 6, 1)
    @test eltype(out) == UInt32
    @test size(out) == (3, 3)
    @test out == ref_out
end

@testset "multiblock LBP" begin
    # Only the center-pixel [4, 4] is well-defined in the reference paper [1].
    # [1] Zhang, Lun, et al. "Face detection based on multi-block lbp representation." _International conference on biometrics_. Springer, Berlin, Heidelberg, 2007.
    #
    # 50  0  0      1 0 0    1 0 0
    #  0  1  0   => 0 x 0 => 0 x 0   => 0b10000001 (129)
    #  0  0  50     0 0 1    0 0 128
    X = [
        50  50  50  0  0  0   0   0   0
        50  50  50  0  0  0   0   0   0
        50  50  50  0  0  0   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  0  0  0  50  50  50
         0   0   0  0  0  0  50  50  50
         0   0   0  0  0  0  50  50  50
    ]
    ref_out = [
        0    0    0  214  214  214  22  22  22
        0    0  128  146  150  150  22  22  22
        0  128  128  130  146  150  22  22  22
      248  200  136  129  145  147  31  31  31
      248  232  200  193  209  209  16  20  20
      248  232  232  201  209  208  16  16  20
      104  104  104  107   64   64   0   0   0
      104  104  104  107   96   64   0   0   0
      104  104  104  107   96   96   0   0   0
    ]

    out = local_binary_pattern(LBP((3, 3)), X)
    @test eltype(out) == UInt8
    @test size(out) == (9, 9)
    @test out == ref_out

    ref_out = [
        0    0    6  214  214  214  214  22  22
        0    0    2  214  214  214  214  22  22
       40    8   11  147   23   22  255  31  31
      248  248  201  129  129  150  255  31  31
      248  248  105  129  129  147  255  31  31
      248  248  104  232  201  208  208  20  22
      248  248  255  255  255  208  208  16  22
      104  104  107  107  107   96   64   0   2
      104  104  107  107  107  104  104   8  11
    ]
    @test ref_out == local_binary_pattern(LBP((2, 2)), X)

    # ensure default parameter values are not changed accidently
    @test out == local_binary_pattern(LBP((3, 3)), X; rotation=false, uniform_degree=nothing)

    # intensity-invariant
    @test_broken local_binary_pattern(LBP((3, 3)), X./3) == out # float-point numerical error
    @test local_binary_pattern(LBP((3, 3)), X./4) == out
    # Gray inputs
    @test local_binary_pattern(LBP((3, 3)), Gray.(X./4)) == out

    # when `block_size == (1, 1)`, it degenerates to the original pixel version
    @test local_binary_pattern(LBP((1, 1)), X) == local_binary_pattern(X)

    # LBP is only defined for scalar values
    @test_throws MethodError local_binary_pattern(LBP((3, 3)), RGB.(Gray.(X./255)))

    # check block_size
    @test_throws ArgumentError local_binary_pattern(LBP((0, 3)), X, (0, 3))
    @test_throws MethodError local_binary_pattern(LBP((1.5, 1.5)), X)
    @test_throws MethodError local_binary_pattern(LBP(3), X)
    @test_throws MethodError local_binary_pattern(LBP((3, )), X)

    # not yet ready for N-dimensional array (although it's doable)
    @test_throws MethodError local_binary_pattern(rand(3, 3, 3), (3, 3))

    @testset "OffsetArrays" begin
        Xo = OffsetArray(X, -1, -1)
        out = local_binary_pattern(LBP((3, 3)), Xo)
        @test axes(out) == axes(Xo)
        @test OffsetArrays.no_offset_view(out) == local_binary_pattern(LBP((3, 3)), X)
    end

    @testset "Rotation Invariant" begin
        ref_out = [
            0   0   0  91  91  91  11  11  11
            0   0   1  37  45  45  11  11  11
            0   1   1   5  37  45  11  11  11
           31  25  17   3  25  39  31  31  31
           31  29  25   7  29  29   1   5   5
           31  29  29  39  29  13   1   1   5
           13  13  13  91   1   1   0   0   0
           13  13  13  91   3   1   0   0   0
           13  13  13  91   3   3   0   0   0
        ]

        out = local_binary_pattern(LBP((3, 3)), X; rotation=true)
        @test eltype(out) == UInt8
        @test size(out) == (9, 9)
        @test out == ref_out
    end

    @testset "Uniform encoding" begin
        ref_out = [
            0    0    0  9   9   9   9   9   9
            0    0  128  9   9   9   9   9   9
            0  128  128  9   9   9   9   9   9
          248    9    9  9   9   9  31  31  31
          248    9    9  9   9   9  16   9   9
          248    9    9  9   9   9  16  16   9
            9    9    9  9  64  64   0   0   0
            9    9    9  9  96  64   0   0   0
            9    9    9  9  96  96   0   0   0
        ]

        out = local_binary_pattern(LBP((3, 3)), X; rotation=false, uniform_degree=2)
        @test eltype(out) == UInt8
        @test size(out) == (9, 9)
        @test out == ref_out
    end

    @testset "Rotation Invariant, Uniform encoding" begin
        ref_out = [
            0  0  0  9  9  9   9   9   9
            0  0  1  9  9  9   9   9   9
            0  1  1  9  9  9   9   9   9
           31  9  9  3  9  9  31  31  31
           31  9  9  7  9  9   1   9   9
           31  9  9  9  9  9   1   1   9
            9  9  9  9  1  1   0   0   0
            9  9  9  9  3  1   0   0   0
            9  9  9  9  3  3   0   0   0
        ]

        out = local_binary_pattern(LBP((3, 3)), X; rotation=true, uniform_degree=2)
        @test eltype(out) == UInt8
        @test size(out) == (9, 9)
        @test out == ref_out
    end

    @testset "Interpolation-based" begin
        @test_throws ArgumentError local_binary_pattern(LBP((3, 3)), X, 4, 1)
        @test_throws ArgumentError local_binary_pattern(LBP((3, 3)), X, 4, 1; rotation=true)
        @test_throws ArgumentError local_binary_pattern(LBP((3, 3)), X, 4, 1; uniform_degree=2)
    end
end

@testset "multiblock LBP + average_mode" begin
    X = [
        50  50  50  0  0  0   0   0   0
        50  50  50  0  0  0   0   0   0
        50  50  50  0  0  0   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  1  1  1   0   0   0
         0   0   0  0  0  0  50  50  50
         0   0   0  0  0  0  50  50  50
         0   0   0  0  0  0  50  50  50
    ]
    ref_out = [
        0    0    0    2    2    2   4   4   4
        0    0    0  130  130  146  16  20  20
        0    0  128  130  146  144  16  20  20
        8  136  136  129  145  145  16  20  20
        8  136  200  193  209  208  16  20  20
        8  200  192  193  208  208  16  16  20
       32   64   64   64   64   64   0   0   0
       32   96   96   96   96   64   0   0   0
       32   96   96   96   96   96   0   0   0
    ]
    out = local_binary_pattern(LBP(average_mode, (3, 3)), X)
    @test out == ref_out

    # the center pixel is well-defined here
    @test out[4, 4] == local_binary_pattern(LBP(average_mode, (1, 1)), [50 0 0; 0 1 0; 0 0 50])[2, 2]
end


text:
n, p = 500, 5
((X, y, θ), (X1, y1, θ1)) = generate_binary(n, p; seed=52)

@testset "Logreg" begin
    # No intercept
    λ = 5.0
    lr = LogisticRegression(λ; fit_intercept=false)
    J  = objective(lr, X, y)
    o  = LogisticLoss() + λ * L2Penalty()
    @test J(θ) == o(y, X*θ, θ)
    θ_newton   = fit(lr, X, y, solver=Newton())
    θ_newtoncg = fit(lr, X, y, solver=NewtonCG())
    θ_lbfgs    = fit(lr, X, y, solver=LBFGS())
    @test isapprox(J(θ),          282.09960, rtol=1e-5)
    @test isapprox(J(θ_newton),   280.37472, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_newtoncg), 280.37472, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    280.37472, rtol=1e-5)

    # With intercept
    lr1 = LogisticRegression(λ; penalize_intercept=true)
    J   = objective(lr1, X, y1)
    θ_newton   = fit(lr1, X, y1, solver=Newton())
    θ_newtoncg = fit(lr1, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(lr1, X, y1, solver=NewtonCG())
    @test isapprox(J(θ1),         213.80970, rtol=1e-5)
    @test isapprox(J(θ_newton),   209.32539, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_newtoncg), 209.32539, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    209.32539, rtol=1e-5)

    # with intercept and not penalizing it
    lr1 = LogisticRegression(λ)
    J   = objective(lr1, X, y1)
    θ_newton   = fit(lr1, X, y1, solver=Newton())
    θ_newtoncg = fit(lr1, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(lr1, X, y1, solver=NewtonCG())
    @test isapprox(J(θ1),         212.29192, rtol=1e-5)
    @test isapprox(J(θ_newton),   207.71595, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_newtoncg), 207.71595, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    207.71595, rtol=1e-5)

    if DO_COMPARISONS
        # This checks that the parameters recovered using Sklearn lead
        # to a similar loss than the one given by our code to verify the
        # correctness of the code.
        lr_sk_ncg = SKLEARN_LM.LogisticRegression(C=1.0/λ, solver="newton-cg")
        lr_sk_ncg.fit(X, y1)
        θ1_sk_ncg = vcat(lr_sk_ncg.coef_[:], lr_sk_ncg.intercept_)
        lr_sk_lbfgs = SKLEARN_LM.LogisticRegression(C=1.0/λ, solver="lbfgs")
        lr_sk_lbfgs.fit(X, y1)
        θ1_sk_lbfgs = vcat(lr_sk_lbfgs.coef_[:], lr_sk_lbfgs.intercept_)
        # Comparing with ours
        @test isapprox(J(θ1_sk_ncg),   207.71595, rtol=1e-5)
        @test isapprox(J(θ1_sk_lbfgs), 207.71595, rtol=1e-5)
    end
end

n, p, c = 500, 5, 4
((X, y, θ), (X1, y1, θ1)) = generate_multiclass(n, p, c; seed=525)

@testset "Multinomial" begin
    # No intercept
    λ = 5.0
    mnr = MultinomialRegression(λ; fit_intercept=false)
    J   = objective(mnr, X, y; c=c)
    θ_newtoncg = fit(mnr, X, y, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y, solver=R.LBFGS())
    @test isapprox(J(θ),          419.67795, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 384.33810, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_lbfgs),    384.33810, rtol=1e-5)

    #  With intercept
    mnr = MultinomialRegression(λ; penalize_intercept=true)
    J   = objective(mnr, X, y1; c=c)
    θ_newtoncg = fit(mnr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y1, solver=R.LBFGS())
    @test isapprox(J(θ),          1244.46404, rtol=1e-5)
    @test isapprox(J(θ_newtoncg),  388.75018, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_lbfgs),     388.75018, rtol=1e-5)

    mnr = MultinomialRegression(λ)
    J   = objective(mnr, X, y1; c=c)
    θ_newtoncg = fit(mnr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y1, solver=R.LBFGS())
    @test isapprox(J(θ),          1242.70990, rtol=1e-5)
    @test isapprox(J(θ_newtoncg),  383.93018, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_lbfgs),     383.93018, rtol=1e-5)

    if DO_COMPARISONS
        lr_sk_ncg = SKLEARN_LM.LogisticRegression(C=1.0/λ, solver="newton-cg",
                                             multi_class="multinomial")
        lr_sk_ncg.fit(X, y1)
        θ1_sk_ncg = vec(vcat(lr_sk_ncg.coef_', lr_sk_ncg.intercept_'))
        lr_sk_lbfgs = SKLEARN_LM.LogisticRegression(C=1.0/λ, solver="lbfgs",
                                               multi_class="multinomial")
        lr_sk_lbfgs.fit(X, y1)
        θ1_sk_lbfgs = vec(vcat(lr_sk_ncg.coef_', lr_sk_ncg.intercept_'))
        # Comparing with ours
        @test isapprox(J(θ1_sk_ncg),   385.67895, rtol=1e-5)
        @test isapprox(J(θ1_sk_lbfgs), 385.67895, rtol=1e-5)
    end
end

n, p = 500, 100
(_, (X, y, θ)) = generate_binary(n, p;  seed=52551, sparse=0.1)

@testset "Logreg/EN" begin
    ρ = 0.8
    α = 0.03
    λ = α * (1 - ρ) * n
    γ = α * ρ * n
    enlr = LogisticRegression(λ, γ; penalize_intercept=true)
    J    = objective(enlr, X, y)
    θ_fista = fit(enlr, X, y)
    θ_ista  = fit(enlr, X, y, solver=ISTA())
    @test isapprox(J(θ),       271.43962, rtol=1e-5)
    @test isapprox(J(θ_fista), 250.46547, rtol=1e-5) # <- ref values
    @test isapprox(J(θ_ista),  250.46547, rtol=1e-5)
    # sparsity
    @test nnz(θ)       == 7  # generating θ
    @test nnz(θ_fista) == 18
    @test nnz(θ_ista)  == 18

    # pure l1 regularization (unclear how sklearn does the mixing of l1/l2 for logreg)
    enlr = LogisticRegression(γ; penalty=:l1)
    J    = objective(enlr, X, y)
    θ_fista = fit(enlr, X, y)
    @test isapprox(J(θ),       263.76834, rtol=1e-5)
    @test isapprox(J(θ_fista), 245.65841, rtol=1e-5) # <- ref value
    @test nnz(θ_fista) == 17

    if DO_COMPARISONS
        # NOTE: this algorithm is stochastic so can't have hard tests
        PY_RND.seed(1531)
        enlr_sk = SKLEARN_LM.LogisticRegression(penalty="elasticnet", C=1.0/γ, l1_ratio=1, solver="saga")
        enlr_sk.fit(X, y)
        θ_sk = vcat(enlr_sk.coef_[:], enlr_sk.intercept_)
        @test isapprox(J(θ_sk), 245.658, rtol=1e-3)
        @test 17 ≤ nnz(θ_sk) ≤ 19
    end
end

n, p, c = 1000, 100, 3
(_, (X, y, θ)) = generate_multiclass(n, p, c;  seed=33, sparse=0.1)

@testset "Multin/EN" begin
    λ = 10
    γ = 50
    enmnr = MultinomialRegression(λ, γ; penalize_intercept=true)
    J     = objective(enmnr, X, y; c=c)
    θ_fista = fit(enmnr, X, y)
    θ_ista  = fit(enmnr, X, y)
    @test isapprox(J(θ),       1574.61741, rtol=1e-5)
    @test isapprox(J(θ_fista),  922.36776, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_ista),   922.36776, rtol=1e-5)
    @test nnz(θ_fista) == 15
    @test nnz(θ_ista)  == 15

    # pure l1 regularization for sklearn comp
    enmnr = MultinomialRegression(γ; penalty=:l1)
    J     = objective(enmnr, X, y; c=c)
    θ_fista = fit(enmnr, X, y)
    @test isapprox(J(θ),       1432.37281, rtol=1e-5)
    @test isapprox(J(θ_fista),  912.25297, rtol=1e-5) # <- ref value
    @test nnz(θ_fista) == 16

    if DO_COMPARISONS
        # NOTE: this algorithm is stochastic
        PY_RND.seed(1531)
        enmnr_sk = SKLEARN_LM.LogisticRegression(penalty="elasticnet", C=1.0/γ, l1_ratio=1,
                                            solver="saga", multi_class="multinomial")
        enmnr_sk.fit(X, y)
        θ_sk = enmnr_sk.coef_'[:]
        @test isapprox(J(θ_sk), 912.39215, rtol=1e-3)
        @test 14 ≤ nnz(θ_sk) ≤ 16
    end
end


text:
@testset "nri" begin

  #example om from NRI paper
  om1 = [0 0 0 0 0; 0 2 0 0 1;
         0 0 0 1 0; 0 0 3 0 0;
         0 1 0 0 0];
  om1 = sparse(om1);

  fullTPs, segTPs = NRI.compute_TPs(om1)

  @test segTPs[1] == 0
  @test segTPs[2] == 1
  @test segTPs[3] == 0
  @test segTPs[4] == 3
  @test segTPs[5] == 0

  fullFPs, segFPs = NRI.compute_FPs(om1)

  @test segFPs[1] == 0
  @test segFPs[2] == 1
  @test segFPs[3] == 0
  @test segFPs[4] == 0
  @test segFPs[5] == 1

  fullFNs, segFNs = NRI.compute_FNs(om1)

  @test segFNs[1] == 0
  @test segFNs[2] == 2
  @test segFNs[3] == 0
  @test segFNs[4] == 0
  @test segFNs[5] == 0

  fullNRI, segNRI, segNRIw = NRI.nri(om1)

  @test segNRI[2] == 0.4
  @test segNRI[4] == 1
  @test segNRI[5] == 0

end


text:
using IntervalWavelets

# Interior
x, y = DaubScaling(4,10);
C = ifilter(4,true)
z = DaubScaling(C,10);
plot(x, y, x, z, "r")
savefig("interior.pdf")

# Left boundary
x, Y = DaubScaling(4, 'L', 10);
Y[1,:] = Y[2,:]
plot(x, Y)
legend( ["k=0"; "k=1"; "k=2"; "k=3"], (0.8,0.9) )
savefig("left.pdf")

# Right boundary
x, Y = DaubScaling(4, 'R', 10);
Y[end,:] = Y[end-1,:]
plot(x, Y)
legend( ["k=0"; "k=1"; "k=2"; "k=3"] )
savefig("right.pdf")


text:
# https://pkgdocs.julialang.org/v1/creating-packages/
#   julia
#   ] generate Waterfall
#   ; cd Waterfall
#   ] activate .
#   ] add DataFrames, Luxor, CSV, Statistics
#   import Waterfall

module Waterfall

WATERFALL_DIR = abspath(joinpath(dirname(Base.find_package("Waterfall")), ".."))
export WATERFALL_DIR

# include("includes.jl")

# export WIDTH
# export HEIGHT
# export BORDER
# export TOP_BORDER
# export BOTTOM_BORDER
# export LEFT_BORDER
# export RIGHT_BORDER
# export SEP

# export convert
# export Cascade
# export Data
# export Geometry
# export Sampling
# export Parallel
# export Vertical
# export Horizontal
# export Axis
# export Violin

# export Shape
# # export Legend
# # export Annotation

# export Coloring
# export Blending

# export Plot

# # include("definitions.jl")
# # include("utils.jl")

# # export scale_height!
# # export scale_width!
# # export pick_random
# # export pick_scenario
# # export lower_triangular
# # export combine_over

end


text:
# The abstract types provided by InteratomicPotentials.jl

export AbstractPotential, NonTrainablePotential, TrainablePotential, EmpiricalPotential, MixedPotential

"""
    AbstractPotential

The abstract supertype of all interatomic potentials.
"""
abstract type AbstractPotential end
include("types/abstract_potential.jl")

"""
    NonTrainablePotential <: AbstractPotential

Abstract type for potentials that aren't trainable.
"""
abstract type NonTrainablePotential <: AbstractPotential end
include("types/non_trainable_potential.jl")

"""
    TrainablePotential{P<:NamedTuple,HP<:NamedTuple} <: AbstractPotential

Abstract type for potentials that are trainable.
`P` is a `NamedTuple` of parameter names and `HP`` is a `NamedTuple`` of hyperparameter names.
"""
abstract type TrainablePotential{P<:NamedTuple,HP<:NamedTuple} <: AbstractPotential end
include("types/trainable_potential.jl")

"""
    EmpiricalPotential{P<:NamedTuple,HP<:NamedTuple} <: TrainablePotential{P,HP}

Defines an empirical potential, a heuristic function used to describe the intermolecular potential energy of a configuration of atoms. Various potentials have been found to agree empirically with the experimentally obtained potential energy of a configuration of atoms for particular atoms in particular situations. This package implements a number of such popular potentials.

`P` is a `NamedTuple` of parameter names and `HP`` is a `NamedTuple`` of hyperparameter names.
"""
abstract type EmpiricalPotential{P<:NamedTuple,HP<:NamedTuple} <: TrainablePotential{P,HP} end
include("types/empirical_potential.jl")

"""
    MixedPotential <: AbstractPotential

Abstract type for potentials that are the combination of multiple sub-potentials.
"""
abstract type MixedPotential <: AbstractPotential end
include("types/mixed_potential.jl")


text:
using DifferentialEquations
using Plots
using LinearAlgebra

function DeQ(du,u,p,t)
du[1]=-u[2]*(1+u[1]^2) / ( (u[1]^2+u[2]^2)*(1-u[2]^2));
du[2]=u[1]*(1-u[2]^2)/ ( (u[1]^2+u[2]^2)*(1+u[1]^2));
du[3]= 2*u[1]*u[2]/(u[1]^2+u[2]^2)^2 *u[3]-1/(u[2]^2+u[1]^2)^2* (1-u[2]^2)*(3*u[1]^4+u[2]^2*u[1]^2-u[2]^2+u[1]^2) /( (1+u[1]^2)^2 )*u[4];
du[4]=-(1/(u[1]^2+u[2]^2)^2* (1+u[1]^2)*(3*u[2]^4+u[1]^2*u[2]^2+u[1]^2-u[2]^2)/( (1-u[2]^2)^2) )*u[3]-(2*u[1]*u[2]/(u[1]^2+u[2]^2)^2)*u[4];
du[5]= 2*u[1]*u[2]/(u[1]^2+u[2]^2)^2 *u[5]-1/(u[2]^2+u[1]^2)^2* (1-u[2]^2)*(3*u[1]^4+u[2]^2*u[1]^2-u[2]^2+u[1]^2) /( (1+u[1]^2)^2 )*u[6];
du[6]=-(1/(u[1]^2+u[2]^2)^2* (1+u[1]^2)*(3*u[2]^4+u[1]^2*u[2]^2+u[1]^2-u[2]^2)/( (1-u[2]^2)^2) )*u[5]-(2*u[1]*u[2]/(u[1]^2+u[2]^2)^2)*u[6];
end
condition(u,t,integrator)= u[1]
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)

#
H=range(.05,stop=.2,length=10)
# for k=1:length(H)
H=.125
    h=1/(2*H)
    Y0=1.0/sqrt(h+1)
    tspan0=(0.,1e4)
    u0=[0.,Y0,0.,0.,0.,0.]
    prob = ODEProblem(DeQ,u0,tspan0)
    t,A1=solve(prob, Vern9(),reltol=1e-12,abstol=1e-20,maxiters=1e15,callback=cb)
    FinalT=t.t[end]

    u0=[0.,Y0,1.,0.,0.,1.]
    tspan=(0.,FinalT)
    prob = ODEProblem(DeQ,u0,tspan)
    t,A=solve(prob, Vern9(),reltol=1e-12,abstol=1e-12,maxiters=1e15)
    N=length(t.t)
    r=zeros(N)
    z=zeros(N)
    theta=zeros(N)
    for j=1:N
        #Create Monodromy matrix
        a=A[3,j]
        b=A[4,j]
        c=A[5,j]
        d=A[6,j]
        M=[a c; b d]
        W,S,V = svd(M)
        P=W*Diagonal(S)*W'
        U=W*V'
        r[j]=P[1,1]
        z[j]=P[1,2]
        theta[j]=asin(U[2,1])
    end
    x=r.*cos.(theta)
    y=r.*sin.(theta)

plot3d(x,y,z)


text:
using Test
using TestSetExtensions
using ReferenceTests
using MetaBalls:MetaBalls

include("generators.jl")
include("c_code.jl")

@test_reference "voxels.txt" map((x) -> edge_lookup(x), generate_voxels(16))

@testset "edge lookup" begin
    @test_reference "voxels.txt" map((x) -> MetaBalls.edge_lookup(x), generate_voxels(16))
end


text:
module ReservoirComputing

using SparseArrays
using LinearAlgebra
using MLJLinearModels
using LIBSVM
using GaussianProcesses
using Optim
using Distributions
using Statistics
using Distances

abstract type AbstractReservoirComputer end
abstract type AbstractEchoStateNetwork <: AbstractReservoirComputer end
abstract type NonLinearAlgorithm end

include("leaky_fixed_rnn.jl")
include("train.jl")
export ESNtrain, Ridge, Lasso, ElastNet, RobustHuber, HESNtrain

include("nla.jl")
export nla, NLADefault, NLAT1, NLAT2, NLAT3

include("esn_input_layers.jl")
export init_input_layer, init_dense_input_layer, init_sparse_input_layer, min_complex_input, irrational_sign_input, physics_informed_input
include("esn_reservoirs.jl")
export init_reservoir_givendeg, init_reservoir_givensp, pseudoSVD, DLR, DLRB, SCR, CRJ

include("echostatenetwork.jl")
export ESN, ESNpredict, ESNpredict_h_steps

include("dafesn.jl")
export dafESN, dafESNpredict, dafESNpredict_h_steps

include("svesm.jl")
export SVESMtrain, SVESM_direct_predict, SVESMpredict, SVESMpredict_h_steps

include("esgp.jl")
export ESGPtrain, ESGPpredict, ESGPpredict_h_steps

include("ECA.jl")
export ECA
#include("reca.jl")
#export RECA, reca_predict
include("reca_discrete.jl")
export RECA_discrete, RECAdirect_predict_discrete
include("gameoflife.jl")
export GameOfLife
include("reca_gol.jl")
export RECA_TwoDim, RECATDdirect_predict_discrete, RECATD_predict_discrete

include("rmm.jl")
export RMM, RMMdirect_predict

include("gruesn.jl")
export GRUESN, GRUESNpredict

include("hesn.jl")
export HESN, HESNpredict

end #module


text:
module JuliaLambdaExample

# Calculate π using Monte Carlo simulation
function monte_carlo_pi(n)
    inside((x, y)) = x * x + y * y < 1.0
    points = rand(2, n)
    hits = count(inside(points[:, v]) for v in 1:n)
    return hits / n * 4
end

# Parsing input (this is more commonly structured as JSON string)
parse_data(event_data) = parse(Int, event_data)

# Formatting result (this is more commonly structured as JSON string)
format_result(value) = "$value"

# Define a handler function that is called by the lambda runtime
function handle_event(event_data, headers)
    @info "Handling request" event_data headers
    simulations = parse_data(event_data)
    my_pi = monte_carlo_pi(simulations)
    return format_result(my_pi)
end

end

max_stars_repo_name: ramboldio/aws-lambda-container-julia

max_issues_repo_name: ramboldio/aws-lambda-container-julia

max_forks_repo_name: ramboldio/aws-lambda-container-julia


text:
const UNINITIALIZED_DATETIME = Dates.DateTime(Dates.Minute(0))
const UNINITIALIZED_PERIOD = Dates.Period(Dates.Minute(0))
const UNINITIALIZED_HORIZON = 0

"""Container for forecasts and their metadata.."""
mutable struct ForecastMetadata
    resolution::Dates.Period
    horizon::Int64
end

function ForecastMetadata()
    resolution = UNINITIALIZED_PERIOD
    horizon = UNINITIALIZED_HORIZON

    return ForecastMetadata(resolution, horizon)
end

function reset_info!(forecasts::ForecastMetadata)
    forecasts.resolution = UNINITIALIZED_PERIOD
    forecasts.horizon = UNINITIALIZED_HORIZON
    @info "Reset system forecast information."
end

function is_uninitialized(forecasts::ForecastMetadata)
    return forecasts.resolution == UNINITIALIZED_PERIOD &&
           forecasts.horizon == UNINITIALIZED_HORIZON
end

function _verify_forecast(metadata::ForecastMetadata, forecast::ForecastInternal)
    if forecast.resolution != metadata.resolution
        throw(DataFormatError(
            "Forecast resolution $(forecast.resolution) does not match system " *
            "resolution $(metadata.resolution)",
        ))
    end

    if get_horizon(forecast) != metadata.horizon
        throw(DataFormatError(
            "Forecast horizon $(get_horizon(forecast)) does not match system horizon " *
            "$(metadata.horizon)",
        ))
    end
end

function check_add_forecast!(metadata::ForecastMetadata, forecast::ForecastInternal)
    if is_uninitialized(metadata)
        # This is the first forecast added.
        metadata.horizon = get_horizon(forecast)
        metadata.resolution = forecast.resolution
    end

    # This will throw if something is invalid.
    _verify_forecast(metadata, forecast)
end

function ForecastMetadata(data::NamedTuple)
    resolution = JSON2.read(JSON2.write(data.resolution), Dates.Period)
    horizon = data.horizon
    return ForecastMetadata(initial_time, resolution, horizon)
end

"""Return the horizon for all forecasts."""
get_forecasts_horizon(forecasts::ForecastMetadata)::Int64 = forecasts.horizon

"""Return the resolution for all forecasts."""
get_forecasts_resolution(forecasts::ForecastMetadata)::Dates.Period = forecasts.resolution


text:
module Landau1P

using NumericalIntegration

"""
    generate_solver(f₀, f, v_, Δt, A)

Generates solver for Landau integral in 1-dimensional velocity space. ``A = 4πe⁴N₀/(m²v₀³ω)`` is a normalizing constant.
"""
function generate_solver(f₀, f, v_, Δt, A)
    v = copy(v_)
    L = 20.0 # Coulomb logarithm
    Tₜ = 1e-3 # transversal temperature needed to estimate |vₛ-vₚ|² - (vₛₓ-vₚₓ)² as ≈2Tₜ

    I = similar(f)
    J = similar(v)

    function solve!()
        for i in eachindex(f)
            if i == 1
                Δfᵢ = (f₀[2]-f₀[1])/(v[2]-v[1])
            elseif i == length(f)
                Δfᵢ = (f₀[end]-f₀[end-1])/(v[end]-v[end-1])
            else
                Δfᵢ = (f₀[i+1]-f₀[i-1])/(v[i+1]-v[i-1])
            end
            for j in eachindex(f₀)
                if i != j
                    if i == 1
                        Δfⱼ = (f₀[2]-f₀[1])/(v[2]-v[1])
                    elseif i == length(f)
                        Δfⱼ = (f₀[end]-f₀[end-1])/(v[end]-v[end-1])
                    else
                        Δfⱼ = (f₀[i+1]-f₀[i-1])/(v[i+1]-v[i-1])
                    end
                    J[j] = (f₀[i]*Δfⱼ - f₀[j]*Δfᵢ)*2Tₜ/abs(v[i]-v[j])^3
                end
            end
            I[i] = L*A*integrate(v, J)
        end

        for i in eachindex(f)
            if i == 1
                f[1] = f₀[1] - Δt*(I[2] - I[1])/(v[2] - v[1])
            elseif i == length(f)
                f[end] = f₀[end] - Δt*(I[end] - I[end-1])/(v[end] - v[end-1])
            else
                f[i] = f₀[i] - Δt*(I[i+1] - I[i-1])/(v[i+1] - v[i-1])
            end
        end
    end

    return solve!
end
end  # module Landau1P


text:
module MH

using Distributions
using Random

function ff(x::T) where T <: AbstractFloat
    x /= 100.0
    pdf(Normal(0, 1), x) + pdf(Normal(3, 1), x) + pdf(Normal(6, 1), x)
end

function mh_serial_naive(x₀::T, N::Integer, burn_N::Integer, f::Function) where T <: AbstractFloat
    x = Vector{T}(undef, N)
    xₜ = x₀::T
    σ = one(T)

    @inline function _gen_candidate(xₜ::Y, σ::Y, f::Function) where Y <: AbstractFloat
        g = Normal(xₜ, σ)
        xc = rand(g)
        f_xc = f(xc)
        α = f_xc / f(xₜ)
        xc, α
    end

    target = 0.3
    accepted = 0
    f_xₜ = f(xₜ)
    for i = 1:burn_N
        xc, α = _gen_candidate(xₜ, σ, f)
        u = rand(Uniform(0, 1))
        if (u <= α)
            xₜ = xc
            accepted += 1
        end
        σ += convert(T, 1 / σ * 1000( (accepted / i) - target))
    end

    for t = 1:N
        xc, α = _gen_candidate(xₜ, σ, f)
        u = rand(Uniform(0, 1))
        if (u <= α)
            xₜ = xc
        end
        x[t] = xₜ
    end

    return x
end

function mh_serial_optimized(
    x₀::T,
    N::Integer,
    burn_N::Integer,
    f::Function,
    seed::Integer
) where {T <: AbstractFloat}
    # optymalizacja algorytmu poprzez zapamiętywanie poprzednich ewaluacji funkcji aproksymującej
    # optymalizacja wykonania przez jawne podanie generatora liczb losowych (globalny bywa wolniejszy)

    x = Vector{T}(undef, N)
    xₜ = x₀::T
    σ = one(T)

    rng = MersenneTwister(seed)

    @inline function _gen_candidate(xₜ::Y, σ::Y, f_xₜ, f::Function, rng::RNG) where {Y <: AbstractFloat,RNG <: AbstractRNG}
        g = Normal(xₜ, σ)
        xc = rand(rng, g)
        f_xc = f(xc)
        α = f_xc / f_xₜ
        xc, α, f_xc
    end

    target = 0.3
    accepted = 0
    f_xₜ = f(xₜ)
    for i = 1:burn_N
        xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
        u = rand(rng, Uniform(0, 1))
        if (u <= α)
            xₜ = xc
            f_xₜ = f_xc
            accepted += 1
        end
        σ += convert(T, 1 / σ * 1000( (accepted / i) - target))
    end

    for t = 1:N
        xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
        u = rand(rng, Uniform(0, 1))
        if (u <= α)
            xₜ = xc
            f_xₜ = f_xc
        end
        x[t] = xₜ
    end

    return x
end

function mh_threaded_naive(
    x₀::T,
    N::Integer,
    burn_N::Integer,
    f::Function, seed
) where {T <: AbstractFloat}

    x = Vector{T}(undef, N)
    xₜ = x₀::T
    σ = one(T)

    rng = MersenneTwister(seed)

    @inline function _gen_candidate(xₜ::Y, σ::Y, f_xₜ, f::Function, rng::RNG) where {Y <: AbstractFloat,RNG <: AbstractRNG}
        g = Normal(xₜ, σ)
        xc = rand(rng, g)
        f_xc = f(xc)
        α = f_xc / f_xₜ
        xc, α, f_xc
    end

    target = 0.3
    accepted = 0
    f_xₜ = f(xₜ)
    for i = 1:burn_N
        xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
        u = rand(rng, Uniform(0, 1))
        if (u <= α)
            xₜ = xc
            f_xₜ = f_xc
            accepted += 1
        end
        σ += convert(T, 1 / σ * 1000( (accepted / i) - target))
    end

    Threads.@threads for t = 1:N
        xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
        u = rand(rng, Uniform(0, 1))
        if (u <= α)
            xₜ = xc
            f_xₜ = f_xc
        end
        x[t] = xₜ
    end

    return x
end

function mh_threaded_optimized(
    x₀::T,
    N::Integer,
    burn_N::Integer,
    f::Function,
    seed::Integer
) where {T <: AbstractFloat}

    @inline function _gen_candidate(xₜ::T, σ::T, f_xₜ, f, rng) where T <: AbstractFloat
        g = Normal(xₜ, σ)
        xc = rand(rng, g)
        f_xc = f(xc)
        α = f_xc / f_xₜ
        xc, α, f_xc
    end

    function __generate_loop(x::Vector{T}, xₜ::T, σ::T, range::R, f, rng::RNG) where {T <: AbstractFloat,R <: AbstractRange,RNG <: AbstractRNG}
        f_xₜ = f(xₜ)
        for t in range
            xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
            u = rand(rng, Uniform(0, 1))
            if (u <= α)
                xₜ = xc
                f_xₜ = f_xc
            end
            @inbounds x[t] = xₜ
        end
    end

    function __burn_loop(xₜ::T, σ::T, burn_N::Integer, target, f, rng::RNG) where  {T <: AbstractFloat,RNG <: AbstractRNG}
        accepted = 0
        f_xₜ = f(xₜ)
        for i = 1:burn_N
            xc, α, f_xc = _gen_candidate(xₜ, σ, f_xₜ, f, rng)
            u = rand(rng, Uniform(0, 1))
            if (u <= α)
                xₜ = xc
                f_xₜ = f_xc
                accepted += 1
            end
            σ += convert(T, 1 / σ * 1000( (accepted / i) - target))
        end
        σ, xₜ
    end

    x = Vector{T}(undef, N)
    xₜ = x₀::T
    σ = one(T)

    rng = MersenneTwister(seed)

    TARGET = 0.3
    σ, xₜ = __burn_loop(xₜ, σ, burn_N, TARGET, f, rng)

    N_THR = Threads.nthreads()

    ranges = collect(Iterators.partition(1:N, (N + N_THR) ÷ N_THR))

    _task = (id) -> Threads.@spawn __generate_loop(x, $xₜ, $σ, $ranges[id], f,  MersenneTwister(rng.seed[1] * id))
    wait.(map(id -> _task(id), 1:N_THR))

    return x
end

export ff, mh_serial_naive, mh_serial_optimized, mh_threaded_naive, mh_threaded_optimized

end


text:
using StanSample, DiffEqBayesStan, MCMCChains
using OrdinaryDiffEq
using ParameterizedFunctions
using RecursiveArrayTools
using Distributions
using StatsPlots

function SIR!(du,u,pars,t)
    β, γ = pars
    S, I, R = u

    dS = -β*S*I
    dI = β*S*I - γ*I
    dR = γ*I
    du .= [dS, dI, dR]
end

t_obs = collect(0:1.0:218)
prob = ODEProblem(SIR!, [0.99, 0.01, 0.0], (t_obs[1], t_obs[end]), [0.20, 0.15])
sol = solve(prob, Tsit5(), reltol = 1e-6, abstol = 1e-6, saveat = t_obs)
y = DataFrame(sol',[:S,:I,:R])
plot(t_obs, y."I")


text:
"""
Compute the stopping times of numbers subjected to the Collatz sequence.
https://en.wikipedia.org/wiki/Collatz_conjecture
"""

include("../ragamuffin.jl")

function collatz_length(x)
    if x == 1
        return 0
    elseif x % 2 == 0
        return 1 + collatz_length(x / 2)
    else
        return 1 + collatz_length(3 * x + 1)
    end
end
N = 100
range = 1:N
expected = @time collect(map(collatz_length, range))

function ragamuffin_collatz()
    T = Int32
    program = initialize(map(T, range) |> collect, T(-1), Frame, 32, T)
    add!(program, :root, (:end, :increment), value_function=x -> x.value == 1)
    add!(program, :increment, :parity, dest_function = x -> 1)
    add!(program, :parity, (:even, :odd), value_function = x -> x.value % 2 == 0)
    add!(program, :even, value_function = x -> x.value > 0 ? Cell{T}(T(x.value / 2), x.depth, x.output_index) : x)
    add!(program, :odd, value_function = x -> Cell{T}(T(x.value * 3 + 1), x.depth, x.output_index))
    @time run(program)
    return program.dest
end

batched = ragamuffin_collatz()
@info batched
@assert expected == batched


text:
# [[file:~/Dropbox/Research/TetradTurbulence/MLCoarseGrainedVGT/docs/Article/article.org::*Problem Struct Definition][Problem Struct Definition:1]]
# unfortunately order does matter here, so be careful when modifying
mutable struct LearningProblem

    #--------------set by operator--------------#
    initialConditions;#::Array{Float64, 2}; #shape=(18,trajectories)
    trainingData;#::Array{Float64, 3}; #shape = (18, length(tsteps), trajectories)

    x;#::Array{Float64, 2}; #shape=(numdims, numTrajs*(numSteps-1)
    y;#::Array{Float64, 2}; #shape=(numdims, numTrajs*(numSteps-1)

    validationData;#::Array{Float64, 3};
    tsteps::Array{Float64, 1};
    miniBatchSize::Int;
    maxiters::Int;
    optimizer; #e.g., GradientDescent(), SGD(), etc.

    trueDrift_; #should take du, u, p, t as input
    trueDiff_;  # and modify du, return nothing

    encode_; #takes data of size(18,numTsteps,numTrajs) and outputs size(reducedDims, numTsteps, numTrajs)
    decode_; #inverse of encode

    modelDrift_; #should take u,p,t as params
            # and return tangent of length(u)

    modelDiff_; #should take u,p,t as params
           # and return noise of length(u)

    fullSystemDrift_; #encodes, applies drift, decodes. Mainly for post-processing.

    loss_;
    lossAndGrads_; #should take x,y,p as inputs
           # and return the loss, and gradient w.r.t. parameters p

    params::Array{Float64, 1};

    #parameter lengths
    driftLength;
    diffLength;

    #functions called by modelDrift, modelDiff_ respectively
    parameterizedDriftClosure_;
    parameterizedPrecisionMat_;

    #--------------set by training--------------#
    lossArray::Array{Float64, 1};

    #--------------set by post-processing--------------#
    plts; #dictionary

end

#--------------constructor for use by operator--------------#
function LearningProblem(initialConditions,
                         trainingData,#::Array{Float64, 3},
                         x,#::Array{Float64, 2},
                         y,#::Array{Float64, 2},
                         validationData,#::Array{Float64, 3},
                         tsteps::Array{Float64, 1},
                         miniBatchSize::Int,
                         maxiters::Int,
                         optimizer,
                         trueDrift_,
                         trueDiff_,
                         encode_,
                         decode_,
                         modelDrift_,
                         modelDiff_,
                         fullSystemDrift_,
                         loss_,
                         lossAndGrads_,
                         initialParams::Array{Float64, 1},
                         driftLength,
                         diffLength,
                         parameterizedDriftClosure_,
                         parameterizedPrecisionMat_)

    return LearningProblem(initialConditions,
                           trainingData,
                           x,
                           y,
                           validationData,
                           tsteps,
                           miniBatchSize,
                           maxiters,
                           optimizer,
                           trueDrift_,
                           trueDiff_,
                           encode_,
                           decode_,
                           modelDrift_,
                           modelDiff_,
                           fullSystemDrift_,
                           loss_,
                           lossAndGrads_,
                           initialParams,
                           driftLength,
                           diffLength,
                           parameterizedDriftClosure_,
                           parameterizedPrecisionMat_,
                           zeros(maxiters), #lossArray
                           Dict{String, Plots.Plot{Plots.GRBackend}}()); #plots
end

#delete large allocations before saving problem to disk.
function myserialize(path, l::LearningProblem)
    l.initialConditions = nothing;
    l.trainingData = nothing;
    l.x = nothing;
    l.y = nothing;
    serialize(path, l);
end

# used inside of problem setup to condense along time dimension to provide training data
# of size (numDims, numSamples=numTsteps*numTrajs) and to generate ground truth tangents
# to compare learning against
function parseTrainingData(data, dt)
    numDims, numSteps, numTrajs = size(data);

    x = zeros(numDims, numTrajs*(numSteps-1)); #the last time-step is not a valid sample
    y = [(data[:,j+1,i] - data[:,j,i])/dt for i in 1:numTrajs for j in 1:(numSteps-1)]; #generate array of array of tangents
    y = hcat(y...); #y is 9900 arrays of 18 each. We want size(y)=(18,9900)

    stride = numSteps-1;
    for i in 1:numTrajs
        for j in 1:(numSteps-1)
            x[:,stride*(i-1)+(j-1)+1] = data[:,j,i]; #the perils of 1-based indexing.
            #I suppose i,j could start at 0, and then add one on RHS.
        end
    end

    return x,y;
end
# Problem Struct Definition:1 ends here

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
using StatsBase
using Flux

function sample_keys(n, dicts...)
    common_keys = intersect((collect(keys(d)) for d in dicts)...)
    sample(common_keys, n, replace=false)
end

function sample_dicts(n, dicts...)
    sampled_keys = sample_keys(n, dicts...)
    Tuple(Dict(k => d[k] for k in sampled_keys) for d in dicts)
end

function sample_dicts!(n, dicts...)
    sampled_keys = Set(sample_keys(n, dicts...))

    for dict in dicts, key in keys(dict)
        if !(key in sampled_keys)
            delete!(dict, key)
        end
    end
end

function roll_window(series, input_size, output_size)
    inputs = Vector()
    outputs = Vector()

    for i in input_size:length(series) - output_size
        push!(inputs, series[i + 1 - input_size:i])
        push!(outputs, series[i + 1:i + output_size])
    end

    reduce(hcat, inputs), reduce(hcat, outputs)
end


text:
addprocs(2)
@everywhere require("Images.jl")

# Issue #287
using Images
import Images.RGB
@everywhere function test287(img)
    return 0;
end
Imgs = Array(Image{RGB{Float64}}, 2);
Imgs[1] = convert(Image{RGB}, rand(100, 100, 3));
Imgs[2] = convert(Image{RGB}, rand(100, 100, 3));

let Imgs = Imgs;
    ret = pmap(i -> test287(Imgs[i]), 1:2);
    @test ret == Any[0,0]
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "gemm" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = sparse(rand(elty,m,k))
        B = sparse(rand(elty,k,n))
        C = A * B
        d_A = CudaSparseMatrixCSR(A)
        d_B = CudaSparseMatrixCSR(B)
        d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
        r_r = to_host(d_C.rowPtr)
        r_c = to_host(d_C.colVal)
        r_v = to_host(d_C.nzVal)
        h_C = to_host(d_C)
        @test C ≈ h_C
        @test_throws DimensionMismatch CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O')
        #=A = sparse(rand(elty,m,k))
        B = sparse(rand(elty,k,n))
        d_A = CudaSparseMatrixCSC(A)
        d_B = CudaSparseMatrixCSC(B)
        C = A * B
        d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
        h_C = to_host(d_C)
        @test_approx_eq(C,h_C)
        @test_throws(DimensionMismatch,CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O'))=#
    end
end


text:
let
    include("../src/KiteViewer.jl")
    Utils.test(true)
    main(false)
end

@info "Precompile script has completed execution."


text:
module bwh

using Parameters
using Statistics
using Plots
using ImageFiltering
using DiffEqOperators, SparseArrays, LinearAlgebra
using DifferentialEquations
using FFTW
using BenchmarkTools
using Printf
using DelimitedFiles

include("Params.jl")
include("laplacian.jl")
include("convapprox.jl")
include("rhs.jl")
include("integrate.jl")
include("plotbwh.jl")
include("main.jl")

export main, Params, integrate, plotbwh

end #module

max_stars_repo_name: jhardenberg/bwh

max_issues_repo_name: jhardenberg/bwh

max_forks_repo_name: jhardenberg/bwh


text:
import CSV

function get_detailled_crs_data(path::String, filebase::String)

    # select all relevant files
    files = path .* filter!(s->occursin(r"^" * "$filebase", s), readdir(path))

    # read files
    df = mapreduce(CSV.read, vcat, files)

    # delete all empty rows
    CSV.dropmissing!(df)

    return sort!(df)
end

function get_detailled_shock_data(path::String, selected_id::Integer=-1)
    if (selected_id == -1)
        return get_detailled_crs_data(path, "CRs_details_shock_injection")
    else
        df = get_detailled_crs_data(path, "CRs_details_shock_injection")
        return sort!(df[df[:, :ID] .== selected_id,:])
    end
end

function get_detailled_Dpp_data(path::String, selected_id::Integer=-1)
    if (selected_id == -1)
        return get_detailled_crs_data(path, "CRs_details_Dpp")
    else
        df = get_detailled_crs_data(path, "CRs_details_Dpp")
        return sort!(df[df[:, :ID] .== selected_id,:])
    end
end

function get_detailled_radiative_data(path::String, selected_id::Integer=-1)
    if (selected_id == -1)
        return get_detailled_crs_data(path, "CRs_details_radiative")
    else
        df = get_detailled_crs_data(path, "CRs_details_radiative")
        return sort!(df[df[:, :ID] .== selected_id,:])
    end
end

function get_detailled_adiabatic_data(path::String, selected_id::Integer=-1)
    if (selected_id == -1)
        return get_detailled_crs_data(path, "CRs_details_adiabatic")
    else
        df = get_detailled_crs_data(path, "CRs_details_adiabatic")
        return sort!(df[df[:, :ID] .== selected_id,:])
    end
end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
"""
SDERGM DGP  for density and GWESP as described in the paper
Heavily relying on R to smaple from the ERGM at each t, and it assumes that
all needed packages have been loaded in R via Rcall
"""
function sdergmDGPpaper_Rcall( Model::DynNets.SdErgmDirBinGlobalPseudo, indTvPar::BitArray{1},T::Int,N::Int;vConstPar ::Array{<:Real,1} = zeros(Real,2), ftot_0::Array{<:Real,1} = zeros(Real,2),UM = [-3 , 0.05],B_Re  = 0.9, A_Re  = 0.01 )

    indTvPar = trues(2)
    ftot_0 = zeros(Real,2)
    indTargPar = falses(2)
    NTargPar = 0

    NergmPar = 2#
    NTvPar   = sum(indTvPar)

    function initialize_pars(vParOptim_0)
        last = 0
        for i=1:NergmPar
            if indTvPar[i]
                if indTargPar[i]
                    vParOptim_0[last+1:last+2] = [ B_Re; A_Re]
                    last+=2
                else
                    vParOptim_0[last+1:last+3] = [UM[i]*(1-B_Re) ; B_Re; A_Re]
                    last+=3
                end
            else
                vParOptim_0[last+1] = UM[i]
                last+=1
            end
        end
        return vParOptim_0
    end
    vResGasPar = initialize_pars(zeros(NergmPar + NTvPar*2 -NTargPar))
    # Organize parameters of the GAS update equation
    Wvec = vResGasPar[1:3:3*NTvPar]
    Bvec = vResGasPar[2:3:3*NTvPar]
    Avec = vResGasPar[3:3:3*NTvPar]

    # start values equal the unconditional mean,and  constant ones remain equal to the unconditional mean, hence initialize as:
    UMallPar = zeros(Real,NergmPar)
    UMallPar[indTvPar] =  Wvec ./ (1 .- Bvec)
    if !prod(indTvPar) # if not all parameters are time varying
        UMallPar[.!indTvPar] = vConstPar
    end
    #println(UMallPar)
    fVecT = ones(Real,NergmPar,T)
    obsT = fill(zeros(2,2),T)

    statsT = ones(Real,NergmPar,T)

    sum(ftot_0)==0  ?    ftot_0 = UMallPar : ()# identify(Model,UMallNodesIO)

    ftot_t = copy(ftot_0)
    if NTvPar==0
        I_tm1 = ones(1,1)
    else
        I_tm1 = Float64.(Diagonal{Real}(I,NTvPar))
    end
    formula_ergm_str = " edges +  gwesp(decay = 0.25,fixed = TRUE)"
    ftot_t
    for t=1:T

        ftot_t_R = ftot_t# cat(,ftot_t,dims = 2)
        # Sample the observations using RCall -----------------------
        @rput ftot_t_R ; @rput N;
        reval("formula_ergm = net ~ "*formula_ergm_str)
         #create an empty network, the formula defining ergm, sample the ensemble and
         # store the sufficient statistics and change statistics in R
         R"
         net <- network.initialize(N)
         net <- simulate(formula_ergm, nsim = 1, seed = sample(1:100000000,1), coef = as.numeric(ftot_t_R), control = control.simulate.formula(MCMC.burnin = 100000))
         sampledMat_t_R <- as.matrix.network( net)
         chStat_t <- ergmMPLE(formula_ergm)
         changeStats_t_R <- cbind(chStat_t$response, chStat_t$predictor,chStat_t$weights)
         stats_t_R <- summary(formula_ergm)
        "
        # import sampled networks in julia
        sampledMat_t = BitArray( @rget(sampledMat_t_R))
        changeStats_t = @rget changeStats_t_R;# tmp = Array{Array{Float64,2}}(T); for t=1:T tmp[t] =  changeStats_T[t];end;changeStats_T = tmp
        #print(size(changeStats_t))
        stats_t = @rget(stats_t_R); #tmp = zeros(Nterms,T); for t=1:T tmp[:,t] = stats_T[t];end;stats_T = tmp

        #-----------------------------------------------------------------
        obs_t = changeStats_t # vector of in and out degrees
        obsT[t] = obs_t
        statsT[:,t] = stats_t
        #print((t,I_tm1))
        ftot_t,loglike_t,I_tm1,~ = DynNets.predict_score_driven_par(Model,obs_t,ftot_t,I_tm1,indTvPar,Wvec,Bvec,Avec)
        fVecT[:,t] = ftot_t #store the filtered parameters from previous iteration
    end

    return obsT,statsT,fVecT
end
export sdergmDGPpaper_Rcall


text:
using Statistics
using Plots
using BenchmarkTools
using Profile
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile, DiffResults
using LinearAlgebra
using Zygote
using Deriv_Utils_New
using Data_Utils
push!(LOAD_PATH, pwd()*"/scratch_NM")
using LossFuncs

#Image Cases
Nx=16
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
img = zeros(Float64, Nx, Nx)
img[6,6]=1.0

#ds1s2 = wst_S1S2_derivfast(img, fhash)
#Create wrappers
function DHC_S2(input)
    return DHC(input, fhash, doS2=true)
end

function proxy(input)
    return [input[1]+input[2], input[1]*input[2], input[1]/input[2]]
end

fjac = Zygote.pullback(proxy, [3, 4])

arr = rand(256, 256)
@time foo1 = arr[:]
@time foo2 = reshape(arr, (256*256)) #100x faster

#0.000214 seconds (2 allocations: 512.078 KiB)
#0.000004 seconds (2 allocations: 80 bytes)

img = readdust(64)
fhash = fink_filter_hash(1, 8, nx=64, pc=1, wd=1, Omega=true)
@benchmark ds20 = Deriv_Utils_New.wst_S20_deriv(img, fhash) #157ms
@benchmark ds12 = Deriv_Utils_New.wst_S12_deriv(img, fhash) #3.28 s

Profile.clear()
@profile ds12 = Deriv_Utils_New.wst_S12_deriv(img, fhash)
Juno.profiler()

pnzlist = Array{Array, 2}

function intersect_p(fhash)
    f_ind   = fhash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = fhash["filt_value"]

    (Nf,) = size(fhash["filt_index"])
    for f1=1:Nf
        f_i1 = f_ind[f1]
        f_v1 =f_val[f1]
        for f2=1:Nf
            f_i2 = f_ind[f2]
            f_v2 =f_val[f2]
            pnz = intersect(f_i1, f_i2)
        end
    end
end

function intersect_p_prealloc(fhash)
    f_ind   = fhash["filt_index"]  # (J, L) array of filters represented as index value pairs
    f_val   = fhash["filt_value"]
    pnz_store = []
    (Nf,) = size(fhash["filt_index"])
    for f1=1:Nf
        f_i1 = f_ind[f1]
        f_v1 =f_val[f1]
        for f2=1:Nf
            f_i2 = f_ind[f2]
            f_v2 =f_val[f2]
            push!(pnz_store, intersect(f_i1, f_i2))
        end
    end
end

@benchmark intersect_p_prealloc(fhash) #42ms, same time  as intersect_p

function foo_wrapperout(var)
    a=var+2
    function foo_wrapperin() #Generic Loss wrapper that can see internal variables but only takes vector as an arg
        LossFuncs.foo(a) #Specific loss function to be used here that has all
    end
    outval = foo_wrapperin()
    return outval
end

varout = foo_wrapperout(10)
#DIFFERENT REVERSEDIFF CONFIGS: Which is the fastest reversediff?
#Zyg.forward_jacobian doesnt work: No method matching DHC.FwdDiff
#Zyg.gradient only works if the output is scalar. Can try to benchmark with the chisq function.
#=
function jacobian(f, x)
    y = f(x)
    n = length(y)
    m = length(x)
    j = Array{Float64, 2}(undef, n, m)
    for i in 1:n
        j[i, :] .= gradient(x -> f(x)[i], x)
    end
    return j
end
=#
#This doesnt work either.
#I wont be surprised if there is no off-the-shelf autodiff for the DHC function--complex intermediaries, Jacobian, FFTs
#and if the way we've implemented this is the only way.

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
# License is MIT: https://github.com/JuliaString/Unicode_Entities/LICENSE.md

using StrTables

const VER = UInt32(1)

const datapath = joinpath(@__DIR__, "..", "data")
const dpath = "https://www.unicode.org/Public/UCD/latest/ucd/"
const inpname = "UnicodeData.txt"
const fname = "unicode.dat"
const disp = [false]

const _empty_string = ""
const _empty_val = (false, ' ', _empty_string, _empty_string)

function process_line(vec::Vector{T}) where {T<:AbstractString}
    length(vec) < 11 && return _empty_val
    num = vec[1]
    str = vec[2]
    alias = vec[11]
    ch = parse_hex(UInt32, num)
    str[1] == '<' &&
        return str == "<control>" ? (alias != "", ch, _empty_string, alias) : _empty_val
    # Don't save names that simply contain hex representation
    len = length(num)
    pos = sizeof(str) - len
    pos > 1 && str[pos] == '-' && str[pos+1:end] == num && return _empty_val
    # Check for some characters we won't represent, all outside of BMP range
    ch <= 0x0ffff && return (true, ch, str, alias)
    # Ignore characters in Linear B range (0x10000-0x100ff)
    0x10000 <= ch < 0x10100 && return _empty_val
    # Ignore characters in Linear A range (0x10600-0x107ff)
    0x10600 <= ch < 0x10800 && return _empty_val
    # Ignore characters in hieroglyph range (0x13000-0x14fff)
    0x13000 <= ch < 0x15000 && return _empty_val
    # Ignore characters in Tangut range (0x17000-0x18fff)
    0x17000 <= ch < 0x19000 && return _empty_val
    # Ignore characters in Greek vocal/instrumental range (0x1d200-0x1d2ff)
    0x1d000 <= ch < 0x1d300 && return _empty_val
    # Don't worry about characters outside of BMP/SMP1
    ch > 0x1FFFF && return _empty_val
    (true, ch, str, alias)
end

function load_unicode_data(datapath, dpath, fname)
    lname = joinpath(datapath, fname)
    if isfile(lname)
        println("Loading Unicode Data: ", lname)
        src = lname
    else
        try
            src = string(dpath, fname)
            println("Downloading Unicode Data: ", src)
            download(src, lname)
            println("Saved to: ", lname)
        catch ex
            println("Unable to download $src and save in $lname")
            rethrow(ex)
        end
    end
    symnam = String[]
    symval = UInt32[]
    aliasnam = String[]
    aliasval = UInt32[]
    count = lines = aliascnt = 0
    open(lname, "r") do f
        while (l = chomp(readline(f))) != ""
            lines += 1
            flg, ch, str, alias = process_line(split(l, ";"))
            disp[] && println('#', lines, '\t', Int(flg), " ", l)
            flg || continue
            if symnam != ""
                count += 1
                push!(symnam, str)
                push!(symval, ch)
            end
            if alias != ""
                aliascnt += 1
                push!(aliasnam, alias)
                push!(aliasval, ch)
            end
        end
    end
    # Check for duplicates
    names = Set{String}(symnam)
    dupcnt = 0
    for (str,ch) in zip(aliasnam, aliasval)
        if str in names
            dupcnt += 1
        else
            push!(symnam, str)
            push!(symval, ch)
        end
    end
    println("Removed ",dupcnt," duplicate aliases")
    println("Finished loading ", count, " + ", aliascnt-dupcnt, " entities on ", lines, " lines")
    symnam, symval, src
end

function split_tables(srtval)
    # BMP characters
    l16 = Tuple{UInt16, UInt16}[]
    # non-BMP characters (in range 0x10000 - 0x1ffff)
    l32 = Tuple{UInt16, UInt16}[]

    for (i, ch) in enumerate(srtval)
        ch > 0x1ffff && error("Character $ch too large")
        push!(ch > 0x0ffff ? l32 : l16, (ch%UInt16, i))
    end

    # We now have 2 vectors, one for single BMP characters, the other for SMP-1 characters
    # each has the value and a index into the name table
    # We need to create a vector the same size as the name table, that gives the index
    # into one of the tables, in order to go from names to the output character
    # We also need, for each of the tables, a sorted vector that goes from the indices
    # in each table to the index into the name table (so that we can find multiple names for
    # each character)

    indvec = create_vector(UInt16, length(srtval))
    vec16, ind16, base32 = sortsplit!(indvec, l16, 0)
    vec32, ind32, base2c = sortsplit!(indvec, l32, base32)

    base32%UInt32, indvec, vec16, ind16, vec32, ind32
end

function update_map!(wrdmap, inp, off, wrd_vec, wrd_dict)
    wrd = wrd_vec[inp]
    srt = sortperm(wrd)
    tab = inp[srt]
    map = wrd_vec[tab]
    for (i, v) in enumerate(map)
        wrdmap[wrd_dict[v]] = (i+off)%UInt16
    end
    map
end

function create_map(wrd_vec, wrd_dict, tab1, tab2)
    wrdmap = zeros(UInt16, length(wrd_vec))
    map1 = update_map!(wrdmap, tab1,  53, wrd_vec, wrd_dict)
    map2 = update_map!(wrdmap, tab2, 255, wrd_vec, wrd_dict)
    wrdmap, map1, map2
end

keepword(str) = !_contains(str, r"^[A-Z0-9\-]+$")

function outseg!(out, str)
    for ch in str
        push!(out, ch == '-' ? 0x01 : (ch%UInt8 - ((ch-'0')<=9 ? 0x2e : 0x35)))
    end
end

function packword(inpvec::Vector, wrdmap, wrd_vec, wrd_dict)
    out = UInt8[]
    hasparts = false
    prevw = 0x0000
    for val16 in inpvec
        w = wrdmap[val16]
        if w > 0x00ff
            push!(out, ((w>>>8)+37)%UInt8, w%UInt8)
        elseif w != 0x0000
            push!(out, w%UInt8)
        else
            str = wrd_vec[val16]
            !isempty(out) && (prevw < 0x26 || str[1] == '-') && push!(out, 0x00)
            if '-' in str
                parts = split(str, '-')
                hasparts = true
                disp[] && print(parts)
                len = length(parts)
                for pos = 1:len
                    seg = parts[pos]
                    if (pwrd = get(wrd_dict, seg, 0)) == 0 || (wp = wrdmap[pwrd]) == 0
                        outseg!(out, seg)
                    elseif wp > 0x00ff
                        push!(out, ((wp>>>8)+37)%UInt8, wp%UInt8)
                    else
                        push!(out, wp%UInt8)
                    end
                    pos != len && push!(out, 0x01)
                end
            else
                outseg!(out, str)
            end
        end
        prevw = w
    end
    hasparts && disp[] && println("\t",out)
    out
end

function split_words(input::Vector{<:AbstractString})
    wrd_dict = Dict{String, Int}()
    wrd_vec  = String[]
    wrd_frq  = Int[]
    wrd_loc  = UInt16[]
    str_vec  = create_vector(Vector{UInt16}, length(input))
    #=
    part_dic = Dict{String, Int}()
    part_vec = String[]
    part_frq = Int[]
    wrd_parts = Vector{UInt16}[]
    =#
    ind = 0
    for (i, wrd) in enumerate(input)
        allwrds = split(wrd, ' ')
        outwrds = UInt16[]
        for onewrd in allwrds
            val = get(wrd_dict, onewrd, 0)
            if val == 0
                val = (ind += 1)
                disp[] && println(val, '\t', i, '\t', onewrd)
                push!(wrd_vec, onewrd)
                push!(wrd_frq, 0)
                push!(wrd_loc, i) # location first found (may be only location)
                wrd_dict[onewrd] = val
                if '-' in onewrd
                    allparts = split(onewrd, '-')
                    for part in allparts
                        part == "" && continue
                        if (vp = get(wrd_dict, part, 0)) == 0
                            vp = (ind += 1)
                            disp[] && println("\tparts:\t", vp, '\t', i, '\t', part)
                            push!(wrd_vec, part)
                            push!(wrd_frq, 0)
                            push!(wrd_loc, i) # location first found (may be only location)
                            wrd_dict[part] = vp
                        end
                        wrd_frq[vp] += 1
                    end
                end
            end
            wrd_frq[val] += 1
            push!(outwrds, val)
        end
        str_vec[i] = outwrds
    end

    # Calculate the savings of each word, i.e. frequency * (length-1), where freq > 1
    # take top 256-(16+38) = 202 words
    wrdsav = sort([((wrd_frq[i]-1)*(sizeof(wrd_vec[i])-1), i)
                   for i = 1:length(wrd_vec) if wrd_frq[i]>1 || keepword(wrd_vec[i])],
                  rev=true)
    # This has indexes into wrd_vec for words that will end up as 1-byte
    table1 = [wrdsav[i][2] for i=1:202]
    # This has indexes into wrd_vec for words that will end up as 2-bytes
    table2 = [wrdsav[i][2] for i=203:length(wrdsav)]
    # Calculate the savings of remaining words, i.e. frequency * (length-2) (some will become 0)
    savfrq = Int[]
    savval = UInt16[]
    for i in table2
        savings = (wrd_frq[i]-1)*(sizeof(wrd_vec[i])-2)
        if savings > 2 || keepword(wrd_vec[i])
            push!(savfrq, savings)
            push!(savval, i)
        end
    end

    # For every word in wrd_vec, create an entry to that has 0-37, 38-53, 54-255, 256-and above
    wrd_map, map1, map2 = create_map(wrd_vec, wrd_dict, table1, savval)

    # Pack words
    ent_map = create_vector(Vector{UInt8}, length(str_vec))
    for (i, vec16) in enumerate(str_vec)
	ent_map[i] = packword(vec16, wrd_map, wrd_vec, wrd_dict)
    end
    PackedTable(ent_map), StrTable(map1), StrTable(map2)
end

function make_tables(datapath, dpath, fname)
    symnam, symval, src = load_unicode_data(datapath, dpath, fname)
    srtind = sortperm(symnam)
    srtnam = symnam[srtind]
    srtval = symval[srtind]
    entmap, map1, map2 = split_words(srtnam)
    base32, indvec, vec16, ind16, vec32, ind32 = split_tables(srtval)
    (VER, string(now()), src, base32, entmap, indvec, map1, map2, vec16, ind16, vec32, ind32)
end

savfile = joinpath(datapath, fname)
if isfile(savfile)
    println("Tables already exist")
else
    tup = nothing
    println("Creating tables")
    try
        global tup
        tup = make_tables(datapath, dpath, inpname)
        println("Saving tables to ", savfile)
        StrTables.save(savfile, tup)
        println("Done")
    catch ex
        println(sprint(showerror, ex, catch_backtrace()))
    end
end


text:
import Base.reverse!

mutable struct BitDataBase
    B            :: Int
    guess        :: Vector{Int}
    n_queries    :: Int
    flipdetect   :: Tuple{Int, Int} # index of the flip detector
    paritydetect :: Tuple{Int, Int} # index of the operation parity detector (either zero/two operations or one operation)
    isknown      :: Int
    function BitDataBase(B::Int)
        new(B, fill(-1, B), 0, (0, -1), (0, -1), false)
    end
end

"""
The error is either
- Unknown  (only at the beginning)
- Z2 flip/noflip (when only similars are available)
- Z2 (id, reverseflip)/(reverse,flip)
- Z2xZ2 fully known
"""
function correcterror!(data::BitDataBase)
    isflipped = false
    if data.flipdetect[1] > 0
        b = query!(data, data.flipdetect[1])
        if b != data.flipdetect[2]
            flip!(data)
            isflipped = true
            data.flipdetect= (data.flipdetect[1], b)
        end
    else
        query!(data, 1)
    end
    if data.paritydetect[1] > 0
        b = query!(data, data.paritydetect[1])
        if b != data.paritydetect[2]
            !isflipped && reverse!(data)
            data.paritydetect = (data.paritydetect[1] ,b)
        else
            isflipped && reverse!(data)
        end
    else
        query!(data, 1)
    end
    data
end

"query at position `n` but don't insert the value"
function query!(data::BitDataBase, n)
    println(n)
    b = parse(Int, readline())
    data.n_queries += 1
    b
end

"query for a single bit at position  `n`"
function queryinsert!(data::BitDataBase, n)
    println(n)
    b = parse(Int, readline())
    data.guess[n] = b
    data.n_queries += 1
    data
end

"query for a single bit at position `n` and the reverse position"
function doublequeryinsert!(data::BitDataBase, n)
    println(n)
    b1 = parse(Int, readline())
    data.guess[n] = b1
    m = data.B-n+1
    println(m)
    b2 = parse(Int, readline())
    data.guess[m] = b2
    data.n_queries += 2
    ## similar pairs detect flip, opposite pairs distinguish between
    ## the parity of applied operators (either 02 or 1) (id,
    ## reverseflip) and (reverse, flip).
    if data.flipdetect[1] == 0 && b1 == b2
        data.flipdetect = (n, b1)
    end
    if data.paritydetect[1] == 0 && b1 != b2
        data.paritydetect = (n, b1)
    end
    data
end

function flip!(data::BitDataBase)
    newguess = fill(-1, data.B)
    for i in eachindex(data.guess)
        if data.guess[i] != -1
            b = data.guess[i]
            newguess[i] = b==0 ? 1 : 0
        end
    end
    data.guess = newguess
end

reverse!(data::BitDataBase) = reverse!(data.guess)

function interact(B)
    data = BitDataBase(B)
    for round = 1:14
        correcterror!(data)
        offset = (round-1)*4
        for i=1:4
            doublequeryinsert!(data, offset+i)
            if offset+i == div(data.B, 2)
                if isodd(data.B)
                    doublequeryinsert!(data, offset+i)
                end
                return data
            end
        end
    end
    data
end

function solution()
    n_cases, B = Tuple(parse(Int, i) for i in split(readline()))
    for x in 1:n_cases
        data = interact(B)
        println(reduce(*, string.(data.guess)))
        flush(stdout)
        answer = readline()
        if answer=="N"
            return
        end
    end
end

solution()


text:
using PkgSRA, Plots
using Turing, AxisArrays, DataFrames
using BSON: @load
pyplot()

#####
##### Load the data
#####
this_dat_name = DAT_FOLDERNAME*"dat_ball_"

# Raw data
fname = this_dat_name*"raw.bson";
@load fname dat numerical_grad dyn_control_kick dyn_control_wall

# Naive model
fname = this_dat_name*"naive_model.bson";
@load fname dat_naive naive_model

# Controlled model
fname = this_dat_name*"controlled_model.bson";
@load fname chain_ctr ctr_final dat_final

#####
##### Produce the plots
#####
plot_opt = Dict(:titlefontsize=>48,
        :xticks=>false, :yticks=>false, :fontfamily=>:serif,
        :legendfontsize=>16, :titlefontsize=>24)

## First figure in two panels: intrinsic dynamics

# First: intrinsic dynamics
plot_data = plot(ts, dat[1,:], label="Height", lw=5,
                legend=:topleft);
    plot!(ts, dat[2,:], label="Velocity", lw=5;
    plot_opt...);
    xlabel!("");
    title!("Observable Data")

# Second: Controller
plot_control = plot(ts, vec(dyn_control_wall[2,:]), label="Ground",
                color=COLOR_DICT["control_true"], lw=3);
    plot!(ts, vec(dyn_control_kick[2,:]), label="Kick",
            color=COLOR_DICT["data"], lw=3,
            legendfontsize=14, xticks=false, yticks=false,
            legend=:topright; plot_opt...);
    xlabel!("Time", fontsize=24);
    title!("Hidden Control Signals")

# Create the layout and plot
my_layout = @layout [p1; p2];
    p_final = plot(plot_data, plot_control, layout = my_layout)

# Save
fname = FIGURE_FOLDERNAME * "fig_ball_intrinsic_and_control.png";
savefig(p_final, fname)

## Second figure in three panels: Model dynamics

# First, naive model dynamics
plot_naive = plot(ts, dat_naive[1,:], lw=5,
                legend=:topright);
    plot!(ts, dat_naive[2,:], legend=false, lw=5;
        plot_opt...);
    xlabel!("");
    title!("Naive Model")

# Second, learned control signals
plot_ctr = plot(ts, ctr_final[2,:], lw=5,
                color=COLOR_DICT["control_time"], legend=false;
                plot_opt...);
    xlabel!("");
    title!("Learned Controller")

# Third, Reconstruction of controlled model
plot_final = plot(ts, dat_final', lw=5,
                # color=COLOR_DICT["model_controlled"],
                legend=false;
                plot_opt...);
    xlabel!("Time", guidefontsize=18);
    title!("Controlled Model")

# Create the layout and plot
my_layout = @layout [plot_naive; plot_ctr; plot_final];
    p_final = plot(plot_naive, plot_ctr, plot_final, layout = my_layout)

# Save
fname = FIGURE_FOLDERNAME * "fig_ball_model.png";
savefig(p_final, fname)


text:
"""
    MHAttention(nheads::Int, qkv_layer, attn_drop, projection)

Multi-head self-attention layer.

# Arguments:
- `nheads`: Number of heads
- `qkv_layer`: layer to be used for getting the query, key and value
- `attn_drop`: dropout rate after the self-attention layer
- `projection`: projection layer to be used after self-attention
"""
struct MHAttention{P, Q, R}
  nheads::Int
  qkv_layer::P
  attn_drop::Q
  projection::R
end

"""
    MHAttention(planes, nheads = 8; qkv_bias = false, attn_drop = 0., proj_drop = 0.)

Multi-head self-attention layer.

# Arguments:
- `planes`: number of input channels
- `nheads`: number of heads
- `qkv_bias`: whether to use bias in the layer to get the query, key and value
- `attn_drop`: dropout rate after the self-attention layer
- `proj_drop`: dropout rate after the projection layer
"""
function MHAttention(planes, nheads = 8; qkv_bias = false, attn_drop = 0., proj_drop = 0.)
  @assert planes % nheads == 0 "planes should be divisible by nheads"
  qkv_layer = Dense(planes, planes * 3; bias = qkv_bias)
  attn_drop = Dropout(attn_drop)
  proj = Chain(Dense(planes, planes), Dropout(proj_drop))

  MHAttention(nheads, qkv_layer, attn_drop, proj)
end

@functor MHAttention

function (m::MHAttention)(x::AbstractArray{T, 3}) where T
  features, len_seq, batch_size = size(x)
  q, k, v = chunk(reshape(m.qkv_layer(x), features ÷ m.nheads, m.nheads, len_seq, 3 * batch_size), 3; dims = 4)
  scale = convert(T, sqrt(size(q, 1) / m.nheads))
  attn = m.attn_drop(softmax(NeuralAttentionlib.matmul(q, permutedims(k, (2, 1, 3, 4))) * scale))
  x = m.projection(reshape(NeuralAttentionlib.matmul(attn, v), (features, len_seq, batch_size)))
end


text:
#=
# 3D Poisson, Dirichlet bc
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("poisson3d");
useLog("poisson3dlog")

n = 10;
ord = 2;

domain(3)
functionSpace(order=ord)

mesh(HEXMESH, elsperdim=n)

u = variable("u")

testSymbol("v")

boundary(u, 1, DIRICHLET, 0)

# Write the weak form
coefficient("f", "-3*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)")
weakForm(u, "-dot(grad(u), grad(v)) - f*v")

solve(u);

outputValues(u, "p3dout", format="vtk", asci=false);

# exact solution is sin(pi*x)*sin(pi*y)*sin(pi*z)
# check error
maxerr = 0;
exact(x,y,z) = sin(pi*x)*sin(pi*y)*sin(pi*z);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    z = Finch.grid_data.allnodes[3,i];
    err = abs(u.values[i] - exact(x,y,z));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

finalize_finch()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
@testset "memory allocator" begin

CuArrays.reclaim(1024)
CuArrays.reclaim()

@test CuArrays.@retry_reclaim(42, return 42) == 42
@test CuArrays.@retry_reclaim(42, return 41) == 41

@test_throws OutOfGPUMemoryError CuArray{Int}(undef, 10^20)
@test_throws OutOfGPUMemoryError CuArrays.extalloc() do
    CuArray{Int}(undef, 10^20)
end

end


text:
# Victoria Park using-server
addprocs(3)
using IncrementalInference
using HDF5, JLD, Gadfly, Colors, Cairo
using KernelDensityEstimate, Distributions
using Caesar, RoME
using RoMEPlotting

## load all the model data
# d = odometry information
# f = laser scanner detections
# MM = multi-modal individual id references
# MMr = reworked to map to only one previous feature
# examplefolder, datafolder

function evalLikelihood(fg::FactorGraph, sym::Symbol, point::Vector{Float64})
  p = getVertKDE(fg, sym)
  Ndim(p) == length(point) ? nothing : error("point (dim=$(length(point))) must have same dimension as belief (dim=$(Ndim(p)))")
  evaluateDualTree(p, (point')')[1]
end

# Evaluate the likelihood of an Array{2} of points on the marginal belief of some variable
# note the dimensions must match
function evalLikelihood(fg::FactorGraph, sym::Symbol, points::Array{Float64,2})
  p = getVertKDE(fg, sym)
  Ndim(p) == size(points,1) ? nothing : error("points (dim=$(size(points,1))) must have same dimension as belief (dim=$(Ndim(p)))")
  evaluateDualTree(p, (points))
end

include(joinpath(Pkg.dir("Caesar"),"examples","wheeled","loadVicPrkData.jl"))

include(joinpath(ENV["HOME"],"Documents","blandauthlocal.jl"))
backend_config, user_config = standardcloudgraphsetup(addrdict=addrdict)
# Start new session
addrdict["sessionId"] = "VICPRK_VID"
addrdict["robotId"] = "Ute"

# Figure export folder
currdirtime = now()
imgdir = joinpath(ENV["HOME"], "Pictures", "vicprkimgs", "$(currdirtime)")
mkdir(imgdir)

fg = Caesar.initfg(sessionname=user_config["sessionId"], robotname=user_config["robotId"], cloudgraph=backend_config)
deleteServerSession!(fg.cg, user_config["sessionId"])

# init pose
prevn = initFactorGraph!(fg, init=d[1][1:3])[1]
Podo=diagm([0.5;0.5;0.005])
N=100
lcmode=:mmodal
lsrNoise=diagm([0.1;1.0])

#build :# poses for the factorGraph
for idx=9:30
  prev, X, nextn = getLastPose2D(fg)
  vp, fp = addOdoFG!(fg, nextn, d[idx][1:3], Podo, N=N)
  # add landmarks
  addLandmarksFactoGraph!(fg, f, idx, prevn, nextn, lcmode=lcmode, lsrNoise=lsrNoise, N=N, MM=MM)
  prevn = nextn
  if (idx%10==0)
     # Solve
     tree = wipeBuildNewTree!(fg, drawpdf=false);
    @time inferOverTree!(fg,tree, N=100);
  end
  pl=drawPosesLandms(fg,window=(nextn,25))
  Gadfly.draw(PNG(joinpath(imgdir,"$(nextn).png"),20cm,20cm),pl)
end

0

# run(`ffmpeg -y -i x%d.png -threads 4 -vcodec libx264 -s 1920x1080 -b:v 2M -filter:v "setpts=10*PTS" /home/dehann/Videos/vicprk_test.mp4`)

# pl=drawPosesLandms(fg);
# Gadfly.draw(PDF("/tmp/before.pdf",20cm,20cm),pl)

# batchSolve(fg)
# on ssh terminal, run slamindb(iterations=1)

# fetch a local copy

fg = Caesar.initfg(sessionname=user_config["sessionId"], cloudgraph=backend_config) #, robotname=user_config["robotId"]
fullLocalGraphCopy!(fg)
# pl1=drawPosesLandms(fg)

pl1=drawPosesLandms(fg,window=(:x9,25))

Gadfly.draw(PDF("/tmp/after.pdf",20cm,20cm),pl1)

@async run(`evince /tmp/after.pdf`)

subLocalGraphCopy!(fg, string.([:x1,:x2,:x3,:x4,:x5]), neighbors=1)

X,L = ls(fg)
d=0
t=[]
z=zeros(length(L)+1,length(L)+1)

j=1;
for l1 in L
  i=1
  for l2 in L
    b=getKDEMax(getVertKDE(fg, l1))
    c=evalLikelihood(fg, l2 , b)
    if ( c>0.002 && l1!=l2 )
      d+=1
      #Replace this with A delete and replace
      #what the hell should i do now
      addSoftEqualityPoint2D(fg, l1, l2 )

      push!(t,(l1,l2))
      z[i,j]=1
    elseif (l1==l2)
      z[i,j]=2
    end
    i+=1
  end
  j+=1
end
#run Server Code
fg = Caesar.initfg(sessionname=user_config["sessionId"], cloudgraph=backend_config)
fullLocalGraphCopy!(fg)
pl2=drawPosesLandms(fg)
draw(PDF("daniel.pdf",20cm,20cm),pl2)

# Remove the new session from the server
deleteServerSession!(fg.cg, user_config["sessionId"])

#