text:
🐙 = mapreduce(line -> parse.(Int8, split(line, "")), hcat, readlines("input.txt"))
💡 = fill(false, size(🐙))
const R = CartesianIndices(🐙)
const If, Il = extrema(R)
const I1 = one(If)

function flash!(🐙, 💡, I)
    💡[I] = true
    for i in max(I-I1, If):min(I+I1, Il)
        🐙[i] += 1
        🐙[i] > 9 && !💡[i] && flash!(🐙, 💡, i)
    end
end

function step!(🐙, 💡)
    fill!(💡, false)
    for I in R; 🐙[I] += 1; end
    for I in R; 🐙[I] > 9 && !💡[I] && flash!(🐙, 💡, I); end
    for I in R; 💡[I] && (🐙[I] = 0); end
end

function solve(🐙, 💡)
    🧮 = 💥 = i = 0
    while !all(💡) || i <= 100
        i += 1
        step!(🐙, 💡)
        🧮 = i <= 100 ? 🧮 + sum(💡) : 🧮
        💥 = all(💡) ? 💥 + i : 💥
    end
    return 🧮, 💥
end

🧮, 💥 = solve(🐙, 💡)
@info "Number of flashes in first 100 steps = $🧮"
@info "First step where all flash = $💥"


text:
include("utils.jl")

@testset "Generalized Grassmann" begin
    @testset "Real" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedGrassmann(3, 2, B)
        x = [1.0 0.0; 0.0 0.5; 0.0 0.0]
        @testset "Basics" begin
            @test repr(M) ==
                  "GeneralizedGrassmann(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℝ)"
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 2
            @test base_manifold(M) === M
            @test_throws DomainError is_manifold_point(M, [1.0, 0.0, 0.0, 0.0], true)
            @test_throws DomainError is_manifold_point(
                M,
                1im * [1.0 0.0; 0.0 1.0; 0.0 0.0],
                true,
            )
            @test !is_tangent_vector(M, x, [0.0, 0.0, 1.0, 0.0])
            @test_throws DomainError is_tangent_vector(M, x, [0.0, 0.0, 1.0, 0.0], true)
            @test_throws DomainError is_tangent_vector(
                M,
                x,
                1 * im * zero_tangent_vector(M, x),
                true,
            )
            @test injectivity_radius(M) == π / 2
            @test injectivity_radius(M, ExponentialRetraction()) == π / 2
            @test injectivity_radius(M, x) == π / 2
            @test injectivity_radius(M, x, ExponentialRetraction()) == π / 2
            @test mean(M, [x, x, x]) == x
        end
        @testset "Embedding and Projection" begin
            y = similar(x)
            z = embed(M, x)
            @test z == x
            embed!(M, y, x)
            @test y == z
            a = [1.0 0.0; 0.0 2.0; 0.0 0.0]
            @test !is_manifold_point(M, a)
            b = similar(a)
            c = project(M, a)
            @test c == x
            project!(M, b, a)
            @test b == x
            X = [0.0 0.0; 0.0 0.0; -1.0 1.0]
            Y = similar(X)
            Z = embed(M, x, X)
            embed!(M, Y, x, X)
            @test Y == X
            @test Z == X
        end

        types = [Matrix{Float64}]
        TEST_STATIC_SIZED && push!(types, MMatrix{3,2,Float64,6})
        X = [0.0 0.0; 1.0 0.0; 0.0 2.0]
        Y = [0.0 0.0; -1.0 0.0; 0.0 2.0]
        @test inner(M, x, X, Y) == 0
        y = retract(M, x, X)
        z = retract(M, x, Y)
        @test is_manifold_point(M, y)
        @test is_manifold_point(M, z)
        @test retract(M, x, X) == exp(M, x, X)

        a = project(M, x + X)
        c = retract(M, x, X, ProjectionRetraction())
        d = retract(M, x, X, PolarRetraction())
        @test a == c
        @test c == d
        e = similar(a)
        retract!(M, e, x, X)
        @test e == exp(M, x, X)
        @test vector_transport_to(M, x, X, y, ProjectionTransport()) == project(M, y, X)
        @testset "Type $T" for T in types
            pts = convert.(T, [x, y, z])
            @test !is_manifold_point(M, 2 * x)
            @test_throws DomainError !is_manifold_point(M, 2 * x, true)
            @test !is_tangent_vector(M, x, y)
            @test_throws DomainError is_tangent_vector(M, x, y, true)
            test_manifold(
                M,
                pts,
                test_exp_log=false,
                default_inverse_retraction_method=LogarithmicInverseRetraction(),
                default_retraction_method=ExponentialRetraction(),
                test_injectivity_radius=false,
                test_is_tangent=true,
                test_project_tangent=true,
                test_default_vector_transport=false,
                test_forward_diff=false,
                test_reverse_diff=false,
                projection_atol_multiplier=15.0,
                retraction_atol_multiplier=10.0,
                is_tangent_atol_multiplier=4 * 10.0^2,
                retraction_methods=[PolarRetraction(), ProjectionRetraction()],
                mid_point12=nothing,
            )
        end
    end

    @testset "Complex" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedGrassmann(3, 2, B, ℂ)
        @testset "Basics" begin
            @test repr(M) ==
                  "GeneralizedGrassmann(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℂ)"
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 4
            @test !is_manifold_point(M, [1.0, 0.0, 0.0, 0.0])
            @test !is_tangent_vector(M, [1.0 0.0; 0.0 1.0; 0.0 0.0], [0.0, 0.0, 1.0, 0.0])
            x = [1.0 0.0; 0.0 0.5; 0.0 0.0]

            x = [1im 0.0; 0.0 0.5im; 0.0 0.0]
            @test is_manifold_point(M, x)
            @test !is_manifold_point(M, 2 * x)
        end
    end

    @testset "Quaternion" begin
        B = [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0]
        M = GeneralizedGrassmann(3, 2, B, ℍ)
        @test repr(M) ==
              "GeneralizedGrassmann(3, 2, [1.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 1.0], ℍ)"
        @testset "Basics" begin
            @test representation_size(M) == (3, 2)
            @test manifold_dimension(M) == 8
        end
    end
end

ext: jl


text:
using Stheno, Test, Random, BlockArrays, TimerOutputs
using BlockArrays: _BlockArray

using Stheno: ew, pw, mean_vector, cov, cov_diag
using Stheno: EQ

const to = TimerOutput()

macro timedtestset(name, code)
    return esc(:(@timeit to $name @testset $name $code))
end

include("test_util.jl")

@testset "Stheno" begin

    println("util:")
    @timedtestset "util" begin
        include(joinpath("util", "zygote_rules.jl"))
        include(joinpath("util", "covariance_matrices.jl"))
        @testset "block_arrays" begin
            include(joinpath("util", "block_arrays", "test_util.jl"))
            include(joinpath("util", "block_arrays", "dense.jl"))
            include(joinpath("util", "block_arrays", "diagonal.jl"))
        end
        include(joinpath("util", "abstract_data_set.jl"))
        include(joinpath("util", "distances.jl"))
    end

    println("gp:")
    @timedtestset "gp" begin
        include(joinpath("gp", "mean.jl"))
        include(joinpath("gp", "kernel.jl"))
        include(joinpath("gp", "gp.jl"))
    end

    println("composite:")
    @timedtestset "composite" begin
        include(joinpath("composite", "test_util.jl"))
        include(joinpath("composite", "indexing.jl"))
        include(joinpath("composite", "cross.jl"))
        include(joinpath("composite", "conditioning.jl"))
        include(joinpath("composite", "product.jl"))
        include(joinpath("composite", "addition.jl"))
        include(joinpath("composite", "compose.jl"))
        include(joinpath("composite", "approximate_conditioning.jl"))
    end

    println("abstract_gp:")
    @timedtestset "abstract_gp" begin
        include("abstract_gp.jl")
    end
end

display(to)


text:
@testset "Metabolite" begin
    m1 = Metabolite()
    m1.id = "met1"
    m1.name = "metabolite 1"
    m1.formula = "C6H12O6N"
    m1.charge = 1
    m1.compartment = "c"
    m1.notes = Dict("notes" => ["blah", "blah"])
    m1.annotations = Dict("sboterm" => ["sbo"], "kegg.compound" => ["ads", "asds"])

    @test all(
        contains.(
            sprint(show, MIME("text/plain"), m1),
            ["met1", "metabolite 1", "C6H12O6N", "blah", "asds"],
        ),
    )

    m2 = Metabolite("met2")

    m2.formula = "C6H12O6N"

    m3 = Metabolite("met3")
    m3.formula = "X"
    m3.annotations = Dict("sboterm" => ["sbo"], "kegg.compound" => ["ad2s", "asds"])

    ats = get_atoms(m1)
    @test ats["C"] == 6 && ats["N"] == 1

    m4 = Metabolite("met4")
    m4.formula = "X"
    m4.annotations = Dict("sboterm" => ["sbo"], "kegg.compound" => ["adxxx2s", "asdxxxs"])

    md = OrderedDict(m.id => m for m in [m1, m2, m3])
    @test issetequal(["met1", "met3"], ambiguously_identified_items(annotation_index(md)))
end


text:
using GaussianDistributions

