text:
using DataInterpolations, Test
using StableRNGs
import ForwardDiff

@testset "Linear Interpolation" begin
    u = 2.0collect(1:10)
    t = 1.0collect(1:10)
    A = LinearInterpolation(u,t)

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(0)        == 0.0
    @test A(5.5)      == 11.0
    @test A(11)       == 22

    u = vcat(2.0collect(1:10)', 3.0collect(1:10)')
    A = LinearInterpolation(u,t)

    for (_t, _u) in zip(t, eachcol(u))
        @test A(_t) == _u
    end
    @test A(0)        == [0.0, 0.0]
    @test A(5.5)      == [11.0, 16.5]
    @test A(11)       == [22, 33]

    x = 1:10
    y = 2:4
    u_= x' .* y
    u = [u_[:,i] for i = 1:size(u_,2)]
    A = LinearInterpolation(u,t)
    @test A(0)        == [0.0, 0.0, 0.0]
    @test A(5.5)      == [11.0, 16.5, 22.0]
    @test A(11)       == [22.0, 33.0, 44.0]

    u = [u_[:,i:i+1] for i = 1:2:10]
    t = 1.0collect(2:2:10)
    A = LinearInterpolation(u,t)

    @test A(0)        == [-2.0 0.0; -3.0 0.0; -4.0 0.0]
    @test A(3)        == [4.0 6.0; 6.0 9.0; 8.0 12.0]
    @test A(5)        == [8.0 10.0; 12.0 15.0; 16.0 20.0]

    # with NaNs (#113)
    u = [NaN, 1.0, 2.0, 3.0]
    t = 1:4
    A = LinearInterpolation(u, t)
    @test isnan(A(1.0))
    @test A(2.0) == 1.0
    @test A(2.5) == 1.5
    @test A(3.0) == 2.0
    @test A(4.0) == 3.0

    u = [0.0, NaN, 2.0, 3.0]
    A = LinearInterpolation(u, t)
    @test A(1.0) == 0.0
    @test isnan(A(2.0))
    @test isnan(A(2.5))
    @test A(3.0) == 2.0
    @test A(4.0) == 3.0

    u = [0.0, NaN, 2.0, 3.0]
    A = LinearInterpolation(u, t)
    @test A(1.0) == 0.0
    @test isnan(A(2.0))
    @test isnan(A(2.5))
    @test A(3.0) == 2.0
    @test A(4.0) == 3.0

    u = [0.0, 1.0, NaN, 3.0]
    A = LinearInterpolation(u, t)
    @test A(1.0) == 0.0
    @test A(2.0) == 1.0
    @test isnan(A(2.5))
    @test isnan(A(3.0))
    @test A(4.0) == 3.0

    u = [0.0, 1.0, 2.0, NaN]
    A = LinearInterpolation(u, t)
    @test A(1.0) == 0.0
    @test A(2.0) == 1.0
    @test A(3.0) == 2.0
    @test isnan(A(3.5))
    @test isnan(A(4.0))

    # Test type stability
    u = Float32.(1:5)
    t = Float32.(1:5)
    A1 = LinearInterpolation(u, t)
    u = 1:5
    t = 1:5
    A2 = LinearInterpolation(u, t)
    u = [1//i for i in 1:5]
    t = (1:5)
    A3 = LinearInterpolation(u, t)
    u = [1//i for i in 1:5]
    t = [1//(6-i) for i in 1:5]
    A4 = LinearInterpolation(u, t)

    F32 = Float32(1)
    F64 = Float64(1)
    I32 = Int32(1)
    I64 = Int64(1)
    R32 = Int32(1)//Int32(1)
    R64 = 1//1
    for A in [A1, A2, A3, A4]
        @test @inferred(A(F32)) === A(F32)
        @test @inferred(A(F64)) === A(F64)
        @test @inferred(A(I32)) === A(I32)
        @test @inferred(A(I64)) === A(I64)
        @test @inferred(A(R32)) === A(R32)
        @test @inferred(A(R64)) === A(R64)
    end

    # Nan time value:
    t = 0.0:3;  # Floats
    u = [0, -2, -1, -2];
    A = LinearInterpolation(u, t);
    dA = t -> ForwardDiff.derivative(A, t)
    @test isnan(dA(NaN))

    t = 0:3;  # Integers
    u = [0, -2, -1, -2];
    A = LinearInterpolation(u, t);
    dA = t -> ForwardDiff.derivative(A, t)
    @test isnan(dA(NaN))

    # Test derivative at point gives derivative to the right (except last is to left):
    ts = t[begin:end-1]
    @test dA.(ts) == dA.(ts .+ 0.5)
    # Test last derivitive is to the left:
    @test dA(last(t)) == dA(last(t) - 0.5)
end

@testset "Quadratic Interpolation" begin
    u = [1.0, 4.0, 9.0, 16.0]
    t = [1.0, 2.0, 3.0, 4.0]
    A = QuadraticInterpolation(u,t)

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(0.0) == 0.0
    @test A(1.5) == 2.25
    @test A(2.5) == 6.25
    @test A(3.5) == 12.25
    @test A(5.0) == 25

    u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
    A = QuadraticInterpolation(u,t)

    for (_t, _u) in zip(t, eachcol(u))
        @test A(_t) == _u
    end
    @test A(0.0) == [0.0  ,  0.0 ]
    @test A(1.5) == [2.25 ,  2.25]
    @test A(2.5) == [6.25 ,  6.25]
    @test A(3.5) == [12.25, 12.25]
    @test A(5.0) == [25.0 , 25.0 ]

    u_= [1.0, 4.0, 9.0, 16.0]' .* ones(5)
    u = [u_[:,i] for i = 1:size(u_,2)]
    A = QuadraticInterpolation(u,t)
    @test A(0)   == zeros(5)
    @test A(1.5) == 2.25 * ones(5)
    @test A(2.5) == 6.25 * ones(5)
    @test A(3.5) == 12.25 * ones(5)
    @test A(5.0) == 25.0 * ones(5)

    u = [repeat(u[i], 1, 3) for i=1:4]
    A = QuadraticInterpolation(u,t)
    @test A(0)   == zeros(5, 3)
    @test A(1.5) == 2.25 * ones(5, 3)
    @test A(2.5) == 6.25 * ones(5, 3)
    @test A(3.5) == 12.25 * ones(5, 3)
    @test A(5.0) == 25.0 * ones(5, 3)
end

@testset "Lagrange Interpolation" begin
    u = [1.0, 4.0, 9.0]
    t = [1.0, 2.0, 3.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 4.0
    @test A(1.5) == 2.25

    u = [1.0, 8.0, 27.0, 64.0]
    t = [1.0, 2.0, 3.0, 4.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 8.0
    @test A(1.5) ≈ 3.375
    @test A(3.5) ≈ 42.875

    u = [1.0 4.0 9.0 16.0; 1.0 4.0 9.0 16.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == [4.0,4.0]
    @test A(1.5) ≈ [2.25,2.25]
    @test A(3.5) ≈ [12.25,12.25]

    u_= [1.0, 4.0, 9.0]' .* ones(4)
    u = [u_[:,i] for i = 1:size(u_,2)]
    t = [1.0, 2.0, 3.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 4.0 * ones(4)
    @test A(1.5) == 2.25 * ones(4)

    u_= [1.0, 8.0, 27.0, 64.0]' .* ones(4)
    u = [u_[:,i] for i = 1:size(u_,2)]
    t = [1.0, 2.0, 3.0, 4.0]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 8.0 * ones(4)
    @test A(1.5) ≈ 3.375 * ones(4)
    @test A(3.5) ≈ 42.875 * ones(4)

    u = [repeat(u[i], 1, 3) for i=1:4]
    A = LagrangeInterpolation(u,t)

    @test A(2.0) == 8.0 * ones(4, 3)
    @test A(1.5) ≈ 3.375 * ones(4, 3)
    @test A(3.5) ≈ 42.875 * ones(4, 3)
end

@testset "Akima Interpolation" begin
    u = [0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0]
    t = collect(0.0:10.0)
    A = AkimaInterpolation(u, t)

    @test A(0.0) ≈ 0.0
    @test A(0.5) ≈ 1.375
    @test A(1.0) ≈ 2.0
    @test A(1.5) ≈ 1.5
    @test A(2.5) ≈ 1.953125
    @test A(3.5) ≈ 2.484375
    @test A(4.5) ≈ 4.1363636363636366866103344
    @test A(5.1) ≈ 5.9803623910336236590978842
    @test A(6.5) ≈ 5.5067291516462386624652936
    @test A(7.2) ≈ 5.2031367459745245795943447
    @test A(8.6) ≈ 4.1796554159017080820603951
    @test A(9.9) ≈ 3.4110386597938129327189927
    @test A(10.0) ≈ 3.0
end

@testset "ConstantInterpolation" begin

    t = [1.0, 2.0, 3.0, 4.0]

    @testset "Vector case" for u in
        [[1.0, 2.0, 0.0, 1.0], ["B", "C", "A", "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[2]
        @test A(2.0) == u[2]
        @test A(2.5) == u[3]
        @test A(3.0) == u[3]
        @test A(3.5) == u[1]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[1]
        @test A(2.0) == u[2]
        @test A(2.5) == u[2]
        @test A(3.0) == u[3]
        @test A(3.5) == u[3]
        @test A(4.0) == u[1]
        @test A(4.5) == u[1]
    end

    @testset "Matrix case" for u in
        [[1.0 2.0 0.0 1.0; 1.0 2.0 0.0 1.0], ["B" "C" "A" "B"; "B" "C" "A" "B"]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,2]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,3]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,1]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[:,1]
        @test A(1.0) == u[:,1]
        @test A(1.5) == u[:,1]
        @test A(2.0) == u[:,2]
        @test A(2.5) == u[:,2]
        @test A(3.0) == u[:,3]
        @test A(3.5) == u[:,3]
        @test A(4.0) == u[:,1]
        @test A(4.5) == u[:,1]
    end

    @testset "Vector of Vectors case" for u in
        [[[1.0, 2.0], [0.0, 1.0], [1.0, 2.0], [0.0, 1.0]],
        [["B", "C"], ["A", "B"], ["B", "C"], ["A", "B"]]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[2]
        @test A(2.0) == u[2]
        @test A(2.5) == u[3]
        @test A(3.0) == u[3]
        @test A(3.5) == u[4]
        @test A(4.0) == u[4]
        @test A(4.5) == u[4]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[1]
        @test A(2.0) == u[2]
        @test A(2.5) == u[2]
        @test A(3.0) == u[3]
        @test A(3.5) == u[3]
        @test A(4.0) == u[4]
        @test A(4.5) == u[4]
    end

    @testset "Vector of Matrices case" for u in
        [[[1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0], [1.0 2.0; 1.0 2.0], [0.0 1.0; 0.0 1.0]],
        [["B" "C"; "B" "C"], ["A" "B"; "A" "B"], ["B" "C"; "B" "C"], ["A" "B"; "A" "B"]]]

        A = ConstantInterpolation(u, t, dir=:right)
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[2]
        @test A(2.0) == u[2]
        @test A(2.5) == u[3]
        @test A(3.0) == u[3]
        @test A(3.5) == u[4]
        @test A(4.0) == u[4]
        @test A(4.5) == u[4]

        A = ConstantInterpolation(u, t) # dir=:left is default
        @test A(0.5) == u[1]
        @test A(1.0) == u[1]
        @test A(1.5) == u[1]
        @test A(2.0) == u[2]
        @test A(2.5) == u[2]
        @test A(3.0) == u[3]
        @test A(3.5) == u[3]
        @test A(4.0) == u[4]
        @test A(4.5) == u[4]
    end
end

@testset "QuadraticSpline Interpolation" begin
    u = [0.0, 1.0, 3.0]
    t = [-1.0, 0.0, 1.0]

    A = QuadraticSpline(u,t)

    # Solution
    P₁ = x -> (x + 1)^2 # for x ∈ [-1, 0]
    P₂ = x -> 2*x + 1   # for x ∈ [ 0, 1]

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    @test A(-2.0) == P₁(-2.0)
    @test A(-0.5) == P₁(-0.5)
    @test A(0.7)  == P₂( 0.7)
    @test A(2.0)  == P₂( 2.0)

    u_= [0.0, 1.0, 3.0]' .* ones(4)
    u = [u_[:,i] for i = 1:size(u_,2)]
    A = QuadraticSpline(u,t)
    @test A(-2.0) == P₁(-2.0) * ones(4)
    @test A(-0.5) == P₁(-0.5) * ones(4)
    @test A(0.7)  == P₂( 0.7) * ones(4)
    @test A(2.0)  == P₂( 2.0) * ones(4)

    u = [repeat(u[i], 1, 3) for i=1:3]
    A = QuadraticSpline(u,t)
    @test A(-2.0) == P₁(-2.0) * ones(4, 3)
    @test A(-0.5) == P₁(-0.5) * ones(4, 3)
    @test A(0.7)  == P₂( 0.7) * ones(4, 3)
    @test A(2.0)  == P₂( 2.0) * ones(4, 3)
end

@testset "CubicSpline Interpolation" begin
    u = [0.0, 1.0, 3.0]
    t = [-1.0, 0.0, 1.0]

    A = CubicSpline(u,t)

    # Solution
    P₁ = x -> 1 + 1.5x + x^2 + 0.5x^3 # for x ∈ [-1.0, 0.0]
    P₂ = x -> 1 + 1.5x + x^2 - 0.5x^3 # for x ∈ [0.0, 1.0]

    for (_t, _u) in zip(t, u)
        @test A(_t) == _u
    end
    for x in (-1.5, -0.5, -0.7)
        @test A(x) ≈ P₁(x)
    end
    for x in (0.3, 0.5, 1.5)
        @test A(x) ≈ P₂(x)
    end

    u_= [0.0, 1.0, 3.0]' .* ones(4)
    u = [u_[:,i] for i = 1:size(u_,2)]
    A = CubicSpline(u,t)
    for x in (-1.5, -0.5, -0.7)
        @test A(x) ≈ P₁(x) * ones(4)
    end
    for x in (0.3, 0.5, 1.5)
        @test A(x) ≈ P₂(x) * ones(4)
    end

    u = [repeat(u[i], 1, 3) for i=1:3]
    A = CubicSpline(u,t)
    for x in (-1.5, -0.5, -0.7)
        @test A(x) ≈ P₁(x) * ones(4, 3)
    end
    for x in (0.3, 0.5, 1.5)
        @test A(x) ≈ P₂(x) * ones(4, 3)
    end
end

# BSpline Interpolation and Approximation
t = [0,62.25,109.66,162.66,205.8,252.3]
u = [14.7,11.51,10.41,14.95,12.24,11.22]

A = BSplineInterpolation(u,t,2,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
@test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]
@test [A(t[1]), A(t[end])] == [u[1], u[end]]

A = BSplineInterpolation(u,t,2,:ArcLen,:Average)

@test [A(25.0), A(80.0)] == [13.363814458968486, 10.685201117692609]
@test [A(190.0), A(225.0)] == [13.437481084762863, 11.367034741256463]
@test [A(t[1]), A(t[end])] == [u[1], u[end]]

A = BSplineApprox(u,t,2,4,:Uniform,:Uniform)

@test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
@test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575]
@test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]

# Curvefit Interpolation
rng = StableRNG(12345)
model(x, p) = @. p[1]/(1 + exp(x - p[2]))
t = range(-10, stop=10, length=40)
u = model(t, [1.0, 2.0]) + 0.01*randn(rng, length(t))
p0 = [0.5, 0.5]

A = Curvefit(u, t, model, p0, LBFGS())

ts = [-7.0, -2.0, 0.0, 2.5, 5.0]
vs = [1.0013468217936277, 0.9836755196317837, 0.8833959853995836, 0.3810348276782708, 0.048062978598861855]
us = A.(ts)

@test vs ≈ us

# missing values handling tests
u = [1.0, 4.0, 9.0, 16.0, 25.0, missing, missing]
t = [1.0, 2.0, 3.0, 4.0, missing, 6.0, missing]
A = QuadraticInterpolation(u,t)

@test A(2.0) == 4.0
@test A(1.5) == 2.25
@test A(3.5) == 12.25
@test A(2.5) == 6.25

u = copy(hcat(u, u)')
A = QuadraticInterpolation(u,t)

@test A(2.0) == [4.0, 4.0]
@test A(1.5) == [2.25, 2.25]
@test A(3.5) == [12.25,12.25]
@test A(2.5) == [6.25, 6.25]


text:
module Magnetic

using ModelingToolkit

# FluxTubes
include("FluxTubes/FluxTubes.jl")

# QuasiStatic

# FundamentalWave

end #module

ext: jl


text:
abstract type Segway end

################################# Mutable version ###########################

mutable struct Segway_mutable<:Segway
    Mb::Float64 # kg # Masse des Körpers
    Mw::Float64 # kg # Masse der Räder
    R::Float64 # m #Radius der Räder
    cz::Float64 # m #Schwerpunkt des Körpers Höhe
    b::Float64 # m #Half the distance between the wheels (OO_w in paper)
    Ixx::Float64 # kg m^2 #Inertia of the body
    Iyy::Float64 # kg m^2 #Inertia of the body
    Izz::Float64 # kg m^2 #Inertia of the body
    Iwa::Float64 # kg m^2 #Inertia of the wheels about their axis
    Iwd::Float64 # kg m^2 #Inertia of the wheels about a diameter
end

function Segway_mutable(;
    Mb::Float64 = 1.76, # kg # Masse des Körpers
    Mw::Float64 = 0.147, # kg # Masse der Räder
    R::Float64 = 0.07, # m #Radius der Räder
    cz::Float64 = 0.077, # m #Schwerpunkt des Körpers Höhe über Radachse
    b::Float64 = 0.1985/2, # m #Half the distance between the wheels (OO_w in paper)
    Ixx::Float64 = (0.166^2+0.21^2)*Mb/12 + Mb*cz^2, # kg m^2 #Inertia of the body
    Iyy::Float64 = (0.072^2+0.21^2)*Mb/12 + Mb*cz^2, # kg m^2 #Inertia of the body
    Izz::Float64 = (0.072^2+0.166^2)*Mb/12, # kg m^2 #Inertia of the body
    Iwa::Float64 = Mw*R^2/2, # kg m^2 #Inertia of the wheels about their axis
    Iwd::Float64 = Mw*b^2) # kg m^2 #Inertia of the wheels about their axis)

    return Segway_mutable(Mb, Mw, R, cz, b, Ixx, Iyy, Izz, Iwa, Iwd)
end

function Segway_mutable(params::Vector{Float64}, x::Vector{Float64}=zeros(7))
    @assert(length(params) == 10)
    return Segway_mutable(params...)
end

################################# Unmutable version ###########################
struct Segway_unmutable<:Segway
    Mb::Float64 # kg # Masse des Körpers
    Mw::Float64 # kg # Masse der Räder
    R::Float64 # m #Radius der Räder
    cz::Float64 # m #Schwerpunkt des Körpers Höhe
    b::Float64 # m #Half the distance between the wheels (OO_w in paper)
    Ixx::Float64 # kg m^2 #Inertia of the body
    Iyy::Float64 # kg m^2 #Inertia of the body
    Izz::Float64 # kg m^2 #Inertia of the body
    Iwa::Float64 # kg m^2 #Inertia of the wheels about their axis
    Iwd::Float64 # kg m^2 #Inertia of the wheels about a diameter
end

function Segway_unmutable(;
    Mb::Float64 = 1.76, # kg # Masse des Körpers
    Mw::Float64 = 0.147, # kg # Masse der Räder
    R::Float64 = 0.07, # m #Radius der Räder
    cz::Float64 = 0.077, # m #Schwerpunkt des Körpers Höhe über Radachse
    b::Float64 = 0.1985/2, # m #Half the distance between the wheels (OO_w in paper)
    Ixx::Float64 = (0.166^2+0.21^2)*Mb/12 + Mb*cz^2, # kg m^2 #Inertia of the body
    Iyy::Float64 = (0.072^2+0.21^2)*Mb/12 + Mb*cz^2, # kg m^2 #Inertia of the body
    Izz::Float64 = (0.072^2+0.166^2)*Mb/12, # kg m^2 #Inertia of the body
    Iwa::Float64 = Mw*R^2/2, # kg m^2 #Inertia of the wheels about their axis
    Iwd::Float64 = Mw*b^2) # kg m^2 #Inertia of the wheels about their axis)

    return Segway_unmutable(Mb, Mw, R, cz, b, Ixx, Iyy, Izz, Iwa, Iwd)
end

function Segway_unmutable(params::Vector{Float64}, x::Vector{Float64}=zeros(7))
    @assert(length(params) == 10)
    return Segway_unmutable(params...)
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using Plots
using DataStructures
using Statistics
using ProgressBars
L=200
function BoundaryCondition(i)
    if i==L+1
        return 1
    elseif i==0
        return L
    else
        return i
    end
end
function height_cal(i,arr) #calculating the height of the layer in a column
    row=1
    height=0
    for j in 1:L
        if arr[row,i]!=0.0
            height+=1
            row+=1
        else
            break #exit the loop for the next column
        end
    end
    return height
end
function ColumnChoosing(arr, col) #choosing which column to fall into
    h_dict=DataStructures.OrderedDict() #an ordered dictionary : height=>column_num
    for i in (col-1):1:(col+1)
        i=BoundaryCondition(i)
        height= height_cal(i,arr)
        h_dict[height]=i
    end
    hArr=collect(keys(h_dict)) #having all the heights in an array
    return get(h_dict,min(hArr...), "ERROR") #returning the column respectively
end
function mean_calculater(arr,L)
    h=[]
    for i in 1:L
        i=BoundaryCondition(i)
        height= height_cal(i,arr)
        push!(h, height)
    end
    return mean(h)
end
function deposing(L,n)
    color=1
    meanList=[]
    arr=zeros((L,L))
    anim = Animation()
    for particle in ProgressBar(1:n)
        column=rand(1:L)
        col=ColumnChoosing(arr, column)
        row=1
        for j in 1:L
            if arr[row,col]==0.0
                arr[row,col]=color
                break
            else
                row+=1
            end
        end
        if particle%(10*L*color)==0
            color+=1
        end
        p = heatmap(hcat(arr))
        frame(anim, p)
        #meanNum=mean_calculater(arr,L)
        #push!(meanList, meanNum)
    end
    gif(anim,"C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\bdw.gif",fps=80)
    return arr, meanList
end
arr,meanList=deposing(L,30000) #for plotting the dynamics and mean: n=30000

heatmap(hcat(arr), c=cgrad(:roma, 10, categorical = true, scale = :exp), xlabel="L")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\RandomBallisticDepositionWithRelaxation.png")
scatter(1:30000,meanList, xlabel="time", ylabel="mean height of the layer")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter3\\Fig\\RandomBallisticDepositionWithRelaxation2.png")


text:
@testset "Conv Layers" begin
    RTOL_LOW = 1e-2
    RTOL_HIGH = 1e-5

    in_channel = 3
    out_channel = 5
    N = 4
    T = Float32

    adj1 =  [0 1 0 1
             1 0 1 0
             0 1 0 1
             1 0 1 0]

    g1 = GNNGraph(adj1,
            ndata=rand(T, in_channel, N),
            graph_type=GRAPH_T)

    adj_single_vertex =  [0 0 0 1
                          0 0 0 0
                          0 0 0 1
                          1 0 1 0]

    g_single_vertex = GNNGraph(adj_single_vertex,
                                ndata=rand(T, in_channel, N),
                                graph_type=GRAPH_T)

    test_graphs = [g1, g_single_vertex]

    @testset "GCNConv" begin
        l = GCNConv(in_channel => out_channel)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        l = GCNConv(in_channel => out_channel, tanh, bias=false)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        l = GCNConv(in_channel => out_channel, add_self_loops=false)
        test_layer(l, g1, rtol=RTOL_HIGH, outsize=(out_channel, g1.num_nodes))

        @testset "edge weights" begin
            s = [2,3,1,3,1,2]
            t = [1,1,2,2,3,3]
            w = T[1,2,3,4,5,6]
            g = GNNGraph((s, t, w), ndata=ones(T, 1, 3), graph_type=GRAPH_T)
            x = g.ndata.x
            l = GCNConv(1 => 1, add_self_loops=false, use_edge_weight=true)
            l.weight .= 1
            d = degree(g, dir=:in, edge_weight=true)
            y = l(g, x)
            @test y[1,1] ≈ w[1] / √(d[1]*d[2]) + w[2] / √(d[1]*d[3])
            @test y[1,2] ≈ w[3] / √(d[2]*d[1]) + w[4] / √(d[2]*d[3])
            @test y ≈ l(g, x, w)

            # test gradient with respect to edge weights
            w = rand(T, 6)
            x = rand(T, 1, 3)
            g = GNNGraph((s, t, w), ndata=x, graph_type=GRAPH_T, edata=w)
            l = GCNConv(1 => 1, add_self_loops=false, use_edge_weight=true)
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(1, g.num_nodes), test_gpu=false)
            @test gradient(w -> sum(l(g, x, w)), w)[1] isa AbstractVector{T}   # redundan test but more esplicit
        end
    end

    @testset "ChebConv" begin
        k = 2
        l = ChebConv(in_channel => out_channel, k)
        @test size(l.weight) == (out_channel, in_channel, k)
        @test size(l.bias) == (out_channel,)
        @test l.k == k
        for g in test_graphs
            g = add_self_loops(g)
            test_layer(l, g, rtol=RTOL_HIGH, test_gpu=false, outsize=(out_channel, g.num_nodes))
            if TEST_GPU
                @test_broken test_layer(l, g, rtol=RTOL_HIGH, test_gpu=true, outsize=(out_channel, g.num_nodes))
            end
        end

        @testset "bias=false" begin
            @test length(Flux.params(ChebConv(2=>3, 3))) == 2
            @test length(Flux.params(ChebConv(2=>3, 3, bias=false))) == 1
        end
    end

    @testset "GraphConv" begin
        l = GraphConv(in_channel => out_channel)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        l = GraphConv(in_channel => out_channel, relu, bias=false, aggr=mean)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        @testset "bias=false" begin
            @test length(Flux.params(GraphConv(2=>3))) == 3
            @test length(Flux.params(GraphConv(2=>3, bias=false))) == 2
        end
    end

    @testset "GATConv" begin

        for heads in (1, 2), concat in (true, false)
            l = GATConv(in_channel => out_channel; heads, concat)
            for g in test_graphs
                test_layer(l, g, rtol=RTOL_LOW,
                    exclude_grad_fields = [:negative_slope],
                    outsize=(concat ? heads*out_channel : out_channel, g.num_nodes))
            end
        end

        @testset "edge features" begin
            ein = 3
            l = GATConv((in_channel, ein) => out_channel, add_self_loops=false)
            g = GNNGraph(g1, edata=rand(T, ein, g1.num_edges))
            test_layer(l, g, rtol=RTOL_LOW,
                exclude_grad_fields = [:negative_slope],
                outsize=(out_channel, g.num_nodes))
        end

        @testset "num params" begin
            l = GATConv(2 => 3, add_self_loops=false)
            @test length(Flux.params(l)) == 3
            l = GATConv((2,4) => 3, add_self_loops=false)
            @test length(Flux.params(l)) == 4
            l = GATConv((2,4) => 3, add_self_loops=false, bias=false)
            @test length(Flux.params(l)) == 3
        end
    end

    @testset "GATv2Conv" begin

        for heads in (1, 2), concat in (true, false)
            l = GATv2Conv(in_channel => out_channel, tanh; heads, concat)
            for g in test_graphs
                test_layer(l, g, rtol=RTOL_LOW,
                    exclude_grad_fields = [:negative_slope],
                    outsize=(concat ? heads*out_channel : out_channel, g.num_nodes))
            end
        end

        @testset "edge features" begin
            ein = 3
            l = GATv2Conv((in_channel, ein) => out_channel, add_self_loops=false)
            g = GNNGraph(g1, edata=rand(T, ein, g1.num_edges))
            test_layer(l, g, rtol=RTOL_LOW,
                exclude_grad_fields = [:negative_slope],
                outsize=(out_channel, g.num_nodes))
        end

        @testset "num params" begin
            l = GATv2Conv(2 => 3, add_self_loops=false)
            @test length(Flux.params(l)) == 5
            l = GATv2Conv((2,4) => 3, add_self_loops=false)
            @test length(Flux.params(l)) == 6
            l = GATv2Conv((2,4) => 3, add_self_loops=false, bias=false)
            @test length(Flux.params(l)) == 4
        end

        @testset "edge features" begin
            ein = 3
            l = GATv2Conv((in_channel, ein) => out_channel, add_self_loops=false)
            g = GNNGraph(g1, edata=rand(T, ein, g1.num_edges))
            test_layer(l, g, rtol=1e-3,
                exclude_grad_fields = [:negative_slope],
                outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "GatedGraphConv" begin
        num_layers = 3
        l = GatedGraphConv(out_channel, num_layers)
        @test size(l.weight) == (out_channel, out_channel, num_layers)

        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "EdgeConv" begin
        l = EdgeConv(Dense(2*in_channel, out_channel), aggr=+)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "GINConv" begin
        nn = Dense(in_channel, out_channel)

        l = GINConv(nn, 0.01f0, aggr=mean)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        @test !in(:eps, Flux.trainable(l))
    end

    @testset "NNConv" begin
        edim = 10
        nn = Dense(edim, out_channel * in_channel)

        l = NNConv(in_channel => out_channel, nn, tanh, bias=true, aggr=+)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(T, edim, g.num_edges))
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "SAGEConv" begin
        l = SAGEConv(in_channel => out_channel)
        @test l.aggr == mean

        l = SAGEConv(in_channel => out_channel, tanh, bias=false, aggr=+)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "ResGatedGraphConv" begin
        l = ResGatedGraphConv(in_channel => out_channel, tanh, bias=true)
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end
    end

    @testset "CGConv" begin
        edim = 10
        l = CGConv((in_channel, edim) => out_channel, tanh, residual=false, bias=true)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(T, edim, g.num_edges))
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(out_channel, g.num_nodes))
        end

        # no edge features
        l1 = CGConv(in_channel => out_channel, tanh, residual=false, bias=true)
        @test l1(g1, g1.ndata.x) == l1(g1).ndata.x
        @test l1(g1, g1.ndata.x, nothing) == l1(g1).ndata.x
    end

    @testset "AGNNConv" begin
        l = AGNNConv()
        @test l.β == [1f0]
        for g in test_graphs
            test_layer(l, g, rtol=RTOL_HIGH, outsize=(in_channel, g.num_nodes))
        end
    end

    @testset "MEGNetConv" begin
        l = MEGNetConv(in_channel => out_channel, aggr=+)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(T, in_channel, g.num_edges))
            test_layer(l, g, rtol=RTOL_LOW,
                outtype=:node_edge,
                outsize=((out_channel, g.num_nodes), (out_channel, g.num_edges)))
        end
    end

    @testset "GMMConv" begin
        ein_channel = 10
        K = 5
        l = GMMConv((in_channel, ein_channel )=> out_channel, K=K)
        for g in test_graphs
            g = GNNGraph(g, edata=rand(Float32, ein_channel, g.num_edges))
            test_layer(l, g, rtol=RTOL_HIGH, outsize = (out_channel, g.num_nodes))
        end
    end
end


text:
using Images

# struct to store the Lab color values
# and center position of each cluster
mutable struct Cluster
    l
    a
    b
    y
    x
end

function slic(img, K, M, iterations=10, connectivity=false)
    img_lab = Lab.(img)
    size_tuple = size_spatial(img)
    image_height = size_tuple[1]
    image_width = size_tuple[2]
    S = round(Int, (sqrt((image_height * image_width) / K)))
    clusters = Cluster[] # The properties of each cluster
    labels = fill(-1, image_height, image_width) # Label of each pixel
    distance = fill(Inf, image_height, image_width) # Distance matrix of each pixel to belonging cluster
    pixel_count = Integer[] # Pixel counts of each cluster

    # Initialize each cluster and its fields
    for x = div(S, 2):S:image_width
        for y = div(S, 2):S:image_height
            push!(clusters, Cluster(img_lab[y, x].l,
                                   img_lab[y, x].a,
                                   img_lab[y, x].b,
                                   y,
                                   x))
            push!(pixel_count, 0)
        end
    end

    # Move the center of each cluster to the local lowgest gradient position
    function get_gradient(y, x)
        if x + 1 > image_width x = image_width - 2 end
        if y + 1 > image_height y = image_height - 2 end

        return img_lab[y + 1, x + 1].l - img_lab[y, x].l +
               img_lab[y + 1, x + 1].a - img_lab[y, x].a +
               img_lab[y + 1, x + 1].b - img_lab[y, x].b
    end
    for i = 1:length(clusters)
        # Get current gradient of this center
        current_gradient = get_gradient(clusters[i].y, clusters[i].x)

        for dh = -1:1
            for dw = -1:1
                _y = clusters[i].y + dh
                _x = clusters[i].x + dw
                new_gradient = get_gradient(_y, _x)
                if new_gradient < current_gradient
                    clusters[i].l = img_lab[_y, _x].l
                    clusters[i].a = img_lab[_y, _x].a
                    clusters[i].b = img_lab[_y, _x].b
                    clusters[i].y = _y
                    clusters[i].x = _x

                    current_gradient = new_gradient
                end
            end
        end
    end

    # SLIC superpixle calculation
    function cluster_pixels()
        for i = 1:length(clusters)
            for x = (clusters[i].x - 2 * S):(clusters[i].x + 2 * S)
                if x <= 0 || x > image_width continue end

                for y = (clusters[i].y - 2 * S):(clusters[i].y + 2 * S)
                    if y <= 0 || y > image_height continue end

                    L = img_lab[y, x].l
                    A = img_lab[y, x].a
                    B = img_lab[y, x].b
                    Dc = sqrt((L - clusters[i].l)^2 +
                              (A - clusters[i].a)^2 +
                              (B - clusters[i].b)^2)
                    Ds = sqrt((y - clusters[i].y)^2 +
                              (x - clusters[i].x)^2)
                    D = sqrt((Dc / M)^2 + (Ds / S)^2)

                    if D < distance[y, x]
                        distance[y, x] = D
                        labels[y, x] = i
                    end
                end
            end
        end
    end
    function update_cluster_position()
        # Clear the position value and pixel counts of each cluster
        for i = 1:length(clusters)
           clusters[i].y = clusters[i].x = pixel_count[i] = 0
        end

        # Compute the new position of new cluster center
        for x in 1:image_width
            for y in 1:image_height
                label_index = labels[y, x]
                if label_index == -1 continue end

                clusters[label_index].y += y
                clusters[label_index].x += x
                pixel_count[label_index] += 1
            end
        end

        # Assign the new position to each cluster
        for i = 1:length(clusters)
            new_y = div(clusters[i].y, pixel_count[i])
            new_x = div(clusters[i].x, pixel_count[i])
            clusters[i].l = img_lab[new_y, new_x].l
            clusters[i].a = img_lab[new_y, new_x].a
            clusters[i].b = img_lab[new_y, new_x].b
            clusters[i].y = new_y
            clusters[i].x = new_x
        end
    end
    @time for i = 1:iterations
        println("SLIC iteration $(i) ...")
        cluster_pixels()
        update_cluster_position()
    end

    # Enforce connectivity
    # Reference: https://github.com/scikit-image/scikit-image/blob/7e4840bd9439d1dfb6beaf549998452c99f97fdd/skimage/segmentation/_slic.pyx#L240-L348
    function enforce_connectivity(labels, min_size, max_size, start_label=1)
        depth = 1 # reserved for supervoxel used, now set to 1
        height = size(labels, 1)
        width = size(labels, 2)

        dx = [1, -1, 0, 0]
        dy = [0, 0, 1, -1]
        #dz = [] # reserved for supervoxel used

        # indicates that the label of this pixel has not been assigned
        mask_label = start_label - 1
        labels_final = fill(mask_label, height, width)

        current_new_label = start_label
        #label = start_label

        # used for BFS
        current_segment_size = 1
        bfs_visited = 0

        # store neighboring pixels
        # now set the dimension to 2 because we are using superpixel
        coord_list = fill(0, max_size, 2)

        for x = 1:width
            for y = 1:height
                if labels[y, x] == mask_label continue end
                if labels_final[y, x] > mask_label continue end

                adjacent = 0
                label = labels[y, x]
                labels_final[y, x] = current_new_label
                current_segment_size = 1
                bfs_visited = 0
                coord_list[bfs_visited + 1, 1] = y
                coord_list[bfs_visited + 1, 2] = x

                # Preform BFS to find the size of superpixel with
                # same lable number
                while bfs_visited < current_segment_size <= max_size
                    for i = 1:4
                        yy = coord_list[bfs_visited + 1, 1] + dy[i]
                        xx = coord_list[bfs_visited + 1, 2] + dx[i]

                        if 1 <= yy <= height &&  1 <= xx <= width
                            if labels[yy, xx] == label && labels_final[yy, xx] == mask_label
                                labels_final[yy, xx] = current_new_label
                                coord_list[current_segment_size + 1, 1] = yy # <-- index problem in the future
                                coord_list[current_segment_size + 1, 2] = xx # <-- index problem in the future
                                current_segment_size += 1

                                if current_segment_size > max_size break end
                            elseif labels_final[yy, xx] > mask_label &&
                                   labels_final[yy, xx] != current_new_label
                                adjacent = labels_final[yy, xx]
                            end
                        end
                    end
                    bfs_visited += 1
                end

                #for i = 1:current_segment_size println(coord_list[i, 1]) end
                #println(max_size, " ", current_segment_size)

                # merge the superpixel to its neighbor if it is too small
                if current_segment_size < min_size
                    @inbouds @simd for i = 1:current_segment_size
                        #println(i, " ", coord_list[i, 1])
                        labels_final[coord_list[i, 1],
                                     coord_list[i, 2]] = adjacent
                    end
                else
                    current_new_label += 1
                end
            end
        end

        return labels_final
    end
    if connectivity
        segment_size = image_height * image_width / size(clusters, 1)
        labels = enforce_connectivity(labels, round(Int, 0.5 * segment_size), round(Int, 3.0 * segment_size))
    end

    print(size(labels))
    print(size(clusters))
    print(labels[1, 1])

    # Create output image
    # The color of each cluster is as same as its center
    # except the center
    out_image = copy(img_lab)
    for x = 1:image_width
        for y = 1:image_height
            out_image[y, x, :] .= Lab(clusters[labels[y, x]].l,
                                   clusters[labels[y, x]].a,
                                   clusters[labels[y, x]].b)
        end
    end
    out_image = RGB.(out_image)

    # Return processed result
    return out_image
end


text:
"""

"""
module QuantumSpatialSearch

    using SparseArrays
    using LinearAlgebra
    using Distributed

    export spatialsearch, generateRTN, superposition_state

    const jumpProb = 0.05;   #jump probability in a timestep

    include("graphs.jl")
    include("spatialsearch.jl")
    include("utilities.jl")
    include("RTN.jl")

end # module

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
struct DescriptorSetLayoutCache
    device::Device
    dset_layouts::Dictionary{DescriptorSetLayoutCreateInfo,DescriptorSetLayout}
end

device(cache::DescriptorSetLayoutCache) = cache.device

DescriptorSetLayoutCache(device) = DescriptorSetLayoutCache(device, Dictionary{DescriptorSetLayoutCreateInfo,DescriptorSetLayout}())

@forward DescriptorSetLayoutCache.dset_layouts Base.haskey, Base.getindex, Base.insert!

function Base.get!(cache::DescriptorSetLayoutCache, info::DescriptorSetLayoutCreateInfo)
    if haskey(cache, info)
        cache[info]
    else
        dset_layout = unwrap(create_descriptor_set_layout(device(cache), info))
        insert!(cache, info, dset_layout)
        dset_layout
    end
end

function allocate_descriptor_sets!(da::DescriptorAllocator, cache::DescriptorSetLayoutCache, shaders)
    allocate_descriptor_sets!(da, create_descriptor_set_layouts!(cache, shaders))
end

function create_descriptor_set_layouts!(cache::DescriptorSetLayoutCache, shaders)
    map(collect_bindings(shaders)) do layout_binding
        info = DescriptorSetLayoutCreateInfo(layout_binding)
        Created(get!(cache, info), info)
    end
end

function create_descriptor_set_layouts(device, shaders)
    map(collect_bindings(shaders)) do layout_binding
        info = DescriptorSetLayoutCreateInfo(layout_binding)
        Created(unwrap(create_descriptor_set_layout(device, info)), info)
    end
end

function collect_bindings(shaders)
    binding_sets = Dictionary{Int,Vector{DescriptorSetLayoutBinding}}()
    for shader ∈ shaders
        for info ∈ shader.descriptor_infos
            push!(get!(binding_sets, info.index, DescriptorSetLayoutBinding[]), DescriptorSetLayoutBinding(info.binding, info.type, shader.source.stage; descriptor_count=1))
        end
    end
    if !all(collect(keys(binding_sets)) .== 0:length(binding_sets) - 1)
        error("Invalid layout description (non-contiguous binding sets from 0) in $binding_sets.")
    end
    collect(values(binding_sets))
end


text:
module ModuleLogging

export ModuleLogger

# using Lazy: @forward
using MacroTools
import Logging: AbstractLogger, SimpleLogger, Info
import Logging: shouldlog, handle_message, min_enabled_level, catch_exceptions

# SimpleLogger wrapper
struct ModuleLogger{T <:AbstractLogger} <: AbstractLogger
    logger::T
    modules::Set{Module}
end

# Constructors

ModuleLogger(stream::IO=stderr, level=Info, modules::Module...; kwargs...) =
    ModuleLogger(SimpleLogger(stream, level; kwargs...), Set(modules))

ModuleLogger(logger::T, modules::Module...) where T <: AbstractLogger =
    ModuleLogger(logger, Set(modules))

# Should log only if the module is being monitored

shouldlog(logger::ModuleLogger, level, _module, group, id) = begin
    _module ∈ logger.modules &&
        shouldlog(logger.logger, level, _module, group, id)
end

# forward macro that supports keyword arguments
# (upstream PR - https://github.com/MikeInnes/Lazy.jl/pull/112)
macro forward(ex, fs)
    @capture(ex, T_.field_) || error("Syntax: @forward T.x f, g, h")
    T = esc(T)
    fs = isexpr(fs, :tuple) ? map(esc, fs.args) : [esc(fs)]
    :($([:($f(x::$T, args...; kwargs...) = (Base.@_inline_meta; $f(x.$field, args...; kwargs...)))
         for f in fs]...);
      nothing)
end

# Delegate to underlying logger
@forward ModuleLogger.logger handle_message, min_enabled_level, catch_exceptions

end # module


text:
install_tips = """
OpenGL/GLFW wasn't loaded correctly or couldn't be initialized.
This likely means, you don't have an OpenGL capable Graphic Card,
or you don't have an OpenGL 3.3 capable video driver installed.
If you're on a headless server, you still need to install x-server and
proper graphics drivers.
If you don't have a GPU, there is also a Cairo software backend
for Makie which you can use:
https://github.com/JuliaPlots/CairoMakie.jl.
Please check the below error and open an issue at:
https://github.com/JuliaPlots/GLMakie.jl.
After you fixed your OpenGL install, please run `]build GLMakie`!
Otherwise, GLMakie would still load, but would be disabled as a
default backend for Makie.
"""
try
    using GLFW
catch e
    open("deps.jl", "w") do io
        println(io, "const WORKING_OPENGL = false")
    end
    # it would be nice to check if this is a GLFW error, but if GLFW doesn't actually load
    # we can't easily use GLFW.GLFWError. Well, GLFW error is the most likely, and
    # we will print the error, to inform the user what happens, so I think this should be fine!
    println(stderr,"Load error for GLFW")
    error(install_tips)
end

try
    using ModernGL
    # Create a windowed mode window and its OpenGL context
    window = GLFW.Window(resolution = (10, 10), major = 3, minor = 3, visible = false, focus = false)
    glversion = unsafe_string(glGetString(GL_VERSION))
    m = match(r"(\d+)\.(\d+)(.\d+)?\s", glversion)
    # I don't really trust that all vendors have a version that matches
    # the above regex, so let's make no match non fatal!
    if m === nothing
        @warn("Unknown OpenGL version format: $glversion. You need to verify if it's above OpenGL 3.3 yourself!")
    else
        v = VersionNumber(parse(Int, m[1]), parse(Int, m[2]))
        if !(v >= v"3.3")
            open("deps.jl", "w") do io
                println(io, "const WORKING_OPENGL = false")
            end
            println(stderr, "Your OpenGL version is too low! Update your driver or GPU! Version found: $v, version required: 3.3")
            error(install_tips)
        end
    end

    open("deps.jl", "w") do io
        println(io, "const WORKING_OPENGL = true")
    end
catch e
    open("deps.jl", "w") do io
        println(io, "const WORKING_OPENGL = false")
    end
    # it would be nice to check if this is a GLFW error, but if GLFW doesn't actually load
    # we can't easily use GLFW.GLFWError. Well, GLFW error is the most likely, and
    # we will print the error, to inform the user what happens, so I think this should be fine!
    println(stderr, "Initialization error of GLFW")
    error(install_tips)
end


text:
"""
    CosmoParticles._applyind(a, ind::Union{AbstractVector,Colon})

Apply indices or mask to a `Number`, `Vector`, or `Matrix`.

Other possible types: `Nothing`

The following indexing is applied:
- `a::Number`: `a` is returned directly.
- `a::Vector`: `a[ind]` is returned.
- `a::Matrix`: `a[:, ind]` is returned.

This is not exported.
"""
_applyind(a::Union{Number,Nothing}, _::Union{AbstractVector,Colon}) = a
_applyind(a::AbstractVector, ind::Union{AbstractVector,Colon}) = a[ind]
_applyind(a::AbstractMatrix, ind::Union{AbstractVector,Colon}) = a[:, ind]
_applyind(a::ApplyVector{<:Any,F}, mask::AbstractVector{Bool}) where {F<:Union{typeof(hcat),typeof(vcat)}} =
    a[mask]
_applyind(a::ApplyMatrix{<:Any,F}, mask::AbstractVector{Bool}) where {F<:Union{typeof(hcat),typeof(vcat)}} =
    a[:, mask]

function _applyind(a::ApplyVector{<:Any,F}, ind::AbstractVector) where {F<:Union{typeof(hcat),typeof(vcat)}}
    anew = similar(a, length(ind))
    @inbounds for (i, indi) in enumerate(ind)
        anew[i] = a[indi]
    end
    return anew
end

function _applyind(a::ApplyMatrix{<:Any,F}, ind::AbstractVector) where {F<:Union{typeof(hcat),typeof(vcat)}}
    n = size(a, 1)
    anew = similar(a, n, length(ind))
    @inbounds for (i, indi) in enumerate(ind)
        for j in 1:n
            anew[j, i] = a[j, indi]
        end
    end
    return anew
end

"""
    CosmoParticles.applyind!(p::AbstractParticles, ind::Union{AbstractVector,Colon})

In-place application of indices or mask to all particle properties.

This is not exported.
"""
function applyind!(p::AbstractParticles, ind::Union{AbstractVector,Colon})
    Threads.@threads for key in keys(p) |> collect # collect used for compatibility with threads
        p[key] = _applyind(p[key], ind)
    end

    return p
end

"""
    CosmoParticles.applyind(p::AbstractParticles, ind::Union{AbstractVector,Colon}; affect=keys(p))

Create new particles with the given indices or mask applied to all particle properties.

This can also be called by the simple syntax `p[ind]`.
If the keyword argument `affect` is a vector of `Symbol`s, only those properties are indexed into and
added to the newly created particles object.

This is not exported.
"""
function applyind(p::AbstractParticles, ind::Union{AbstractVector,Colon}; affect=keys(p))
    affected_keys = intersect(affect, keys(p)) |> collect # collect used for compatibility with threads

    # preallocate particles (to avoid threading issues)
    pnew = empty(p)
    for key in affected_keys
        pnew[key] = p[key]
    end

    Threads.@threads for key in affected_keys # collect used for compatibility with threads
        pnew[key] = _applyind(p[key], ind)
    end

    return pnew
end

"""
    CosmoParticles._removeind(a, ind::AbstractVector)

Remove indices of a `Number`, `Vector`, or `Matrix`.

The vector `ind` has to contain sorted and unique indices: `sort!(unique!(ind))`.
Additionally, all elements in `ind` have to be valid indices in `a`.

Other possible types: `Nothing`

The following indexing is applied:
- `a::Number`: `a` is returned directly.
- `a::Vector`: `a[Not(ind)]` is returned.
- `a::Matrix`: `a[:, Not(ind)]` is returned.

This is not exported.
"""
_removeind(a::Union{Number,Nothing}, _::AbstractVector) = a

function _removeind(a::AbstractVector, ind::AbstractVector)
    na = length(a)
    nind = length(ind)
    nanew = na - nind
    anew = similar(a, nanew)

    # catch empty array (otherwise throws BoundsError in loop)
    na == nind && return anew

    ia = 1
    iind = 1
    ianew = 1
    @inbounds while true
        while iind > nind || ia < ind[iind]
            anew[ianew] = a[ia]
            ia += 1
            ianew += 1
            ianew > nanew && return anew
        end
        ia += 1
        iind += 1
    end
end

function _removeind(a::AbstractMatrix, ind::AbstractVector)
    n = size(a, 1)
    na = size(a, 2)
    nind = length(ind)
    nanew = na - nind
    anew = similar(a, n, nanew)

    # catch empty array (otherwise throws BoundsError in loop)
    na == nind && return anew

    ia = 1
    iind = 1
    ianew = 1
    @inbounds while true
        while iind > nind || ia < ind[iind]
            for j in 1:n
                anew[j, ianew] = a[j, ia]
            end
            ia += 1
            ianew += 1
            ianew > nanew && return anew
        end
        ia += 1
        iind += 1
    end
end

"""
    CosmoParticles.removeind!(p::AbstractParticles, ind::AbstractVector)

In-place removal of indices to all particle properties.

This is not exported.
"""
function removeind!(p::AbstractParticles, ind::AbstractVector)
    Threads.@threads for key in keys(p) |> collect # collect used for compatibility with threads
        p[key] = _removeind(p[key], ind)
    end

    return p
end

"""
    CosmoParticles.removeind(p::AbstractParticles, ind::AbstractVector; affect=keys(p))

Create new particles without the given indices or mask applied to all particle properties.

This can also be called by the simple syntax `p[ind]`.
If the keyword argument `affect` is a vector of `Symbol`s, only those properties are indexed into and
added to the newly created particles object.

This is not exported.
"""
function removeind(p::AbstractParticles, ind::AbstractVector; affect=keys(p))
    affected_keys = intersect(affect, keys(p)) |> collect # collect used for compatibility with threads

    # preallocate particles (to avoid threading issues)
    pnew = empty(p)
    for key in affected_keys
        pnew[key] = p[key]
    end

    Threads.@threads for key in affected_keys
        pnew[key] = _removeind(p[key], ind)
    end

    return pnew
end

"""
    CosmoParticles.findall_in(a::AbstractVector, set)

Return all indices of `a` that are in `set`.

If both `a` and `set` are sorted `AbstractVector`s, then the optimized [`findall_in_sorted`](@ref) is called.
Otherwise, a `Set` is constructed from the `Vector` to perform the checks with `in`.

This is not exported.
"""
function findall_in(a::AbstractVector, set::AbstractVector)
    if issorted(a) && issorted(set)
        return findall_in_sorted(a, set)
    else
        return findall_in(a, Set(set))
    end
end

findall_in(a::AbstractVector, set::AbstractSet) = findall(in.(a, (set,)))

"""
    CosmoParticles.findall_in_sorted(a::AbstractVector, set::AbstractVector)

Return all indices of `a` that are in `set`, where both `a` and `set` are assumed to be sorted.

This uses an optimized algorithm that is faster than creating a `Set` from `set` and performing checks with `in`.

This is not exported.
"""
function findall_in_sorted(a::AbstractVector, set::AbstractVector)
    if isempty(a) || isempty(set)
        return Int64[]
    end

    na = length(a)
    nset = length(set)

    ind_all = Vector{Int64}(undef, min(na, nset))

    # NOTE: this is an alternative algorithm, which performs slightly worse than the below
    # iind = 0
    # iset = 1
    # for ia in 1:na
    # while a[ia] > set[iset] && iset < nset
    # iset += 1
    # end
    # if a[ia] == set[iset]
    # iind += 1
    # iset += 1
    # ind_all[iind] = ia
    # end
    # end

    iind = 0
    iset = 1
    ia = 1

    # for a large number of values, search for the first indices in a faster way
    if na > 1_000_000
        if first(a) < first(set)
            ia = searchsortedfirst(a, first(set))
        else
            iset = searchsortedfirst(set, first(a))
            ia = searchsortedfirst(a, set[iset])
        end
    end

    if ia > na || iset > nset
        return Int64[]
    end

    while true
        if a[ia] == set[iset]
            iind += 1
            ind_all[iind] = ia
            iset += 1
            ia += 1
            ((ia > na) || (iset > nset)) && break
        else
            if a[ia] < set[iset]
                ia += 1
                ia > na && break
            else
                iset += 1
                iset > nset && break
            end
        end
    end

    return resize!(ind_all, iind)
end

"""
    CosmoParticles.product_preserve_type(arr::AbstractArray{T}, b::Real) where {T}

Multiply an array with a scalar while preserving the element type of the array.

For unitful arrays, the scalar factor is converted to the number type of the quantity before multiplying.

This is not exported.
"""
product_preserve_type(arr::AbstractArray{T}, b::Real) where {T<:Real} = arr .* convert(T, b)
product_preserve_type(arr::AbstractArray{<:Quantity{T}}, b::Real) where {T<:Real} = arr .* convert(T, b)

"""
    CosmoParticles product_preserve_type!(arr::AbstractArray{T}, b::Real) where {T}

Multiply an array with a scalar in-place while preserving the element type of the array.

By converting the scalar factor to the array element type, this can be more performant than a normal
broadcasted product.
For unitful arrays, the scalar factor is converted to the number type of the quantity before multiplying.

This is not exported.
"""
function product_preserve_type!(arr::AbstractArray{T}, b::Real) where {T<:Real}
    arr .*= convert(T, b)
end
function product_preserve_type!(arr::AbstractArray{<:Quantity{T}}, b::Real) where {T<:Real}
    arr .*= convert(T, b)
end
function product_preserve_type!(
    arr::ApplyArray{<:Any,N,F},
    b::Real,
) where {N,F<:Union{typeof(hcat),typeof(vcat)}}
    product_preserve_type!.(arr.args, b)
    return arr
end

"""
    CosmoParticles.ustrip_lazy([unit,] a::Number)
    CosmoParticles.ustrip_lazy([unit,] a::AbstractArray)
    CosmoParticles.ustrip_lazy!(unit, a::AbstractArray)

Strips off the units of unitful arrays in a performant way while keeping the number type conserved (not always the case for `ustrip`).

Like [`ustrip`](https://painterqubits.github.io/Unitful.jl/stable/manipulations/#Unitful.ustrip) for normal
unitful arrays, but without reallocating non-unitful or lazy arrays.

Note that after calling this method with a different conversion unit, the original array still has the old units,
but the numerical values correspond to the new units. This method should always be called as
`a = uconvert_lazy!(u, a)`.

This is not exported.
"""
ustrip_lazy(a::Unitful.AbstractQuantity) = ustrip(a)
ustrip_lazy(a::Number) = a
ustrip_lazy(a::AbstractArray) = a
ustrip_lazy(a::Array{<:Quantity}) = ustrip(a)
ustrip_lazy(a::Diagonal{<:Quantity}) = ustrip(a)
ustrip_lazy(a::Bidiagonal{<:Quantity}) = ustrip(a)
ustrip_lazy(a::Tridiagonal{<:Quantity}) = ustrip(a)
ustrip_lazy(a::SymTridiagonal{<:Quantity}) = ustrip(a)

ustrip_lazy(a::AbstractArray{<:Quantity{T}}) where {T} = reinterpret(T, a) #ustrip.(a) # alternative implementation
ustrip_lazy(a::ApplyArray{<:Quantity{T}}) where {T} = reinterpret(T, a)
ustrip_lazy(a::Fill{<:Quantity{T}}) where {T} = reinterpret(T, a)

ustrip_lazy(u::Unitful.Units, a) = ustrip.(u, a)
ustrip_lazy(u::Unitful.Units, a::AbstractVector{<:Quantity{T}}) where {T<:AbstractFloat} = ustrip_lazy.(u, a)
ustrip_lazy(u::Unitful.Units, a::Quantity{T}) where {T<:AbstractFloat} = ustrip_lazy(uconvert_lazy(u, a))
ustrip_lazy!(u::Unitful.Units, a) = ustrip_lazy(uconvert_lazy!(u, a))

uconvert_lazy(u::Unitful.Units, a::AbstractArray{<:Quantity}) = uconvert_lazy!(u, copy(a))

"""
    CosmoParticles.uconvert_lazy(u::Unitful.Units, a::AbstractArray{<:Quantity})

Converts the unitful array to the given units while preserving the number type passed (not always the case for `uconvert`).

This is not exported.
"""
function uconvert_lazy(u::Unitful.Units, a::Quantity{T,D,U}) where {T,D,U}
    return convert(T, ustrip(u, a)) * u
end

"""
    CosmoParticles.uconvert_lazy!(u::Unitful.Units, a::AbstractArray{<:Quantity{T,D,U}}) where {T,D,U}

Converts the unitful array to the specified units in-place.

Note that after calling this method the original array still has the old units, but the numerical values
correspond to the new units. This method should always be called as `a = uconvert_lazy!(u, a)`.

This is not exported.
"""
function uconvert_lazy!(u::Unitful.Units, a::AbstractArray{<:Quantity{T,D,U}}) where {T,D,U}
    Q = Quantity{T,dimension(u),typeof(u)}
    if typeof(u) == U
        return reinterpret(Q, a)
    elseif (u isa Unitful.AffineUnits) || (U <: Unitful.AffineUnits)
        a_unitless = ustrip_lazy(a)
        a_unitless .= ustrip_lazy(Unitful.uconvert_affine.(u, a))
        return reinterpret(Q, a_unitless)
    else
        product_preserve_type!(a, Unitful.convfact(u, U()))
        return reinterpret(Q, a)
    end
end


text:
using SpinMC
using Test
using MPI

#init MPI
MPI.Init()
comm = MPI.COMM_WORLD
commRank = MPI.Comm_rank(comm)
commSize = MPI.Comm_size(comm)
@test commSize == 2

#run tests
if commRank == 0
    rec = SpinMC.MPISendrecvFloat(42.0, 1, comm)
    @test rec == 43.0
elseif commRank == 1
    rec = SpinMC.MPISendrecvFloat(43.0, 0, comm)
    @test rec == 42.0
end

if commRank == 0
    SpinMC.MPISendBool(true, 1, comm)
    SpinMC.MPISendBool(false, 1, comm)
elseif commRank == 1
    rec = SpinMC.MPIRecvBool(0, comm)
    @test rec == true
    rec = SpinMC.MPIRecvBool(0, comm)
    @test rec == false
end

rec = SpinMC.MPIBcastBool(true, 0, comm)
@test rec == true
rec = SpinMC.MPIBcastBool(false, 0, comm)
@test rec == false

#finalize MPI
MPI.Finalize()
@test MPI.Finalized()


text:
function boundary_traction_component_operator(
    basis,
    quad,
    component,
    normal,
    stiffness,
    dim,
    facedetjac,
    jac,
    vectosymmconverter,
)

    numqp = length(quad)
    normals = repeat(normal, inner = (1, numqp))
    components = repeat(component, inner = (1, numqp))
    scalearea = repeat([facedetjac], numqp)
    return -1.0 * surface_traction_component_operator(
        basis,
        quad,
        quad,
        components,
        normals,
        stiffness,
        dim,
        scalearea,
        jac,
        vectosymmconverter,
    )
end

function assemble_boundary_traction_component_operator!(
    sysmatrix,
    basis,
    facequads,
    stiffness,
    mesh,
    onboundary,
    component,
)

    facemidpoints = reference_face_midpoints()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    jac = jacobian(mesh)
    dim = dimension(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    normals = reference_face_normals()
    uniformquads = uniform_face_quadratures(facequads)

    faceids = 1:nfaces

    uniformtop1 = [
        vec(
            boundary_traction_component_operator(
                basis,
                q,
                component,
                n,
                stiffness[+1],
                dim,
                d,
                jac,
                vectosymmconverter,
            ),
        ) for (q, n, d) in zip(uniformquads, normals, facedetjac)
    ]
    uniformtop2 = [
        vec(
            boundary_traction_component_operator(
                basis,
                q,
                component,
                n,
                stiffness[-1],
                dim,
                d,
                jac,
                vectosymmconverter,
            ),
        ) for (q, n, d) in zip(uniformquads, normals, facedetjac)
    ]

    uniformtop = [uniformtop1, uniformtop2]

    ncells = number_of_cells(mesh)

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )

                if cellsign == +1 || cellsign == -1
                    row = cell_sign_to_row(cellsign)
                    nodeids = nodal_connectivity(mesh, cellsign, cellid)
                    assemble_cell_matrix!(
                        sysmatrix,
                        nodeids,
                        dim,
                        uniformtop[row][faceid],
                    )
                else
                    assemble_cut_cell_boundary_traction_component_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        component,
                        normals,
                        stiffness,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        jac,
                        vectosymmconverter,
                    )
                    assemble_cut_cell_boundary_traction_component_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        component,
                        normals,
                        stiffness,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        jac,
                        vectosymmconverter,
                    )
                end

            end
        end
    end
end

function assemble_cut_cell_boundary_traction_component_operator!(
    sysmatrix,
    basis,
    facequads,
    component,
    normals,
    stiffness,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    jac,
    vectosymmconverter,
)

    dim = dimension(mesh)
    operator = vec(
        boundary_traction_component_operator(
            basis,
            facequads[cellsign, faceid, cellid],
            component,
            normals[faceid],
            stiffness[cellsign],
            dim,
            facedetjac[faceid],
            jac,
            vectosymmconverter,
        ),
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_matrix!(sysmatrix, nodeids, dim, operator)
end

function boundary_mass_component_operator(basis, quad, component, dim, scale)
    numqp = length(quad)
    scalearea = repeat([scale], numqp)
    components = repeat(component, inner = (1, numqp))
    return component_mass_matrix(basis, quad, quad, components, dim, scalearea)
end

function assemble_boundary_mass_component_operator!(
    sysmatrix,
    basis,
    facequads,
    mesh,
    onboundary,
    component,
    penalty,
)

    uniformquads = uniform_face_quadratures(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    nfaces = number_of_faces_per_cell(facequads)
    dim = dimension(mesh)
    facemidpoints = reference_face_midpoints()

    faceids = 1:nfaces

    uniformmassops = [
        vec(
            boundary_mass_component_operator(
                basis,
                q,
                component,
                dim,
                penalty * d,
            ),
        ) for (q, d) in zip(uniformquads, facedetjac)
    ]

    ncells = number_of_cells(mesh)

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )

                if cellsign == +1 || cellsign == -1
                    nodeids = nodal_connectivity(mesh, cellsign, cellid)
                    assemble_cell_matrix!(
                        sysmatrix,
                        nodeids,
                        dim,
                        uniformmassops[faceid],
                    )
                else
                    assemble_cut_cell_boundary_mass_component_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        component,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                    assemble_cut_cell_boundary_mass_component_operator!(
                        sysmatrix,
                        basis,
                        facequads,
                        component,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
            end
        end
    end
end

function assemble_cut_cell_boundary_mass_component_operator!(
    sysmatrix,
    basis,
    facequads,
    component,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    penalty,
)

    dim = dimension(mesh)
    operator = vec(
        boundary_mass_component_operator(
            basis,
            facequads[cellsign, faceid, cellid],
            component,
            dim,
            penalty * facedetjac[faceid],
        ),
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_matrix!(sysmatrix, nodeids, dim, operator)
end

function assemble_boundary_displacement_component_rhs!(
    sysrhs,
    bcfunc,
    basis,
    facequads,
    mesh,
    onboundary,
    component,
    penalty,
)

    facemidpoints = reference_face_midpoints()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    ncells = number_of_cells(mesh)
    faceids = 1:nfaces

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        cellmap1 = cell_map(mesh, +1, cellid)
        cellmap2 = cell_map(mesh, -1, cellid)

        check_cellsign(cellsign)
        for faceid in faceids
            if cell_connectivity(mesh, faceid, cellid) == 0 && (
                onboundary(cellmap1(facemidpoints[faceid])) ||
                onboundary(cellmap2(facemidpoints[faceid]))
            )
                if cellsign == +1 || cellsign == 0
                    assemble_face_boundary_displacement_component_rhs!(
                        sysrhs,
                        bcfunc,
                        basis,
                        facequads,
                        component,
                        mesh,
                        +1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
                if cellsign == -1 || cellsign == 0
                    assemble_face_boundary_displacement_component_rhs!(
                        sysrhs,
                        bcfunc,
                        basis,
                        facequads,
                        component,
                        mesh,
                        -1,
                        faceid,
                        cellid,
                        facedetjac,
                        penalty,
                    )
                end
            end
        end
    end
end

function assemble_face_boundary_displacement_component_rhs!(
    sysrhs,
    bcfunc,
    basis,
    facequads,
    component,
    mesh,
    cellsign,
    faceid,
    cellid,
    facedetjac,
    penalty,
)

    dim = dimension(mesh)
    cellmap = cell_map(mesh, cellsign, cellid)
    rhsvals =
        penalty * component_linear_form(
            bcfunc,
            basis,
            facequads[cellsign, faceid, cellid],
            component,
            cellmap,
            dim,
            facedetjac[faceid],
        )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_rhs!(sysrhs, nodeids, dim, rhsvals)
end

function assemble_penalty_displacement_component_bc!(
    sysmatrix,
    sysrhs,
    bcfunc,
    basis,
    facequads,
    stiffness,
    mesh,
    onboundary,
    component,
    penalty,
)

    assemble_boundary_traction_component_operator!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        onboundary,
        component,
    )
    assemble_boundary_mass_component_operator!(
        sysmatrix,
        basis,
        facequads,
        mesh,
        onboundary,
        component,
        penalty,
    )
    assemble_boundary_displacement_component_rhs!(
        sysrhs,
        bcfunc,
        basis,
        facequads,
        mesh,
        onboundary,
        component,
        penalty,
    )
end


text:
"""
$(TYPEDEF)

Options for using SymbolicRegression.jl within the `solve` function.
Automatically creates [`Options`](https://astroautomata.com/SymbolicRegression.jl/stable/api/#Options) with the given specification.
Sorts the operators stored in `functions` into unary and binary operators on conversion.

# Fields
$(FIELDS)

"""
struct EQSearch <: AbstractSymbolicRegression
    "Operators used for symbolic regression"
    functions::AbstractVector{Function}
    "Additionally keyworded arguments passed to SymbolicRegression.Options"
    kwargs
end

function EQSearch(functions = Function[/, +, *, exp, cos]; kwargs...)
    return EQSearch(functions, kwargs)
end

function to_options(x::EQSearch)
    binops = Tuple([fi for fi in x.functions if is_binary(fi)])
    unaops = Tuple([fi for fi in x.functions if is_unary(fi)])
    fnames = fieldnames(SymbolicRegression.Options)
    ks = Dict(
        [k => v for (k,v) in x.kwargs if k ∈ fnames]
    )
    return SymbolicRegression.Options(;
        binary_operators = binops, unary_operators = unaops,
        ks...
    )
end

function DiffEqBase.solve(prob::AbstractDataDrivenProblem, alg::EQSearch;
    max_iter::Int = 10,
    weights = nothing,
    numprocs = nothing, procs = nothing,
    multithreading = false,
    runtests::Bool = true,
    eval_expression = false,
    kwargs...
    )

    opt = to_options(alg)

    Y = get_target(prob)

    # Inputs
    X̂, _, t, U = get_oop_args(prob)

    t = iszero(t) ? [] : t
    # Cat the inputs
    X = vcat([x for x in (X̂, U, permutedims(t)) if !isempty(x)]...)
    # Use the eqssearch of symbolic regression
    hof = SymbolicRegression.EquationSearch(X, Y, niterations = max_iter, weights = weights, options = opt,
            numprocs = numprocs, procs = procs, multithreading = multithreading,
            runtests = runtests)

    build_solution(prob, alg, hof; eval_expression = eval_expression)
end

function pareto_optimal_equations(hof::HallOfFame, prob, alg)
    return pareto_optimal_equations([hof], prob, alg)
end

function pareto_optimal_equations(hof::Vector{HallOfFame}, prob, alg)

    opts = DataDrivenDiffEq.to_options(alg)
    y = DataDrivenDiffEq.get_target(prob)
    x, _, t, c = DataDrivenDiffEq.get_oop_args(prob)
    X =  vcat([x for x in (x, c, permutedims(t)) if !isempty(x)]...)

    @parameters t
    @variables x[1:size(prob.X, 1)](t) u[1:size(prob.U,1)](t)
    x = collect(x)
    u = collect(u)
    x_ = Num[x;u;t]

    # Build a dict

    subs = Dict([SymbolicUtils.Sym{LiteralReal}(Symbol("x$(i)")) => x_[i] for i in 1:size(x_, 1)]...)

    eqs = map(1:size(hof, 1)) do i
        d = calculateParetoFrontier(X, y[i,:], hof[i], opts)
        isempty(d) && return Num(0)
        eq_ = node_to_symbolic(last(d).tree, opts)
        substitute(eq_, subs)
    end

    return Num.(eqs), x, u, t
end

function build_solution(prob::AbstractDataDrivenProblem, alg::EQSearch, hof; eval_expression = false)

    eqs, x, u, t = pareto_optimal_equations(hof, prob, alg)

    lhs, dt = assert_lhs(prob)

    # Build the lhs
    if (length(eqs) == length(x)) && (lhs != :direct)
            if lhs == :continuous
                d = Differential(t)
            elseif lhs == :discrete
                d = Difference(t, dt = dt)
            end
            eqs = [d(x[i]) ~ eq for (i,eq) in enumerate(eqs)]
    end

    res_ = Basis(
        eqs, x, iv = t, controls = u, eval_expression = eval_expression
    )

    X = get_target(prob)
    Y = res_(get_oop_args(prob)...)

    retcode = :converged

    return DataDrivenSolution(
        res_, [], retcode, alg, hof, prob, false
    )
end


text:
"""
Package: Forecast

    quakes(query = "count")

Return data about earthquakes with a magnitude higher or equal to six from 1638-06-11 to 2021-10-02. The data has been collected from https://earthquake.usgs.gov/

# Arguments
- `query`: Defaults to "count" and accepts the following values:
    - "full" : Full original dataset
    - "count" : Number of earthquakes by year

# Returns
Dataframe containing information based on the parameter `query`

# Examples
```julia-repl
julia> quakes()
71×2 DataFrame
 Row │ year        quakes
     │ Date        Int64
─────┼────────────────────
   1 │ 1950-01-01     138
   2 │ 1951-01-01     151
   3 │ 1952-01-01     181
   [...]
```
"""
function quakes(query::String = "count")

    @assert query in ["full","count"]

    path = art_path("quakes")

    qk_df = open(joinpath(path, "quakes.csv")) do file
        CSV.read(file,DataFrame,
                 dateformat = "yyyy-mm-ddTHH:MM:SS.sssZ",
                 types = Dict(:time => Date))
    end

    if query == "full"
        @info "Full dataset from 1638-06-11 to 2021-10-02"
        return qk_df
    end

    if query == "count"
        @info "Number of earthquakes per year"
        qk_df = DataFrame(reshape(Dates.year.(qk_df.time),:,1),:auto)

        db = SQLite.DB();
        Q = q -> DBInterface.execute(db, q);

        qk_df |> SQLite.load!(db, "qk_df", temp=false);

        query_count = """select x1 as year, count(*) as total_quakes
                                 from qk_df
                                 group by x1
                                 order by x1"""

        return Q(query_count)  |> DataFrame

    end

end


text:
include("src/JuliaTwitter.jl")
using JuliaTwitter

#Parses a string into a number if it can be represented as one
parse_if_num(s::AbstractString) = isa(parse(s), Number)  ? parse(s) : s
parse_if_num(s::Number) = s

#Evaluates a function with the given args as strings
function eval_str(fn::AbstractString, args...)
  #Parse the arguments
  parsed_args = map(parse_if_num, [args...])

  #Evaluate the parsed expression
  eval(Expr(:call, parse(fn), parsed_args...))
end

VALID_ARG_NAMES = Dict(
  #Maximum number of tweets to search for
  "-tweets" => "maxtweets",
  "-n" => "maxtweets",
  #Query
  "-query" => "query",
  "-q" => "query",
  #Max time to continue searching
  "-time" => "time",
  "-t" => "time",
  #Databse Handle
  "-dbh" => "dbh",
  #Alert Every
  "-alert-every" => "alertevery",
  "-a" => "alertevery",
  #Max consecutive errors
  "-max-errors" => "maxerrors",
  "-e" => "maxerrors"

)

function is_valid_arg(s::AbstractString)
  haskey(VALID_ARG_NAMES, s)
end

proper_argname(s::AbstractString) = VALID_ARG_NAMES[s]

function parse_args(args)#::Array{AbstractString, 1})
  dict = Dict{AbstractString, Union{Int, AbstractString}}(
    "maxtweets" => 100000,
    "time" => 60,
    "alertevery" => 20,
    "maxerrors" => 3
  )

  ix = 1
  while ix < length(args)
    try
      if is_valid_arg(args[ix])
        dict[proper_argname(args[ix])] = args[ix+1]
      else
        print_with_color(:red, "$(args[ix]) is not a valid argument\n")
        throw(error("Invalid argument"))
      end
    catch y
      println(y)
      print_with_color(:red, "Invalid argument near $(args[ix])\n")
      throw(error("Invalid argument"))
    end
    ix += 2
  end

  #Check for dbh
  if !haskey(dict, "dbh")
    println("[ERROR] Database not provided. Use -dbh "" to provide.")
    error("database not provided")
  end

  #Check for q
  if !haskey(dict, "q")
    println("[ERROR] Query not provided. Use -q "" to provide.")
    error("query not provided")
  end

  #Convert numbers
  dict["maxtweets"] = parse_if_num(dict["maxtweets"])
  dict["time"] = parse_if_num(dict["time"])
  dict["alertevery"] = parse_if_num(dict["alertevery"])
  dict["maxerrors"] = parse_if_num(dict["maxerrors"])

  #Check numbers
  typeof(dict["maxtweets"]) != Int && throw(ArgumentException(string("Cannot parse input for maxtweets: ", dict["maxtweets"])))
  typeof(dict["time"]) != Int && throw(ArgumentException(string("Cannot parse input for time: ", dict["time"])))
  typeof(dict["alertevery"]) != Int && throw(ArgumentException(string("Cannot parse input for alertevery: ", dict["alertevery"])))
  typeof(dict["maxerrors"]) != Int && throw(ArgumentException(string("Cannot parse input for maxerrors: ", dict["maxerrors"])))

  return dict
end

function create_dbh(in_str::AbstractString)
  in_str = string("JuliaTwitter.", in_str)
  eval(parse(in_str))
end

function main()
  println("Julia Loaded")
  #Parse and validate all arguments
  args = parse_args(ARGS)
  println("Using arguments:")
  println("\tMax number of tweets:\t", args["maxtweets"])
  println("\tMax time to search:\t", args["time"])
  println("\tQuery:\t\t\t", args["query"])
  println("\tDatabase:\t\t", args["dbh"])
  println("\tMax consecutive errors:\t", args["maxerrors"])
  println("\tAlert every N tweets:\t", args["alertevery"])

  dbh = 0
  try
    dbh = create_dbh(args["dbh"])
  catch y
    println("[ERROR] Could not initialize database: ", args["dbh"])
    println("        See error message below")
    throw(y)
  end
  println("Successfully initialized database: ", dbh)

  options = Dict("track" => args["query"])

  println("Starting stream...")
  stream_tweets(JuliaTwitter.auth, dbh, options;
        max_tweets = args["maxtweets"],
        max_time = args["time"],
        alert_every = args["alertevery"],
        max_err_attemps = args["maxerrors"])

  println("Complete.")
end

main()

#=
Load 10 trump tweets into a file

julia stream_tweets.jl -dbh TwTXTFile(\"path/to/txtfile.txt\") -maxtweets 10 -q Trump

Possible arguments
-max-tweets (-m) <number>
  the maximum number of tweets to return
-query (-q) <string>
  the query
-time (-t) <number>
  the max number of seconds to search
  -time-<<unit>> (-t-<<unit>>) <number>
  the max number of <<unit>> to search
  will be converted to seconds before entering into dict
  example:
    -t-hour 4

-db-handle (-d) <string>
=#

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm),  Float64, (Float64,), x), x)
gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm),  Float32, (Float32,), x), x)

