text:
@testset "Static graphs" begin
    g = @inferred(CompleteDiGraph(5))
    @test nv(g) == 5 && ne(g) == 20
    g = @inferred(CompleteGraph(5))
    @test nv(g) == 5 && ne(g) == 10

    g = @inferred(CompleteBipartiteGraph(5, 8))
    @test nv(g) == 13 && ne(g) == 40

    g = @inferred(StarDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(StarGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(StarGraph(1))
    @test nv(g) == 1 && ne(g) == 0

    g = @inferred(PathDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(PathGraph(5))
    @test nv(g) == 5 && ne(g) == 4

    g = @inferred(CycleDiGraph(5))
    @test nv(g) == 5 && ne(g) == 5
    g = @inferred(CycleGraph(5))
    @test nv(g) == 5 && ne(g) == 5

    g = @inferred(WheelDiGraph(5))
    @test nv(g) == 5 && ne(g) == 8
    g = @inferred(WheelGraph(5))
    @test nv(g) == 5 && ne(g) == 8

    g = @inferred(Grid([3,3,4]))
    @test nv(g) == 3*3*4
    @test ne(g) == 75
    @test Δ(g) == 6
    @test δ(g) == 3

    g = @inferred(Grid([3,3,4], periodic=true))
    @test nv(g) == 3*3*4
    @test ne(g) == 108
    @test Δ(g) == 6
    @test δ(g) == 6

    g = @inferred(CliqueGraph(3,5))
    @test nv(g) == 15 && ne(g) == 20
    @test g[1:3] == CompleteGraph(3)

    g = @inferred(crosspath(3, BinaryTree(2)))
    # f = Vector{Vector{Int}}[[2 3 4];
    # [1 5];
    # [1 6];
    # [1 5 6 7];
    # [2 4 8];
    # [3 4 9];
    # [4 8 9];
    # [5 7];
    # [6 7]
    # ]
    I = [1,1,1,2,2,3,3,4,4,4,4,5,5,5,6,6,6,7,7,7,8,8,9,9]
    J = [2,3,4,1,5,1,6,1,5,6,7,2,4,8,3,4,9,4,8,9,5,7,6,7]
    V = ones(Int, length(I))
    Adj = sparse(I,J,V)
    @test Adj == sparse(g)

    g = @inferred(DoubleBinaryTree(3))
    # [[3,2,8]
    # [4,1,5]
    # [1,6,7]
    # [2]
    # [2]
    # [3]
    # [3]
    # [10,9,1]
    # [11,8,12]
    # [8,13,14]
    # [9]
    # [9]
    # [10]
    # [10]]
    I = [1,1,1,2,2,2,3,3,3,4,5,6,7,8,8,8,9,9,9,10,10,10,11,12,13,14]
    J = [3,2,8,4,1,5,1,6,7,2,2,3,3,10,9,1,11,8,12,8,13,14,9,9,10,10]
    V = ones(Int, length(I))
    Adj = sparse(I,J,V)
    @test Adj == sparse(g)

    rg3 = @inferred(RoachGraph(3))
    # [3]
    # [4]
    # [1,5]
    # [2,6]
    # [3,7]
    # [4,8]
    # [9,8,5]
    # [10,7,6]
    # [11,10,7]
    # [8,9,12]
    # [9,12]
    # [10,11]
    I = [ 1,2,3,3,4,4,5,5,6,6,7,7,7,8,8,8,9,9,9,10,10,10,11,11,12,12 ]
    J = [ 3,4,1,5,2,6,3,7,4,8,9,8,5,10,7,6,11,10,7,8,9,12,9,12,10,11 ]
    V = ones(Int, length(I))
    Adj = sparse(I,J,V)
    @test Adj == sparse(rg3)
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
include("..//src//utils//logplots.jl")

using Scruff
using Scruff.SFuncs
using Scruff.Utils
using Scruff.RTUtils
using Scruff.Models
using Scruff.Algorithms

pyplot()

x1m = Cat([1,2], [0.1, 0.9])()
x2m = Cat([1,2,3], [0.2, 0.3, 0.5])()
cpd2 = Dict((1,1) => [0.3, 0.7], (1,2) => [0.6, 0.4], (2,1) =>[0.4, 0.6],
            (2,2) => [0.7, 0.3], (3,1) => [0.5, 0.5], (3,2) => [0.8, 0.2])
x3m = DiscreteCPT([1,2], cpd2)()
x4m = DiscreteCPT([1,2], Dict((1,) => [0.15, 0.85], (2,) => [0.25, 0.75]))()
x5m = DiscreteCPT([1,2], Dict((1,) => [0.35, 0.65], (2,) => [0.45, 0.55]))()
x6m = DiscreteCPT([1,2], Dict((1,) => [0.65, 0.35], (2,) => [0.75, 0.25]))()

x1 = x1m(:x1)
x2 = x2m(:x2)
x3 = x3m(:x3)
x4 = x4m(:x4)
x5 = x5m(:x5)
x6 = x6m(:x6)

fivecpdnet = InstantNetwork(Variable[x1,x2,x3,x4,x5], VariableGraph(x3=>[x2,x1], x4=>[x3], x5=>[x3]))

run = Runtime(fivecpdnet)
# default_initializer(run)

logger = BPLogger([:x1,:x2])
@info "logger=$logger"
# run em
with_logger(logger) do
    Scruff.Algorithms.three_pass_BP(run)
end

ext: jl


text:
"""
Angle bar, transient vibration.

The bar is excited by a Hann-windowed force pulse.

Reference:
Ostachowicz W, Kudela P, Krawczuk M, Zak A. Guided Waves in Structures for SHM: The Time - domain Spectral Element Method. A John Wiley & Sons, Ltd., 2011.
"""
module angle_bar_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 barside = 250*phun("mm")
const barlength = 1000*phun("mm")
const barthickness = 10*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 totalforce = 1*phun("N")
const forcepatchradius = 2*phun("mm")
const forcedensity = 2*totalforce/(pi*forcepatchradius^2)
const tend = 0.3*phun("milli*s")
const visualize = true
const color = "magenta"
const distributedforce = 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)

    tolerance = barthickness/nref/10
    X = [0.0 0.0 0.0;
    0.0 barside 0.0;
    0.0 0.0 barside;
    ]
    l2fens = FENodeSet(X);
    conn = [2 1; 1 3; ]
    l2fes = FESetL2(conn);
    nLayers = 4
    ex = function (x, k)
        x[1] += k*barlength/nLayers
        x
    end
    fens,fes = Q4extrudeL2(l2fens, l2fes, nLayers, ex); # Mesh
    fens, fes = Q4toT3(fens, fes)
    for r in 1:nref
        fens, fes = T3refine(fens, fes)
    end
    @show count(fens)
    vtkwrite("angle_bar.vtu", fens, fes)

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

    # 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), barthickness), 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=[barlength/2 barside 0.0])[1]
    cpoint = selectnode(fens; nearestto=[barlength/2 0 0])[1]

    points = Dict(
    "Load" => selectnode(fens; nearestto=[barlength/2 barside 0.0])[1],
    "Crease" => selectnode(fens; nearestto=[barlength/2 0 0])[1],
    )

    pointdofs = Dict()
    for k in keys(points)
        pointdofs[k] = dchi.dofnums[points[k], 3]
    end

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[points["Load"], 1]
        dy = XYZ[2] - fens.xyz[points["Load"], 2]
        dz = XYZ[3] - fens.xyz[points["Load"], 3]
        d = sqrt(dx^2+dy^2+dz^2)
        forceout .= 0.0
        if d < forcepatchradius
            forceout[3] = -forcedensity
        end
        return forceout
    end

    if distributedforce
        lfemm = FEMMBase(IntegDomain(fes, TriRule(6)))
        fi = ForceIntensity(FFlt, 6, computetrac!);
        Fmag = distribloads(lfemm, geom0, dchi, fi, 2);
    else
        Fmag = fill(0.0, dchi.nfreedofs)
        Fmag[pointdofs["Load"]] = -totalforce
    end

    @show sum(Fmag), -totalforce

    function force!(F, t)
        mul = 0.0
        if t <= 1/modulation_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))
    pointdeflections = Dict()
    for k in keys(points)
        pointdeflections[k] = fill(0.0, nsteps+1)
    end
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        for k in keys(points)
            b = pointdeflections[k]
            b[step+1] = U[pointdofs[k]]
        end
        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' "  :-
        for k in sort(string.(keys(points)))
            @gp  :- collect(0.0:dt:(nsteps*dt)) pointdeflections[k] " lw 2 lc rgb '$color' with lines title 'Deflection at $(string(k))' "  :-
        end

        @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("angle_bar_$nref", fens, fes, times; vectors = vectors)
    end
end

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

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

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

end # module
nothing


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 8bc61658-b880-11eb-0b6c-a5b7ff4dcf3e
using SampledSignals

# ╔═╡ 821a6565-7046-46ff-b874-521fd846d703
using WAV

# ╔═╡ 06853df6-44a0-42cb-8e5a-95bbc8962b81
using Base, UnicodePlots

# ╔═╡ aa8e6100-5919-461a-8394-dea283f15b00
using Unitful

# ╔═╡ bc1888ec-9c33-45e7-91cb-4a7024f1acd7
begin

# some examples using signal operators
using SignalOperators

# a pure tone 20 dB below a power 1 signal, with on and off ramps (for
# a smooth onset/offset)
sound1 = Signal(sin,ω=1kHz) |> Until(5s) |> Ramp |> Normpower |> Amplify(-20dB)

# a sound defined by a file, matching the overall power to that of sound1
sound2 = "example.wav" |> Normpower |> Amplify(-20dB)

# a 1kHz sawtooth wave
sound3 = Signal(ϕ -> ϕ-π,ω=1kHz) |> Ramp |> Normpower |> Amplify(-20dB)

# a 5 Hz amplitude modulated noise
sound4 = randn |>
    Amplify(Signal(ϕ -> 0.5sin(ϕ) + 0.5,ω=5Hz)) |>
    Until(5s) |> Normpower |> Amplify(-20dB)

# a 1kHz tone surrounded by a notch noise
SNR = 5dB
x = Signal(sin,ω=1kHz) |> Until(1s) |> Ramp |> Normpower |> Amplify(-20dB + SNR)
z = Signal(randn) |> Until(1s) |> Filt(Bandstop,0.5kHz,2kHz) |> Normpower |> Amplify(-20dB)
scene = Mix(x,z)

# # write all of the signals to a single file, at 44.1 kHz
# Append(sound1,sound2,sound3,sound4,scene) |> ToFramerate(44.1kHz) |> sink("examples.wav")
end

# ╔═╡ 770f4a6b-2fe3-4493-a324-08a54131ea88
md"""
# Oscillators
It is the most important component of a synthesizer, it is a component that
that generates a sequence of numbers that repeats after a certain interval.

The most simplest oscillator is one that produces which is generated by
unsurprisingly yet delightfully, the sine function.

https://www.reddit.com/r/Python/comments/lw50ne/making_a_synthesizer_using_python/
https://python.plainenglish.io/making-a-synth-with-python-oscillators-2cb8e68e9c3b
"""

# ╔═╡ fe3abf37-bcad-4717-9ebd-702838c5ccbf
begin    # Example audio
    audio_one_channel = SampleBuf(rand(100000), 10000)
    audio_two_channel = SampleBuf(rand(100000,2), 10000)
    audio_multi_channel = SampleBuf(rand(100000,2), 10000)
end

# ╔═╡ 961a4dfa-0cbb-46aa-8af1-a056ade120ad
function get_sin_oscillator(freq=1, sample_rate=512)
	increment = ( 2 * pi * freq) / sample_rate
	return(sin(v) for v in Iterators.countfrom(0, increment))
end

# ╔═╡ dfe31fdd-6be0-41cf-8fae-77c0e8692627
begin
osc =  get_sin_oscillator(1, 512)
samples = [osc for i in 1:512]
# samples = [next(osc,3) for i in 1:512]
end

# ╔═╡ d1890135-cbcf-4b23-a32c-418be77debc0
samples

# ╔═╡ ea2c8ea0-7f9d-4687-b657-21e4d2a0c253
sizeof(osc)

# ╔═╡ 764cb127-bcec-450f-8ed5-933e3bc1db6c
begin

# some examples from WaveForm.jl
squarewave(x::Real) = ifelse(mod2pi(x) < π, 1.0, -1.0)

md"""
    squarewave(x, θ)
Compute ``2\pi``-periodic square wave of `x` with a duty cycle `θ`
and a peak amplitude ``1``.
"""
function squarewave(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod2pi(x) < 2π * θ, 1.0, -1.0)
end

md"""
    squarewave1(x)
Compute ``1``-periodic square wave of `x` with a peak amplitude ``1``.
"""
squarewave1(x::Real) = ifelse(mod(x, 1) < 1/2, 1.0, -1.0)

md"""
    squarewave1(x, θ)
Compute ``1``-periodic square wave of `x` with a duty cycle `θ`
and a peak amplitude ``1``.
"""
function squarewave1(x::Real, θ::Real)
    0 ≤ θ ≤ 1 || throw(DomainError(θ, "squwarewave1(x, θ) is only defined for 0 ≤ θ ≤ 1."))
    ifelse(mod(x, 1) < θ, 1.0, -1.0)
end

md"""
    trianglewave(x)
Compute ``2\pi``-periodic triangle wave of `x` with a peak amplitude ``1``.
"""
function trianglewave(x::Real)
    modx = mod2pi(x + π/2)
    ifelse(modx < π, 2modx/π - 1, -2modx/π + 3)