"""
    φ(::Val{T}, args...)

Compute the φ function appearing in the Girsanov formula and needed for
sampling from the full conditional distribution of the parameters (whose
indices are specified by the `Val`) conditional on the path,
observations and other parameters.
"""
@generated function φ(::Val{T}, args...) where T
    z = Expr(:tuple, (:(phi(Val($i), args...)) for i in 1:length(T) if T[i])...)
    return z
end

"""
    φᶜ(::Val{T}, args...)

Compute the φᶜ function appearing in the Girsanov formula. This function
complements φ.
"""
@generated function φᶜ(::Val{T}, args...) where T
    z = Expr(:tuple, (:(phi(Val($i), args...)) for i in 0:length(T) if i==0 || !T[i])...)
    return z
end

phi(::Val{0}, t, x, P::FitzhughDiffusion) = -x[2]
phi(::Val{1}, t, x, P::FitzhughDiffusion) = x[1]-x[1]^3+(1-3*x[1]^2)*x[2]
phi(::Val{2}, t, x, P::FitzhughDiffusion) = one(x[1])
phi(::Val{3}, t, x, P::FitzhughDiffusion) = -x[1]
phi(::Val{4}, t, x, P::FitzhughDiffusion) = zero(x[1])
phi(::Val{5}, t, x, P::FitzhughDiffusion) = zero(x[1])

"""
    conjugateDraw(θ, XX, PT, prior, ::updtIdx)

Draw from the full conditional distribution of the parameters whose indices are
specified by the object `updtIdx`, conditionally on the path given in container
`XX`, and conditionally on all other parameter values given in vector `θ`.
"""
function conjugateDraw(θ, XX, PT, prior, updtIdx)
    μ = mustart(updtIdx)
    𝓦 = μ*μ'
    ϑ = SVector(thetaex(updtIdx, θ))
    μ, 𝓦 = _conjugateDraw(ϑ, μ, 𝓦, XX, PT, updtIdx)

    Σ = inv(𝓦 + inv(Matrix(prior.Σ)))
    Σ = (Σ + Σ')/2 # eliminates numerical inconsistencies
    μₚₒₛₜ = Σ * (μ + Vector(prior.Σ\prior.μ))
    rand(Gaussian(μₚₒₛₜ, Σ))
end
mustart(::Val{T}) where {T} = @SVector zeros(sum(T))
@generated function thetaex(::Val{T}, θ) where T
    z = Expr(:tuple, 1.0, (:(θ[$i]) for i in 1:length(T) if  !T[i])...)
    return z
end

function _conjugateDraw(ϑ, μ, 𝓦, XX, PT, updtIdx)
    for X in XX
        for i in 1:length(X)-1
            φₜ = SVector(φ(updtIdx, X.tt[i], X.yy[i], PT))
            φᶜₜ = SVector(φᶜ(updtIdx, X.tt[i], X.yy[i], PT))
            dt = X.tt[i+1] - X.tt[i]
            dy = X.yy[i+1][2]-X.yy[i][2]
            μ = μ + φₜ*dy - φₜ*dot(ϑ, φᶜₜ)*dt
            𝓦 = 𝓦 + φₜ*φₜ'*dt
        end
    end
    μ = μ/PT.σ^2
    𝓦 = 𝓦/PT.σ^2
    μ, 𝓦
end


text:
for f in base_functions
    @testset "Testing $f" begin
        F = (f)()
        for T in FloatingPointTypes
            s = MLK.base_initiate(F, T)
            s_tmp = convert(T, get(base_functions_initiate, f, 0))

            @test typeof(s) == T
            @test s == s_tmp

            f_tmp = get(base_functions_aggregate, f, (s,x,y) -> convert(T, NaN))
            for x in (zero(T), one(T), convert(T,2)), y in (zero(T), one(T), convert(T,2))
                s = MLK.base_aggregate(F, one(T), x, y)
                s_tmp = f_tmp(one(T), x, y)

                @test typeof(s) == T
                @test s == s_tmp
            end

            f_tmp = get(base_functions_return, f, s -> s)
            for x in (zero(T), one(T), convert(T,2))
                s = MLK.base_return(F, x)
                s_tmp = f_tmp(x)

                @test typeof(s) == T
                @test s == s_tmp
            end

            test_properties = get(base_functions_properties, f, (false,false))
            @test MLK.isstationary(F) == test_properties[1]
            @test MLK.isisotropic(F)  == test_properties[2]
        end
    end
end


text:
function insertsort!(A)
    for i in 2:length(A)
        key = A[i]
        j = i - 1
        while j >= 1 && A[j] > key
            A[j + 1] = A[j]
            j -= 1
        end
        A[j + 1] = key
    end
    return A
end

function insertsort(A)
    insertsort!(deepcopy(A))
end

function insertsortperm!(A)
    numbers = collect(1:length(A))
    for i in 2:length(A)
        key = A[i]
        j = i - 1
        while j >= 1 && A[j] > key
            A[j + 1] = A[j]
            numbers[j + 1],numbers[j] =numbers[j],numbers[j + 1]
            j -= 1
        end
        A[j + 1] = key
        numbers[j + 1],numbers[i] = numbers[i],numbers[j + 1]
    end
    return numbers
end

function insertsortperm(A)
    insertsortperm!(deepcopy(A))
end


text:
################################################################################
# 3-D model

"""
    dett(::Type{DeforModelRed3D},  C::FMat{T}) where {T}

Compute the determinant of a general square matrix.
"""
function dett(::Type{DeforModelRed3D},  C::FMat{T}) where {T}
	return (C[1,1] * C[2,2] * C[3,3] +
		C[1,2] * C[2,3] * C[3,1] +
		C[1,3] * C[2,1] * C[3,2] -
		C[1,3] * C[2,2] * C[3,1] -
		C[1,2] * C[2,1] * C[3,3] -
		C[1,1] * C[2,3] * C[3,2])
end

"""
    strainvdet(::Type{DeforModelRed3D},  Cv::FVec{T}) where {T}

Compute the determinant of a symmetric strain-like square matrix represented
as a vector. Remember that the shear strain components are twice the entries
of the matrix representation.
"""
function strainvdet(::Type{DeforModelRed3D},  Cv::FVec{T}) where {T}
	return (Cv[1] * Cv[2] * Cv[3] +
		Cv[4]/2 * Cv[6]/2 * Cv[5]/2 +
		Cv[5]/2 * Cv[4]/2 * Cv[6]/2 -
		Cv[5]/2 * Cv[2] * Cv[5]/2 -
		Cv[4]/2 * Cv[4]/2 * Cv[3] -
		Cv[1] * Cv[6]/2 * Cv[6]/2)
end

"""
    strainvtr(::Type{DeforModelRed3D},  Cv::FVec{T}) where {T}

Compute the trace of a symmetric strain-like square matrix represented as a
vector.
"""
function strainvtr(::Type{DeforModelRed3D},  Cv::FVec{T}) where {T}
	return (Cv[1] + Cv[2] + Cv[3])
end

"""
    strainttov!(::Type{DeforModelRed3D}, v::FVec{T}, t::FMat{T}) where {T}

Convert a symmetric matrix of 3x3 strain components  into a 6-component vector.
"""
function strainttov!(::Type{DeforModelRed3D}, v::FVec{T}, t::FMat{T}) where {T}
	v[1] = t[1,1];
	v[2] = t[2,2];
	v[3] = t[3,3];
	v[4] = t[1,2] + t[2,1];
	v[5] = t[1,3] + t[3,1];
	v[6] = t[3,2] + t[2,3];
	return v
end

"""
    strainvtot!(::Type{DeforModelRed3D}, t::FMat{T}, v::FVec{T}) where {T}

Convert a strain 3-vector to a  matrix of 2x2 strain components (symmetric tensor).
"""
function strainvtot!(::Type{DeforModelRed3D}, t::FMat{T}, v::FVec{T}) where {T}
	t[1,1] = v[1];
	t[2,2] = v[2];
	t[3,3] = v[3];
	t[1,2] = v[4]/2.;
	t[2,1] = v[4]/2.;
	t[1,3] = v[5]/2.;
	t[3,1] = v[5]/2.;
	t[3,2] = v[6]/2.;
	t[2,3] = v[6]/2.;
    return t
end

"""
    stressttov!(::Type{DeforModelRed3D}, v::FVec{T}, t::FMat{T}) where {T}

Convert a symmetric matrix of 3x3 stress components to a 6-component vector.
"""
function stressttov!(::Type{DeforModelRed3D}, v::FVec{T}, t::FMat{T}) where {T}
	v[1] = t[1,1];
	v[2] = t[2,2];
	v[3] = t[3,3];
	v[4] = 1/2.0*(t[1,2] + t[2,1]);
	v[5] = 1/2.0*(t[1,3] + t[3,1]);
	v[6] = 1/2.0*(t[3,2] + t[2,3]);
	return v
end

"""
    stressvtot!(::Type{DeforModelRed3D}, t::FMat{T}, v::FVec{T}) where {T}

Convert a 6-vector to a  matrix of 3x3 stress components (symmetric tensor).
"""
function stressvtot!(::Type{DeforModelRed3D}, t::FMat{T}, v::FVec{T}) where {T}
	t[1,1] = v[1];
	t[2,2] = v[2];
	t[3,3] = v[3];
	t[1,2] = v[4];
	t[2,1] = v[4];
	t[1,3] = v[5];
	t[3,1] = v[5];
	t[3,2] = v[6];
	t[2,3] = v[6];
	return t
end

################################################################################
# 2-D plane strain model

"""
    dett(::Type{DeforModelRed2DStrain},  C::FMat{T}) where {T}

Compute the determinant of a general square matrix.
"""
function dett(::Type{DeforModelRed2DStrain},  C::FMat{T}) where {T}
	return (C[1,1] * C[2,2] - C[1,2] * C[2,1])
end

"""
    strainvdet(::Type{DeforModelRed2DStrain},  Cv::FVec{T}) where {T}

Compute the determinant of a symmetric strain-like square matrix represented
as a vector. Remember that the shear strain components are twice the entries
of the matrix representation.
"""
function strainvdet(::Type{DeforModelRed2DStrain},  Cv::FVec{T}) where {T}
	return (Cv[1] * Cv[2] - Cv[3]/2 * Cv[3]/2)
end

"""
    strainvtr(::Type{DeforModelRed2DStrain},  Cv::FVec{T}) where {T}

Compute the trace of a symmetric strain-like square matrix represented as a
vector.
"""
function strainvtr(::Type{DeforModelRed2DStrain},  Cv::FVec{T}) where {T}
	return (Cv[1] + Cv[2])
end

"""
    strainttov!(::Type{DeforModelRed2DStrain}, v::FVec{T}, t::FMat{T}) where {T}

Convert a symmetric matrix of 2x2 strain components  into a 3-component vector.
"""
function strainttov!(::Type{DeforModelRed2DStrain}, v::FVec{T}, t::FMat{T}) where {T}
    v[1] = t[1,1];
    v[2] = t[2,2];
    v[3] = t[1,2] + t[2,1];
    return v
end

"""
    strainvtot!(::Type{DeforModelRed2DStrain}, t::FMat{T}, v::FVec{T}) where {T}

Convert a strain 3-vector to a  matrix of 2x2 strain components (symmetric tensor).
"""
function strainvtot!(::Type{DeforModelRed2DStrain}, t::FMat{T}, v::FVec{T}) where {T}
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3]/2.;
    t[2,1] = v[3]/2.;
    return t