"""
    gamma(x)

Compute the gamma function of `x`.
"""
gamma(x::Real) = gamma(float(x))

function lgamma_r(x::Float64)
    signp = Ref{Int32}()
    y = ccall((:lgamma_r,libm),  Float64, (Float64, Ptr{Int32}), x, signp)
    return y, signp[]
end
function lgamma_r(x::Float32)
    signp = Ref{Int32}()
    y = ccall((:lgammaf_r,libm),  Float32, (Float32, Ptr{Int32}), x, signp)
    return y, signp[]
end
lgamma_r(x::Real) = lgamma_r(float(x))
lgamma_r(x::Number) = lgamma(x), 1 # lgamma does not take abs for non-real x
"`lgamma_r(x)`: return L,s such that `gamma(x) = s * exp(L)`" lgamma_r

"""
    lfact(x)

Compute the logarithmic factorial of a nonnegative integer `x`.
Equivalent to [`lgamma`](@ref) of `x + 1`, but `lgamma` extends this function
to non-integer `x`.
"""
lfact(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : lgamma(x + oneunit(x))

"""
    lgamma(x)

Compute the logarithm of the absolute value of [`gamma`](@ref) for
[`Real`](@ref) `x`, while for [`Complex`](@ref) `x` compute the
principal branch cut of the logarithm of `gamma(x)` (defined for negative `real(x)`
by analytic continuation from positive `real(x)`).
"""
function lgamma end

# asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)|
@inline function lgamma_asymptotic(z::Complex{Float64})
    zinv = inv(z)
    t = zinv*zinv
    # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1))
    return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + # <-- log(2pi)/2
       zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03,
                         7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04,
                         8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03,
                         6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02)
