text:
@testset "rectifiedlinearunit" begin
    @testset "one" begin
        X = [0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0]
        W = [1.0; 1.0]
        bias = 0.0
        rlu = RectifiedLinearUnit(W, bias)
        @test forward(rlu, X) == [0.0, 1.0, 1.0, 2.0]
    end

end

max_stars_repo_name: helgith/ExampleNeuralNetworks

max_issues_repo_name: helgith/ExampleNeuralNetworks

max_forks_repo_name: helgith/ExampleNeuralNetworks


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

module StudentData

export Student, students, sortedStudents, studentCompareFn

# Complex type
struct Student
    name::String
    rollno::String
    cgpa::Float64
end
# instances of the complex type Student
students = [
    Student("jack", "s2020123", 8.9),
    Student("jill", "s2020124", 4.9),
    Student("tom", "s2020125", 6.9),
    Student("pip", "s2020126", 9.9),
    Student("mark", "s2020127", 7.9)
]

# expected result after sorting students (descending) by CGPA
sortedStudents = [
    Student("pip", "s2020126", 9.9),
    Student("jack", "s2020123", 8.9),
    Student("mark", "s2020127", 7.9),
    Student("tom", "s2020125", 6.9),
    Student("jill", "s2020124", 4.9)
]
# comparison function for student
function studentCompareFn(x::Student, y::Student)
    x.cgpa > y.cgpa
end
end

ext: jl

max_stars_repo_name: codekomali/JuliaPractice

max_issues_repo_name: codekomali/JuliaPractice

max_forks_repo_name: codekomali/JuliaPractice


text:
module DS19Presentation

export load, slide4, slide5, slide6, slide8_9, slide10_11, slide13, slide15,
    slide17_18_19_20, slide22_23, slide24_25, slide27, slide28, slide30_31,
    slide_a2, slide_a3, slide_a4_5, slide_a6_7, slide_a8, slide_a9,
    pgfplots, save_animation, animate

# Fix for https://github.com/JuliaIO/ImageMagick.jl/issues/140
using ImageMagick

using DataDeps
using MD5
using Serialization
using StorageGraphs

# ENV["DATADEPS_ALWAYS_ACCEPT"] = true

function __init__()
    include("$(@__DIR__)/data.jl")
end

using NuclearSurfaceVibrations
using .Classical
using .Visualizations
using .InitialConditions: depchain, initial_conditions!
using .Classical.ParallelTrajectories
using StaticArrays
using OrdinaryDiffEq
using LaTeXStrings
using IntervalArithmetic
using StatsBase
using AbstractPlotting: Node, save, hbox, vbox

include("benchmarks.jl")

# open("src/data.jl", "w") do fh
#     registration = generate(Figshare(), "https://figshare.com/articles/DS19_presentation_materials/8146106")
#     print(fh, registration)
# end

const bg = colorant"#FAFAFA"
const ac = colorant"#EB811B"

function load()
    deserialize(datadep"DS19 presentation materials/graph.jls")
end

function ics_separate(g, E=0.1, ic_alg=PoincareRand(n=5), B=0.15)
    p = PhysicalParameters(B=B)
    q0, p0 = initial_conditions!(g, E, alg=ic_alg, params=p)
    p₀ = [SVector{2}(p0[i, :]) for i ∈ axes(p0, 1)]
    q₀ = [SVector{2}(q0[i, :]) for i ∈ axes(q0, 1)]

    return p₀, q₀
end

function ics(g, E=0.1, ic_alg=PoincareRand(n=5), B=0.15)
    p₀, q₀ = ics_separate(g, E, ic_alg, B)
    z0 = [vcat(p₀[i], q₀[i]) for i ∈ axes(q₀, 1)]
end

function get_sol(g, i=1, ic_alg=PoincareRand(n=5); B=0.15, E=0.1)
    p = PhysicalParameters(B=B)
    z0 = ics(g, E, ic_alg, B)
    prob = ODEProblem(ż, z0[i], 40., p)
    sol = solve(prob, Vern9(), abstol=1e-14, reltol=1e-14, maxiters=1e9)
end

function get_sim(g, ic_alg=PoincareRand(n=5); B=0.15, E=0.1)
    p = PhysicalParameters(B=B)
    q0, p0 = initial_conditions!(g, E, alg=ic_alg, params=p)
    poincaremap(q0, p0, params=p, t=40)
end

function get_psol(g, i=1, ic_alg=PoincareRand(n=5); B=0.15, E=0.1)
    p = PhysicalParameters(B=B)
    z0 = ics(g, E, ic_alg, B)
    d0 = 1e-3
    pprob = parallel_problem(ż, (z0[i], z0[i].+d0/√4), 40., p)
    psol = solve(pprob, Vern9(), abstol=1e-14, reltol=1e-14, maxiters=1e9)
end

function slide4(g; saveimage=false)
    sol = get_sol(g)
    t = Node(0.)

    surface_sc = animate_solution(sol, t)
    if saveimage
        path = joinpath("assets", "nucleus.png")
        save(path, surface_sc)
    end
    return surface_sc
end

"""
    slide5(g; saveimage=false, savevideo=false)

This function returns the time node and the scene for the
simulation of the nuclear surface and its sections.
Call as `t, sc = slide5(g);` and then display the scene
using `display(sc)` and use `animate(t, (0, 40))` to animate
it.
"""
function slide5(g; saveimage=false, savevideo=false)
    sol = get_sol(g)
    t = Node(0.)

    surface_sc = animate_solution(sol, t)
    section_sc = θϕ_sections(sol, t, surface_sc.limits[])
    sc = vbox(surface_sc, section_sc, sizes=[0.7, 0.3])
    if saveimage
        path = joinpath("assets", "nucleus-with-sections.png")
        save(path, sc)
    end
    if savevideo
        path = joinpath("assets", "nucleus-with-sections.mp4")
        save_animation(sc, t, (0, 10), path)
    end
    return t, sc
end

function nucleus_poincare(g, B, E; i=1)
    sol = get_sol(g, i, B=B, E=E)
    sim = get_sim(g, B=B, E=E)
    t = Node(0.)

    surface_sc = animate_solution(sol, t)
    line_sc = path_animation3D(sol, t)
    plot_slice!(line_sc, sim[i])

    sc = vbox(surface_sc, line_sc, sizes=[0.5,0.5])
    return t, sc
end

"""
    slide6(g; saveimage=false, savevideo=false)

This function returns the time node and the scene for the
simulation of the nuclear surface and the motion in phase space.
Call as `t, sc = slide6(g);` and then display the scene
using `display(sc)` and use `animate(t, (0, 40))` to animate
it.
"""
function slide6(g; saveimage=false, savevideo=false)
    t, sc = nucleus_poincare(g, 0.5, 0.3, i=3)
    if saveimage
        path = joinpath("assets", "nucleus-with-poincare.png")
        save(path, sc)
    end
    if savevideo
        path = joinpath("assets", "nucleus-with-poincare.mp4")
        save_animation(sc, t, (0, 40), path)
    end
    return t, sc
end

"""
    slide8_9(g)

This function returns the two plots for short time benchmarks.
Call as `p1, p2 = slide8_9(g)`.
"""
function slide8_9(g)
    p1, p2 = short_benchmark(g)
    savefig(p1, "assets/short-benchmark-E.tex")
    savefig(p2, "assets/short-benchmark-t.tex")

    return p1, p2
end

"""
    slide10_11(g)

This function returns the two plots for long time benchmarks.
Run twice in a row for meaningfull results, as the first time
will include precompilation times. See the julia documentation for
details. Call as `p1, p2 = slide10_11(g)`.
"""
function slide10_11(g)
    p1, p2 = long_benchmark(g)
    savefig(p1, "assets/long-benchmark-E.tex")
    savefig(p2, "assets/long-benchmark-t.tex")

    return p1, p2