end

"""
    stressttov!(::Type{DeforModelRed2DStrain}, v::FVec{T}, t::FMat{T}) where {T}

Convert a symmetric matrix of 2x2 stress components to a 3-component vector.
"""
function stressttov!(::Type{DeforModelRed2DStrain}, v::FVec{T}, t::FMat{T}) where {T}
	v[1] = t[1,1];
	v[2] = t[2,2];
	v[3] = 0.5*(t[1,2] + t[2,1]);
	return v
end

"""
    stressvtot!(::Type{DeforModelRed2DStrain}, t::FMat{T}, v::FVec{T}) where {T}

Convert a vector to a  matrix of 2x2 stress components (symmetric tensor).

If `v` has 4 entries, also the `t[3,3]` matrix entry is set.

The stress vector components need to be ordered as:
    sigmax, sigmay, tauxy, sigmaz,
which is the ordering used for the plane-strain model reduction.
"""
function stressvtot!(::Type{DeforModelRed2DStrain}, t::FMat{T}, v::FVec{T}) where {T}
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3];
    t[2,1] = v[3];
    (length(v) == 4)  && (t[3,3] = v[4]);
    return t
end

################################################################################
# 2-D plane stress model

"""
    stressttov!(::Type{DeforModelRed2DStress}, v::FVec{T}, t::FMat{T}) where {T}

Convert a symmetric matrix of 2x2 stress components to a 3-component vector.
"""
function stressttov!(::Type{DeforModelRed2DStress}, v::FVec{T}, t::FMat{T}) where {T}
	v[1] = t[1,1];
	v[2] = t[2,2];
	v[3] = 0.5*(t[1,2] + t[2,1]);
	return v
end

"""
    stressvtot!(::Type{DeforModelRed2DStress}, t::FMat{T}, v::FVec{T}) where {T}

Convert a 3-vector to a  matrix of 2x2 stress components (symmetric tensor).
"""
function stressvtot!(::Type{DeforModelRed2DStress}, t::FMat{T}, v::FVec{T}) where {T}
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[3];
    t[2,1] = v[3];
    return t
end

################################################################################
# 2-D axially symmetric stress model

"""
    stressvtot!(::Type{DeforModelRed2DAxisymm}, t::FMat{T}, v::FVec{T}) where {T}

Convert a 4-vector to a  matrix of 3x3 stress components (tensor).

Convert a 4-vector to a *symmetric* matrix of 3x3 stress components (tensor).

The stress vector components need to be ordered as:
    sigmax, sigmay, sigmaz, tauxy.
"""
function stressvtot!(::Type{DeforModelRed2DAxisymm}, t::FMat{T}, v::FVec{T}) where {T}
    t[1,1] = v[1];
    t[2,2] = v[2];
    t[1,2] = v[4];
    t[2,1] = v[4];
    t[3,3] = v[3];
    return t
end


text:
import Test: @test, @testset, @test_throws
import GSL: sf_bessel_zero_Jnu
import SpecialFunctions: besselj
import LinearAlgebra: norm
import FFTW
import Luna: Modes, Maths, Capillary, Grid, PhysData, Hankel, NonlinearRHS, Fields
import Luna.PhysData: wlfreq

@testset "delegation" begin
    m = Modes.delegated(Capillary.MarcatiliMode(75e-6, :He, 5.9))
    @test Modes.Aeff(m) ≈ 8.42157534886545e-09
    @test abs(1e9*Modes.zdw(m) - 562) < 1
    n(m, ω; z=0) = real(Modes.neff(m, ω; z=z))
    m2 = Modes.delegated(Capillary.MarcatiliMode(75e-6, :He, 1.0), neff=n)
    @test Modes.α(m2, 2e15) == 0
    @test Modes.α(m2, wlfreq(800e-9)) == 0

    cm = Capillary.MarcatiliMode(75e-6, :He, 5.9)
    dm = Modes.delegated(cm)
    @test Modes.Aeff(dm) ≈ 8.42157534886545e-09

    dm = Modes.delegated(cm, Aeff=(m; z=0) -> 0.5)
    @test Modes.Aeff(dm) == 0.5
    @test Modes.N(dm) == Modes.N(m)
    @test Modes.β(dm, 2e15) == Modes.β(m, 2e15)

    # fully delegated test
    n(ω; z=0) = 1.0 + 3.0im
    m3 = Modes.arbitrary(neff=n,
                         dimlimits=(;z=0)->(:polar, (0.0, 0.0), (100e-6, 2π)),
                         Aeff=(;z=0) -> 1e-5,
                         N=(;z=0) -> 1)
    @test Modes.α(m3, 1.5e15) == 2*3*1.5e15/PhysData.c
    @test Modes.β(m3, 1.5e15) == 1.5e15/PhysData.c
    @test Modes.Aeff(m3) == 1e-5
    @test Modes.N(m3) == 1
    @test_throws ErrorException Modes.field(m3, (0.0, 0.0))

    m4 = Modes.arbitrary(neff=Modes.neff_from_αβ((ω; z=0) -> 6*ω/PhysData.c,
                                                 (ω; z=0) -> ω/PhysData.c))
    @test Modes.α(m4, 1.5e15) ≈ 2*3*1.5e15/PhysData.c
    @test Modes.β(m4, 1.5e15) ≈ 1.5e15/PhysData.c
    @test Modes.neff(m4, 1.5e15) ≈ 1.0 + 3.0im
    @test Modes.neff(m4, 1.5e15) ≈ 1.0 + 3.0im
end

@testset "overlap" begin
a = 100e-6
m = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced)
r = collect(range(0, a, length=2^16))
unm = sf_bessel_zero_Jnu(0, 1)
Er = besselj.(0, unm*r/a) # spatial profile of the HE11 mode - overlap should be perfect
η = Modes.overlap(m, r, Er; dim=1)
@test abs2(η[1]) ≈ 1
unm = sf_bessel_zero_Jnu(0, 2)
Er = besselj.(0, unm*r/a) # spatial profile of HE12 - overlap should be 0
η = Modes.overlap(m, r, Er; dim=1)
@test isapprox(abs2(η[1]), 0, atol=1e-18)
# Check that we reproduce the optimal coupling for Gaussian beams at w₀ = 0.64a
fac = collect(range(0.3, stop=0.9, length=128))
ηn = zero(fac)
r = collect(range(0, 4a, length=2^16))
m = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=1)
for i in eachindex(fac)
    w0 = fac[i]*a
    Er = Maths.gauss.(r, w0/sqrt(2))
    ηn[i] = abs2.(Modes.overlap(m, r, Er, dim=1)[1])
end
@test 0.63 < fac[argmax(ηn)] < 0.65

#= Testing that the sum of all modes is 1, i.e. the normalised overlap preserves the
    total energy =#
fac = 0.45
w0 = fac*a
Er = Maths.gauss.(r, w0/sqrt(2))
s = 0
for mi = 1:10
    m = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=mi)
    η = Modes.overlap(m, r, Er, dim=1)[1]
    s += abs2(η[1])
end
@test isapprox(s, 1, rtol=1e-4)

#= Testing non-normalised overlap integrals. First we create a spatial field that consists
    of two pulses: one 30 fs pulse at 800 nm with the spatial shape of the HE11 mode,
    and one 15 fs pulse at 400 nm with the spatial shape of the HE12 mode.
    We then calculate the non-normalised overlap integral of this spatial field with
    the HE11 and HE12 modes of a capillary and check that we reproduce the pulses that we
    put in at the beginning, as well as their energy. =#