end

# Compute the logΓ(z) function using a combination of the asymptotic series,
# the Taylor series around z=1 and z=2, the reflection formula, and the shift formula.
# Many details of these techniques are discussed in D. E. G. Hare,
# "Computing the principal branch of log-Gamma," J. Algorithms 25, pp. 221-236 (1997),
# and similar techniques are used (in a somewhat different way) by the
# SciPy loggamma function.  The key identities are also described
# at http://functions.wolfram.com/GammaBetaErf/LogGamma/
function lgamma(z::Complex{Float64})
    x = real(z)
    y = imag(z)
    yabs = abs(y)
    if !isfinite(x) || !isfinite(y) # Inf or NaN
        if isinf(x) && isfinite(y)
            return Complex(x, x > 0 ? (y == 0 ? y : copysign(Inf, y)) : copysign(Inf, -y))
        elseif isfinite(x) && isinf(y)
            return Complex(-Inf, y)
        else
            return Complex(NaN, NaN)
        end
    elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y|
        return lgamma_asymptotic(z)
    elseif x < 0.1 # use reflection formula to transform to x > 0
        if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0
            return Complex(Inf, signbit(x) ? copysign(oftype(x, pi), -y) : -y)
        end
        # the 2pi * floor(...) stuff is to choose the correct branch cut for log(sinpi(z))
        return Complex(1.1447298858494001741434262, # log(pi)
                       copysign(6.2831853071795864769252842, y) # 2pi
                       * floor(0.5*x+0.25)) -
               log(sinpi(z)) - lgamma(1-z)
    elseif abs(x - 1) + yabs < 0.1
        # taylor series around zero at z=1
        # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]]
        w = Complex(x - 1, y)
        return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,
                                -4.0068563438653142846657956e-01,2.705808084277845478790009e-01,
                                -2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,
                                -1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,
                                -1.1133426586956469049087244e-01,1.000994575127818085337147e-01,
                                -9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,
                                -7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,
                                -6.6668705882420468032903454e-02)
    elseif abs(x - 2) + yabs < 0.1
        # taylor series around zero at z=2
        # ... coefficients are [1-eulergamma; [(-1)^k * (zeta(k)-1)/k for k in 2:12]]
        w = Complex(x - 2, y)
        return w * @evalpoly(w, 4.2278433509846713939348812e-01,3.2246703342411321823620794e-01,
                               -6.7352301053198095133246196e-02,2.0580808427784547879000897e-02,
                               -7.3855510286739852662729527e-03,2.8905103307415232857531201e-03,
                               -1.1927539117032609771139825e-03,5.0966952474304242233558822e-04,
                               -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,
                               -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05)
    end
    # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series
    shiftprod = Complex(x,yabs)
    x += 1
    sb = false # == signbit(imag(shiftprod)) == signbit(yabs)
    # To use log(product of shifts) rather than sum(logs of shifts),
    # we need to keep track of the number of + to - sign flips in
    # imag(shiftprod), as described in Hare (1997), proposition 2.2.
    signflips = 0
    while x <= 7
        shiftprod *= Complex(x,yabs)
        sb′ = signbit(imag(shiftprod))
        signflips += sb′ & (sb′ != sb)
        sb = sb′
        x += 1
    end
    shift = log(shiftprod)
    if signbit(y) # if y is negative, conjugate the shift
        shift = Complex(real(shift), signflips*-6.2831853071795864769252842 - imag(shift))
    else
        shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842)
    end
    return lgamma_asymptotic(Complex(x,y)) - shift