end

"""
    slide13(g; saveimage=false, savevideo=false)

This function returns the time node and the scene for the
simulation of the distance between two close trajectories.
Call as `t, sc = slide13(g);` and then display the scene
using `display(sc)` and use `animate(t, (0, 40))` to animate
it.
"""
function slide13(g; saveimage=false, savevideo=false)
    psol = get_psol(g, 3, B=0.5, E=0.3)
    t = Node(0.)

    parallel_sc = parallel_paths(psol, t, linecolor=ac)
    psc = paths_distance(psol, t)
    lpsc = paths_distance_log(psol, t)
    dsc = hbox(psc, lpsc)
    sc = vbox(parallel_sc, dsc, sizes=[0.5, 0.5])

    if saveimage
        path = joinpath("assets", "dist-with-log.png")
        save(path, sc)
    end
    if savevideo
        path = joinpath("assets", "dist-with-log.mp4")
        save_animation(sc, t, (0, 40), path)
    end
    return t, sc
end

"""
    slide15(g; saveimage=false)

This function returns the interactive Poincaré section.
Call as `slide15(g)`.
"""
function slide15(g; saveimage=false)
    sc = poincare_explorer(g, 120., DynSys(T=1e5), PoincareRand(n=500), t=1e4,
        params=PhysicalParameters(B=0.5), markersize=0.08)
    if saveimage
        path = joinpath("assets", "poincare.png")
        save(path, sc)
    end
    return sc
end

"""
    slide17_18_19_20(g)

This function produces the .tex files for the histogram evolution
of λ with respect to the integration time in the `assets` folder.
"""
function slide17_18_19_20(g)
    E = 120.
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    ic_dep = InitialConditions.depchain(p,E,ic_alg)
    for T in 10. .^(4:7)
        λs = g[:λ, ic_dep..., (λ_alg=DynSys(T=T),)][1]
        tpow = Int(log10(T))
        plt = histogram(λs, nbins = 50, xlabel=L"\lambda", ylabel=L"N",
            label=latexstring("T=10^$tpow"), lw=0, framestyle=:grid,
            color=colorant"#6699CC", background_color=bg, markerstrokealpha=0,
            tex_output_standalone=true)
        savefig(plt, "assets/hist-lambda-$tpow.tex")
    end
end

"""
    slide22_23(g)

This functuion returns the histogram with the selected chaotic part.
"""
function slide22_23(g)
    E = 120.
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    λhist, shist = selected_hist(g, E, DynSys(T=1e5), ic_alg, params=p)
    plt = plot(λhist, xlabel=L"\lambda", ylabel=L"N", label=L"T=10^5",
        tex_output_standalone=true, lw=0, framestyle=:grid, background_color=bg,
        color=colorant"#6699CC")
    savefig(plt, "assets/hist-lambda-selected-before.tex")

    plot!(plt, shist, label="chaotic", lw=0, color=colorant"#DDCC77")
    savefig(plt, "assets/hist-lambda-selected-after.tex")
    return plt
end