a = 100e-6 # capillary radius
# spatial grid, with a bigger aperture than the capillary - as we would have in a simulation
q = Hankel.QDHT(2a, 512)
grid = Grid.RealGrid(1, 800e-9, (200e-9, 2000e-9), 0.5e-12)
# First pulse
It1 = Maths.gauss.(grid.t, fwhm=30e-15)
Et1 = @. sqrt(It1)*cos(2π*PhysData.c/800e-9*grid.t)
Eω1 = FFTW.rfft(Et1)
# Spatial profile of the first pulse
unm = sf_bessel_zero_Jnu(0, 1)
Er1 = besselj.(0, unm*q.r/a)'
Er1[q.r .> a] .= 0
Etr1 = Et1 .* Er1 # create spatio-temporal pulse profile
# Second pulse
It2 = 4*Maths.gauss.(grid.t, fwhm=15e-15)
Et2 = @. sqrt(It2)*cos(2π*PhysData.c/400e-9*grid.t)
Eω2 = FFTW.rfft(Et2)
# Spatial profile of the second pulse
unm = sf_bessel_zero_Jnu(0, 2)
Er2 = besselj.(0, unm*q.r/a)'
Er2[q.r .> a] .= 0
Etr2 = Et2 .* Er2 # create spatio-temporal pulse profile

# The total spatio-temporal field is the sum of the two pulses
Etr = Etr1 .+ Etr2
Eωr = FFTW.rfft(Etr, 1)

ert, ekω = Fields.energyfuncs(grid, q)
energy1 = ert(Etr1)
energy2 = ert(Etr2)
@test ert(Etr) ≈ energy1 + energy2

m1 = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=1)
m2 = Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=2)
Eωm1 = Modes.overlap(m1, q.r, Eωr; dim=2, norm=false)
Eωm2 = Modes.overlap(m2, q.r, Eωr; dim=2, norm=false)

Etm1 = FFTW.irfft(Eωm1[:, 1], length(grid.t))
Etm2 = FFTW.irfft(Eωm2[:, 1], length(grid.t))

et, eω = Fields.energyfuncs(grid)
# check that the non-normalised overlap integral preserves the total energy
@test isapprox(et(Etm1), energy1, rtol=1e-3)
@test isapprox(et(Etm2), energy2, rtol=1e-3)

# check that the non-normalised overlap integral preserves the spectral shape
En1 = Eω1/norm(Eω1)
Emn1 = Eωm1/norm(Eωm1)
@test all(isapprox.(En1, Emn1, atol=1e-3*maximum(abs.(Emn1))))

En2 = Eω2/norm(Eω2)
Emn2 = Eωm2/norm(Eωm2)
@test all(isapprox.(En2, Emn2, atol=1e-3*maximum(abs.(Emn2))))
end

@testset "overlap with interpolation" begin
#= Here we test the complete decomposition of a field into two modes while also re-gridding
    onto a new time grid. We create two pulses in different spatial modes with different
    arrival times and durations on one grid and check that we reproduce them on a new grid.
    =#
a = 100e-6 # capillary radius
# spatial grid, with a bigger aperture than the capillary - as we would have in a simulation
q = Hankel.QDHT(2a, 512)
grid = Grid.RealGrid(1, 800e-9, (400e-9, 1000e-9), 0.5e-12)

fwhm1 = 30e-15
τ1 = -5e-15
fwhm2 = 15e-15
τ2 = 10e-15

# First pulse
It1 = Maths.gauss.(grid.t, x0=τ1, fwhm=fwhm1)
Et1 = @. sqrt(It1)*cos(2π*PhysData.c/800e-9*grid.t)
Eω1 = FFTW.rfft(Et1)
# Spatial profile of the first pulse
unm = sf_bessel_zero_Jnu(0, 1)
Er1 = besselj.(0, unm*q.r/a)'
Er1[q.r .> a] .= 0
Etr1 = Et1 .* Er1 # create spatio-temporal pulse profile
# Second pulse
It2 = 4*Maths.gauss.(grid.t, x0=τ2, fwhm=fwhm2)
Et2 = @. sqrt(It2)*cos(2π*PhysData.c/800e-9*grid.t)
Eω2 = FFTW.rfft(Et2)
# Spatial profile of the second pulse
unm = sf_bessel_zero_Jnu(0, 2)
Er2 = besselj.(0, unm*q.r/a)'
Er2[q.r .> a] .= 0
Etr2 = Et2 .* Er2 # create spatio-temporal pulse profile

# The total spatio-temporal field is the sum of the two pulses
Etr = Etr1 .+ Etr2
Eωr = FFTW.rfft(Etr, 1)

ert, ekω = Fields.energyfuncs(grid, q)
energy1 = ert(Etr1)
energy2 = ert(Etr2)
@test ert(Etr) ≈ energy1 + energy2

modes = (Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=1),
         Capillary.MarcatiliMode(a, :He, 1.0, model=:reduced, m=2))
newgrid = Grid.RealGrid(1, 800e-9, (160e-9, 3000e-9), 1e-12)

Eωm = Modes.overlap(modes, newgrid, grid, q.r, Eωr)

et, eω = Fields.energyfuncs(newgrid)
# do we preserve the energy?
@test isapprox(eω(Eωm[:, 1]), energy1, rtol=1e-3)
@test isapprox(eω(Eωm[:, 2]), energy2, rtol=1e-3)

Etm = FFTW.irfft(Eωm, length(newgrid.t), 1)
Itm = abs2.(Maths.hilbert(Etm))
# is fwhm the same?
@test isapprox(Maths.fwhm(newgrid.t, Itm[:, 1]), fwhm1, rtol=1e-4)
@test isapprox(Maths.fwhm(newgrid.t, Itm[:, 2]), fwhm2, rtol=1e-4)
# do we preserve the temporal location?
@test abs(newgrid.t[argmax(Itm[:, 1])] - τ1) < (newgrid.t[2]-newgrid.t[1])
@test isapprox(Maths.moment(newgrid.t, Itm[:, 1]), τ1, rtol=1e-3)
@test abs(newgrid.t[argmax(Itm[:, 2])] - τ2) < (newgrid.t[2]-newgrid.t[1])
@test isapprox(Maths.moment(newgrid.t, Itm[:, 2]), τ2, rtol=1e-3)

# are the temporal shapes reproduced?
It1new = Maths.gauss.(newgrid.t, x0=τ1, fwhm=fwhm1)
It1new /= norm(It1new)
Itm1 = Itm[:, 1]/norm(Itm[:, 1])
@test all(isapprox.(Itm1, It1new, atol=1e-3*maximum(abs.(It1new))))

It2new = Maths.gauss.(newgrid.t, x0=τ2, fwhm=fwhm2)
It2new /= norm(It2new)
Itm2 = Itm[:, 2]/norm(Itm[:, 2])
@test all(isapprox.(Itm2, It2new, atol=1e-3*maximum(abs.(It2new))))
end

##
@testset "Orthogonality" begin
a = 100e-6
m1 = Capillary.MarcatiliMode(a; n=1, m=1)
m2 = Capillary.MarcatiliMode(a; n=1, m=2)
@test Modes.overlap(m1, m1) ≈ 1
@test isapprox(Modes.overlap(m1, m2), 0; atol=1e-10)
@test Modes.orthonormal((m1, m2))
# Capillary EH modes are all orthonormal
mm = collect([Capillary.MarcatiliMode(a; n=1, m=mi) for mi=1:5])
@test Modes.orthonormal(mm)
# Two identical modes in the list - not a valid basis set
@test !Modes.orthonormal([mm..., m1])

m1 = Capillary.MarcatiliMode(a)
m2 = Capillary.MarcatiliMode(a; ϕ=π/2)
@test isapprox(Modes.overlap(m1, m2), 0; atol=1e-10)

m2 = Capillary.MarcatiliMode(2a; ϕ=π/2)
@test_throws ErrorException Modes.overlap(m1, m2)
end


text:
function ex13(r::Robot)
  x = 0
  y = 0
  while !isborder(r, Sud)
    move!(r, Sud)
    y += 1
  end
  while !isborder(r, West)
    move!(r, West)
    x += 1
  end

  putmarkerss(r)

  for i in 1:x
    move!(r, Ost)
  end

  for i in 1:y
    move!(r, Nord)
  end

end

function putmarkerss(r::Robot)
  putmarker!(r)
  funcname(r, Nord)
  while !isborder(r, West)
    move!(r, West)
  end
  putmarker!(r)
  funcname(r, Sud)
  while !isborder(r, West)
    move!(r, West)
  end
end

function funcname(r::Robot, side::HorizonSide)
  while !isborder(r, side) && !isborder(r, Ost)
    move!(r, side)
    move!(r, Ost)
    putmarker!(r)
  end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
import Pkg
println(pwd())
if !isfile("Project.toml")
    println("Project.toml が見つからない。")
    exit(1)
end
Pkg.activate(".")
using HtmlTable, DataFrames
A = [1000 200 430; 2050 560 870]
html = fromMatrix(A, ["A", "B", "C"], rownum=true)
write("./test/result1.html", html)
df = DataFrame(A, [:A, :B, :C])
html = fromDataFrame(df, rownum=true)
write("./test/result2.html", html)
html = fromCSV("./test/data.csv", rownum=true)
write("./test/result3.html", html)
println("test/result1.html, test/result2.html, test/result3.htmlなどが更新されました。")