end
lgamma(z::Complex{T}) where {T<:Union{Integer,Rational}} = lgamma(float(z))
lgamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(lgamma(Complex{Float64}(z)))

gamma(z::Complex) = exp(lgamma(z))

"""
    beta(x, y)

Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``.
"""
function beta(x::Number, w::Number)
    yx, sx = lgamma_r(x)
    yw, sw = lgamma_r(w)
    yxw, sxw = lgamma_r(x+w)
    return exp(yx + yw - yxw) * (sx*sw*sxw)
end

"""
    lbeta(x, y)

Natural logarithm of the absolute value of the [`beta`](@ref)
function ``\\log(|\\operatorname{B}(x,y)|)``.
"""
lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
module TestNAs
    using Base.Test
    using DataArrays

    # anyna(a::AbstractArray)
    @test anyna(Any[NA, 1])
    @test !anyna([1, 2])
    anyna(repeat([1, 2], outer = [1, 2]))
    @test !anyna(repeat([1, 2], outer = [1, 2, 2]))

    # anyna(da::DataArray)
    anyna(DataArray([1, 2], falses(2)))
    anyna(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
    da = DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2))
    @test !anyna(da)
    da[2] = NA
    @test anyna(da)

    # anyna(pda::PooledDataArray)
    anyna(PooledDataArray([1, 2], falses(2)))
    anyna(PooledDataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
    pda = PooledDataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2))
    @test !anyna(pda)
    pda[2] = NA
    @test anyna(pda)

    # allna(a::AbstractArray)
    @test allna(Any[NA, NA])
    @test !allna(Any[NA, 1])
    @test !allna([1, 2])
    allna(repeat([1, 2], outer = [1, 2]))
    @test !allna(repeat([1, 2], outer = [1, 2, 2]))

    # allna(da::DataArray)
    allna(DataArray([1, 2], falses(2)))
    allna(DataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
    da = DataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2))
    da[1] = NA
    @test !allna(da)
    da[:] = NA
    @test allna(da)

    # allna(da::PooledDataArray)
    allna(PooledDataArray([1, 2], falses(2)))
    allna(PooledDataArray(repeat([1, 2], outer = [1, 2]), falses(2, 2)))
    pda = PooledDataArray(repeat([1, 2], outer = [1, 2, 2]), falses(2, 2, 2))
    pda[1] = NA
    @test !allna(pda)
    pda[:] = NA
    @test allna(pda)

    dv = DataArray([1, 2, 3], BitArray([false, false, false]))

    dv = DataArray(collect(1:6), fill(false, 6))

    a = dropna(dv)
    for v in each_failna(dv); end
    @test collect(each_failna(dv)) == a
    @test collect(each_dropna(dv)) == a
    @test collect(each_replacena(dv, 4)) == a

    dv[[1, 2, end]] = NA

    a = dropna(dv)
    @test_throws NAException for v in each_failna(dv); end
    @test collect(each_dropna(dv)) == a
    @test collect(each_replacena(dv, 4)) == [4, 4, a..., 4]