end

md"""
    trianglewave1(x)
Compute ``1``-periodic triangle wave of `x` with a peak amplitude ``1``.
"""
function trianglewave1(x::Real)
    modx = mod(x + 1/4, 1.0)
    ifelse(modx < 1/2, 4modx - 1, -4modx + 3)
end

md"""
    sawtoothwave(x)
Compute ``2\pi``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
"""
sawtoothwave(x::Real) = rem2pi(x, RoundNearest) / π

md"""
    sawtoothwave1(x)
Compute ``1``-periodic sawtooth wave of `x` with a peak amplitude ``1``.
"""
sawtoothwave1(x::Real) = rem(x, 1.0, RoundNearest) * 2

end

# ╔═╡ 29ef8525-5c0e-4dfd-8313-c3232ae362ab
squarewave(1)

# ╔═╡ a98ee93f-aac1-4b80-b337-f96f17744561
lineplot([squarewave], -π/2, 2π)

# ╔═╡ 466a36e8-c298-4ff9-ae0a-87c92905cf2a
lineplot([sawtoothwave], -π/2, 2π)

# ╔═╡ 3122c57a-f46f-49d7-9770-b51b1347b8f9
lineplot([trianglewave], -π/2, 2π)

# ╔═╡ 11b6cb26-da27-4222-94fe-40ef41f1c954
  audio_multi1_channel = SampleBuf(rand(100000,2), 10000)

# ╔═╡ acde3a7e-b314-4e3d-96d7-81b63821e060
typeof(SampleBuf(rand(100000,2), 10000))

# ╔═╡ 1dd3020a-9d13-4ab5-9967-14644d72f6df
audio_multi1_channel.data

# ╔═╡ 7f0b5cd9-24c7-4c40-a62c-0d2e909c1ae9
# Using Dan Casmiro WAV.jl
begin
fs = 8e3
t = 0.0:1/fs:prevfloat(1.0)
f = 1e3
y = sin.(2pi * f * t) * 0.1
wavwrite(y, "example.wav", Fs=fs)

y, fs = wavread("example.wav")
y = sin.(2pi * 2f * t) * 0.1
wavappend(y, "example.wav")

y, fs = wavread("example.wav")
wavplay(y, fs)
end

# ╔═╡ e7a0111f-1df0-4a36-8a24-d68c059b8cad
begin
nsamples = []
for i in 0:1/10000:10
	push!(nsamples, sinpi(f * i))
end
end

# ╔═╡ 3b19dcd5-37ef-4676-924b-e0b0778b1dd5
nsamples1 = hcat(nsamples, nsamples)

# ╔═╡ 4ff67c7e-73c8-4cfa-91a0-392ab3af2318
# using sampledsignals
audio_samples = SampleBuf(nsamples1, 10000)

# ╔═╡ 9b3bbcb7-544c-43b6-957a-495291525736
# using wavplay
wavplay(audio_samples.data,10000)

# ╔═╡ b92e9a63-f08e-47fe-baab-f5cc309bcdbc
# using SignalOperators.Units # allows the use of dB, Hz, s etc... as unitful values

# ╔═╡ b271de0b-3566-4712-bcdc-06a66f6f0b86
sound1

# ╔═╡ Cell order:
# ╠═8bc61658-b880-11eb-0b6c-a5b7ff4dcf3e
# ╠═821a6565-7046-46ff-b874-521fd846d703
# ╠═06853df6-44a0-42cb-8e5a-95bbc8962b81
# ╠═770f4a6b-2fe3-4493-a324-08a54131ea88
# ╠═fe3abf37-bcad-4717-9ebd-702838c5ccbf
# ╠═961a4dfa-0cbb-46aa-8af1-a056ade120ad
# ╠═dfe31fdd-6be0-41cf-8fae-77c0e8692627
# ╠═d1890135-cbcf-4b23-a32c-418be77debc0
# ╠═ea2c8ea0-7f9d-4687-b657-21e4d2a0c253
# ╠═764cb127-bcec-450f-8ed5-933e3bc1db6c
# ╠═29ef8525-5c0e-4dfd-8313-c3232ae362ab
# ╠═a98ee93f-aac1-4b80-b337-f96f17744561
# ╠═466a36e8-c298-4ff9-ae0a-87c92905cf2a
# ╠═3122c57a-f46f-49d7-9770-b51b1347b8f9
# ╠═e7a0111f-1df0-4a36-8a24-d68c059b8cad
# ╠═3b19dcd5-37ef-4676-924b-e0b0778b1dd5
# ╠═9b3bbcb7-544c-43b6-957a-495291525736
# ╠═4ff67c7e-73c8-4cfa-91a0-392ab3af2318
# ╠═11b6cb26-da27-4222-94fe-40ef41f1c954
# ╠═acde3a7e-b314-4e3d-96d7-81b63821e060
# ╠═1dd3020a-9d13-4ab5-9967-14644d72f6df
# ╠═7f0b5cd9-24c7-4c40-a62c-0d2e909c1ae9
# ╠═b92e9a63-f08e-47fe-baab-f5cc309bcdbc
# ╠═aa8e6100-5919-461a-8394-dea283f15b00
# ╠═bc1888ec-9c33-45e7-91cb-4a7024f1acd7
# ╠═b271de0b-3566-4712-bcdc-06a66f6f0b86

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
# Included in benchmark.jl
using Distributed
using Dates
"Adds process 2, same LOAD_PATH as process 1"
function prepareworker()
    # Prepare worker
    FULLLOADPATH = LOAD_PATH
    if nworkers() < 2
        addprocs(1)
    end
    futur = @spawnat 2  for p in FULLLOADPATH
                            p ∉ LOAD_PATH && push!(LOAD_PATH, p)
                        end
    # waits for process 2 to get ready, trigger errors if any
    fetch(futur)
end

"Start and wait for async hts server"
function start_hts(timeout)
    hts_task = @async ws_hts.listen_hts()
    t1 = now() + timeout
    while now() < t1
        sleep(0.5)
        isdefined(ws_hts.TCPREF, :x) && break
    end
    if now()>= t1
        msg = " did not establish server before timeout "
        clog("start_hts", :red, msg, timeout)
        error(msg)
    end
    hts_task
end

"""
Make connection and get a reference to the HTS side of HTS-JCE connection.
Assumes a server is already running.
"""
function get_hts_jce()
    id = "get_hts_jce"
    # Make the parallel worker initiate a websocket connection to this process.
    # It will finish execution when either peer closes the connection.
    @spawnat 2  ws_jce.echowithdelay_jce()
    # async HTS will not accept the connection before we yield to it
    zlog(id, "JCE will be connecting to HTS. We ask HTS for a reference...")
    zflush()
    hts = Union{WebSocket, String}("")
    t1 = now() + TIMEOUT
    while now() < t1
        yield()
        hts = ws_hts.getws_hts()
        isa(hts, WebSocket) && break
        sleep(0.5)
    end
    return hts
end
"""
Make and get a reference to the server side of HTS-BCE websocket connection.
Will launch a different web browser every time. Returns a string
if all available browsers have been launched one time.

Assumes a server is already running.
"""
function get_hts_bce()
    id = "get_hts_jce"
    hts = Union{WebSocket, String}("")
    browser =""
    opened, browser = open_a_browser()
    # Launch the next browser in line.
    if opened
        zlog(id, "BCE in, ", browser, " will be connecting to HTS. Getting reference...")
        zflush()
        hts = Union{WebSocket, String}("")
        t1 = now() + TIMEOUT
        while now() < t1
            yield()
            hts = ws_hts.getws_hts()
            isa(hts, WebSocket) && break
            sleep(0.5)
        end
        return hts, browser
    end
    return hts, browser
end
"""
Send n messages of length messagesize HTS-JCE and back.
    -> id, serverlatencies, clientlatencies
Starts and closes a new HTS-JCE connection every function call.
In the terminology of BenchmarkTools, a call is a sample
consisting of n evaluations.
"""
function HTS_JCE(n, messagesize)
    id = "HTS_JCE"
    zlog(id, "Warming up, compiling")
    zflush()
    hts = get_hts_jce()
    if !isa(hts, WebSocket)
        msg = " could not get websocket reference"
        clog(id, msg)
        error(id * msg)
    end
    clog(id, hts)
    # Random seeding, same for all samples
    srand(1)
    clog(id, "Sending ", n, " messages of ", messagesize , " random bytes")
    sendtime = Int64(0)
    receivereplytime = Int64(0)
    sendtimes = Vector{Int64}()
    receivereplytimes =  Vector{Int64}()
    for i = 1:n
        msg = rand(0x20:0x7f, messagesize)
        sendtime = time_ns()
        write(hts, msg)
        # throw away replies
        readguarded(hts)
        receivereplytime = time_ns()
        push!(receivereplytimes, Int64(receivereplytime < typemax(Int64) ? receivereplytime : 0 ))
        push!(sendtimes,  Int64(sendtime < typemax(Int64) ? sendtime : 0 ))
    end
    #
    zlog(id, "Sending 'exit', JCE will send its time log and then exit.")
    zflush()
    write(hts, "exit")
    # We deserialize JCE's time records from this sample
    bs = Base.BufferStream()
    write(bs, read(hts))
    close(bs)
    if bytesavailable(bs) > 0
        receivetimes, replytimes = deserialize(bs)
    else
        error("Did not receive receivetimes, replytimes")
    end
    # We must read from the websocket in order for it to respond to
    # a closing message from JCE. It's also nice to yield to the async server
    # so it can exit from it's handler and release the websocket reference.
    isopen(hts) && readguarded(hts)
    yield()
    serverlatencies = receivetimes - sendtimes
    clientlatencies = receivereplytimes - replytimes
    return id, serverlatencies, clientlatencies
end
"""
Use the next browser in line, start a new HTS-BCE connection and collect n evaluations.
This is one sample.
Returns browser name and vectors with n rows:
# t1  Send 0, receive 0, measure time interval
# t2  Send 0, receive 0 twice, measure time interval
# t3  Send x, receive 0, measure time interval
# t4  Send x, receive 0, receive x, measure time interval
"""
function HTS_BCE(n, x)
    id = "HTS_BCE"
    zlog(id)
    zflush()
    msg = ""
    (hts, browser) = get_hts_bce()
    if browser == ""
        msg = "Could not find and open more browser types"
    elseif !isa(hts, WebSocket)
        msg = " could not get ws reference from " * browser * " via HTS"
    end
    if msg != ""
        clog(id, msg)
        return "", Vector{Int}(), Vector{Int}(), Vector{Int}(), Vector{Int}()
    end
    clog(id, hts)
    # Random seeding, same for all samples
    srand(1)
    clog(id, "Sending ", n, " messages of ", x , " random bytes")
    st1 = UInt64(0)
    st2 = UInt64(0)
    rt1 = UInt64(0)
    rt2 = UInt64(0)
    rt3 = UInt64(0)
    rt4 = UInt64(0)

    t1 =  Vector{Int}()
    t2 =  Vector{Int}()
    t3 =  Vector{Int}()
    t4 =  Vector{Int}()

    for i = 1:n
            # Send empty message, measure time after two empty replies
            msg = Vector{UInt8}()
            st1 = time_ns()
            write(hts, msg)
            read(hts)
            rt1 = time_ns()
            read(hts)
            rt2 = time_ns()

            # Send message, measure time after one empty and one echoed replies
            msg = rand(0x20:0x7f, x)
            st2 = time_ns()
            write(hts, msg)
            # throw away replies. We expect a sequence per sent message:
            # two empty messages, then the original content.
            read(hts)
            rt3 = time_ns()
            read(hts)
            rt4 = time_ns()
            # Store time intervals
            push!(t1, Int(rt1 - st1))
            push!(t2, Int(rt2 - st1))
            push!(t3, Int(rt3 - st2))
            push!(t4, Int(rt4 - st2))
    end
    #
    zlog(id, "Close websocket")
    close(hts)
    # Also yield to the server task so it can release it's reference.
    sleep(1)
    # Return browser name and measured time interval vectors [ns]
    browser, t1, t2, t3, t4
end

"
Constant message size [b], measured time interval vectors [ns]
    -> server and client speeds, server and client bandwidth [ns/b]
"
function serverandclientspeeds(messagesize, serverlatencies, clientlatencies)
    serverspeeds =  serverlatencies / messagesize
    clientspeeds =  clientlatencies / messagesize
    n = length(serverspeeds)
    serverbandwidth = sum(serverlatencies) / (n * messagesize)
    clientbandwidth = sum(clientlatencies) / (n * messagesize)
    serverspeeds, clientspeeds, serverbandwidth, clientbandwidth
end