max_stars_repo_name: makandat/julia

max_issues_repo_name: makandat/julia

max_forks_repo_name: makandat/julia


text:
# ---
# title: 1000. Minimum Cost to Merge Stones
# id: problem1000
# author: Indigo
# date: 2021-02-18
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-cost-to-merge-stones/description/>
# hidden: true
# ---
#
# There are `N` piles of stones arranged in a row.  The `i`-th pile has
# `stones[i]` stones.
#
# A _move_ consists of merging **exactly  `K` consecutive** piles into one pile,
# and the cost of this move is equal to the total number of stones in these `K`
# piles.
#
# Find the minimum cost to merge all piles of stones into one pile.  If it is
# impossible, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: stones = [3,2,4,1], K = 2
#     Output: 20
#     Explanation:
#     We start with [3, 2, 4, 1].
#     We merge [3, 2] for a cost of 5, and we are left with [5, 4, 1].
#     We merge [4, 1] for a cost of 5, and we are left with [5, 5].
#     We merge [5, 5] for a cost of 10, and we are left with [10].
#     The total cost was 20, and this is the minimum possible.
#
#
# **Example 2:**
#
#
#
#     Input: stones = [3,2,4,1], K = 3
#     Output: -1
#     Explanation: After any merge operation, there are 2 piles left, and we can't merge anymore.  So the task is impossible.
#
#
# **Example 3:**
#
#
#
#     Input: stones = [3,5,1,2,6], K = 3
#     Output: 25
#     Explanation:
#     We start with [3, 5, 1, 2, 6].
#     We merge [5, 1, 2] for a cost of 8, and we are left with [3, 8, 6].
#     We merge [3, 8, 6] for a cost of 17, and we are left with [17].
#     The total cost was 25, and this is the minimum possible.
#
#
#
#
# **Note:**
#
#   * `1 <= stones.length <= 30`
#   * `2 <= K <= 30`
#   * `1 <= stones[i] <= 100`
#
#
## @lc code=start
using LeetCode

function merge_stones(stones::Vector{Int}, K::Int)
    len = length(stones)
    (len - 1) % (K - 1) != 0 && return -1
    pre = fill(0, len + 1)
    cumsum!(@view(pre[2:len + 1]), stones)
    dp = fill(0, len, len)
    for m in K:len, i in 1:(len - m + 1)
        dp[i, i + m - 1] = minimum(dp[i, k] + dp[k + 1, i + m - 1]
                               for k in i:(K - 1):(i + m - 2)) +
                           ((m - 1) % (K - 1) == 0 ? pre[i + m] - pre[i] : 0)
    end
    return dp[1, len]
end
## @lc code=end


text:
include("../src/utils/pasquill_gifford.jl")
include("../src/models/plume_rise.jl")

@testset "Gaussian plume tests" begin
    # Gaussian plume example, *Guidelines for Consequence Analysis of Chemical
    # Releases* CCPS, 1999, pg 97
    ex = Scenario( Dict([
        :mass_emission_rate => 0.1,
        :jet_diameter => 10.0,
        :jet_velocity => 1.0,
        :release_height => 0.0,
        :windspeed => 2.0,
        :pasquill_gifford => "F"
    ]))
    # known answers
    x₁ = 500.0

    # missing model params
    @test_throws MissingException plume(ambient, model=:gaussian)

    # missing params exclusively for plume rise
    s = Scenario(ex; release_temperature=missing)
    @test_throws MissingException plume(s, model=:gaussian, plumerise=true)

    # invalid stability class for plume rise
    s = Scenario(ex; release_temperature=325, ambient_temperature=298,
                     pasquill_gifford="not a real one")
    @test_throws ErrorException plume_rise(s, true)

    @testset "Gaussian plume tests for class $class" for class in ["A", "B", "C", "D", "E", "F"]
        s = Scenario(ex; pasquill_gifford=class)

        m, u = s.mass_emission_rate, s.windspeed
        σy, σz = crosswind_dispersion(class)(x₁), vertical_dispersion(class)(x₁)
        c₁ = m / (π*σy*σz*u)

        # basic model, no plume rise no stack downwash
        no_plume_rise = plume(s, model=:gaussian)
        @test no_plume_rise(x₁,0,0) ≈ c₁

        # basic model, with stack downwash included
        Dⱼ, uⱼ = s.jet_diameter, s.jet_velocity
        u, h = 2*uⱼ, 1.0
        hₑ = h + 2*Dⱼ*( (uⱼ/u) - 1.5 )
        c₂ = m / (2π*σy*σz*u)*( exp(-0.5*((h-hₑ)/σz)^2) + exp(-0.5*((h+hₑ)/σz)^2) )
        downwash = plume(Scenario(s; release_height=h, windspeed=u), model=:gaussian, downwash=true)
        @test downwash(x₁, 0, h) ≈ c₂

        # models with plume rise
        # buoyancy driven, Fb<55
        u = s.windspeed
        Tᵣ, Tₐ = 315, 298
        bs1 = Scenario(s; release_height=h, release_temperature=Tᵣ,
                                            ambient_temperature=Tₐ)
        Δh = plume_rise(bs1, true)(x₁)
        hₑ  = h + Δh
        σyₑ = √( (Δh/3.5)^2 + σy^2 )
        σzₑ = √( (Δh/3.5)^2 + σz^2 )
        c₃ = m / (2π*σyₑ*σzₑ*u)*( exp(-0.5*((h-hₑ)/σzₑ)^2) + exp(-0.5*((h+hₑ)/σzₑ)^2) )
        b1 = plume(Scenario(bs1), model=:gaussian, plumerise=true)
        @test b1(x₁, 0, h) ≈ c₃

        # buoyancy driven, Fb>55
        Tᵣ = 400
        bs2 = Scenario(bs1; release_temperature=Tᵣ)
        Δh = plume_rise(bs2, true)(x₁)
        hₑ  = h + Δh
        σyₑ = √( (Δh/3.5)^2 + σy^2 )
        σzₑ = √( (Δh/3.5)^2 + σz^2 )
        c₄ = m / (2π*σyₑ*σzₑ*u)*( exp(-0.5*((h-hₑ)/σzₑ)^2) + exp(-0.5*((h+hₑ)/σzₑ)^2) )
        b2 = plume(Scenario(bs2), model=:gaussian, plumerise=true)
        @test b2(x₁, 0, h) ≈ c₄

        # momentum driven, Fb<55
        Tᵣ = Tₐ
        ms1 = Scenario(bs1; release_temperature=Tᵣ)
        Δh = plume_rise(ms1, true)(x₁)
        hₑ  = h + Δh
        σyₑ = √( (Δh/3.5)^2 + σy^2 )
        σzₑ = √( (Δh/3.5)^2 + σz^2 )
        c₅ = m / (2π*σyₑ*σzₑ*u)*( exp(-0.5*((h-hₑ)/σzₑ)^2) + exp(-0.5*((h+hₑ)/σzₑ)^2) )
        m1 = plume(Scenario(ms1), model=:gaussian, plumerise=true)
        @test m1(x₁, 0, h) ≈ c₅

        # momentum driven, Fb>55
        g = 9.80616
        uⱼ = ( (4/0.00575)*55*g )^(3/5) * (1/Dⱼ)
        ms2 = Scenario(bs1; jet_velocity=uⱼ)
        Δh = plume_rise(ms2, true)(x₁)
        hₑ  = h + Δh
        σyₑ = √( (Δh/3.5)^2 + σy^2 )
        σzₑ = √( (Δh/3.5)^2 + σz^2 )
        c₆ = m / (2π*σyₑ*σzₑ*u)*( exp(-0.5*((h-hₑ)/σzₑ)^2) + exp(-0.5*((h+hₑ)/σzₑ)^2) )
        m2 = plume(Scenario(ms2), model=:gaussian, plumerise=true)
        @test m2(x₁, 0, h) ≈ c₆

    end

end

@testset "Gaussian puff tests" begin
    # Gaussian plume example, *Guidelines for Consequence Analysis of Chemical
    # Releases* CCPS, 1999, pg 97
    ex = Scenario( Dict([
        :mass_emission_rate => 1.0,
        :release_duration => 1.0,
        :jet_diameter => 10.0,
        :jet_velocity => 1.0,
        :release_height => 0.0,
        :windspeed => 2.0,
        :pasquill_gifford => "F"
    ]))
    # knowns
    m = ex.mass_emission_rate*ex.release_duration
    h = ex.release_height
    u = ex.windspeed
    x₁ = 500.0
    t₁ = x₁/u

    # missing model params
    @test_throws MissingException puff(ambient, model=:gaussian)

    @testset "Gaussian puff tests for class $class" for class in ["A", "B", "C", "D", "E", "F"]
        s = Scenario( ex; pasquill_gifford=class )

        σx = crosswind_dispersion(class, plume=false)(x₁)
        σy = σx
        σz = vertical_dispersion(class, plume=false)(x₁)

        c₁ = m / (√(2π)*π*σx*σy*σz)

        test_puff = puff(s, model=:gaussian)
        @test test_puff(x₁,0,h,t₁) ≈ c₁

    end
end