end


text:
# #TODO: finish
# function hilbert_transform(x::Array{T}; n::Int=10)
# end

function estimate_rsrange(x::Array{T})::T where T<:Real
    n = size(x,1)
    @assert n>2 "need more than two elements, have $x"
    m = sum(x)/n
    y = x .- m
    z = cumsum(y)
    r = maximum(z) - minimum(z)
    s = sqrt(sum((x.-m).^2)/n)
    return r/s
end

function npieces(x::Array{T}; minsize::Int=8)::Int where {T<:Real}
    i = 0
    n = size(x,1)
    while n > minsize
        i += 1
        n = fld(n, 2)
    end
    return i
end

function genrsdata(x::Array{T})::Matrix{T} where {T<:Real}
    depth = npieces(x)
    if depth == 0
        return [size(x,1) estimate_rsrange(x)]
    end
    rsdata = zeros(0,2)
    rsdata = [rsdata; [size(x,1) estimate_rsrange(x)]]
    a, b = divide(x)
    rsdata = [rsdata; genrsdata(a)]
    rsdata = [rsdata; genrsdata(b)]
    return rsdata
end

function divide(x::A)::Tuple{A, A} where {A<:AbstractArray}
    n = size(x,1)
    @assert n>=2
    h = n/2
    a = x[1:floor(Int,h)]
    if floor(h) == h
        b = x[(ceil(Int,h)+1):end]
    else
        b = x[ceil(Int,h):end]
    end
    return a, b