"""
    slide24_25(g)

This function returns the plot of the maximal Lyapunov exponent averaged
over initial conditions.
Call as `p1, p2 = slide24_25(g)`.
"""
function slide24_25(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    p1 = mean_over_ic(g, DynSys(T=1e5), ic_alg, params=p, Einterval=10:10:3000,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC", lw=2.5,
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true)
    savefig(p1, "assets/mean-over-ic.tex")
    p2 = mean_over_ic(g, DynSys(T=1e5), ic_alg, params=p, Einterval=0.01:0.01:10,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC", lw=2.5,
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true)
    savefig(p2, "assets/mean-over-ic-low.tex")
    return p1, p2
end

"""
    slide27(g)

This function returns the plot of d∞ averaged over initial conditions.
Call as `slide27(g)`.
"""
function slide27(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    plt = mean_over_ic(g, DInftyAlgorithm(T=1e5), ic_alg, params=p,
        Einterval=0.01:0.01:10, framestyle=:grid, background_color=bg,
        color=colorant"#6699CC", lw=2.5, markerstrokewidth=0, markerstrokealpha=0,
        tex_output_standalone=true)
    savefig(plt, "assets/mean-over-ic-dinf.tex")
    return plt
end

"""
    slide28(g)

This function returns the plot of Γ averaged over initial conditions.
Call as `slide28(g)`.
"""
function slide28(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    plt = mean_over_ic(g, DynSys(T=1e5), DInftyAlgorithm(T=1e5), ic_alg, params=p,
        Einterval=0.01:0.02:10, framestyle=:grid, background_color=bg,
        color=colorant"#6699CC", lw=2.5, markerstrokewidth=0, markerstrokealpha=0,
        tex_output_standalone=true)
    savefig(plt, "assets/mean-over-ic-Gamma.tex")
    return plt
end

"""
    slide30_31(g)

This function returns the plot of the maximal Lyapunov exponent averaged
over the energy.
Call as `p1, p2 = slide30_31(g)`.
"""
function slide30_31(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    p1 = mean_over_E(g, DynSys(T=1e5), 10:10:3000, ic_alg=ic_alg,
        Binterval=0.1:0.02:0.6, legend=false, markersize=9,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC",
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true)
    savefig(p1, "assets/mean-over-E.tex")
    p2 = mean_over_E(g, DynSys(T=1e5), 0.01:0.01:10, ic_alg=ic_alg,
        Binterval=0.1:0.02:0.6, legend=false, markersize=9,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC",
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true)
    savefig(p2, "assets/mean-over-E-low.tex")
    return p1, p2
end

function slide_a2()
    x = y = range(-5, 5, length=100)
    plt = contour(x,y,
        (x,y)->V((x,y),PhysicalParameters(B=0.5)),
        xlabel=L"q_0", ylabel=L"q_2", colorbar_title=L"V",
        levels=50, framestyle=:grid, tex_output_standalone=true,
        background_color=bg)
    savefig(plt, "assets/potential-contour.tex")

    return plt
end

function slide_a3(g; saveimage=false, savevideo=false)
    t, sc = nucleus_poincare(g, 0.4, 0.1, i=1)
    if saveimage
        path = joinpath("assets", "nucleus-with-poincare-regular.png")
        save(path, sc)
    end
    if savevideo
        path = joinpath("assets", "nucleus-with-poincare-regular.mp4")
        save_animation(sc, t, (0, 40), path)
    end
    return t, sc
end

function slide_a4_5(g)
    p1, p2 = short_benchmark(g, rescaling=true)
    savefig(p1, "assets/short-benchmark-rescaling-E.tex")
    savefig(p2, "assets/short-benchmark-rescaling-t.tex")

    return p1, p2
end

function slide_a6_7(g)
    p1, p2 = long_benchmark(g, rescaling=true)
    savefig(p1, "assets/long-benchmark-rescaling-E.tex")
    savefig(p2, "assets/long-benchmark-rescaling-t.tex")

    return p1, p2
end

function slide_a8(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    plt = mean_over_ic(g, DynSys(T=1e5), ic_alg, params=p, Einterval=10:10:3000,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC", lw=2.5,
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true,
        label=latexstring("T=10^5"), legend=:topright)
    mean_over_ic(g, DynSys(T=1e6), ic_alg, params=p, Einterval=10:10:3000,
        plt=plt, color=colorant"#DDCC77", lw=2.5, markerstrokewidth=0,
        markerstrokealpha=0, label=latexstring("T=10^6"), legend=:topright)
    savefig(plt, "assets/mean-over-ic-T-comparison.tex")

    return plt
end

function slide_a9(g)
    p = PhysicalParameters(B=0.5)
    ic_alg = PoincareRand(n=500)
    plt = mean_over_ic(g, DynSys(T=1e5), ic_alg, params=p, Einterval=10:10:3000,
        framestyle=:grid, background_color=bg, color=colorant"#6699CC", lw=2.5,
        markerstrokewidth=0, markerstrokealpha=0, tex_output_standalone=true,
        label=latexstring("B=$(p.B)"), legend=:topright)

    p = PhysicalParameters(B=0.6)
    mean_over_ic(g, DynSys(T=1e6), ic_alg, params=p, Einterval=10:10:3000,
        plt=plt, color=colorant"#DDCC77", lw=2.5, markerstrokewidth=0,
        markerstrokealpha=0, label=latexstring("B=$(p.B)"), legend=:topright)
    savefig(plt, "assets/mean-over-ic-B-difference.tex")

    return plt
end

end # module


text:
#######
# Tag #
#######

struct Tag{F,V}
end

const TAGCOUNT = Threads.Atomic{UInt}(0)

# each tag is assigned a unique number
# tags which depend on other tags will be larger
@generated function tagcount(::Type{Tag{F,V}}) where {F,V}
    :($(Threads.atomic_add!(TAGCOUNT, UInt(1))))
end

function Tag(f::F, ::Type{V}) where {F,V}
    tagcount(Tag{F,V}) # trigger generated function
    Tag{F,V}()
end

Tag(::Nothing, ::Type{V}) where {V} = nothing

@inline function ≺(::Type{Tag{F1,V1}}, ::Type{Tag{F2,V2}}) where {F1,V1,F2,V2}
    tagcount(Tag{F1,V1}) < tagcount(Tag{F2,V2})
end

struct InvalidTagException{E,O} <: Exception
end

Base.showerror(io::IO, e::InvalidTagException{E,O}) where {E,O} =
    print(io, "Invalid Tag object:\n  Expected $E,\n  Observed $O.")

checktag(::Type{Tag{FT,VT}}, f::F, x::AbstractArray{V}) where {FT,VT,F,V} =
    throw(InvalidTagException{Tag{F,V},Tag{FT,VT}}())

checktag(::Type{Tag{F,V}}, f::F, x::AbstractArray{V}) where {F,V} = true

# no easy way to check Jacobian tag used with Hessians as multiple functions may be used
checktag(::Type{Tag{FT,VT}}, f::F, x::AbstractArray{V}) where {FT<:Tuple,VT,F,V} = true

# custom tag: you're on your own.
checktag(z, f, x) = true

##################
# AbstractConfig #
##################

abstract type AbstractConfig{N} end

Base.copy(cfg::AbstractConfig) = deepcopy(cfg)

Base.eltype(cfg::AbstractConfig) = eltype(typeof(cfg))

@inline chunksize(::AbstractConfig{N}) where {N} = N

####################
# DerivativeConfig #
####################

struct DerivativeConfig{T,D} <: AbstractConfig{1}
    duals::D
end

"""
    ForwardDiff.DerivativeConfig(f!, y::AbstractArray, x::AbstractArray)

Return a `DerivativeConfig` instance based on the type of `f!`, and the types/shapes of the
output vector `y` and the input vector `x`.

The returned `DerivativeConfig` instance contains all the work buffers required by
`ForwardDiff.derivative` and `ForwardDiff.derivative!` when the target function takes the form
`f!(y, x)`.

If `f!` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `y` or `x`.
"""
function DerivativeConfig(f::F,
                          y::AbstractArray{Y},
                          x::X,
                          tag::T = Tag(f, X)) where {F,X<:Real,Y<:Real,T}
    duals = similar(y, Dual{T,Y,1})
    return DerivativeConfig{T,typeof(duals)}(duals)
end

checktag(::DerivativeConfig{T},f,x) where {T} = checktag(T,f,x)
Base.eltype(::Type{DerivativeConfig{T,D}}) where {T,D} = eltype(D)

##################
# GradientConfig #
##################

struct GradientConfig{T,V,N,D} <: AbstractConfig{N}
    seeds::NTuple{N,Partials{N,V}}
    duals::D
end

"""
    ForwardDiff.GradientConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `GradientConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `GradientConfig` instance contains all the work buffers required by
`ForwardDiff.gradient` and `ForwardDiff.gradient!`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function GradientConfig(f::F,
                        x::AbstractArray{V},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(f, V)) where {F,V,N,T}
    seeds = construct_seeds(Partials{N,V})
    duals = similar(x, Dual{T,V,N})
    return GradientConfig{T,V,N,typeof(duals)}(seeds, duals)
end

checktag(::GradientConfig{T},f,x) where {T} = checktag(T,f,x)
Base.eltype(::Type{GradientConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N}

##################
# JacobianConfig #
##################

struct JacobianConfig{T,V,N,D} <: AbstractConfig{N}
    seeds::NTuple{N,Partials{N,V}}
    duals::D
end

"""
    ForwardDiff.JacobianConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `JacobianConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `JacobianConfig` instance contains all the work buffers required by
`ForwardDiff.jacobian` and `ForwardDiff.jacobian!` when the target function takes the form
`f(x)`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function JacobianConfig(f::F,
                        x::AbstractArray{V},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(f, V)) where {F,V,N,T}
    seeds = construct_seeds(Partials{N,V})
    duals = similar(x, Dual{T,V,N})
    return JacobianConfig{T,V,N,typeof(duals)}(seeds, duals)
end

"""
    ForwardDiff.JacobianConfig(f!, y::AbstractArray, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `JacobianConfig` instance based on the type of `f!`, and the types/shapes of the
output vector `y` and the input vector `x`.

The returned `JacobianConfig` instance contains all the work buffers required by
`ForwardDiff.jacobian` and `ForwardDiff.jacobian!` when the target function takes the form
`f!(y, x)`.

If `f!` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `y` or `x`.
"""
function JacobianConfig(f::F,
                        y::AbstractArray{Y},
                        x::AbstractArray{X},
                        ::Chunk{N} = Chunk(x),
                        ::T = Tag(f, X)) where {F,Y,X,N,T}
    seeds = construct_seeds(Partials{N,X})
    yduals = similar(y, Dual{T,Y,N})
    xduals = similar(x, Dual{T,X,N})
    duals = (yduals, xduals)
    return JacobianConfig{T,X,N,typeof(duals)}(seeds, duals)
end

checktag(::JacobianConfig{T},f,x) where {T} = checktag(T,f,x)
Base.eltype(::Type{JacobianConfig{T,V,N,D}}) where {T,V,N,D} = Dual{T,V,N}

#################
# HessianConfig #
#################

struct HessianConfig{T,V,N,DG,DJ} <: AbstractConfig{N}
    jacobian_config::JacobianConfig{T,V,N,DJ}
    gradient_config::GradientConfig{T,Dual{T,V,N},N,DG}
end

"""
    ForwardDiff.HessianConfig(f, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `HessianConfig` instance based on the type of `f` and type/shape of the input
vector `x`.

The returned `HessianConfig` instance contains all the work buffers required by
`ForwardDiff.hessian` and `ForwardDiff.hessian!`. For the latter, the buffers are
configured for the case where the `result` argument is an `AbstractArray`. If
it is a `DiffResult`, the `HessianConfig` should instead be constructed via
`ForwardDiff.HessianConfig(f, result, x, chunk)`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function HessianConfig(f::F,
                       x::AbstractArray{V},
                       chunk::Chunk = Chunk(x),
                       tag = Tag(f, V)) where {F,V}
    jacobian_config = JacobianConfig(f, x, chunk, tag)
    gradient_config = GradientConfig(f, jacobian_config.duals, chunk, tag)
    return HessianConfig(jacobian_config, gradient_config)
end

"""
    ForwardDiff.HessianConfig(f, result::DiffResult, x::AbstractArray, chunk::Chunk = Chunk(x))

Return a `HessianConfig` instance based on the type of `f`, types/storage in `result`, and
type/shape of the input vector `x`.

The returned `HessianConfig` instance contains all the work buffers required by
`ForwardDiff.hessian!` for the case where the `result` argument is an `DiffResult`.

If `f` is `nothing` instead of the actual target function, then the returned instance can
be used with any target function. However, this will reduce ForwardDiff's ability to catch
and prevent perturbation confusion (see https://github.com/JuliaDiff/ForwardDiff.jl/issues/83).

This constructor does not store/modify `x`.
"""
function HessianConfig(f::F,
                       result::DiffResult,
                       x::AbstractArray{V},
                       chunk::Chunk = Chunk(x),
                       tag = Tag(f, V)) where {F,V}
    jacobian_config = JacobianConfig((f,gradient), DiffResults.gradient(result), x, chunk, tag)
    gradient_config = GradientConfig(f, jacobian_config.duals[2], chunk, tag)
    return HessianConfig(jacobian_config, gradient_config)
end

checktag(::HessianConfig{T},f,x) where {T} = checktag(T,f,x)
Base.eltype(::Type{HessianConfig{T,V,N,DG,DJ}}) where {T,V,N,DG,DJ} =
    Dual{T,Dual{T,V,N},N}


text:
SEED_Char(io::IO, BUF::SeisIOBuf, nb::UInt16) = replace(String(read(io, nb)),
                                            ['\r', '\0'] =>"")

# This gets special handling as Float64 arrays are assumed in S.x[c]
function SEED_Float64!(io::IO, S::SeisData, c::Int64, xi::Int64, nb::UInt16)
  buf = getfield(BUF, :buf)
  x = getindex(getfield(S, :x), c)
  readbytes!(io, buf, nb)
  nx = div(nb, 8)
  if getfield(BUF, :xs) == true
    j = xi
    @inbounds for i = 1:nx
      j += 1
      y  = UInt32(buf[8*i-7]) << 56
      y |= UInt32(buf[8*i-6]) << 48
      y |= UInt32(buf[8*i-5]) << 40
      y |= UInt32(buf[8*i-4]) << 32
      y |= UInt32(buf[8*i-3]) << 24
      y |= UInt32(buf[8*i-2]) << 16
      y |= UInt32(buf[8*i-1]) << 8
      y |= UInt32(buf[8*i])
      x[j] = y
    end
  else
    xr = reinterpret(Float64, buf)
    copyto!(x, xi+1, xr, 1, nx)
  end
  setfield!(BUF, :k, Int64(nx))
  return nothing
end

function SEED_Unenc!(io::IO, S::SeisData, c::Int64, xi::Int64, nb::UInt16, nx::UInt16)
  buf = getfield(BUF, :buf)
  checkbuf_8!(buf, xi)
  x = getindex(getfield(S, :x), c)
  readbytes!(io, buf, nb)
  T::Type = if BUF.fmt == 0x01
      Int16
    elseif BUF.fmt == 0x03
      Int32
    elseif BUF.fmt == 0x04
      Float32
    end
  if getfield(BUF, :swap) == true
    xr = bswap.(reinterpret(T, buf))
  else
    xr = reinterpret(T, buf)
  end
  copyto!(x, xi+1, xr, 1, nx)
  setfield!(BUF, :k, Int64(nx))
  return nothing
end

function SEED_Geoscope!(io::IO, BUF::SeisIOBuf)
  mm = 0x0fff
  gm = BUF.fmt == 0x0d ? 0x7000 : 0xf000
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(read(io, UInt16)) : read(io, UInt16)
    m = Int32(x & mm)
    g = Int32((x & gm) >> 12)
    ex = -1*g
    setindex!(BUF.x, ldexp(Float64(m-2048), ex), i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_CDSN!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(read(io, UInt16)) : read(io, UInt16)
    m = Int32(x & 0x3fff)
    g = Int32((x & 0xc000) >> 14)
    mult = 4^g * g==3 ? 2 : 1
    m -= 0x1fff
    setindex!(BUF.x, m*mult, i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_SRO!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = BUF.swap ? bswap(read(io, UInt16)) : read(io, UInt16)
    m = Int32(x & 0x0fff)
    g = Int32((x & 0xf000) >> 12)
    if m > 0x07ff
      m -= 0x1000
    end
    ex = -1*g + 10
    setindex!(BUF.x, ldexp(Float64(m), ex), i)
  end
  BUF.k = BUF.n
  return nothing
end

function SEED_DWWSSN!(io::IO, BUF::SeisIOBuf)
  for i = 0x0001:BUF.n
    x = signed(UInt32(BUF.swap ? bswap(read(io, UInt16)) : read(io, UInt16)))
    BUF.x[i] = x > 32767 ? x - 65536 : x
  end
  BUF.k = BUF.n
  return nothing
end

# Steim1 or Steim2
function SEED_Steim!(io::IO, BUF::SeisIOBuf, nb::UInt16)
  x = getfield(BUF, :x)
  buf = getfield(BUF, :buf)
  ff = getfield(BUF, :x32)
  nc = Int64(div(nb, 0x0040))
  ni = div(nb, 0x0004)
  readbytes!(io, buf, nb)

  # Parse buf as UInt32s
  if ni > lastindex(ff)
    resize!(ff, ni)
  end
  yy = zero(UInt32)
  if getfield(BUF, :xs) == true
    @inbounds for ib = 1:ni
      yy  = UInt32(buf[4*ib-3]) << 24
      yy |= UInt32(buf[4*ib-2]) << 16
      yy |= UInt32(buf[4*ib-1]) << 8
      yy |= UInt32(buf[4*ib])
      ff[ib] = yy
    end
  else
    @inbounds for il = 1:ni
      yy  = UInt32(buf[4*il]) << 24
      yy |= UInt32(buf[4*il-1]) << 16
      yy |= UInt32(buf[4*il-2]) << 8
      yy |= UInt32(buf[4*il-3])
      ff[il] = yy
    end
  end

  k = zero(Int64)
  x0 = zero(Float32)
  xn = zero(Float32)
  a = zero(UInt8)
  b = zero(UInt8)
  c = zero(UInt8)
  d = zero(UInt8)
  fq = zero(Float32)
  m = zero(UInt8)
  p = zero(UInt32)
  q = zero(Int32)
  u = zero(UInt32)
  y = zero(UInt32)
  z = zero(UInt32)
  χ = zero(UInt32)
  r = zero(Int64)
  for i = 1:nc
    z = getindex(ff, 1+r)
    for j = 1:16
      χ = getindex(ff, j+r)
      y = (z >> steim[j]) & 0x00000003
      if y == 0x00000001
        a = 0x00
        b = 0x08
        c = 0x04
      elseif BUF.fmt == 0x0a
        a = 0x00
        if y == 0x00000002
          b = 0x10
          c = 0x02
        elseif y == 0x00000003
          b = 0x20
          c = 0x01
        end
      else
        p = χ >> 0x0000001e
        if y == 0x00000002
          a = 0x02
          if p == 0x00000001
            b = 0x1e
            c = 0x01
          elseif p == 0x00000002
            b = 0x0f
            c = 0x02
          elseif p == 0x00000003
            b = 0x0a
            c = 0x03
          end
        elseif y == 0x00000003
          if p == 0x00000000
            a = 0x02
            b = 0x06
            c = 0x05
          elseif p == 0x00000001
            a = 0x02
            b = 0x05
            c = 0x06
          else
            a = 0x04
            b = 0x04
            c = 0x07
          end
        end
      end
      if y != 0x00000000
        u = χ << a
        m = zero(UInt8)
        d = 0x20 - b
        while m < c
          k = k + 1
          q = signed(u)
          q >>= d
          fq = Float32(q)
          setindex!(x, fq, k)
          m = m + 0x01
          u <<= b
        end
      end
      if i == 1
        if j == 2
          x0 = Float32(signed(χ))
        elseif j == 3
          xn = Float32(signed(χ))
        end
      end
    end
    r = r+16
  end

  if BUF.wo != 0x01
    vx = view(getfield(BUF, :x), 1:k)
    reverse!(vx)
  end
  setindex!(x, x0, 1)

  # Cumsum by hand
  xa = copy(x0)
  @inbounds for i1 = 2:k
    xa = xa + getindex(x, i1)
    setindex!(x, xa, i1)
  end

  # Check data values
  if isapprox(getindex(x, k), xn) == false
    println(stdout, string("RDMSEED: data integrity -- Steim-",
                            getfield(BUF, :fmt) - 0x09, " sequence #",
                            String(copy(getfield(BUF, :seq))),
                            " integrity check failed, last_data=",
                            getindex(getfield(BUF, :x), k),
                            ", should be xn=", xn))
  end
  setfield!(BUF, :k, k)
  return nothing
end


text:
#=
  Test of Bernoulli

  See
     ~/blog/bernoulli_test.blog
     ~/webppl/bernoulli_test.wppl
=#

# using Memoization
using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)
include("jl_utils.jl")

@model function bernoulli_test(y)
    prob ~ Beta(4,4)
    n = length(y)
    for i in 1:n
        y[i] ~ Bernoulli(prob)
    end

end

# Generated by R:
# > rbinom(120,1,0.3)
# I.e. prob should be around 0.3
#=
y = [0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,1,1,0,1,0,
    1,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,1,1,0,1,0,1,0,0,1,1,0,0,0,1,0,1,0,0,0,
    0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,1,0,1,1,0,0,1,0,0,1,0,1,1,
    0,1,0,1,0,0,0,0,0,1,0,0]
=#
# Generate data
y = rand(Bernoulli(0.3),1_000)
println("length(y):$(length(y))")
# println("y:$y")
model = bernoulli_test(y)
num_chns = 4

# chns = sample(model, Prior(), MCMCThreads(), 1000, num_chns)

# chns = sample(model, MH(), 10_000)
chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, PG(20), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, PG(20), 1_000)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns) # Nope: does not handle the observation
# chns = sample(model, IS(), 10_000)

# chns = sample(model, SMC(), 10_000)

# chns = sample(model, NUTS(1000,0.65), MCMCThreads(), 40_000, num_chns)
# chns = sample(model, Gibbs(MH(:zlabels),NUTS(1000,0.65,:m,:b,:sigma)), MCMCThreads(), 40_000, num_chns)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "GPS L5" begin
    gpsl5 = GPSL5()
    @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x35, 50)) == true

    @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x35, 5)) == false

    @test @inferred(Tracking.is_upcoming_integration_new_bit(gpsl5, 0x3ca, 10)) == true # 0x3ca == 1111001010

    @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1))
    @test @inferred(Tracking.get_default_correlator(gpsl5, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3))
end


text:
using NRES
using Test

include("test_NRES.jl")
include("test_tile_map.jl")
include("test_composite_map.jl")

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
module TestPreservingFuncs

using Base.Test
using DataArrays
using DataFrames
using Dates

using TimeData

println("Running type preserving function tests")

allTypes = (:Timedata, :Timenum, :Timematr)

################
## hcat tests ##
################

tm = Timematr(rand(2, 3))
hcat(tm[:, 1], tm[:, 2], tm[:, 3])

td = Timedata(rand(2, 3))
hcat(td[:, 1], td[:, 2])
hcat(td[:, 1], td[:, 2], td[:, 1])
hcat(tm[:, 1], tm[:, 2], tm[:, 1])
@test_throws ErrorException hcat(tm[:, 1], td[:, 2])
@test_throws ErrorException hcat(tm[:, 1], tm[:, 2], td[:, 2])

################
## vcat tests ##
################

## test information retrieval functions
for t in allTypes
    eval(quote
        ## unequal column names
        td = $(t)(rand(10, 2))
        td2 = $(t)(rand(10,2), [:a, :b])
        @test_throws ErrorException vcat(td, td2)

        ## unequal indices
        td = $(t)(rand(10,2), [:a, :b])
        td2 = $(t)(rand(10,2), [:a, :b],
                   [Date(2010,1,1):Date(2010,1,10)])
        @test_throws ErrorException vcat(td, td2)

        ## vcat
        df = DataFrame()
        df[:a] = @data([0.4, 0.3])
        df[:b] = @data([0.3, 0.8])

        dats1 = [Date(2014,1,1):Date(2014,1,2)]
        dats2 = [Date(2014,1,3):Date(2014,1,4)]
        td = $(t)(df, dats1)
        td2 = $(t)(df, dats2)
        td3 = $(t)(vcat(df, df), [dats1, dats2])
        td4 = $(t)(vcat(df, df, df), [dats1, dats2, dats1])

        ## multiple arguments
        @test isequal(vcat(td, td2), td3)
        @test isequal(vcat(td, td2, td), td4)

        ## throw error for unequal types
        td = Timedata(rand(4, 3))
        tm = Timematr(rand(4, 3))
        @test_throws ErrorException vcat(tm, td)
        @test_throws ErrorException vcat(tm, td, td)

    end)
end

##################
## flipud tests ##
##################

for t in allTypes
    eval(quote

        nams = [:a, :b]
        vals = [0.2 0.4; 0.6 0.8]
        df = composeDataFrame(vals, nams)
        df2 = composeDataFrame(flipud(vals), nams)
        td = $(t)(df)
        expFlipped = $(t)(df2, flipud(idx(td)))

        @test flipud(td) == expFlipped
        @test isequal(expFlipped, flipud(td))

        @test flipud(flipud(td)) == td
        @test isequal(flipud(flipud(td)), td)

    end)
end

##########
## narm ##
##########

df = DataFrame()
df[:a] = @data([3, NA])
df[:b] = @data([4, NA])
td = Timedata(df)
tn = Timenum(df)
tm = Timematr(rand(2, 3))

td = Timedata(df)
dfShort = DataFrame()
dfShort[:a] = @data([3])
dfShort[:b] = @data([4])
tdShort = Timedata(dfShort)
tnShort = Timenum(dfShort)
tm = Timematr(rand(2, 3))

@test narm(td) == tdShort
@test narm(tn) == tnShort
@test narm(tm) == tm

########################
## asTd / asTm / asTn ##
########################

## test isaRowVector
@test isaRowVector([1 2 3])
@test !isaRowVector([1, 2, 3])

kk = Array(Float64, 1, 1)
@test isaRowVector(kk)

@test !isaRowVector([1])
@test_throws ErrorException isaRowVector(1)

## test asArrayOfEqualDimensions
td = Timedata(rand(2, 2))
@test asArrayOfEqualDimensions([1 2], td) == [1 2; 1 2]
@test asArrayOfEqualDimensions([1, 2], td) == [1 1; 2 2]

kk = Array(Int, 1, 1)
kk[1] = 1
@test asArrayOfEqualDimensions([1], td[1, 1]) == kk
@test asArrayOfEqualDimensions(1, td[1, 1]) == kk
@test asArrayOfEqualDimensions(1, td) == [1 1; 1 1]
@test_throws ErrorException asArrayOfEqualDimensions([3 4; 5 6], td[1, 1])

df = DataFrame(x1 = [1, 2], x2 = [1, 2])
@test TimeData.asTd([1, 2], td) == TimeData.Timedata(df)

df = DataFrame(x1 = [1, 1], x2 = [2, 2])
@test TimeData.asTd([1 2], td) == TimeData.Timedata(df)

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Abaqus
using JuliaFEM.Testing
using DataFrames

# to turn on automatic file download, set
# ENV["ABAQUS_DOWNLOAD_URL"] = "http://<domain>:2080/v2016/books/eif"
# if don't want to download all stuff to current directory, set also
# ENV["ABAQUS_DOWNLOAD_DIR"] = "/tmp"

""" Run test, return true if simulation is succesfull, i.e. no errors raise
during parsing .inp file or execution of model. This doesn't mean that results
are meaningful; they must be checked in separately. Running model only verifies
that no catastrophic failures happen during file parsing. """
function abaqus_run_test(name)
    return_code = abaqus_run_model(name; fetch=true, verbose=true)
    return_code == 0 && return true
    return false
end

@testset "JuliaFEM-ABAQUS interface" begin
  @testset "1 Element Verification" begin
    @testset "1.2 Eigenvalue tests" begin
      @testset "1.2.1 Eigenvalue extraction for single unconstrained elements" begin
        @testset "Acoustic elements" begin
          @testset "AC1D2 elements." begin
#           abaqus_run_test("ec12afe1") || return
          end
        end
        @testset "Three-dimensional continuum elements" begin
          @testset "C3D10 elements." begin
#           abaqus_run_test("ec3asfe1") || return
          end
        end
      end
    end
    @testset "1.3 Simple load tests" begin
      @testset "1.3.1 Membrane loading of plane stress, plane strain, membrane, and shell elements" begin
        @testset "CPS4 elements." begin
#         abaqus_run_test("ecs4sfs1") || return
        end
      end
      @testset "1.3.3 Three-dimensional solid elements" begin
        @testset "C3D8 elements." begin
          abaqus_run_test("ec38sfs2") || return
          res = abaqus_open_results("ec38sfs2")

node_output1 = wsv"""
NODE           U1          U2          U3          COOR1       COOR2       COOR3
1     -2.0000E-33 -2.0000E-33 -2.0000E-33          0.000       0.000       0.000
2     -2.6667E-05 -1.0000E-33 -1.7333E-04          2.000       0.000       0.000
3     -2.0000E-04 -2.6667E-05 -1.7333E-04          2.000       2.000       0.000
4     -1.7333E-04 -2.6667E-05 -1.0000E-33          0.000       2.000       0.000
5     -3.6777E-48 -8.6667E-05 -1.3333E-05          0.000       0.000       1.000
6     -2.6667E-05 -8.6667E-05 -1.8667E-04          2.000       0.000       1.000
7     -2.0000E-04 -1.1333E-04 -1.8667E-04          2.000       2.000       1.000
8     -1.7333E-04 -1.1333E-04 -1.3333E-05          0.000       2.000       1.000
"""

node_output_2 = wsv"""
NODE          RF1         RF2         RF3            CF1         CF2         CF3
1        1500.000    1500.000    1000.000          0.000       0.000       0.000
2           0.000     500.000       0.000       1500.000       0.000       0.000
3           0.000       0.000       0.000        500.000     500.000   -1000.000
4           0.000       0.000       0.000        500.000    1500.000       0.000
5        -500.000       0.000       0.000          0.000    -500.000    1000.000
6           0.000       0.000       0.000       -500.000   -1500.000       0.000
7           0.000       0.000       0.000      -1500.000   -1500.000   -1000.000
8           0.000       0.000       0.000      -1500.000    -500.000       0.000
"""
          #= to check also results:
          side, opts = read_result(xdmf, "SECTION/side")
          @test isapprox(side["SOFM"], 3464.0)
          @test isapprox(side["SOF1"], 2000.0)
          @test isapprox(side["SOF2"], 2000.0)
          @test isapprox(side["SOF3"], 2000.0)
          @test isapprox(side["SOMM"], 2828.0)
          @test isapprox(side["SOM1"], 0.0)
          @test isapprox(side["SOM2"], 2000.0)
          @test isapprox(side["SOM3"], -2000.0)
          @test isapprox(side["SOAREA"], 2.000)
          @test isapprox(side["SOCF1"], 2/3)
          @test isapprox(side["SOCF2"], 2/3)
          @test isapprox(side["SOCF3"], 1/6)
          =#
        end
      end
    end
  end
end


text:
G = [0 1 0 0 1 0 ;
         1 0 1 0 1 0 ;
         0 1 0 1 0 0 ;
         0 0 1 0 1 1 ;
         1 1 0 1 0 0 ;
         0 0 0 1 0 0 ;
        ]

L = [2 -1 0 0 -1 0;
        -1 3 -1 0 -1 0 ;
        0 -1 2 -1 0 0;
        0 0 -1 3 -1 -1;
        -1 -1 0 -1 3 0;
        0 0 0 -1 0 1]

D = Diagonal([2,3,2,3,3,1])

L_norm = (D^-0.5)*L*(D^-0.5)

@testset "graphs" begin
        @test degree(G) == D
        @test laplacian(G) == L
        @test norm_laplacian(G) ≈ L_norm
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [289*x1^2*x3^2+1088*x1^2*x3*x7+3249*x1^2*x5^2+8550*x1^2*x5*x6+10450*x1^2*x6^2+1024*x1^2*x7^2+1326*x1*x2*x3^2-4212*x1*x2*x3*x6+2496*x1*x2*x3*x7+8094*x1*x2*x4*x5+10650*x1*x2*x4*x6+2550*x1*x3*x4*x7-2788*x1*x3*x4*x8-1350*x1*x3*x5*x6+3162*x1*x3*x5*x8+9600*x1*x4*x6^2-11520*x1*x4*x6*x8+4800*x1*x4*x7^2-5248*x1*x4*x7*x8-2394*x1*x5^3-3150*x1*x5^2*x6+5952*x1*x5*x7*x8+10240*x1*x6*x7*x8+7605*x2^2*x3^2+5041*x2^2*x4^2+6084*x2^2*x8^2+3900*x2*x3^2*x5+5850*x2*x3*x4*x7-6396*x2*x3*x4*x8+7254*x2*x3*x5*x8-2982*x2*x4*x5^2-15288*x2*x4*x8^2-13416*x2*x6^2*x8+2184*x2*x7^2*x8-4212*x2*x8^3+625*x3^2*x5^2+5625*x4^2*x6^2-13500*x4^2*x6*x8+5625*x4^2*x7^2-12300*x4^2*x7*x8+24428*x4^2*x8^2+13950*x4*x5*x7*x8-15252*x4*x5*x8^2+16856*x4*x6^2*x8+12000*x4*x6*x7*x8-2744*x4*x7^2*x8-14400*x4*x7*x8^2+5292*x4*x8^3+441*x5^4+8649*x5^2*x8^2+7396*x6^4-2408*x6^2*x7^2+4644*x6^2*x8^2+196*x7^4+5644*x7^2*x8^2+729*x8^4-4510*x1^2*x6-5460*x1*x2*x3-3120*x1*x2*x8-1750*x1*x3*x5+7040*x1*x3*x6+1368*x1*x4*x5-5700*x1*x4*x6+12920*x1*x4*x8+10032*x1*x5*x7+3440*x1*x6^2+13200*x1*x6*x7-1566*x1*x6*x8-560*x1*x7^2-8000*x1*x7*x8+1080*x1*x8^2+4524*x2*x3*x8+1704*x2*x4^2+12496*x2*x4*x7+8250*x3*x4*x6-9900*x3*x4*x8+1450*x3*x5*x8+8800*x3*x7*x8-504*x4*x5^2-3696*x5^2*x7+4125*x1^2-5500*x1*x3+4158*x1*x6-2030*x1*x8-12012*x2*x3+3025*x3^2-3850*x3*x5+144*x4^2+2112*x4*x7+7744*x7^2+841*x8^2+5390*x1-4466*x8+5929, 1326*x1*x3^2-4212*x1*x3*x6+2496*x1*x3*x7+8094*x1*x4*x5+10650*x1*x4*x6+15210*x2*x3^2+10082*x2*x4^2+12168*x2*x8^2+3900*x3^2*x5+5850*x3*x4*x7-6396*x3*x4*x8+7254*x3*x5*x8-2982*x4*x5^2-15288*x4*x8^2-13416*x6^2*x8+2184*x7^2*x8-4212*x8^3-5460*x1*x3-3120*x1*x8+4524*x3*x8+1704*x4^2+12496*x4*x7-12012*x3, 578*x1^2*x3+1088*x1^2*x7+2652*x1*x2*x3-4212*x1*x2*x6+2496*x1*x2*x7+2550*x1*x4*x7-2788*x1*x4*x8-1350*x1*x5*x6+3162*x1*x5*x8+15210*x2^2*x3+7800*x2*x3*x5+5850*x2*x4*x7-6396*x2*x4*x8+7254*x2*x5*x8+1250*x3*x5^2-5460*x1*x2-1750*x1*x5+7040*x1*x6+4524*x2*x8+8250*x4*x6-9900*x4*x8+1450*x5*x8+8800*x7*x8-5500*x1-12012*x2+6050*x3-3850*x5, 8094*x1*x2*x5+10650*x1*x2*x6+2550*x1*x3*x7-2788*x1*x3*x8+9600*x1*x6^2-11520*x1*x6*x8+4800*x1*x7^2-5248*x1*x7*x8+10082*x2^2*x4+5850*x2*x3*x7-6396*x2*x3*x8-2982*x2*x5^2-15288*x2*x8^2+11250*x4*x6^2-27000*x4*x6*x8+11250*x4*x7^2-24600*x4*x7*x8+48856*x4*x8^2+13950*x5*x7*x8-15252*x5*x8^2+16856*x6^2*x8+12000*x6*x7*x8-2744*x7^2*x8-14400*x7*x8^2+5292*x8^3+1368*x1*x5-5700*x1*x6+12920*x1*x8+3408*x2*x4+12496*x2*x7+8250*x3*x6-9900*x3*x8-504*x5^2+288*x4+2112*x7, 6498*x1^2*x5+8550*x1^2*x6+8094*x1*x2*x4-1350*x1*x3*x6+3162*x1*x3*x8-7182*x1*x5^2-6300*x1*x5*x6+5952*x1*x7*x8+3900*x2*x3^2+7254*x2*x3*x8-5964*x2*x4*x5+1250*x3^2*x5+13950*x4*x7*x8-15252*x4*x8^2+1764*x5^3+17298*x5*x8^2-1750*x1*x3+1368*x1*x4+10032*x1*x7+1450*x3*x8-1008*x4*x5-7392*x5*x7-3850*x3, 8550*x1^2*x5+20900*x1^2*x6-4212*x1*x2*x3+10650*x1*x2*x4-1350*x1*x3*x5+19200*x1*x4*x6-11520*x1*x4*x8-3150*x1*x5^2+10240*x1*x7*x8-26832*x2*x6*x8+11250*x4^2*x6-13500*x4^2*x8+33712*x4*x6*x8+12000*x4*x7*x8+29584*x6^3-4816*x6*x7^2+9288*x6*x8^2-4510*x1^2+7040*x1*x3-5700*x1*x4+6880*x1*x6+13200*x1*x7-1566*x1*x8+8250*x3*x4+4158*x1, 1088*x1^2*x3+2048*x1^2*x7+2496*x1*x2*x3+2550*x1*x3*x4+9600*x1*x4*x7-5248*x1*x4*x8+5952*x1*x5*x8+10240*x1*x6*x8+5850*x2*x3*x4+4368*x2*x7*x8+11250*x4^2*x7-12300*x4^2*x8+13950*x4*x5*x8+12000*x4*x6*x8-5488*x4*x7*x8-14400*x4*x8^2-4816*x6^2*x7+784*x7^3+11288*x7*x8^2+10032*x1*x5+13200*x1*x6-1120*x1*x7-8000*x1*x8+12496*x2*x4+8800*x3*x8-3696*x5^2+2112*x4+15488*x7, -2788*x1*x3*x4+3162*x1*x3*x5-11520*x1*x4*x6-5248*x1*x4*x7+5952*x1*x5*x7+10240*x1*x6*x7+12168*x2^2*x8-6396*x2*x3*x4+7254*x2*x3*x5-30576*x2*x4*x8-13416*x2*x6^2+2184*x2*x7^2-12636*x2*x8^2-13500*x4^2*x6-12300*x4^2*x7+48856*x4^2*x8+13950*x4*x5*x7-30504*x4*x5*x8+16856*x4*x6^2+12000*x4*x6*x7-2744*x4*x7^2-28800*x4*x7*x8+15876*x4*x8^2+17298*x5^2*x8+9288*x6^2*x8+11288*x7^2*x8+2916*x8^3-3120*x1*x2+12920*x1*x4-1566*x1*x6-8000*x1*x7+2160*x1*x8+4524*x2*x3-9900*x3*x4+1450*x3*x5+8800*x3*x7-2030*x1+1682*x8-4466];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
@muladd function perform_step!(integrator,
                               cache::Union{ImplicitEMConstantCache,
                                            ImplicitEulerHeunConstantCache,
                                            ImplicitRKMilConstantCache},
                                            f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack nlsolver = cache
  @unpack uf = nlsolver
  alg = unwrap_alg(integrator, true)
  theta = alg.theta
  alg.symplectic ? a = dt/2 : a = theta*dt
  repeat_step = false
  if isnewton(nlsolver)
    uf.t = t
    J = update_W!(integrator, cache, a, repeat_step)
  end

  # TODO: Stochastic extrapolants?
  u = uprev

  L = integrator.g(uprev,p,t)
  ftmp = integrator.f(uprev,p,t)
  gtmp = L.*integrator.W.dW

  if typeof(cache) <: ImplicitEulerHeunConstantCache
    utilde = uprev + gtmp
    gtmp = ((integrator.g(utilde,p,t) + L)/2)*integrator.W.dW
  end

  if typeof(cache) <: ImplicitRKMilConstantCache || integrator.opts.adaptive == true
    if alg_interpretation(alg) == :Ito ||
       typeof(cache) <: ImplicitEMConstantCache
      K = @.. uprev + dt * ftmp
      utilde =  K + L*integrator.sqdt
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      mil_correction = ggprime .* (integrator.W.dW.^2 .- dt)./2
      gtmp += mil_correction
    elseif alg_interpretation(alg) == :Stratonovich ||
           typeof(cache) <: ImplicitEulerHeunConstantCache
      utilde = uprev + L*integrator.sqdt
      ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
      mil_correction = ggprime.*(integrator.W.dW.^2)./2
      gtmp += mil_correction
    end
  end

  if alg.symplectic
    z = zero(u) # constant extrapolation, justified by ODE IM
  else
    z = dt*ftmp # linear extrapolation
  end
  nlsolver.z = z

  nlsolver.c = a
  if alg.symplectic
    # u = uprev + z then  u = (uprev+u)/2 = (uprev+uprev+z)/2 = uprev + z/2
    #u = uprev + z/2 + gtmp/2
    tmp = uprev + gtmp/2
  else
    #u = uprev + dt*(1-theta)*ftmp + theta*z + gtmp
    tmp = uprev + dt*(1-theta)*ftmp + gtmp
  end
  nlsolver.tmp = tmp
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return nothing

  if alg.symplectic
    u = uprev + z + gtmp
  else
    #u = uprev + dt*(1-theta)*ftmp + theta*z + gtmp
    u = tmp + theta*z
  end

  if integrator.opts.adaptive

    if !isnewton(nlsolver)
      is_compos = isa(integrator.alg, StochasticDiffEqCompositeAlgorithm)
      J = calc_J(nlsolver, integrator,cache,is_compos)
    end

    Ed = _reshape(dt*(J*_vec(ftmp))/2, axes(ftmp))
    if typeof(cache) <: Union{ImplicitEMConstantCache,ImplicitEulerHeunConstantCache}
        En = mil_correction
    else
        En = integrator.opts.internalnorm.(integrator.W.dW.^3, t) .*
             integrator.opts.internalnorm.(ggprime, t).^2 ./ 6
    end

    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end

  integrator.u = u
  return nothing
end

@muladd function perform_step!(integrator,
                               cache::Union{ImplicitEMCache,
                                            ImplicitEulerHeunCache,
                                            ImplicitRKMilCache},
                               f=integrator.f)
  @unpack t,dt,uprev,u,p = integrator
  @unpack gtmp,gtmp2,nlsolver = cache
  @unpack du1,dz,z,k,tmp = nlsolver
  J = (isnewton(nlsolver) ? nlsolver.cache.J : nothing)
  alg = unwrap_alg(integrator, true)
  alg.symplectic ? a = dt/2 : a = alg.theta*dt
  dW = integrator.W.dW
  mass_matrix = integrator.f.mass_matrix
  theta = alg.theta

  repeat_step = false

  if integrator.success_iter > 0 && !integrator.u_modified && alg.extrapolant == :interpolant
    current_extrapolant!(u,t+dt,integrator)
  elseif alg.extrapolant == :linear
    @.. u = uprev + integrator.fsalfirst*dt
  else # :constant
    copyto!(u,uprev)
  end

  update_W!(integrator, cache, a, repeat_step)

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

  # Handle noise computations

  integrator.g(gtmp,uprev,p,t)
  integrator.f(tmp,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    @.. gtmp2 = gtmp*dW
  else
    mul!(gtmp2,gtmp,dW)
  end

  if typeof(cache) <: ImplicitEulerHeunCache
    gtmp3 = cache.gtmp3
    @.. z = uprev + gtmp2
    integrator.g(gtmp3,z,p,t)
    @.. gtmp = (gtmp3 + gtmp)/2
    if is_diagonal_noise(integrator.sol.prob)
      @.. gtmp2 = gtmp*dW
    else
      mul!(gtmp2,gtmp,dW)
    end
  end

  if typeof(cache) <: ImplicitRKMilCache
    gtmp3 = cache.gtmp3
    if alg_interpretation(alg) == :Ito
      @.. z = uprev + dt * tmp + integrator.sqdt * gtmp
      integrator.g(gtmp3,z,p,t)
      @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation
      @.. gtmp2 += gtmp3*(dW.^2 - dt)/2
    elseif alg_interpretation(alg) == :Stratonovich
      @.. z = uprev + integrator.sqdt * gtmp
      integrator.g(gtmp3,z,p,t)
      @.. gtmp3 = (gtmp3-gtmp)/(integrator.sqdt) # ggprime approximation
      @.. gtmp2 += gtmp3*(dW.^2)/2
    end
  end

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

  if alg.symplectic
    @.. z = zero(eltype(u)) # Justified by ODE solvers, constrant extrapolation when IM
  else
    @.. z = dt*tmp # linear extrapolation
  end

  nlsolver.c = a
  if alg.symplectic
    #@.. u = uprev + z/2 + gtmp2/2
    @.. tmp = uprev + gtmp2/2
  else
    #@.. u = uprev + dt*(1-theta)*tmp + theta*z + gtmp2
    @.. tmp = uprev + dt*(1-theta)*tmp + gtmp2
  end
  z = nlsolve!(integrator, cache)

  if alg.symplectic
    @.. u = uprev + z + gtmp2
  else
    @.. u = tmp + theta*z
  end

  if integrator.opts.adaptive

    if has_Wfact(f)
      # This means the Jacobian was never computed!
      f.jac(J,uprev,p,t)
    end

    mul!(vec(z),J,vec(tmp))
    @.. k = dt*dt*z/2

    # k is Ed
    # dz is En
    if typeof(cache) <: Union{ImplicitEMCache,ImplicitEulerHeunCache}
      dW_cache = cache.dW_cache
      if !is_diagonal_noise(integrator.sol.prob)
        g_sized = norm(gtmp,2)
      else
        g_sized = gtmp
      end

      if typeof(cache) <: ImplicitEMCache
        @.. z = uprev + dt * tmp + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @.. dW_cache = dW.^2 - dt
          diff_tmp = integrator.opts.internalnorm(dW_cache, t)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @.. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @.. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2 - dt)
        end

      elseif typeof(cache) <: ImplicitEulerHeunCache
        @.. z = uprev + integrator.sqdt * g_sized

        if !is_diagonal_noise(integrator.sol.prob)
          integrator.g(gtmp,z,p,t)
          g_sized2 = norm(gtmp,2)
          @.. dW_cache = dW.^2
          diff_tmp = integrator.opts.internalnorm(dW_cache, t)
          En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
          @.. dz = En
        else
          integrator.g(gtmp2,z,p,t)
          g_sized2 = gtmp2
          @.. dz = (g_sized2-g_sized)/(2integrator.sqdt)*(dW.^2)
        end

      end

    elseif typeof(cache) <: ImplicitRKMilCache
      # gtmp3 is ggprime
      @.. dz = integrator.opts.internalnorm(dW^3, t)*integrator.opts.internalnorm(gtmp3, t)^2 / 6
    end

    calculate_residuals!(tmp, k, dz, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.delta,
                         integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)

  end
end


text:
# Starting point for tuples example for Learning Julia

# a tuple is a sequence of values and is immutable once defined
# like an array, it can contain different types
tup1 = (1, 3.0, "five", 7, '9')
println("The third element of tup1 is ", tup1[3])
println("The second to fourth elements of tup1 is ", tup1[2:4], "\n")

# Tuples can be iterated over like arrays
for i in tup1
    println(i)
end
println()

# The "in" operator can be used to see if a tuple contains a value
println("Is 3 in tup1? ", 3 in tup1)
println("Is 4 in tup1? ", 4 in tup1)
println()

# tuples can have multiple dimensions
tup2 = ((1, 2), (3, 4), (5, 6))
println("The third element in tup2 is ", tup2[3])
println("The second element of the first element of tup2 is ", tup2[1][2], "\n")

# attempting to modify a tuple will produce an error
# ex. tup1[2] = 11 will result in an error

# Named Tuples provide a way to associate names with tuple values
namedTup = (a = 1, b = "two", c = '3')
println("The element associated with name a in namedTup is ", namedTup.a)

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course