text:
using RDatasets, DataArrays, DataFrames, Gadfly

set_default_plot_size(10inch, 10inch)

barley = dataset("lattice", "barley")
levels!(barley[:Year], ["1931", "1932"])

idx = [startswith(x,"No.") for x in barley[:Variety]]
plot(barley[idx,:],
     xgroup="Variety", ygroup="Site", x="Year", y="Yield",
     Geom.subplot_grid(Geom.line, Geom.point,
                       free_x_axis=true, free_y_axis=true))


text:
function draw_line(vent)
    x1, y1, x2, y2 = vent
    xs = x1 == x2 ? fill(x1, abs(y2 - y1) + 1) : x1:sign(x2 - x1):x2
    ys = y1 == y2 ? fill(y1, abs(x2 - x1) + 1) : y1:sign(y2 - y1):y2
    [zip(xs, ys)...]
end

function overlap_lines(lines)
    overlaps = Dict{Tuple{Int, Int}, Int}()
    points = reduce(vcat, lines)
    for point in points
        overlaps[point] = get(overlaps, point, 0) + 1
    end
    overlaps
end

function is_hv(vent)
    x1, y1, x2, y2 = vent
    x1 == x2 || y1 == y2
end

vents = [
    parse.(Int, split(replace(vent, " -> " => ','), ','))
    for vent in readlines("day5/vents.txt")
]

hv_lines = draw_line.(filter(is_hv, vents))
hv_overlaps = overlap_lines(hv_lines)

p1 = sum(values(hv_overlaps) .>= 2)

@show p1

lines = draw_line.(vents)
overlaps = overlap_lines(lines)

p2 = sum(values(overlaps) .>= 2)

@show p2


text:
# the need to define python functions is a bit annoying, and against
# the spirit of the package, however, the amount of work necessary to implement
# a pure-julia interface is quite high.
#
# for now I am satisfied that likat 95% of kwant constructions so far
# have a pure-julia interface (i.e. no need to use PyCall directly)

function make_system(a=1, t=1.0, W=10, r1=10, r2=20)
    # Start with an empty tight-binding system and a single square lattice.
    # `a` is the lattice constant (by default set to 1 for simplicity).

    lat = kwant.lattice.square(a)

    syst = kwant.Builder()

    #### Define the scattering region. ####
    # Now, we aim for a more complex shape, namely a ring (or annulus)
    function ring(pos)
        (x, y) = pos
        rsq = x^2 + y^2
        return (r1^2 < rsq < r2^2)
    end

    # and add the corresponding lattice points using the `shape`-function
    syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
    syst[lat.neighbors()] = -t

    # In order to introduce a flux through the ring, we introduce a phase on
    # the hoppings on the line cut through one of the arms.  Since we want to
    # change the flux without modifying the Builder instance repeatedly, we
    # define the modified hoppings as a function that takes the flux as its
    # parameter phi.
    py"""
    def hopping_phase_py(site1, site2, phi):
        return -$t * $exp(1j * phi)
    """

    py"""
    def crosses_branchcut_py(hop):
        ix0, iy0 = hop[0].tag
        # builder.HoppingKind with the argument (1, 0) below
        # returns hoppings ordered as ((i+1, j), (i, j))
        return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
    """
    # Modify only those hopings in x-direction that cross the branch cut
    py"""
    import kwant
    def hops_across_cut_py(syst):
        for hop in kwant.builder.HoppingKind((1, 0), $(lat.o), $(lat.o))(syst):
            if crosses_branchcut_py(hop):
                yield hop
    """
    syst[py"hops_across_cut_py"] = py"hopping_phase_py"

    #### Define the leads. ####
    # left lead
    sym_lead = kwant.TranslationalSymmetry((-a, 0))
    lead = kwant.Builder(sym_lead)

    function lead_shape(pos)
        (x,y) = pos
        return (-W / 2 < y < W / 2)
    end

    lead[lat.shape(lead_shape, (0, 0))] = 4 * t
    lead[lat.neighbors()] = -t

    #### Attach the leads and return the system. ####
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    return syst
end

function plot_conductance(syst; energy, fluxes)
    # compute conductance

    normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
    data = []
    for flux in fluxes
        smatrix = kwant.smatrix(syst, energy, params=Dict(:phi=>flux))
        append!(data,smatrix.transmission(1, 0))
    end

    pyplot.figure()
    pyplot.plot(normalized_fluxes, data)
    pyplot.xlabel("flux [flux quantum]")
    pyplot.ylabel("conductance [e^2/h]")
    pyplot.gcf()
end


text:
using CSV, HTTP, JSON

mutable struct Preference
    name                            ::String
    preferences                     ::Dict{String, Int}
    currentLT                       ::String
    alreadyAssigned                 ::String
    enableAssignmentToCurrentLT     ::Bool
end

Preference(name) = Preference(
    name,
    Dict{String, Int}(),
    "",
    "",
    0
)

function load_data(args)
    env = Dict{String, Any}()

    data = JSON.parse(open("..\\data\\preferences.json", "r"))
    # data = HTTP.get("https://firestore.googleapis.com/v1/projects/shadow-5f1c5123/databases/(default)/documents/preferences")
    # println(data)

    MC = Dict{String, Any}()
    for fields in data["documents"]
        name = fields["fields"]["name"]["stringValue"]
        MC[name] = Preference(name)

        MC[name].preferences[fields["fields"]["preference1"]["stringValue"]] = 1
        MC[name].preferences[fields["fields"]["preference2"]["stringValue"]] = 2
        MC[name].preferences[fields["fields"]["preference3"]["stringValue"]] = 3
        MC[name].preferences[fields["fields"]["preference4"]["stringValue"]] = 4
        MC[name].preferences[fields["fields"]["preference5"]["stringValue"]] = 5
        MC[name].preferences[fields["fields"]["preference6"]["stringValue"]] = 6
        MC[name].preferences[fields["fields"]["preference7"]["stringValue"]] = 7
        MC[name].preferences[fields["fields"]["preference8"]["stringValue"]] = 8
        MC[name].preferences[fields["fields"]["preference9"]["stringValue"]] = 9
        MC[name].preferences[fields["fields"]["preference10"]["stringValue"]] = 10
        MC[name].preferences[fields["fields"]["preference11"]["stringValue"]] = 11
        MC[name].preferences[fields["fields"]["preference12"]["stringValue"]] = 12

        MC[name].currentLT = fields["fields"]["currentLT"]["stringValue"]

        MC[name].alreadyAssigned = fields["fields"]["alreadyAssigned"]["stringValue"]
        MC[name].preferences[MC[name].alreadyAssigned] = 1000

        MC[name].enableAssignmentToCurrentLT = fields["fields"]["enableAssignmentToCurrentLT"]["booleanValue"]
    end

    println(MC)

    #data = HTTP.get("https://firestore.googleapis.com/v1/projects/shadow-5f1c5123/databases/(default)/documents/preferences")

    env["data"] = MC

    env["LT"] = [
        "Ben Schurr",
        "Brendan Ryan",
        "Brent Mealings",
        "Chris Greenough",
        "Chris Rowe",
        "Eric Hansen",
        "Gordon Rettkowicz",
        "Jo Martell",
        "Mark Kennerley",
        "Paul Campbell",
        "Piers Shore",
        "Priya Singh",
        "Simon Till"]
    return env
end


text:
x_(i::Int) = :(x[$i])

"""
    batch_arithmetic_ops(op, operands)

Adds additional parenthes due to the following julia bug:
https://github.com/JuliaLang/julia/issues/14821
"""
function batch_arithmetic_ops(op::Symbol, ops)
    batch_size = 3
    l = 1
    if length(ops) == 1
        return :($(ops[1]))
    elseif length(ops) < batch_size + 2
        return Expr(:call, op, ops...)
    end
    batches = []
    while l ≤ length(ops)
        batch_end = min(length(ops), l + batch_size)
        if l == batch_end
            push!(batches, :($(ops[l])))
        else
            push!(batches, Expr(:call, op, ops[l:batch_end]...))
        end
        l = batch_end + 1
    end
    return length(batches) == 1 ? batches[1] : batch_arithmetic_ops(op, batches)
end

"""
    normalized_poly_coefficents(::Type{T}, degrees, coefficients)

Return the cofficients of an univariate polynomial in standard form.
Assume the degrees are in ascending form.

## Example
```julia
julia> normalized_poly_coefficents(Float64, [1, 3], [:c_3, :c_5])
[:(zero(Float64)), :c_3, :(zero(Float64)), :c_5]
```
"""
function normalized_coefficients(::Type{T}, degrees::AbstractVector, coefficients::AbstractVector) where T
    ops = []
    d = degrees[end]
    k = 0
    i = 1
    while k ≤ d
        if i ≤ length(degrees) && degrees[i] == k
            push!(ops, coefficients[i])
            i += 1
        else
            push!(ops, :(zero($T)))
        end
        k += 1
    end
    ops
end