end

function estimate_hurst(x::Array{T}; intercept::Bool=false)::T where {T<:Real}
    RS = genrsdata(x)
    RS = RS[sortperm(RS[:,1]),:]
    xx = log2.(RS[:,1])
    yy = log2.(RS[:,2])
    if intercept
        _, beta = [ones(size(xx)) xx]\yy
    else
        beta = xx\yy
    end
    return beta
end

"""
```
rsrange(x::Array{T}; n::Int=100, cumulative::Bool=false, intercept::Bool=true)
```

Compute the rescaled range of a time series
"""
function rsrange(x::Array{T}; n::Int=100, cumulative::Bool=false, intercept::Bool=true) where {T<:Real}
    @assert size(x,1) >= n
    return runfun(x, estimate_rsrange; n=n, cumulative=cumulative)
end

"""
```
hurst(x::Array{T}; n::Int=100, cumulative::Bool=false, intercept::Bool=false)
```

Compute the Hurst exponent of a time series
"""
function hurst(x::Array{T}; n::Int=100, cumulative::Bool=false, intercept::Bool=false) where {T<:Real}
    @assert size(x,1) >= n
    return runfun(x, estimate_hurst; n=n, cumulative=cumulative, intercept=intercept)
end


text:
function johnson_shortest_paths(g::AbstractGraph{U},
distmx::AbstractMatrix{T}=weights(g)) where T <: Real where U <: Integer

    nvg = nv(g)
    type_distmx = typeof(distmx)