"
Measured time interval vectors [ns] -> client and server speeds, bandwidth [ns/b]
x is message size [b].
"
function serverandclientspeeds_indirect(x, t1, t2, t3, t4)
    # assuming a measured time interval consists of
    # t(x) = t0s + t0c + a*x + b*x
    # where
    # t(x)      measured time at the server
    # t0s       initial serverlatency, for a "zero length" message
    # t0c       initial clientlatency, for a "zero length" message
    # x         message length [b]
    # a         marginal server speed [ns/b]
    # b         marginal client speed [ns/b]
    #
    #
    # t1 = t0s +  t0c              Send 0, receive 0, measure t1
    # t2 = t0s + 2t0c              Send 0, receive 0 twice, measure t2
    # t3 = t0s +  t0c + a*x        Send x, receive 0, measure t3
    # t4 = t0s + 2t0c + a*x + b*x  Send x, receive 0, receive x, measure t4
    #
    # hence,
    t0s =  2t1 - t2
    t0c = -t1 + t2
    a = (- t1   + t3) / x
    b =   (t1 - t2  - t3 + t4) / x
    # Drop the marginal speed info, return effective speeds
    # serverspeeds = (t0s + a*x) / x
    # clientspeeds = (t0s + a*x) / x
    serverspeeds = t0s / x + a
    clientspeeds = t0c / x + b
    n = length(serverspeeds)
    serverbandwidth = sum(serverspeeds) / n
    clientbandwidth = sum(clientspeeds) / n
    return serverspeeds, clientspeeds, serverbandwidth, clientbandwidth
end

## Note that the shorthand plot functions below require input symbols
## that are defined at module-level (not in a local scope)

"Generate a time series lineplot"
lp(sy::Symbol) = lineplot(collect(1:length(eval(sy))), eval(sy), title = String(sy), width = displaysize(stdout)[2]-20, canvas = AsciiCanvas)

"Generate a vector of time series lineplots with a common title prefix"
function lp(symbs::Vector{Symbol}, titleprefix)
    map(symbs) do sy
        pl = lp(sy)
        title!(pl, titleprefix * " " * title(pl))
    end
end

"Generate an x-y lineplot in REPL"
function lp(syx::Symbol, syy::Symbol)
    lpl = lineplot(eval(syx), eval(syy), title = String(syy), width = displaysize(stdout)[2]-20, canvas = AsciiCanvas)
    xlabel!(lpl, String(syx))
end

"Generate a vector of x-y lineplots with a common title prefix"
function lp(syxs::Vector{Symbol}, syys::Vector{Symbol}, titleprefix)
    map(zip(syxs, syys)) do pair
        pl = lp(pair[1], pair[2])
        title!(pl, titleprefix * " " * title(pl))
    end
end

"Generate an x-y scatterplot"
function sp(syx::Symbol, syy::Symbol)
    spl = scatterplot(eval(syx), eval(syy), title = String(syy), width = displaysize(stdout)[2]-15, canvas = DotCanvas)
    xlabel!(spl, String(syx))
end


text:
# Create CompositeMean, CompositeKernel and CompositeCrossKernel.
for (composite_type, parent_type) in [(:CompositeMean, :MeanFunction),
                                      (:CompositeKernel, :Kernel),
                                      (:CompositeCrossKernel, :CrossKernel),]
@eval struct $composite_type{Tf, N} <: $parent_type
    f::Tf
    x::Tuple{Vararg{Any, N}}
    $composite_type(f::Tf, x::Vararg{Any, N}) where {Tf, N} = new{Tf, N}(f, x)
end
end

# CompositeMean definitions.
eachindex(c::CompositeMean) = eachindex(c.x[1])
length(c::CompositeMean) = length(c.x[1])
(μ::CompositeMean)(x) = map(μ.f, map(f->f(x), μ.x)...)
_map(f::CompositeMean, X::AV) = map(f.f, map(f->map(f, X), f.x)...)
map(f::CompositeMean, ::Colon) = map(f.f, map(f->map(f, :), f.x)...)

# CompositeKernel definitions.
(k::CompositeKernel)(x, x′) = map(k.f, map(f->f(x, x′), k.x)...)
(k::CompositeKernel)(x) = map(k.f, map(f->f(x), k.x)...)
length(k::CompositeKernel) = size(k.x[1], 1)
isstationary(k::CompositeKernel) = all(map(isstationary, k.x))
@noinline eachindex(k::CompositeKernel) = eachindex(k.x[1], 1)

_map(f::CompositeKernel, X::AV) = map(f.f, map(f->map(f, X), f.x)...)
_pairwise(f::CompositeKernel, X::AV) = LazyPDMat(map(f.f, map(f->pairwise(f, X), f.x)...))
_map(f::CompositeKernel, X::AV, X′::AV) = map(f.f, map(f->map(f, X, X′), f.x)...)
_pairwise(f::CompositeKernel, X::AV, X′::AV) = map(f.f, map(f->pairwise(f, X, X′), f.x)...)

map(f::CompositeKernel, ::Colon) = map(f.f, map(f->map(f, :), f.x)...)
pairwise(f::ConstantKernel, ::Colon) = LazyPDMat(map(f.f, map(f->pairwise(f, :), f.x)...))

# CompositeCrossKernel definitions.
(k::CompositeCrossKernel)(x, x′) = map(k.f, map(f->f(x, x′), k.x)...)
size(k::CompositeCrossKernel, N::Int) = size(k.x[1], N)
isstationary(k::CompositeCrossKernel) = all(map(isstationary, k.x))
eachindex(k::CompositeCrossKernel, dim::Int) = eachindex(k.x[1], dim)

_map(f::CompositeCrossKernel, X::AV, X′::AV) = map(f.f, map(f->map(f, X, X′), f.x)...)
function _pairwise(f::CompositeCrossKernel, X::AV, X′::AV)
    return map(f.f, map(f->pairwise(f, X, X′), f.x)...)
end

map(f::CompositeCrossKernel, ::Colon) = map(f.f, map(f->map(f, :), f.x)...)
pairwise(f::CompositeCrossKernel, ::Colon) = map(f.f, map(f->pairwise(f, :), f.x)...)

############################## Multiply-by-Function Kernels ################################

"""
    LhsCross <: CrossKernel

A cross-kernel given by `k(x, x′) = f(x) * k(x, x′)`.
"""
struct LhsCross <: CrossKernel
    f::Any
    k::CrossKernel
end
(k::LhsCross)(x, x′) = k.f(x) * k.k(x, x′)
size(k::LhsCross, N::Int) = size(k.k, N)
_pairwise(k::LhsCross, X::AV, X′::AV) = map(k.f, X) .* pairwise(k.k, X, X′)
eachindex(k::LhsCross, dim::Int) = eachindex(k.k, dim)

# map(k::LhsCross, ::Colon) = map(k.f, :) .* map(k.k, :)
# pairwise(k::LhsCross, ::Colon) = map(k.f, :) .* pairwise(k.k, :)

"""
    RhsCross <: CrossKernel

A cross-kernel given by `k(x, x′) = k(x, x′) * f(x′)`.
"""
struct RhsCross <: CrossKernel
    k::CrossKernel
    f
end
(k::RhsCross)(x, x′) = k.k(x, x′) * k.f(x′)
size(k::RhsCross, N::Int) = size(k.k, N)
_pairwise(k::RhsCross, X::AV, X′::AV) = pairwise(k.k, X, X′) .* map(k.f, X′)'
eachindex(k::RhsCross, dim::Int) = eachindex(k.k, dim)

# map(k::RhsCross, ::Colon) = map(k.k, :) .* map(k.f, :)
# pairwise(k::RhsCross, ::Colon) = pairwise(k.k, :) .* map(k.f, :)'

"""
    OuterCross <: CrossKernel

A kernel given by `k(x, x′) = f(x) * k(x, x′) * f(x′)`.
"""
struct OuterCross <: CrossKernel
    f
    k::CrossKernel
end
(k::OuterCross)(x, x′) = k.f(x) * k.k(x, x′) * k.f(x′)
size(k::OuterCross, N::Int) = size(k.k, N)
function _pairwise(k::OuterCross, X::AV, X′::AV)
    return map(k.f, X) .* pairwise(k.k, X, X′) .* map(k.f, X′)'
end
eachindex(k::OuterCross, dim::Int) = eachindex(k.k, dim)

# map(k::OuterCross, ::Colon) = map(k.f, :) .* pairwise(k.k, :) .* map(k.f, :)'

"""
    OuterKernel <: Kernel

A kernel given by `k(x, x′) = f(x) * k(x, x′) * f(x′)`.
"""
struct OuterKernel <: Kernel
    f
    k::Kernel
end
(k::OuterKernel)(x, x′) = k.f(x) * k.k(x, x′) * k.f(x′)
(k::OuterKernel)(x) = k.f(x)^2 * k.k(x)
length(k::OuterKernel) = length(k.k)
_pairwise(k::OuterKernel, X::AV) = Xt_A_X(pairwise(k.k, X), Diagonal(map(k.f, X)))
function _pairwise(k::OuterKernel, X::AV, X′::AV)
    return map(k.f, X) .* pairwise(k.k, X, X′) .* map(k.f, X′)'
end
eachindex(k::OuterKernel) = eachindex(k.k)

############################## Convenience functionality ##############################

import Base: +, *, promote_rule, convert

function *(f, k::CrossKernel)
    return iszero(k) ? k : LhsCross(f, k)
end
function *(k::CrossKernel, f)
    return iszero(k) ? k : RhsCross(k, f)
end
function *(f, k::RhsCross)
    if k.k isa Kernel && f == k.f
        return OuterKernel(f, k.k)
    else
        return LhsCross(f, k)
    end
end
function *(k::LhsCross, f)
    if k.k isa Kernel && f == k.f
        return OuterKernel(f, k.k)
    else
        return RhsCross(k, f)
    end
end

promote_rule(::Type{<:MeanFunction}, ::Type{<:Real}) = MeanFunction
convert(::Type{MeanFunction}, x::Real) = ConstantMean(x)

promote_rule(::Type{<:Kernel}, ::Type{<:Real}) = Kernel
convert(::Type{<:CrossKernel}, x::Real) = ConstantKernel(x)

# Composing mean functions with Reals.
+(μ::MeanFunction, μ′::Real) = +(promote(μ, μ′)...)
+(μ::Real, μ′::MeanFunction) = +(promote(μ, μ′)...)

*(μ::MeanFunction, μ′::Real) = *(promote(μ, μ′)...)
*(μ::Real, μ′::MeanFunction) = *(promote(μ, μ′)...)

# Composing kernels with Reals.
+(k::CrossKernel, k′::Real) = +(promote(k, k′)...)
+(k::Real, k′::CrossKernel) = +(promote(k, k′)...)

*(k::CrossKernel, k′::Real) = *(promote(k, k′)...)
*(k::Real, k′::CrossKernel) = *(promote(k, k′)...)


text:
# Test whether the numerically computed recursion coefficients are correct
using PolyChaos, Test
import LinearAlgebra: norm
# load config
myfile = open("dataRecCoeffs/config.txt")
ns = parse.(Int,readlines(myfile))

nodes = ns[1]:ns[2]:ns[3]
mus = 0.:0.1:1
albe = 0.:0.2:2
tol = 1e-9

@time for n in nodes
        @testset "Hermite $n nodes" begin
            for m in mus
                myfile = open("dataRecCoeffs/hermite$(n)mu$m.txt")
                αβref = parse.(Float64, readlines(myfile))
                αβcom = rm_hermite(n, m)
                @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                αβcom = coeffs(genHermiteOrthoPoly(n-1, m; addQuadrature=false))
                # αβcom = coeffs(OrthoPoly("genhermite",n-1,Dict(:mu=>m)))
                @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
                close(myfile)
            end
        end
end