"""
    monomial_product(::Type{T}, exponent, coefficient, i::Union{Nothing, Int}=nothing)

Generate the monomial product defined by `exponent` with `coefficient.`
If `i` is an `Int` the partial derivative will be generated.
"""
function monomial_product(::Type{T}, exponent, coefficient, i::Union{Nothing, Int}=nothing; access_input=x_) where T
    if i !== nothing && exponent[i] == 0
        return (:(zero($T)), true)
    end
    ops = []
    push!(ops, coefficient)
    for (k, e) in enumerate(exponent)
        if k == i && e == 1
            continue
        elseif k == i && e > 1
            pushfirst!(ops, :($e))
            push!(ops, pow(access_input(k), e - 1))
        elseif e > 0
            push!(ops, pow(access_input(k), e))
        end
    end
    (batch_arithmetic_ops(:*, ops), false)
end

function monomial_product_with_derivatives(::Type{T}, exponent, coefficient; access_input=x_) where T
    val, _ = monomial_product(T, exponent, coefficient, access_input=access_input)
    dvals = map(1:length(exponent)) do i
        monomial_product(T, exponent, coefficient, i, access_input=access_input)
    end

    val, dvals
end

function pow(expr::Union{Expr,Symbol}, k::Integer)
    if k == 2
        :($expr * $expr)
    elseif k == 3
        :(($expr * $expr * $expr))
    elseif k == 4
        @gensym p
        quote
            $p = $expr * $expr
            $p * $p
        end
    elseif k == 1
        :($expr)
    elseif k == 0
        :(one($expr))
    else
        :(pow($expr, $k))
    end
end


text:
@testset "Interpretations" begin
    ifd = TiffImages.IFD(UInt32)
    ifd[TiffImages.PHOTOMETRIC] = TiffImages.PHOTOMETRIC_MINISBLACK
    ifd[TiffImages.SAMPLESPERPIXEL] = 1

    @test TiffImages.interpretation(ifd) == (Gray, false)

    ifd[TiffImages.SAMPLESPERPIXEL] = 2
    @test TiffImages.interpretation(ifd) == (Gray, true)

    # handle case where SamplesPerPixel is missing, issue #56
    delete!(ifd, TiffImages.SAMPLESPERPIXEL)
    @test TiffImages.interpretation(ifd) == (Gray, false)

    ifd[TiffImages.SAMPLESPERPIXEL] = 2
    ifd[TiffImages.EXTRASAMPLES] = TiffImages.EXTRASAMPLE_ASSOCALPHA

    @test TiffImages.interpretation(ifd) == (GrayA, false)

    ifd[TiffImages.EXTRASAMPLES] = TiffImages.EXTRASAMPLE_UNSPECIFIED

    @test TiffImages.interpretation(ifd) == (Gray, true)

    ifd[TiffImages.SAMPLESPERPIXEL] = 3

    @test TiffImages.interpretation(ifd) == (Gray, true)

    ifd[TiffImages.PHOTOMETRIC] = TiffImages.PHOTOMETRIC_RGB

    @test TiffImages.interpretation(ifd) == (RGB, false)

    ifd[TiffImages.SAMPLESPERPIXEL] = 4

    @test TiffImages.interpretation(ifd) == (RGBX, false)

    ifd[TiffImages.EXTRASAMPLES] = TiffImages.EXTRASAMPLE_ASSOCALPHA

    @test TiffImages.interpretation(ifd) == (RGBA, false)
end


text:
using JSON

include("latex_symbols.jl")

result = [
	Dict(
		:label => chop(k, head = 1, tail = 0),
		:body => v,
		:description => v,
	) for (k, v) in latex_symbols
]

result = sort(result, by = e -> e[:label])

JSON.print(result, 4)


text:
## Packages and settings.
include("../src/StencilBuilder.jl")

using BenchmarkTools
using LoopVectorization
using .StencilBuilder

float_type = Float32

## CPU dynamics kernel.
macro make_kernel(do_advec::Bool, do_diff::Bool, do_buoy::Bool)
    ex_rhs_l = []
    if do_advec
        push!(ex_rhs_l, :( -gradx(interpz(u) * interpx(w)) ))
        push!(ex_rhs_l, :( -grady(interpz(v) * interpy(w)) ))
        push!(ex_rhs_l, :( -gradz(interpz(w) * interpz(w)) ))
    end
    if do_diff
        push!(ex_rhs_l, :( visc * (gradx(gradx(w))) ))
        push!(ex_rhs_l, :( visc * (grady(grady(w))) ))
        push!(ex_rhs_l, :( visc * (gradz(gradz(w))) ))
    end
    if do_buoy
        push!(ex_rhs_l, :( alpha*interpz(s) ))
    end

    if length(ex_rhs_l) > 0
        ex_rhs = Expr(:call, :+, ex_rhs_l...)

        ex = quote
            function dynamics_w_kernel!(
                    wt, u, v, w, s,
                    visc, alpha,
                    dxi, dyi, dzi, dzhi,
                    is, ie, js, je, ks, ke)

                @fast3d begin
                    @fd (wt, u, v, w, s) wt += $ex_rhs
                end
            end
        end
    else
        ex = quote
            function dynamics_w_kernel!(
                    wt, u, v, w, s,
                    visc, alpha,
                    dxi, dyi, dzi, dzhi,
                    is, ie, js, je, ks, ke)
            end
        end
    end

    return esc(ex)
end

## User input
itot = 384; jtot = 384; ktot = 384;
igc = 1; jgc = 1; kgc = 1;

do_advec = true; do_diff = true; do_buoy = true

icells = itot + 2igc; jcells = jtot + 2jgc; kcells = ktot + 2kgc
is = igc+1; ie = igc+itot; js = jgc+1; je = jgc+jtot; ks = kgc+1; ke = kgc+ktot

wt = rand(float_type, icells, jcells, kcells)
u = rand(float_type, icells, jcells, kcells)
v = rand(float_type, icells, jcells, kcells)
w = rand(float_type, icells, jcells, kcells)
s = rand(float_type, icells, jcells, kcells)
dxi = rand(float_type)
dyi = rand(float_type)
dzi = rand(float_type, kcells)
dzhi = rand(float_type, kcells)
visc = convert(float_type, 1.)
alpha = convert(float_type, 9.81/300)

## Run CPU kernel
@eval @make_kernel($do_advec, $do_diff, $do_buoy)

@btime dynamics_w_kernel!(
    $wt, $u, $v, $w, $s,
    $visc, $alpha,
    $dxi, $dyi, $dzi, $dzhi,
    $is, $ie, $js, $je, $ks, $ke)


text:
	#
	#	Script for calculation of far-field pattern in the liquid case (sphere)
	#

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# User configurable parameters
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	# Physical parameters
	rho10 = 5.00 ; # Density ratio
	k_0 = 4 ; # Wave number in media 0
	k_1 = 6 ; # Wave number in media 1
	a = 1 ; # Sphere radius
	theta_inc = pi/4 ; # Incidence angle in radians

	# Software parameters
	grid_size = 200 ; # Angular grid size
	N = 150 ; # Number of terms in the series solution

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Calculation
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	# Far-Field pattern calculation
	Pattern = Pattern_LiquidSphere( rho10, k_0, k_1, a, grid_size, N, theta_inc ) ;

	# Saving to disk
	writedlm( "Out.Pattern.dat", Pattern , '\t' ) ;


text:
using Gtk, Cairo, ColorTypes

function generategrays(n, screenwidth)
    verts = Vector{RGB}()
    hwidth = Int(ceil(screenwidth/n))
    for x in 00:Int(floor(0xff/(n-1))):0xff
        rgbgray = RGB(x/255, x/255, x/255)
        for i in 1:hwidth
            push!(verts, rgbgray)
        end
    end
    verts
end

function drawline(ctx, p1, p2, color, width)
    move_to(ctx, p1.x, p1.y)
    set_source(ctx, color)
    line_to(ctx, p2.x, p2.y)
    set_line_width(ctx, width)
    stroke(ctx)
end

const can = @GtkCanvas()
const win = GtkWindow(can, "Grayscale bars/Display", 400, 400)
fullscreen(win)  # start full screen, then reduce to regular window in 5 seconds.

draw(can) do widget
    ctx = getgc(can)
    h = height(can)
    w = width(can)
    gpoints = generategrays(8, w)
    for (i, x) in enumerate(0:w-1)
        drawline(ctx, Point(x, 0.25*h), Point(x, 0), gpoints[i], 1)
    end
    gpoints = reverse(generategrays(16, w))
    for (i, x) in enumerate(0:w-1)
        drawline(ctx, Point(x, 0.5*h), Point(x, 0.25*h), gpoints[i], 1)
    end
    gpoints = generategrays(32, w)
    for (i, x) in enumerate(0:w-1)
        drawline(ctx, Point(x, 0.75*h), Point(x, 0.5*h), gpoints[i], 1)
    end
    gpoints = reverse(generategrays(64, w))
    for (i, x) in enumerate(0:w-1)
        drawline(ctx, Point(x, h), Point(x, 0.75*h), gpoints[i], 1)
    end
end

show(can)
sleep(5)
unfullscreen(win)
const cond = Condition()
endit(w) = notify(cond)
signal_connect(endit, win, :destroy)
wait(cond)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""