#Change when parallel implementation of Bellman Ford available
    wt_transform = bellman_ford_shortest_paths(g, vertices(g), distmx).dists

    if !type_distmx.mutable && type_distmx !=  LightGraphs.DefaultDistance
        distmx = sparse(distmx) #Change reference, not value
    end

#Weight transform not needed if all weights are positive.
    if type_distmx !=  LightGraphs.DefaultDistance
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[src(e)] - wt_transform[dst(e)]
        end
    end

    dijk_state = Parallel.dijkstra_shortest_paths(g, vertices(g), distmx)
    dists = dijk_state.dists
    parents = dijk_state.parents

    broadcast!(-, dists, dists, wt_transform)
    for v in vertices(g)
        dists[:, v] .+= wt_transform[v] #Vertical traversal prefered
    end

    if type_distmx.mutable
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[dst(e)] - wt_transform[src(e)]
        end
    end

    return JohnsonState(dists, parents)
end


text: println("hey world")

max_stars_repo_path: Julia/hey world.jl

max_stars_repo_name: kennethsequeira/Hello-world

max_issues_repo_path: Julia/hey world.jl

max_issues_repo_name: kennethsequeira/Hello-world

max_forks_repo_path: Julia/hey world.jl

max_forks_repo_name: kennethsequeira/Hello-world