@time for n in nodes
        @testset "Logistic $n nodes" begin
            myfile = open("dataRecCoeffs/log$(n).txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = rm_logistic(n)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            αβcom = coeffs(LogisticOrthoPoly(n-1, addQuadrature=false))
            @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
            close(myfile)
    end
end

@time for n in nodes
        @testset "Jacobi $n nodes" begin
            for al in albe, be in albe
                    myfile = open("dataRecCoeffs/jac$(n)al$(al)be$(be).txt")
                    αβref = parse.(Float64,readlines(myfile))
                    αβcom = rm_jacobi(n,al,be)
                    @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                    αβcom = coeffs(JacobiOrthoPoly(n-1, al, be, addQuadrature=false))
                    @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
                    @test rm_jacobi(n,al) == rm_jacobi(n,al,al)
                    close(myfile)
            end
    end
end

@time for n in nodes
        @testset "Jacobi01 $n nodes" begin
            for al in albe, be in albe
                    myfile = open("dataRecCoeffs/jac01$(n)al$(al)be$(be).txt")
                    αβref = parse.(Float64,readlines(myfile))
                    αβcom = rm_jacobi01(n,al,be)
                    @test rm_jacobi01(n,al) == rm_jacobi01(n,al,al)
                    @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                    close(myfile)
            end
    end
end

@time for n in nodes
        @testset "Laguerre $n nodes" begin
            for m in mus
            myfile = open("dataRecCoeffs/laguerre$(n)a$m.txt")
            αβref = parse.(Float64,readlines(myfile))
            αβcom = rm_laguerre(n,m)
            @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
            αβcom = coeffs(genLaguerreOrthoPoly(n-1, m, addQuadrature=false))
            @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
            close(myfile)
        end
    end
end

meipol = 0.1:0.2:2
@time for n in nodes
        @testset "meixner_pollaczek $n nodes" begin
            for lambda in meipol, phi in meipol
                myfile = open("dataRecCoeffs/meixpol$(n)la$(lambda)phi$(phi).txt")
                αβref = parse.(Float64,readlines(myfile))
                αβcom = rm_meixner_pollaczek(n,lambda,phi)
                @test isapprox(norm(αβref-[αβcom[1];αβcom[2]],Inf),0.;atol=tol)
                αβcom = coeffs(MeixnerPollaczekOrthoPoly(n-1, lambda, phi, addQuadrature=false))
                @test isapprox(norm(αβref-[αβcom[:,1];αβcom[:,2]],Inf),0.;atol=tol)
                close(myfile)
                @test rm_meixner_pollaczek(n,lambda) == rm_meixner_pollaczek(n,lambda,pi/2)
                @test rm_meixner_pollaczek(n,lambda,pi/2)[1] == zeros(n)
            end

    end
end


text:
"""
```julia
lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3) -> ds
```
```math
\\begin{aligned}
\\dot{X} &= \\sigma(Y-X) \\\\
\\dot{Y} &= -XZ + \\rho X -Y \\\\
\\dot{Z} &= XY - \\beta Z
\\end{aligned}
```
The famous three dimensional system due to Lorenz [1], shown to exhibit
so-called "deterministic nonperiodic flow". It was originally invented to study a
simplified form of atmospheric convection.

Currently, it is most famous for its strange attractor (occuring at the default
parameters), which resembles a butterfly. For the same reason it is
also associated with the term "butterfly effect" (a term which Lorenz himself disliked)
even though the effect applies generally to dynamical systems.
Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : E. N. Lorenz, J. atmos. Sci. **20**, pp 130 (1963)
"""
function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(loop, u0, [σ, ρ, β], loop_jac)
end
function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function loop_jac(u, p, t)
    @inbounds begin
        σ, ρ, β = p
        J = @SMatrix [-σ  σ  0;
        ρ - u[3]  (-1)  (-u[1]);
        u[2]   u[1]  -β]
        return J
    end
end

function lorenz_iip(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(liip, u0, [σ, ρ, β], liip_jac)
end
function liip(du, u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
        return nothing
    end
end
function liip_jac(J, u, p, t)
    @inbounds begin
    σ, ρ, β = p
    J[1,1] = -σ; J[1, 2] = σ; J[1,3] = 0
    J[2,1] = ρ - u[3]; J[2,2] = -1; J[2,3] = -u[1]
    J[3,1] = u[2]; J[3,2] = u[1]; J[3,3] = -β
    return nothing
    end
end

"""
```julia
chua(u0 = [0.7, 0.0, 0.0]; a = 15.6, b = 25.58, m0 = -8/7, m1 = -5/7)
```
```math
\\begin{aligned}
\\dot{x} &= \\alpha (y - h(x))\\\\
\\dot{y} &= x - y+z \\\\
\\dot{z} &= \\beta y
\\end{aligned}
```
where h(x) is defined by
```math
h(x) = m_1 x + \\frac 1 2 (m_0 - m_1)(|x + 1| - |x - 1|)
```
This is a 3D continuous system that exhibits chaos.

Chua designed an electronic circuit with the expressed goal of exhibiting
chaotic motion, and this system is obtained by rescaling the circuit units
to simplify the form of the equation. [1]

The parameters are a, b, m0 and m1. Setting a = 15.6, m0 = -8/7 and m1 = -5/7,
and varying the parameter b from b = 25 to b = 51, one observes a classic
period-doubling bifurcation route to chaos. [2]

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : Chua, Leon O. "The genesis of Chua's circuit". Berkeley, CA, USA: Electronics Research Laboratory, College of Engineering, University of California, 1992.

[2] : [Leon O. Chua (2007) "Chua circuit", Scholarpedia, 2(10):1488.](http://www.scholarpedia.org/article/Chua_circuit)

"""
function chua(u0 = [0.7, 0.0, 0.0]; a = 15.6, b = 25.58, m0 = -8/7, m1 = -5/7)
    return CDS(chua_eom, u0, [a, b, m0, m1], chua_jacob)
end
function chua_eom(u, p, t)
    @inbounds begin
    a, b, m0, m1 = p
    du1 = a * (u[2] - u[1] - chua_element(u[1], m0, m1))
    du2 = u[1] - u[2] + u[3]
    du3 = -b * u[2]
    return SVector{3}(du1, du2, du3)
    end
end
function chua_jacob(u, p, t)
    a, b, m0, m1 = p
    return @SMatrix[-a*(1 + chua_element_derivative(u[1], m0, m1)) a 0;
                    1 -1 1;
                    0 -b 0]
end
# Helper functions for Chua's circuit.
function chua_element(x, m0, m1)
    return m1 * x + 0.5 * (m0 - m1) * (abs(x + 1.0) - abs(x - 1.0))
end
function chua_element_derivative(x, m0, m1)
    return m1 + 0.5 * (m0 - m1) * (-1 < x < 1 ? 2 : 0)
end

"""
```julia
roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
```
```math
\\begin{aligned}
\\dot{x} &= -y-z \\\\
\\dot{y} &= x+ay \\\\
\\dot{z} &= b + z(x-c)
\\end{aligned}
```
This three-dimensional continuous system is due to Rössler [1].
It is a system that by design behaves similarly
to the `lorenz` system and displays a (fractal)
strange attractor. However, it is easier to analyze qualitatively, as for example
the attractor is composed of a single manifold.
Default values are the same as the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : O. E. Rössler, Phys. Lett. **57A**, pp 397 (1976)
"""
function roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
    return CDS(roessler_eom, u0, [a, b, c], roessler_jacob)
end
function roessler_eom(u, p, t)
    @inbounds begin
    a, b, c = p
    du1 = -u[2]-u[3]
    du2 = u[1] + a*u[2]
    du3 = b + u[3]*(u[1] - c)
    return SVector{3}(du1, du2, du3)
    end
end
function roessler_jacob(u, p, t)
    a, b, c = p
    return @SMatrix [0.0 (-1.0) (-1.0);
                     1.0 a 0.0;
                     u[3] 0.0 (u[1]-c)]
end

"""
    double_pendulum(u0 = [π/2, 0, 0, rand()];
                    G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
Famous chaotic double pendulum system (also used for our logo!). Keywords
are gravity (G), lengths of each rod and mass of each ball (all assumed SI units).

The variables order is [θ1, dθ1/dt, θ2, dθ2/dt].

Jacobian is created automatically (thus methods that use the Jacobian will be slower)!

(please contribute the Jacobian and the e.o.m. in LaTeX :smile:)

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function double_pendulum(u0=[π/2, 0, 0, rand()]; G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
    return CDS(doublependulum_eom, u0, [G, L1, L2, M1, M2])
end
@inbounds function doublependulum_eom(state, p, t)
    G, L1, L2, M1, M2 = p

    du1 = state[2]
    del_ = state[3] - state[1]
    den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
    du2 = (M2*L1*state[2]*state[2]*sin(del_)*cos(del_) +
               M2*G*sin(state[3])*cos(del_) +
               M2*L2*state[4]*state[4]*sin(del_) -
               (M1 + M2)*G*sin(state[1]))/den1

    du3 = state[4]

    den2 = (L2/L1)*den1
    du4 = (-M2*L2*state[4]*state[4]*sin(del_)*cos(del_) +
               (M1 + M2)*G*sin(state[1])*cos(del_) -
               (M1 + M2)*L1*state[2]*state[2]*sin(del_) -
               (M1 + M2)*G*sin(state[3]))/den2
    return SVector{4}(du1, du2, du3, du4)
end

"""
    henonheiles(u0=[0, -0.25, 0.42081,0])
```math
\\begin{aligned}
\\dot{x} &= p_x \\\\
\\dot{y} &= p_y \\\\
\\dot{p}_x &= -x -2 xy \\\\
\\dot{p}_y &= -y - (x^2 - y^2)
\\end{aligned}
```

The Hénon–Heiles system [1] was introduced as a simplification of the motion
of a star around a galactic center. It was originally intended to study the
existence of a "third integral of motion" (which would make this 4D system integrable).
In that search, the authors encountered chaos, as the third integral existed
for only but a few initial conditions.

The default initial condition is a typical chaotic orbit.
The function `Systems.henonheiles_ics(E, n)` generates a grid of
`n×n` initial conditions, all having the same energy `E`.

[1] : Hénon, M. & Heiles, C., The Astronomical Journal **69**, pp 73–79 (1964)
"""
function henonheiles(u0=[0, -0.25, 0.42081, 0]#=; conserveE::Bool = true=#)
    i = one(eltype(u0))
    o = zero(eltype(u0))
    J = zeros(eltype(u0), 4, 4)
    return CDS(hheom!, u0, nothing, hhjacob!, J)
end
function hheom!(du, u, p, t)
    @inbounds begin
        du[1] = u[3]
        du[2] = u[4]
        du[3] = -u[1] - 2u[1]*u[2]
        du[4] = -u[2] - (u[1]^2 - u[2]^2)
        return nothing
    end
end
function hhjacob!(J, u, p, t)
    @inbounds begin
    o = 0.0; i = 1.0
    J[1,:] .= (o,    o,     i,    o)
    J[2,:] .= (o,    o,     o,    i)
    J[3,:] .= (-i - 2*u[2],   -2*u[1],   o,   o)
    J[4,:] .= (-2*u[1],  -1 + 2*u[2],  o,   o)
    return nothing
    end
end
_hhenergy(x,y,px,py) = 0.5(px^2 + py^2) + _hhpotential(x,y)
_hhpotential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
function henonheiles_ics(E, n=10)
    ys = range(-0.4, stop = 1.0, length = n)
    pys = range(-0.5, stop = 0.5, length = n)
    ics = Vector{Vector{Float64}}()
    for y in ys
        V = _hhpotential(0.0, y)
        V ≥ E && continue
        for py in pys
            Ky = 0.5*(py^2)
            Ky + V ≥ E && continue
            px = sqrt(2(E - V - Ky))
            ic = [0.0, y, px, py]
            push!(ics, [0.0, y, px, py])
        end
    end
    return ics
end

"""
    qbh([u0]; A=1.0, B=0.55, D=0.4)

```math
\\begin{aligned}
\\dot{q}_0 &= A p_0 \\\\
\\dot{q}_2 &= A p_2 \\\\
\\dot{p}_0 &= -A q_0 -3 \\frac{B}{\\sqrt{2}} (q_2^2 - q_1^2) - D q_1 (q_1^2 + q_2^2) \\\\
\\dot{p}_2 &= -q_2 (A + 3\\sqrt{2} B q_1 + D (q_1^2 + q_2^2)) (x^2 - y^2)
\\end{aligned}
```

These equations of motion correspond to a Hamiltonian used in nuclear
physics to study the quadrupole vibrations of the nuclear surface [1,2].

```math
H(p_0, p_2, q_0, q_2) = \\frac{A}{2}\\left(p_0^2+p_2^2\\right)+\\frac{A}{2}\\left(q_0^2+q_2^2\\right)
			 +\\frac{B}{\\sqrt{2}}q_0\\left(3q_2^2-q_0^2\\right) +\\frac{D}{4}\\left(q_0^2+q_2^2\\right)^2
```

The Hamiltonian has a similar structure with the Henon-Heiles one, but it has an added fourth order term
and presents a nontrivial dependence of chaoticity with the increase of energy [3].
The default initial condition is chaotic.

[1]: Eisenberg, J.M., & Greiner, W., Nuclear theory 2 rev ed. Netherlands: North-Holland pp 80 (1975)

[2]: Baran V. and Raduta A. A., International Journal of Modern Physics E, **7**, pp 527--551 (1998)

[3]: Micluta-Campeanu S., Raportaru M.C., Nicolin A.I., Baran V., Rom. Rep. Phys. **70**, pp 105 (2018)
"""
function qbh(u0=[0., -2.5830294658973876, 1.3873470962626937, -4.743416490252585];  A=1., B=0.55, D=0.4)
    return CDS(qeom, u0, [A, B, D])
end
function qeom(z, p, t)
    @inbounds begin
        A, B, D = p
        p₀, p₂ = z[1], z[2]
        q₀, q₂ = z[3], z[4]

        return SVector{4}(
            -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
            -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
            A * p₀,
            A * p₂
        )
    end
end

"""
    lorenz96(N::Int, u0 = rand(M); F=0.01)

```math
\\frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F
```

`N` is the chain length, `F` the forcing. Jacobian is created automatically.
(parameter container only contains `F`)
"""
function lorenz96(N::Int, u0 = rand(N); F=0.01)
    @assert N ≥ 4 "`N` must be at least 4"
    lor96 = Lorenz96{N}() # create struct
    return CDS(lor96, u0, [F])
end
struct Lorenz96{N} end # Structure for size type
function (obj::Lorenz96{N})(dx, x, p, t) where {N}
    F = p[1]
    # 3 edge cases
    @inbounds dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F
    @inbounds dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
    @inbounds dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F
    # then the general case
    for n in 3:(N - 1)
      @inbounds dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F
    end
    return nothing
end

"""
    duffing(u0 = [rand(), rand(), 0]; ω = 2.2, f = 27.0, d = 0.2, β = 1)
The (forced) duffing oscillator, that satisfies the equation
```math
\\ddot{x} + d\\cdot\\dot{x} + β*x + x^3 = f\\cos(\\omega t)
```
with `f, ω` the forcing strength and frequency and `d` the dampening.

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function duffing(u0 = [rand(), rand()]; ω = 2.2, f = 27.0, d = 0.2, β = 1)

    J = zeros(eltype(u0), 2, 2)
    J[1,2] = 1
    return CDS(duffing_eom, u0, [ω, f, d, β], duffing_jacob)
end
@inbounds function duffing_eom(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector{2, Float64}(dx1, dx2)
end
@inbounds function duffing_jacob(u, p, t)
    ω, f, d, β = p
    return @SMatrix [0 1 ;
    (-β - 3u[1]^2) -d]
end

"""
    shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
Shinriki oscillator with all other parameters (besides `R1`) set to constants.
*This is a stiff problem, be careful when choosing solvers and tolerances*.
"""
function shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
    # # Jacobian caller for Shinriki:
    # shinriki_eom(::Type{Val{:jac}}, J, u, p, t) = (shi::Shinriki)(t, u, J)
    return CDS(shinriki_eom, u0, [R1])
end
shinriki_voltage(V) = 2.295e-5*(exp(3.0038*V) - exp(-3.0038*V))
function shinriki_eom(u, p, t)
    R1 = p[1]

    du1 = (1/0.01)*(
    u[1]*(1/6.9 - 1/R1) - shinriki_voltage(u[1] - u[2]) - (u[1] - u[2])/14.5
    )

    du2 = (1/0.1)*(
    shinriki_voltage(u[1] - u[2]) + (u[1] - u[2])/14.5 - u[3]
    )

    du3 = (1/0.32)*(-u[3]*0.1 + u[2])
    return SVector{3}(du1, du2, du3)
end

"""
```julia
gissinger(u0 = 3rand(3); μ = 0.119, ν = 0.1, Γ = 0.9)
```
```math
\\begin{aligned}
\\dot{Q} &= \\mu Q - VD \\\\
\\dot{D} &= -\\nu D + VQ \\\\
\\dot{V} &= \\Gamma -V + QD
\\end{aligned}
```
A continuous system that models chaotic reversals due to Gissinger [1], applied
to study the reversals of the magnetic field of the Earth.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : C. Gissinger, Eur. Phys. J. B **85**, 4, pp 1-12 (2012)
"""
function gissinger(u0 = 3rand(3); μ = 0.119, ν = 0.1, Γ = 0.9)
    return CDS(gissinger_eom, u0, [μ, ν, Γ], gissinger_jacob)
end
function gissinger_eom(u, p, t)
    μ, ν, Γ = p
    du1 = μ*u[1] - u[2]*u[3]
    du2 = -ν*u[2] + u[1]*u[3]
    du3 = Γ - u[3] + u[1]*u[2]
    return SVector{3}(du1, du2, du3)
end
function gissinger_jacob(u, p, t)
    μ, ν, Γ = p
    return @SMatrix [μ -u[3] -u[2];
                     u[3] -ν u[1];
                     u[2] u[1] -1]
end

"""
```julia
rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
```
```math
\\begin{aligned}
\\dot{x} &= -\\mu x +yz \\\\
\\dot{y} &= -\\mu y +x(z-\\alpha) \\\\
\\dot{z} &= 1 - xz
\\end{aligned}
```
Rikitake's dynamo is a system that tries to model the magnetic reversal events
by means of a double-disk dynamo system.

[1] : T. Rikitake Math. Proc. Camb. Phil. Soc. **54**, pp 89–105, (1958)
"""
function rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
    return CDS(rikitake_eom, u0, [μ, α], rikitake_jacob)
end
function rikitake_eom(u, p, t)
    μ, α = p
    x,y,z = u
    xdot = -μ*x + y*z
    ydot = -μ*y + x*(z - α)
    zdot = 1 - x*y
    return SVector{3}(xdot, ydot, zdot)
end
function rikitake_jacob(u, p, t)
    μ, α = p
    x,y,z = u
    xdot = -μ*x + y*z
    ydot = -μ*y + x*(z - α)
    zdot = 1 - x*y
    return @SMatrix [-μ  z  y;
                     z-α -μ x;
                     -y  -x 0]
end

"""
```julia
nosehoover(u0 = [0, 0.1, 0])
```
```math
\\begin{aligned}
\\dot{x} &= y \\\\
\\dot{y} &= yz - x \\\\
\\dot{z} &= 1 - y^2
\\end{aligned}
```
Three dimensional conservative continuous system, discovered in 1984 during
investigations in thermodynamical chemistry by Nosé and Hoover, then
rediscovered by Sprott during an exhaustive search as an extremely simple
chaotic system. [1]

See Chapter 4 of "Elegant Chaos" by J. C. Sprott. [2]

[1] : Hoover, W. G. (1995). Remark on ‘‘Some simple chaotic flows’’. *Physical Review E*, *51*(1), 759.

[2] : Sprott, J. C. (2010). *Elegant chaos: algebraically simple chaotic flows*. World Scientific.
"""
nosehoover(u0 = [0, 0.1, 0]) = CDS(nosehoover_eom, u0, nothing, nosehoover_jacob)
function nosehoover_eom(u, p, t)
    x,y,z = u
    xdot = y
    ydot = y*z - x
    zdot  = 1.0 - y*y
    return SVector{3}(xdot, ydot, zdot)
end
function nosehoover_jacob(u, p, t)
    x,y,z = u
    return @SMatrix [0 1 0;
                     -1 z y;
                     0 2y 0]
end

"""
```julia
labyrinth(u0 = [1.0, 0, 0])
```
```math
\\begin{aligned}
\\dot{x} &= \\sin(y) \\\\
\\dot{y} &= \\sin(z) \\\\
\\dot{V} &= \\sin(x)
\\end{aligned}
```
Three dimensional conservative continuous system, whose evolution in 3D space
looks like a speudo-random walk, the orbit moving around like in a labyrinth.

First proposed by René Thomas (1999). [1] See discussion in Section 4.4.3 of
"Elegant Chaos" by J. C. Sprott. [2]

[1] : Thomas, R. (1999). Deterministic chaos seen in terms of feedback circuits: Analysis, synthesis," labyrinth chaos". *International Journal of Bifurcation and Chaos*, *9*(10), 1889-1905.

[2] : Sprott, J. C. (2010). *Elegant chaos: algebraically simple chaotic flows*. World Scientific.
"""
labyrinth(u0 = [1.0, 0, 0]) = CDS(labyrinth_eom, u0, nothing, labyrinth_jacob)
function labyrinth_eom(u, p, t)
    x,y,z = u
    xdot = sin(y)
    ydot = sin(z)
    zdot = sin(x)
    return SVector{3}(xdot, ydot, zdot)
end
function labyrinth_jacob(u, p, t)
    x,y,z = u
    return @SMatrix [0 cos(y) 0;
                     0 0 cos(z);
                     cos(x) 0 0]
end

"""
    antidots(u; B = 1.0, d0 = 0.3, c = 0.2)
An antidot "superlattice" is a Hamiltonian system that corresponds to a
smoothened periodic Sinai billiard with disk diameter `d0` and smooth
factor `c` [1].

This version is the two dimensional
classical form of the system, with quadratic equations of motion and
a perpendicular magnetic field. Notice that the equations of motion
are with respect to the velocity instead of momentum, i.e.:
```math
\\begin{aligned}
\\dot{x} &= v_x \\\\
\\dot{y} &= v_y \\\\
\\dot{v_x} &= B*v_y - U_x \\\\
\\dot{v_y} &= -B*v_x - U_X \\\\
\\end{aligned}
```
with ``U`` the potential energy:
```math
U = \\left(\\tfrac{1}{c^4}\\right) \\left[\\tfrac{d_0}{2} + c - r_a\\right]^4
```
if ``r_a = \\sqrt{(x%1)^2 + (y%1)^2} < \\frac{d_0}{2} + c`` and 0
otherwise. I.e. the potential is periodic with period 1 in both ``x, y`` and
normalized such that for energy value of 1 it is a circle of diameter ``d0``.
The magnetic field is also normalized such that for value `B=1` the cyclotron
diameter is 1.

Fo more details see [1].

[1] : G. Datseris *et al*, [arXiv:1711.05833v3](https://arxiv.org/abs/1711.05833v3)
"""
function antidots(u0 = [0.5, 0.5, rand(2)...];
    d0 = 0.5, c = 0.2, B = 1.0)
    return CDS(antidot_eom, u0, [B, d0, c], antidot_jacob)
end

function antidot_eom(u, p, t)
    B, d0, c = p
    x, y, vx, vy = u
    # Calculate quadrant of (x,y):
    U = Uy = Ux = 0.0
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Calculate derivatives and potential:
    if ρ < 0.5*d0 + c
        sharedfactor = -(4*(c + d0/2 - ρ)^(3))/(c^4*ρ)
        Ux = sharedfactor*xtilde # derivatives
        Uy = sharedfactor*ytilde
    end
    Br = 2*√(2)*B # magnetic field with prefactor
    return SVector{4}(vx, vy, Br*vy - Ux, -vx*Br - Uy)
end
function antidot_jacob(u, p, t)
    B, d0, c = p
    x, y, vx, vy = u
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Calculate derivatives and potential:
    if ρ < 0.5*d0 + c
        Uxx, Uyy, Uxy = antidot_secondderv(xtilde, ytilde, p)
    else
        Uxx, Uyy, Uxy = 0.0, 0.0, 0.0
    end
    Br = 2*√(2)*B # magnetic field with prefactor
    return SMatrix{4, 4}(0.0, 0, -Uxx, -Uxy, 0, 0, -Uxy, -Uyy,
                         1, 0, 0, -Br,        0, 1, Br, 0)
end

function antidot_potential(x::Real, y::Real, p)
    B, d0, c = p
    δ = 4
    # Calculate quadrant of (x,y):
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Check distance:
    pot = ρ > 0.5*d0 + c ? 0.0 : (1/(c^δ))*(0.5*d0 + c - ρ)^δ
end

function antidot_secondderv(x, y, p)
    B, d0, c = p
    r = sqrt(x^2 + y^2)
    Uxx = _Uxx(x, y, c, d0, r)
    Uyy = _Uxx(y, x, c, d0, r)
    Uxy =  (2c + d0 - 2r)^2 * (2c + d0 + 4r)*x*y / (2c^4*r^3)
    return Uxx, Uyy, Uxy
end

function _Uxx(x, y, c, d0, r)
    Uxx =  (2c + d0 - 2r)^2 * (-2c*y^2 - d0*y^2 + 2*r*(3x^2 + y^2)) /
    (2 * (c^4) * r^3)
end

"""
```julia
ueda(u0 = [3.0, 0]; k = 0.1, B = 12.0)
```
```math
\\ddot{x} + k \\dot{x} + x^3 = B\\cos{t}
```
Nonautonomous Duffing-like forced oscillation system, discovered by Ueda in
1961. It is one of the first chaotic systems to be discovered.

The stroboscopic plot in the (x, ̇x) plane with period 2π creates a "broken-egg
attractor" for k = 0.1 and B = 12. Figure 5 of [1] is reproduced by

```julia
using Plots
ds = Systems.ueda()
a = trajectory(ds, 2π*5e3, dt = 2π)
scatter(a[:, 1], a[:, 2], markersize = 0.5, markercolor=:black, leg=false, title="Ueda attractor")
```

For more forced oscillation systems, see Chapter 2 of "Elegant Chaos" by
J. C. Sprott. [2]

[1] : [Ruelle, David, ‘Strange Attractors’, The Mathematical Intelligencer, 2.3 (1980), 126–37](https://doi.org/10/dkfd3n)

[2] : Sprott, J. C. (2010). *Elegant chaos: algebraically simple chaotic flows*. World Scientific.
"""
function ueda(u0 = [3.0, 0]; k = 0.1, B = 12.0)
    return CDS(ueda_eom, u0, [k, B], ueda_jacob)
end
function ueda_eom(u, p, t)
    x,y = u
    k, B = p
    xdot = y
    ydot = B*cos(t) - k*y - x^3
    return SVector{2}(xdot, ydot)
end
function ueda_jacob(u, p, t)
    x,y = u
    k, B = p
    return @SMatrix [0      1;
                     -3*x^2 -k]
end


text:
module SynthControl

using Dates, DataFrames, Optim, LinearAlgebra, RecipesBase, Parameters
import CSV

export TreatmentPanel, SynthControlModel, fit!, isfitted, load_brexit, load_germany, load_basque, load_smoking

include("treatmentpanel.jl")
include("synthcontrolmodels.jl")
include("plotrecipes.jl")
include("load_data.jl")

end # module


text:
using ApproxFun, Test
import ApproxFun: expα, expβ, expγ

#
# This tests the numerical stability of the formulæ appearing in Eq. (2.5) of:
#
# A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput., 26:1214--1233, 2005.
#

@testset "ETDRK4" begin
    x = 10 .^ range(-15, stop=2, length=18)

    @test norm(expα.(x)./expα.(big.(x)).-1,Inf) < eps()
    @test norm(expβ.(x)./expβ.(big.(x)).-1,Inf) < eps()
    @test norm(expγ.(x)./expγ.(big.(x)).-1,Inf) < eps()
end


text:
@testset "Explicit Scitype tests" begin
    STB = ScientificTypesBase
    ST = ScientificTypes
    @test STB.Scitype(Int, ST.CONV) == Count
    @test STB.Scitype(Float64, ST.CONV) == Continuous
    @test STB.Scitype(Float32, ST.CONV) == Continuous
    @test STB.Scitype(Float16, ST.CONV) == Continuous
    @test STB.Scitype(SubString, ST.CONV) == Textual
end

@testset "scitype of (In)Finite elements" begin
    cv = categorical(['x', 'y'])
    c  = cv[1]
    uv = categorical(['x', 'y'], ordered=true)
    u  = uv[1]
    @test scitype((4, 4.5, c, u, "X")) == Tuple{
                       Count,Continuous,Multiclass{2},OrderedFactor{2},Textual}
end

@testset "elscitype" begin
    X = randn(5, 5)
    @test scitype(X) == Arr{Continuous, 2}
    @test elscitype(X) == Continuous
    X = [1, 2, missing, 5]
    @test elscitype(X) == Union{Missing, Count}
    @test elscitype(X[1:2]) == Union{Missing, Count}
    X = categorical([1, 2, 3, missing, 2, 1, missing])
    @test elscitype(X) == Union{Missing,Multiclass{3}}
    X = categorical("lksfjalksjdflkjsdlkfjasldkfj" |> collect)
    @test elscitype(X) <: Multiclass
    X = coerce(categorical([1,2,3,1,2,3]), OrderedFactor)
    @test elscitype(X) <: Union{Missing, OrderedFactor}
end

A = Any[2 4.5;
        4 4.5;
        6 4.5]

@testset "scitype_union" begin
    @test scitype_union(A) == Union{Count,Continuous}
    @test scitype_union(randn(1000000)) == Continuous
    @test scitype_union(1) == Count
    @test scitype_union([1]) == Count
    @test scitype_union(Any[1]) == Count
    @test scitype_union([1, 2.0, "3"]) == Union{Continuous, Count, Textual}
end

@testset "scitype of Arrays" begin
    # COUNT
    @test scitype(A)               == Arr{Union{Count,Continuous}, 2}
    @test scitype([1, 2, 3])       == Vec{Count}
    @test scitype([1, missing, 3]) == Vec{Union{Missing,Count}}
    @test scitype(Any[1, 2, 3])    == Vec{Count}
    @test scitype(Any[1, missing]) == Vec{Union{Missing,Count}}

    # CONTINUOUS
    @test scitype([1.0, 2.0, 3.0])    == Vec{Continuous}
    @test scitype(Any[1.0, missing])  == Vec{Union{Missing,Continuous}}
    @test scitype(Any[1.0, 2.0, 3.0]) == Vec{Continuous}
    @test scitype(Any[1.0, missing])  == Vec{Union{Missing,Continuous}}

    # MULTICLASS
    @test scitype(categorical(1:4))         == Vec{Multiclass{4}}
    @test scitype(view(categorical(1:4), 1:3)) == Vec{Multiclass{4}}
    @test scitype(Any[categorical(1:4)...]) == Vec{Multiclass{4}}
    @test scitype(categorical([1, missing, 3])) ==
        Vec{Union{Multiclass{2},Missing}}

        a = ["aa", "bb", "aa", "bb"] |> categorical
    @test scitype(a[1]) == Multiclass{2}

    # NOTE: the slice here does not contain missings but the machine type
    # still contains a missing so the scitype remains with a missing
    @test scitype(categorical([1, missing, 3])[1:1]) ==
        Vec{Union{Multiclass{2},Missing}}

    # ORDEREDFACTOR
    @test scitype(categorical(1:4, ordered=true)) ==
        Vec{OrderedFactor{4}}
    @test scitype(Any[categorical(1:4, ordered=true)...]) ==
        Vec{OrderedFactor{4}}
    @test scitype(categorical([1, missing, 3], ordered=true)) ==
        Vec{Union{OrderedFactor{2},Missing}}
    # NOTE: see note for multiclass
    @test scitype(categorical([1, missing, 3], ordered=true)[1:1]) ==
        Vec{Union{OrderedFactor{2},Missing}}
end

@testset "scitype of Images" begin
    black = RGB(0, 0, 0)
    color_image = fill(black, (10, 20))
    @test scitype(color_image) == ColorImage{10,20}

    color_image2 = fill(black, (5, 3))
    v = [color_image, color_image2, color_image2]
    @test scitype(v) ==
        Vec{Union{ColorImage{10,20},ColorImage{5,3}}}

    white = Gray(1.0)
    gray_image = fill(white, (10, 20))
    @test scitype(gray_image) == GrayImage{10,20}
end

# @testset "PersistenceDiagrams" begin
#     diagram = PersistenceDiagram([(1, Inf), (2, 3)], dim=0)
#     @test scitype(diagram) == PersistenceDiagram

#     diagrams = [diagram, diagram, diagram]
#     @test scitype(diagrams) == Vec{PersistenceDiagram}
# end

@testset "scitype for objects with temporal types" begin
    d = Date(2020, 4, 21)
    t = Time(8, 15, 42)
    dt = now()
    @test scitype(d) == ScientificDate
    @test scitype(t) == ScientificTime
    @test scitype(dt) == ScientificDateTime
    @test scitype(vcat([missing, ], fill(d, 3))) ==
        AbstractVector{Union{Missing,ScientificDate}}
    @test scitype(vcat([missing, ], fill(t, 3))) ==
        AbstractVector{Union{Missing,ScientificTime}}
    @test scitype(vcat([missing, ], fill(dt, 3))) ==
        AbstractVector{Union{Missing,ScientificDateTime}}
    @test scitype([missing, d, t]) ==
        AbstractVector{Union{Missing,ScientificTimeType}}
end

struct FooSampleable <: Dist.Sampleable{Dist.ArrayLikeVariate{0},
    Dist.Discrete}
end

@testset "scitypes of distributions" begin
    @test scitype(Dist.Normal()) == Density{Continuous}
    @test scitype(Dist.Poisson()) == Density{Count}
    @test scitype(Dist.Categorical(3)) == Density{Count}
    X = [1 2; 3 4]
    @test scitype(Dist.MultivariateNormal(X'X)) ==
    Density{AbstractVector{Continuous}}
    @test scitype(FooSampleable()) == Sampleable{Count}
end

@testset "Tables scitype" begin
    ST = ScientificTypes
    ## Test scitype for Column oriented tables
    X = (
        x = rand(5),
        y = rand(Int, 5),
        z = categorical(collect("asdfa")),
        w = rand(5)
    )
    t = scitype(X)
    @test t <: Table(Continuous, Finite, Count)
    @test t <: Table(Infinite, Multiclass)
    @test !(t <: Table(Continuous, Union{Missing, Count}))
    @test ST._cols_scitype(
        Tables.columns(X), Tables.Schema(Tables.columnnames(X), nothing)
    ) == t

    ## scitype tests for row oriented tables.
    Y = Tables.rowtable(X)
    r = scitype(Y)
    @test r <: Table(Continuous, Finite, Count)
    @test r <: Table(Infinite, Multiclass)
    @test !(r <: Table(Continuous, Union{Missing, Count}))
    rows = Tables.rows(Y)
    @test ST._rows_scitype(
        rows, Tables.Schema(Tables.columnnames(iterate(rows, 1)[1]), nothing)
    ) == r
    # ExtremelyWide row oriented table
    @test ST._rows_scitype(
        rows,
        Tables.Schema(
            Tables.columnnames(iterate(rows, 1)[1]),
            (Int, Int, CategoricalValue{Char, UInt32}, Float64);
            stored = true
        )
    ) == r

    # test schema for column oreinted tables with number of columns
    # exceeding COLS_SPECIALIZATION_THRESHOLD.
    nt = NamedTuple{
            Tuple(Symbol("x$i") for i in Base.OneTo(ST.COLS_SPECIALIZATION_THRESHOLD + 1))
        }(
            Tuple(rand(2) for i in Base.OneTo(ST.COLS_SPECIALIZATION_THRESHOLD + 1))
    )
    @test scitype(nt) <: Table(Continuous)

    #issue 146
    X = Tables.table(coerce(rand("abc", 5, 3), Multiclass))
    @test scitype(X) === Table{AbstractVector{Multiclass{3}}}
end

# TODO: re-instate when julia 1.0 is no longer LTS release:

# @testset "csvfile" begin
#     X = (x = rand(4), )
#     CSV.write("test.csv", X)
#     file = CSV.File("test.csv")
#     @test scitype(file) == scitype(X)
#     rm("test.csv")
# end


text:
using Test

prec=12

@testset "faststorage" begin
    for i in 1:100000
        Kf = rand(1)[1]
        Storage = rand(0.1:0.1:20)[1]
        Overland = rand(0:0.1:10)[1]
        # test that output 0, if input 0
        @test faststorage(0, 0, Kf) == (0,0)
        # test that storage decreases if no input
        Discharge, Faststorage = faststorage(0, Storage, Kf)
        @test -eps(Float64)*10 <= Storage - Faststorage - Discharge <= eps(Float64)*10
        # test that discharge should be the Kf-ratio of sum of storage and overlandflow
        Discharge, Faststorage = faststorage(Overland, Storage, Kf)
        @test (Storage + Overland) * Kf == Discharge
        @test Storage + Overland - (Storage + Overland) * Kf == Faststorage
        @test Storage + Overland - Discharge == Faststorage
        # test right behavior if storage = 0
        Discharge, Faststorage = faststorage(Overland, 0, Kf)
        @test Overland * Kf == Discharge
        @test round(Discharge + Faststorage, digits=prec) == round(Overland, digits=prec)
    end
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
#=
   Basic model as an template for other models...
=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function alphametic(problem_str="SEND+MORE=MONEY", base=10, print_solutions=true,all_solutions=false,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[], # :true,
                                                            # "logging"=> :true,

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    problem = split(problem_str,r"[+=-]").|>string
    println("problem: $problem base:$base")
    p_len = length(problem)

    # create the lookup table of (digit -> ix)
    a = join(problem).|>unique
    n = length(a)
    lookup = Dict(zip(a,1:n))
    # println("lookup:",lookup)

    # length of each number
    lens = problem.|>length

    @variable(model, 0 <= x[1:n] <= base-1, Int)
    sums = [@variable(model, integer=true,lower_bound=base^(lens[i]-1),upper_bound=(base^lens[i])-1) for i in 1:length(lens)]
    # @variable(model, 1 <= sums[1:length(lens)] <= base^(maximum(lens))-1, Int)

    @constraint(model, x in CS.AllDifferent())

    ix = 1
    for prob in problem

      # sum all the digits with proper exponents to a number
      this_len = length(prob)
      ss = [(base^i) * x[lookup[prob[this_len - i]]] for i in (this_len-1):-1:0]
      @constraint(model,sums[ix] == sum(ss))

      # leading digits must be > 0
      @constraint(model,x[lookup[prob[1]]] > 0)

      ix += 1

    end

    # the last number (RHS) is the sum of the previous numbers (LHS)
    @constraint(model, sum([sums[i] for i in 1:p_len-1]) == sums[end])

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")
                for c in a
                  println("$c: $(x_val[lookup[c]])")
                end
            end
        end
        println()
    else
        println("status:$status")
    end

    return status, num_sols
end

problems = [
  "SEND+MORE=MONEY",
  "SEND+MOST=MONEY", "VINGT+CINQ+CINQ=TRENTE",
  "EIN+EIN+EIN+EIN=VIER",
  "DONALD+GERALD=ROBERT",
  "SATURN+URANUS+NEPTUNE+PLUTO+PLANETS",  # this is hard
  "WRONG+WRONG=RIGHT"
]

function run_problems(base=10)
  for problem in problems
    # @time alphametic(problem,10,true,true,6)
    @time alphametic(problem,base,true,false,116)
    println()
  end

end

#=
if length(ARGS) > 0
  problem = ARGS[1]
  base = 10
  if length(ARGS) > 1
    base = parse(Int,ARGS[2])
  end
  if problem == "TEST" || problem == "test"
    run_problems(base)
  else
    alphametic(problem)
  end
else
  run_problems()
end
=#

run_problems()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export i_leapfrog!, i_leapfrog_reuse!

@i @inline function acceleration_reuse(y!::V3{T}, ra::V3{T}, rb::V3{T}, mb::Real, G) where T
    @routine @invcheckoff begin
        @zeros T d anc1 anc2 anc3 anc4
        d += sqdistance(ra, rb)
        anc1 += sqrt(d)
        anc2 += anc1 ^ 3
        anc3 += G * mb
        anc4 += anc3 / anc2
        rb -= ra
    end
    y! += anc4 * rb
    ~@routine
end

@i @inline function :(+=)(acceleration)(y!::V3{T}, ra::V3{T}, rb::V3{T}, mb::Real, G) where T
    @routine @invcheckoff begin
        @zeros T d anc1 anc2 anc3 anc4
        rc ← zero(V3{T})
        d += sqdistance(ra, rb)
        anc1 += sqrt(d)
        anc2 += anc1 ^ 3
        anc3 += G * mb
        anc4 += anc3 / anc2
        rc += rb - ra
    end
    y! += anc4 * rc
    ~@routine
end

@i function i_update_acceleration!(a, r, m, G)
    @safe @assert length(a) == length(m) == length(r)
    @inbounds for j=1:length(a)
        for k=1:length(a)
            if j != k
                a[j] += acceleration(r[j], r[k], m[k], G)
            end
        end
    end
end

@i function i_leapfrog!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector{T}; G=G_year_AU, n, dt) where T
    @invcheckoff for i=1:n
        i_leapfrog_step!(r, v, m; G=G, dt=dt)
    end
end

@i function i_leapfrog_step!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector; G=G_year_AU, dt) where T
    @routine @invcheckoff begin
        nplanets ← length(m)
        halfdt ← zero(dt)
        halfdt += dt/2
    end
    @inbounds @invcheckoff begin
        for j=1:nplanets
            r[j] += halfdt*v[j]
        end
        # compute acceleration and update velocity
        @routine begin
            a ← zeros(V3{T}, nplanets)
            i_update_acceleration!(a, r, m, G)
        end
        for j=1:nplanets
            v[j] += dt*a[j]
        end
        ~@routine

        for j=1:nplanets
            r[j] += halfdt*v[j]
        end
    end
    ~@routine
end

@i function i_leapfrog_reuse!(r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}}, m::AbstractVector{T}; G=G_year_AU, n, dt) where T
    @routine @invcheckoff begin
        nplanets ← length(m)
        halfdt ← zero(dt)
        halfdt += dt/2
        a ← zeros(V3{T}, nplanets)
    end
    @inbounds @invcheckoff for i=1:n
        # update position
        for j=1:nplanets
            r[j] += halfdt*v[j]
        end
        # compute acceleration
        @routine for j=1:nplanets
            for k=1:nplanets
                if j != k
                    acceleration_reuse(a[j], r[j], r[k], m[k], G)
                end
            end
        end

        # update velocity
        for j=1:nplanets
            v[j] += dt*a[j]
        end
        ~@routine

        # update position
        for j=1:nplanets
            r[j] += halfdt*v[j]
        end
    end
    ~@routine
end


text:
module TestDoctest

import FoldsThreads
using Documenter: doctest
using Test

@testset "doctest" begin
    doctest(FoldsThreads; manual = false)
end

end  # module


text:
### Compute empirical covariance matrix recursively

function covariance!(
  C::RealMatrix,
  lastC::RealMatrix,
  k::Integer,
  x::RealVector,
  lastmean::RealVector,
  secondlastmean::RealVector
)
  C[:, :] = (k-1)*lastC
  BLAS.ger!(1.0, x, x, C)
  BLAS.ger!(1.0, -(k+1)*lastmean, lastmean, C)
  BLAS.ger!(1.0, k*secondlastmean, secondlastmean, C)
  C[:, :] /= k
end


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

"""
Traverse the entries in a tree and its subtrees in post or pre order.

Function parameter should have following signature:

    (Cstring, Ptr{Void}, Ptr{Void}) -> Cint
"""
function treewalk(f::Function, tree::GitTree, payload=Any[], post::Bool = false)
    cbf = cfunction(f, Cint, Tuple{Cstring, Ptr{Void}, Ptr{Void}})
    cbf_payload = Ref{typeof(payload)}(payload)
    @check ccall((:git_tree_walk, :libgit2), Cint,
                  (Ptr{Void}, Cint, Ptr{Void}, Ptr{Void}),
                   tree.ptr, post, cbf, cbf_payload)
    return cbf_payload
end

repository(tree::GitTree) = tree.owner
repository(te::GitTreeEntry) = repository(te.owner)

function filename(te::GitTreeEntry)
    str = ccall((:git_tree_entry_name, :libgit2), Cstring, (Ptr{Void},), te.ptr)
    str != C_NULL && return unsafe_string(str)
    return nothing
end

function filemode(te::GitTreeEntry)
    return ccall((:git_tree_entry_filemode, :libgit2), Cint, (Ptr{Void},), te.ptr)
end

function entrytype(te::GitTreeEntry)
    otype = ccall((:git_tree_entry_type, :libgit2), Cint, (Ptr{Void},), te.ptr)
    return objtype(Consts.OBJECT(otype))
end

function entryid(te::GitTreeEntry)
    oid_ptr = ccall((:git_tree_entry_id, :libgit2), Ptr{UInt8}, (Ptr{Void},), te.ptr)
    return GitHash(oid_ptr)
end

function Base.count(tree::GitTree)
    return ccall((:git_tree_entrycount, :libgit2), Csize_t, (Ptr{Void},), tree.ptr)
end

function Base.getindex(tree::GitTree, i::Integer)
    if i < 1 || i > count(tree)
        throw(BoundsError(tree, i))
    end
    te_ptr = ccall((:git_tree_entry_byindex, :libgit2),
                   Ptr{Void},
                   (Ptr{Void}, Csize_t), tree.ptr, i-1)
    return GitTreeEntry(tree, te_ptr, false)
end

function (::Type{T})(te::GitTreeEntry) where T<:GitObject
    repo = repository(te)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_tree_entry_to_object, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ref{Void}),
                   obj_ptr_ptr, repo.ptr, te.ptr)
    return T(repo, obj_ptr_ptr[])
end

function Base.show(io::IO, te::GitTreeEntry)
    println(io, "GitTreeEntry:")
    println(io, "Entry name: ", filename(te))
    println(io, "Entry type: ", entrytype(te))
    println(io, "Entry OID: ", entryid(te))
end

function Base.show(io::IO, tree::GitTree)
    println(io, "GitTree:")
    println(io, "Owner: ", repository(tree))
    println(io, "Number of entries: ", count(tree))
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
const TYPE_0 = zero(UInt8)
const TYPE_1 = one(UInt8) << 5
const TYPE_2 = UInt8(2) << 5
const TYPE_3 = UInt8(3) << 5
const TYPE_4 = UInt8(4) << 5
const TYPE_5 = UInt8(5) << 5
const TYPE_6 = UInt8(6) << 5
const TYPE_7 = UInt8(7) << 5

const BITS_PER_BYTE = UInt8(8)
const HEX_BASE = Int(16)
const LOWEST_ORDER_BYTE_MASK = 0xFF

const TYPE_BITS_MASK = UInt8(0b1110_0000)
const ADDNTL_INFO_MASK = UInt8(0b0001_1111)

const ADDNTL_INFO_UINT8 = UInt8(24)
const ADDNTL_INFO_UINT16 = UInt8(25)
const ADDNTL_INFO_UINT32 = UInt8(26)
const ADDNTL_INFO_UINT64 = UInt8(27)

const SINGLE_BYTE_SIMPLE_PLUS_ONE = UInt8(24)
const SIMPLE_FALSE = UInt8(20)
const SIMPLE_TRUE = UInt8(21)
const SIMPLE_NULL = UInt8(22)
const SIMPLE_UNDEF = UInt8(23)

const ADDNTL_INFO_FLOAT16 = UInt8(25)
const ADDNTL_INFO_FLOAT32 = UInt8(26)
const ADDNTL_INFO_FLOAT64 = UInt8(27)

const ADDNTL_INFO_INDEF = UInt8(31)
const BREAK_INDEF = TYPE_7 | 31

const SINGLE_BYTE_UINT_PLUS_ONE = 24
const UINT8_MAX_PLUS_ONE = 0x100
const UINT16_MAX_PLUS_ONE = 0x10000
const UINT32_MAX_PLUS_ONE = 0x100000000
const UINT64_MAX_PLUS_ONE = 0x10000000000000000

const SIZE_OF_FLOAT64 = sizeof(Float64)
const SIZE_OF_FLOAT32 = sizeof(Float32)
const SIZE_OF_FLOAT16 = sizeof(Float16)

const POS_BIG_INT_TAG = UInt8(2)
const NEG_BIG_INT_TAG = UInt8(3)

const CBOR_FALSE_BYTE = UInt8(TYPE_7 | 20)
const CBOR_TRUE_BYTE = UInt8(TYPE_7 | 21)
const CBOR_NULL_BYTE = UInt8(TYPE_7 | 22)
const CBOR_UNDEF_BYTE = UInt8(TYPE_7 | 23)


text:
# Test gausslaguerre.jl

# Test integration
tol = 4.0e-10

# Evaluate the exact integral of x^α * p(x) *exp(-x) on the positive halfline,
# where p(x) = sum( coef[i+1]*x^i, i=0..degree(p)) is a polynomial given by its
# monomial coefficients
function exact_integral_poly(α, coef)
    # do calculations in BigFloat precision
    z = big(0.0)
    for i in 1:length(coef)
        z += gamma(big(α)+big(i)) * big(coef[i])
    end
    Float64(z)
end

function polyval(x, coef)
    z = zeros(size(x))
    d = length(coef)-1
    for i in 1:length(x)
        z[i] = sum(coef .* x[i].^(0:d))
    end
    z
end

# Evaluate the exact integral of cos(a*x) * x^α * exp(-x) on the positive halfline
# - For α = 0
exact_integral_cos1(a) = 1/(a^2+1)
# - For α = 1/2
exact_integral_cos2(a) = sqrt(π)*cos(3/2*atan(a))/2 / (a^2+1)^(3/4)

Random.seed!(0)

@testset "Gauss–Laguerre" begin
    # Check error
    @test_throws DomainError gausslaguerre(1, -1.4)
    @test_throws DomainError gausslaguerre(-1)

    # Check optional argument
    for n in 0:20
        @test gausslaguerre(n) == gausslaguerre(n,0) == gausslaguerre(n,0.0)
    end
    for n in 30:10:200
        @test gausslaguerre(n) == gausslaguerre(n,0) == gausslaguerre(n,0.0)
    end

    ##########
    # Test the special cases
    ##########

    x, w = gausslaguerre(0)
    @test x == Float64[]
    @test w == Float64[]

    x, w = gausslaguerre(1, 0.4)
    @test x == [1.4]
    @test w == [gamma(1.4)]

    x, w = gausslaguerre(2, -0.37)
    @test x ≈ [0.35328546651962944; 2.90671453348037]
    @test w ≈ [1.269857265034167; 0.15433992954023468]

    ##########
    n = 42
    ##########
    α = 0.0
    coef = rand(17)
    Z = exact_integral_poly(α, coef)

    x, w = gausslaguerre(n, α)
    @test isa(x, Vector{Float64})
    @test isa(w, Vector{Float64})
    @test abs(sum(w) - 1) < tol
    @test abs(dot(w, x) - 1) < tol
    @test abs(x[37] - 98.388267163326702) < tol
    @test abs(w[7] - 0.055372813167092) < tol
    @test abs(dot(w, polyval(x, coef)) - Z)/abs(Z) < tol

    x_gw, w_gw = FastGaussQuadrature.gausslaguerre_GW(n, α)
    @test abs(x[37] - 98.388267163326702) < tol
    @test abs(w[7] - 0.055372813167092) < tol
    @test abs(dot(w, polyval(x, coef)) - Z)/abs(Z) < tol

    x_rec, w_rec = FastGaussQuadrature.gausslaguerre_rec(n, α)
    @test abs(x[37] - 98.388267163326702) < tol
    @test abs(w[7] - 0.055372813167092) < tol
    @test abs(dot(w, polyval(x, coef)) - Z)/abs(Z) < tol

    α = 0.5
    coef = rand(17)
    Z = exact_integral_poly(α, coef)
    x, w = gausslaguerre(n, α)
    @test abs(dot(w, polyval(x, coef)) - Z)/abs(Z) < tol

    ##########
    n = 251
    ##########
    α = 0.0
    x, w = gausslaguerre(n, α)
    a = 4
    Z = exact_integral_cos1(a)
    @test isa(x, Vector{Float64})
    @test isa(w, Vector{Float64})
    @test abs(x[37] - 13.309000189442097) < tol
    @test abs(w[3] - 0.050091759039996) < tol
    @test abs(dot(w, cos.(a*x)) - Z) < tol

    α = 0.5
    x, w = gausslaguerre(n, α)
    @test abs(x[44] -19.095577327730616 ) < tol
    @test abs(w[18] - 0.026245779174690266) < tol
    a = 4
    Z = exact_integral_cos2(a)
    @test abs(dot(w, cos.(a*x)) - Z) < tol

    ############
    n = 350000
    ############
    α = 0.0
    a = 50
    Z = exact_integral_cos1(a)
    x, w = gausslaguerre(n, α)
    @test isa(x, Vector{Float64})
    @test isa(w, Vector{Float64})
    @test abs(dot(w, cos.(a*x)) - Z) < tol

    α = 0.44
    x, w = gausslaguerre(n, α; reduced = true)
    @test w[end] < 100*floatmin(Float64)
end


text:
module GPC

include("gpcbase.jl")

end


text:
const Rect2D = HyperRectangle{2, Float32}

iter_or_array(x) = repeated(x)
iter_or_array(x::Repeated) = x
iter_or_array(x::AbstractArray) = x
# We treat staticarrays as scalar
iter_or_array(x::StaticArray) = repeated(x)

function metrics_bb(char::Char, font::FTFont, scale)
    extent = get_extent(font, char) .* Vec2f0(scale)
    mini = bearing(extent)
    return Rect2D(mini, Vec2f0(extent.scale)), extent
end

"""
    iterate_extents(f, line::AbstractString, fonts, scales)
Iterates over the extends of the characters (glyphs) in line!
Newlines will be drawn like any other character.
`fonts` can be a vector of fonts, or a single font.
`scales` can be a single float or a Vec2, or a vector of any of those.

`f` will get called with `(char::Char, glyph_box::Rec2D, glyph_advance::Point2f0)`.

`char` is the currently iterated char.

`glyph_box` is the boundingbox of the glyph.
widths(box) will be the size of the bitmap, while minimum(box) is where one starts drawing the glyph.
For the minimum at y position, 0 is the where e.g. `m` starts, so `g` will start in the negative, while `^` will start positive.

`glyph_advance` The amount one advances after glyph, before drawing next glyph.
"""
function iterate_extents(f, line::AbstractString, fonts, scales)
    iterator = zip(line, iter_or_array(scales), iter_or_array(fonts))
    lastpos = 0.0
    for (char, scale, font) in iterator
        extent, glyph_box = metrics_bb(char, font, scale)
        mini = minimimum(glyph_box) .+ Vec2f0(lastpos, 0.0)
        glyph_box = Rect2D(mini, widths(glyph_box))
        glyph_advance = Point2f0(extent.advance)
        lastpos += glyph_advance[1]
        f(char, glyph_box, glyph_advance)
    end
end

function glyph_rects(line::AbstractString, fonts, scales)
    rects = Rect2D[]
    iterate_extents(line, fonts, scales) do char, box, advance
        push!(rects, box)
    end
    return rects
end

function boundingbox(line::AbstractString, fonts, scales)
    reduce(union, glyph_rects(lines, fonts, scales))
end


text:
function ArrayToCudaArray!(net,inds)
    for n in inds
        net.value[n]=CudaArray(net.value[n])
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# run_alg.jl
using SPECTrecon
include("eml_em.jl")


text:
#!/usr/bin/env julia

# rfweights create weights for RainFARM downscaling

# RainFARM
# Stochastic downscaling following
# D'Onofrio et al. 2014, J of Hydrometeorology 15 , 830-843 and
# Rebora et. al 2006, JHM 7, 724
# Includes orographic corrections

# Implementation in Julia language
# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy

using RainFARM
using ArgParse
using Compat, Compat.Printf

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--nf", "-n"
            help = "Subdivisions for downscaling"
            arg_type = Int
            default = 2
        "--weights", "-w", "--weight", "--outfile", "-o", "--out"
            help = "Output weights filename"
            arg_type = AbstractString
            default = "weights.nc"
        "--varname", "-v"
            help = "Input variable name (in orofile)"
            arg_type = AbstractString
            default = ""
	"--pass2", "-p"
            help = "Previous rainfarm output for 2-pass weighting"
            arg_type = AbstractString
            default = ""
	"--inweight", "-i"
            help = "Previous weights file to correct for 2-pass"
            arg_type = AbstractString
            default = ""
        "orofile"
            help = "The input file to use for orography"
            arg_type = AbstractString
            required = true
        "reffile"
            help = "A reference file (e.g. the file to downscale)"
            arg_type = AbstractString
            required = true
        "--conv", "-c"
            action = :store_true
            help = "conserve precipitation using convolution"
    end

    s.description="Create weights for RainFARM downscaling"
    s.version="0.1"
    s.add_version=true

    return parse_args(s)
end

args = parse_commandline()
nf=args["nf"]
reffile=args["reffile"]
orofile=args["orofile"]
weightsfn=args["weights"]
varname=args["varname"]
pass1fn=args["pass2"]
inweight=args["inweight"]
fsmooth=args["conv"]

println("Creating weights from file ",orofile)

# Create a reference gridrf.nc file (same grid as rainfarm output files)
(pr,lon_mat,lat_mat)=read_netcdf2d(reffile, varname);
# Creo la griglia fine
nss=size(pr)
if (length(nss)>=3)
    pr=pr[:,:,1]
end
ns=nss[1];
(lon_f, lat_f)=lon_lat_fine(lon_mat, lat_mat,nf);

rr=round.(Int,rand(1)*100000)

println("Output size: ",size(lon_f))
if(varname=="")
   varname="pr"
   run(`cdo -s  setname,pr $reffile reffile_rr.nc`)
   reffile="reffile_rr.nc"
end

# The rest is done in CDO

if(inweight=="")
println("Computing weights")
write_netcdf2d("gridrf.nc",reshape(pr,ns,ns,1),lon_f,lat_f,varname,reffile)
run(`cdo -s timmean $orofile pr_orofile_$rr.nc`)
run(`cdo -s -f nc copy gridrf.nc gridrf_2_$rr.nc`)
run(`cdo -s -f nc remapbil,gridrf_2_$rr.nc pr_orofile_$rr.nc pr_remap_rr.nc`)
if(fsmooth)
  (prr,lon,lat)=read_netcdf2d("pr_remap_rr.nc","")
  ww=prr./smoothconv(prr,ns);
  write_netcdf2d(weightsfn,ww,lon_f,lat_f,varname,reffile)
  run(`rm -f pr_remap_rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
else
  run(`cdo -s gridboxmean,$nf,$nf pr_remap_rr.nc pr_remap_gbm_$rr.nc`)
  run(`cdo -s remapnn,pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc`)
  run(`cdo -s div pr_remap_rr.nc pr_remap_nn_$rr.nc $weightsfn`)
  run(`rm -f pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
  inweight=weightsfn
end
else
println("Correcting weights in ",inweight)
end

if(pass1fn!="")
println("2-pass weighting using: ",pass1fn)

run(`cdo -s timmean $pass1fn pr_remap_$rr.nc`)
run(`cdo -s gridboxmean,$nf,$nf pr_remap_$rr.nc pr_remap_gbm_$rr.nc`)
run(`cdo -s -f nc copy pr_remap_$rr.nc pr_remap_2_$rr.nc`)
run(`cdo -s remapnn,pr_remap_2_$rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc`)
run(`cdo -s div pr_remap_2_$rr.nc pr_remap_nn_$rr.nc neweights_$rr.nc`)

run(`cdo -s div $inweight neweights_$rr.nc weightsn_$rr.nc`)
run(`mv weightsn_$rr.nc $weightsfn`)
run(`rm -f pr_remap_$rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc neweights_$rr.nc pr_remap_2_$rr.nc`)
end


text:
using PreprocessData
#to use dataset for functions, name_of_dataset.jl must be present in PreprocessData/src/datasets

#split dataset into train and test data, also download dataset and create a simply preprocessed csv file
#by default values are returned as `DataFrame`(toarray=true)
#two DataFrames returned
#targer values and attributes are not separated
train, test = split_traintest(
    PreprocessData.Iris(),
    trainSize=0.7,
    seed=123
)
#split dataset into train and test data
#toarray=false two `Tuple` variables are returned
#attributes - first index of tuple, labels - second index of tuple
train, test = split_traintest(
    PreprocessData.Iris(),
    trainSize=0.7,
    seed=123,
    toarray=true
)
train_attributes = train[1]
train_labels = train[2]
test_attributes = test[1]
test_labels = test[2]

#load column names from file header.csv
train, test = split_traintest(
    PreprocessData.Wine(),
    header = true
)

#split dataset into train, test data and data for validation
#by default `DataFrame` is returned
#load column names from function header (then file header.csv)
train, valid, test = split_trainvalidtest(
    PreprocessData.Iris(),
    trainSize=0.7,
    validSize=0.15,
    seed=123,
    header=true
)

#split dataset that already has valid and train data separated
train, valid, test = split_trainvalidtest(
    PreprocessData.Gisette(),
    trainSize=0.7,
    validSize=0.15,
    seed=123,
)

#split dataset into train, test data and data for validation
#arrays (Tuple) are returned
#attributes - first index of tuple, labels - second index of tuple
train, valid, test = split_trainvalidtest(
    PreprocessData.Iris(),
    toarray = true
)
train_attributes = train[1]
train_labels = train[2]
valid_attributes = valid[1]
valid_labels = valid[2]
test_attributes = test[1]
test_labels = test[2]

#split image data from MLDatasets
train, test = split_traintest(PreprocessData.CIFAR10())
train, valid, test = split_trainvalidtest(PreprocessData.CIFAR100())

#download multiple files where labels and data separated, downloading merges labels and attributes
call(PreprocessData.Gisette())

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
versioninfo()

#gc_enable(false)

using BandedMatrices

A = brand(10000,10000,2,3)
B = brand(1000,1000,200,300)
E = brand(10000,10000,-2,2)
C = rand(1000, 1000)
D = rand(10000, 1000)
V = view(brand(10000,10000,2,3), :, :)
tD = D'
v=rand(10000)
w=rand(1000)

print("A*A:")
A*A;
@time for k=1:100
    A*A
end

print("V*V:")
V*V;
@time for k=1:100
    V*V
end

print("A*D:")
A*D;
@time for k=1:10
    A*D
end

print("D*A:")
tD*A;
@time for k=1:10
    tD*A
end

print("A*E:")
A*E
@time for k=1:10
    A*E
end

print("E*D:")
E*D
@time for k=1:10
    E*D
end

print("B*B:")
B*B;
@time for k=1:10
    B*B
end

print("B*C:")
B*C;
@time for k=1:10
    B*C
end

print("C*B:")
C*B;
@time for k=1:10
    C*B
end

print("A*v:")
A*v;
@time for k=1:1000
    A*v
end

print("B*w:")
B*w;
@time for k=1:100
    B*w
end

#gc_enable(true)

# benchmark pure julia routine vs BLAS routine

# Banded * Dense

julia(A, D) = BandedMatrices._banded_generic_matmatmul!(Matrix{Float64}(size(A, 1), size(D, 2)), 'N', 'N', A, D)
blas(A, D) = BandedMatrices.gbmm!('N', 'N', one(Float64), A, D, zero(Float64), Matrix{Float64}(size(A, 1), size(D, 2)))
dense(A, D) = Array(A)*Array(D)

for n in [100, 1000, 10000, 100000]
    for b in [1, 50, 200, 500]
        A, D = brand(n, n, b, b), rand(n, 2)
        println("n = $n, b = $b")
        if b == 1 && n ≤ 10000
            print("dense")
            @time dense(A, D)
        end
        print("julia")
        @time julia(A, D);
        print("blas ")
        @time blas(A, D);
    end
end

# Banded * Banded

julia(A) = BandedMatrices._banded_generic_matmatmul!(
    BandedMatrix{Float64}(undef, size(A), (2*bandwidth(A, 1), 2*bandwidth(A, 2))), 'N', 'N', A, A)
blas(A) = BandedMatrices.gbmm!(
    'N', 'N', one(Float64), A, A, zero(Float64), BandedMatrix{Float64}(undef, size(A), (2*bandwidth(A, 1), 2*bandwidth(A, 2))))
dense(A) = (B=Array(A); B*B)

for n in [100, 1000, 10000]
    for b in [1, 50, 200, 500]
        A = brand(n, n, b, b)
        println("n = $n, b = $b")
        if b == 1 && n ≤ 1000
            print("dense")
            @time dense(A)
        end
        print("julia")
        @time julia(A);
        print("blas ")
        @time blas(A);
    end
end

# println("Large dimension QR solve")
# A=brand(100000,100000,3,4)
# b=rand(size(A,1))
# @time Q,R=qr(A)
# @time R\(Q'*b)

# println("Time should be   0.030149 seconds (23 allocations: 17.548 MB)")
# println("Time should be   0.008210 seconds (29 allocations: 8.012 MB)")

# A=brand(2000,2000,400,400)
# b=rand(size(A,1))
# @time Q,R=qr(A)
# @time R\(Q'*b)

# println("Time should be   0.401671 seconds (24 allocations: 36.665 MB)")
# println("Time should be   0.009276 seconds (30 allocations: 12.262 MB

nothing


text:
using NPZ, Compat

@static if VERSION >= v"0.7.0-DEV.2575"
    using Test
else
    using Base.Test
end

Debug = false

tmp = mktempdir()
if Debug
    println("temporary directory: $tmp")
end

TestArrays = Any[
    true,
    false,
    Int8(-42),
    Int16(-42),
    Int32(-42),
    Int64(-42),
    UInt8(42),
    UInt16(42),
    UInt32(42),
    UInt64(42),
    Float16(3.1415),
    Float32(3.1415),
    Float64(3.1415),
    Complex{Float32}(1, 7),
    Complex{Float64}(1, 7),
    Bool[0, 1, 0, 1, 1, 0],
    Int8[-42, 0, 1, 2, 3, 4],
    Int16[-42, 0, 1, 2, 3, 4],
    Int32[-42, 0, 1, 2, 3, 4],
    Int64[-42, 0, 1, 2, 3, 4],
    UInt8[0, 1, 2, 3, 4],
    UInt16[0, 1, 2, 3, 4],
    UInt32[0, 1, 2, 3, 4],
    UInt64[0, 1, 2, 3, 4],
    Float64[-42, 0, 1, 2, 3.14, 4],
    Float32[-42, 0, 1, 2, 3.14, 4],
    Float16[-42, 0, 1, 2, 3.14, 4],
    Complex{Float32}[1+2im, 3, 4+5im, 6im, 7+8im],
    Complex{Float64}[1+2im, 3, 4+5im, 6im, 7+8im],
    [i-j for i in 1:3, j in 1:5],
    [i-j+k for i in 1:3, j in 1:4, k in 1:5],
    rand(512, 512, 16)
]

# Write a NPZ file with all the test arrays and numbers,
# and read it back in.
old = Dict{String, Any}()
for (i, x) in enumerate(TestArrays)
    old["testvar_" * string(i)] = x
end
filename = "$tmp/big.npz"
npzwrite(filename, old)
output = npzread(filename)
for (k, v) in old
    @test v == output[k]
end
@test old == output

testdict = Dict("one"=>21, "two"=>34)
dictfilename = "$tmp/dict.npz"
npzwrite(dictfilename, testdict)
dictoutput = npzread(dictfilename)
for (k, v) in testdict
    @test v == dictoutput[k]
end
# Write and then read NPY files for each array
for x in TestArrays
    let filename = "$tmp/test.npy"
        npzwrite(filename, x)
        y = npzread(filename)
        @test typeof(x) == typeof(y)
        @test size(x) == size(y)
        @test x == y
    end
end

#v2 tests
v2_solution = Array{Float64,1}([0.1*i for i in 0:10])
v2_read = npzread("linspace.npy")
@test(typeof(v2_read) == typeof(v2_solution))
@test(size(  v2_read) == size(v2_solution))
@test(       v2_read  == v2_solution)

if !Debug
    run(`rm -rf $tmp`)
end

@test NPZ.parseinteger("123Hello") == (123, "Hello")
@test NPZ.parseinteger("123LHello") == (123, "Hello")