The cylinder is banged at the initial time (i.e. it is given an initial velocity).
"""
module plate_with_crack_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using SymRCM
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot; @gp "clear"
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection
using ThreadedSparseCSR
using UnicodePlots
using InteractiveUtils
using BenchmarkTools
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection, vtkwrite

const E = 72.7*phun("GPa");
const nu = 0.33;
const rho = 2700.0*phun("kg/m^3")
const d1 = 500*phun("mm")
const d2 = 50*phun("mm")
const d3 = 750*phun("mm")
const d4 = 250*phun("mm")
const thickness = 5*phun("mm")
const ksi = 0.0
const omegad = 1000*phun("Hz")
const carrier_frequency = 75*phun("kilo*Hz")
const modulation_frequency = carrier_frequency/4
const forcepatchradius = 20*phun("mm")
const forcedensity = 1*phun("N")
const color = "red"
const tend = 0.9*phun("milli*s")
const visualize = true

function parloop_csr!(M, K, ksi, U0, V0, tend, dt, force!, peek, nthr)
    U = deepcopy(U0)
    V = deepcopy(U0)
    A = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    C = deepcopy(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = deepcopy(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invMC .* force!(F, t);
    peek(0, U, t)
    @time for step in 1:nsteps
        # Displacement update
        @. U += dt*V + ((dt^2)/2)*A;
        # External loading
        force!(F, t);
        # Add elastic restoring forces
        F .-= ThreadedSparseCSR.bmul!(E, K, U)
        @inbounds @simd for i in eachindex(U)
            _Fi = F[i]; _Ai = A[i]; _Vi = V[i]
            # Add damping forces
            _Fi -= C[i] * (_Vi + (dt/2) * _Ai)
            # Compute the new acceleration.
            _A1i = invMC[i] * _Fi
            A[i] = _A1i
            # Update the velocity
            V[i] += (dt/2)* (_Ai + _A1i);
        end
        t = t + dt
        peek(step, U, t)
    end
end

function _execute_parallel_csr(nref = 2, nthr = 0, color = "red")
    tolerance = d2/nref/100

    d5 = (d1 - d2) / 2
    d6 = (d3 - d4)
    nd5 = 5
    xs = sort(unique(vcat(
        linearspace(0.0, d5, nd5),
        linearspace(d5, d5+d2, 2),
        linearspace(d5+d2, d1, nd5)
        )))
    ys = collect(linearspace(0.0, d4, nd5))

    fens1, fes1 = T3blockx(xs, ys); # Mesh
    for r in 1:nref
        fens1, fes1 = T3refine(fens1, fes1)
    end
    @show count(fens1), count(fes1)

    ys = collect(linearspace(d4, d3, nd5*2))

    fens2, fes2 = T3blockx(xs, ys); # Mesh
    for r in 1:nref
        fens2, fes2 = T3refine(fens2, fes2)
    end
    @show count(fens2), count(fes2)

    fens, fes1, fes2 = mergemeshes(fens1, fes1,  fens2, fes2, 0.0)
    fes = cat(fes1, fes2)
    @show count(fens), count(fes)

    offset = max(d5 / 4 / nref / 10, 2*tolerance)
    l1 = selectnode(fens; box = [0 d5+offset d4 d4], inflate = tolerance)
    l2 = selectnode(fens; box = [d5+d2-offset d1 d4 d4], inflate = tolerance)

    fens.xyz = xyz3(fens)

    candidates = vcat(l1, l2)
    fens, fes = mergenodes(fens, fes, tolerance, candidates)
    bfes = meshboundary(fes)
    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    @show count(fens)
    vtkwrite("plate_with_crack.vtu", fens, fes)
    vtkwrite("plate_with_crack-boundary.vtu", fens, bfes)

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # No EBC's
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            ThreadedSparseCSR.bmul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * (K * v)) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max = max(omega_max, 20*2*pi*carrier_frequency)
    @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    mpoint = selectnode(fens; nearestto=[d1/2 0 0.0])[1]
    cpoint = selectnode(fens; nearestto=[d1/2 d3 0])[1]

    mpointdof = dchi.dofnums[mpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]

    # Four cycles of the carrier frequency

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[mpoint, 1]
        dy = XYZ[2] - fens.xyz[mpoint, 2]
        dz = XYZ[3] - fens.xyz[mpoint, 3]
        forceout[1:2] .= 0.0
        forceout[3] = forcedensity*exp(-20*sqrt(dx^2+dy^2+dz^2)/forcepatchradius)
        forceout[4:6] .= 0.0
        return forceout
    end

    # Sinusoidal loading on the surface of the shell
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        mul = 0.0
        if t <= 4/carrier_frequency
            mul = 0.5 * (1 - cos(2*pi*modulation_frequency*t)) * sin(2*pi*carrier_frequency*t)
        end
        F .= mul .* Fmag
        return F
    end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    mdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        mdeflections[step+1] = U[mpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"
    parloop_csr!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek, nthr)

    if visualize
        # @gp  "set terminal windows 0 "  :-
        # @gp "clear"
        @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-
        @gp  :- collect(0.0:dt:(nsteps*dt)) mdeflections " lw 4 lc rgb '$color' with lines title 'Deflection at the source' "  :-
        @gp  :- "set xlabel 'Time'" :-
        @gp  :- "set ylabel 'Deflection'" :-
        @gp  :- "set title 'Free-floating plate'"

        # Visualization
        @info "Dumping visualization"
        times = Float64[]
        vectors = []
        for i in 1:length(displacements)
            scattersysvec!(dchi, displacements[i])
            push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
            push!(times, i*dt*nbtw)
        end
        vtkwritecollection("plate_with_crack_$nref", fens, fes, times; vectors = vectors)
    end
end

function test_parallel_csr(nrefs = [4], nthr = 0, color = "red")
    @info "Cracked Plate, nrefs = $nrefs: parallel CSR"
    for nref in nrefs
        _execute_parallel_csr(nref, nthr, color)
    end
    return true
end

function allrun(nrefs = [4], nthr = 0)
    println("#####################################################")
    println("# test_parallel_csr ")
    test_parallel_csr(nrefs, nthr, "blue")
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
# DensityEstimation.jl

export getEyPy!;

function getEyPy!(sims::Array, X_grid::Matrix{T},
    X::Matrix{T}, D::Vector{Matrix{T}}; densout::Bool=false, y_grid::Vector{T}=[0.0],
    rng::MersenneTwister=MersenneTwister(), tol=1.0e-9) where T <: Real

    ## modifies the rng passed to it

    n_star, L = size(X_grid)
    n_y = length(y_grid)

    L == (length(sims[1][:lλ])-1) == length(D) == size(X,2) || throw("L is inconsistent in one of the args.")

    nsim = length(sims)

    Dstar = [ pairDistMat( X_grid[:,[j]] ) for j = 1:L ]
    Dstdat = [ pairDistMat( X_grid[:,[j]],  X[:,[j]] ) for j = 1:L ]

    Fstar = Array{Float64}(undef, L, n_star, nsim)

    for j = 1:L
        Fstar[j,:,:] = hcat([ sims[ii][:mixcomps][j][:μ] .+ rfullcond_fstar!(sims[ii][:mixcomps][j][:fx],
            covMat(D[j], sims[ii][:mixcomps][j][:κ]*sims[ii][:mixcomps][j][:σ2], sims[ii][:mixcomps][j][:corParams]),
            covMat(Dstar[j], sims[ii][:mixcomps][j][:κ]*sims[ii][:mixcomps][j][:σ2], sims[ii][:mixcomps][j][:corParams]),
            covMat(Dstdat[j], sims[ii][:mixcomps][j][:κ]*sims[ii][:mixcomps][j][:σ2], sims[ii][:mixcomps][j][:corParams]),
            rng, tol) for ii = 1:nsim ]...)
    end

    Ey = Matrix{Float64}(undef, n_star, nsim)
    if densout
        Py = Array{Float64}(undef, n_y, n_star, nsim)
    end

    for ii = 1:nsim
        for i = 1:n_star

            λ = exp.(sims[ii][:lλ])
            Ey[i,ii] = λ[1] * sims[ii][:intercept].μ + sum( λ[2:(L+1)] .* Fstar[:,i,ii] )

            if densout
                lp = Vector{Float64}(undef, L+1)
                for iii = 1:n_y
                    lp[1] = sims[ii][:lλ][1] - 0.5*log(2π * sims[ii][:intercept].σ2) - 0.5*( y_grid[iii] - sims[ii][:intercept].μ )^2 / sims[ii][:intercept].σ2
                    for j = 1:L
                        lp[j+1] = sims[ii][:lλ][j+1] - 0.5*log(2π * sims[ii][:mixcomps][j][:σ2]) - 0.5*( y_grid[iii] - Fstar[j,i,ii] )^2 / sims[ii][:mixcomps][j][:σ2]
                    end
                    Py[iii,i,ii] = exp( logsumexp(lp) )
                end
            end
        end
    end

    if densout
        return Ey, Py
    else
        return Ey
    end

end


text:
function copy_self(obj::SubHyperplane)
    return jcall(obj, "copySelf", SubHyperplane, ())
end

function get_hyperplane(obj::SubHyperplane)
    return jcall(obj, "getHyperplane", Hyperplane, ())
end

function get_size(obj::SubHyperplane)
    return jcall(obj, "getSize", jdouble, ())
end

function is_empty(obj::SubHyperplane)
    return jcall(obj, "isEmpty", jboolean, ())
end

function reunite(obj::SubHyperplane, arg0::SubHyperplane)
    return jcall(obj, "reunite", SubHyperplane, (SubHyperplane,), arg0)
end

function split(obj::SubHyperplane, arg0::Hyperplane)
    return jcall(obj, "split", SubHyperplane_SplitSubHyperplane, (Hyperplane,), arg0)
end