text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()
# using Revise
using LazySets, Optim, Plots
import LinearAlgebra: I

# STRUCTURE AND CONSTRUCTORS:

struct NormalForm
    nPlayers::Int64
    nMoves::Tuple{Vararg{Int64}}
    payoffMatList::Tuple{ Vararg{ Array{<:Real} } }

    function NormalForm(players, nmoves, payofflist)
        if size(payofflist[1])==nmoves && players==ndims(payofflist[1])
            new(players, nmoves, payofflist)
        else
            error("Violates dimension requirements for NormalForm.")
        end
    end

end

function NormalForm(normalFormTable::Array{<:Tuple})
    nplayers=ndims(normalFormTable)
    nmoves=size(normalFormTable)

    payoffmatlist=ntuple(
        i -> [normalTable[k][i] for k in CartesianIndices(normalFormTable)],
        nplayers
    )

    return NormalForm(nplayers, nmoves, payoffmatlist)
end

#CALCULATING VALUES

function pureMinMax2(normalform::NormalForm)
    pure_mm1= minimum(
        maximum(
            normalform.payoffMatList[1];
            dims=1
        ),
        dims=2
    )
    pure_mm2= minimum(
        maximum(
            normalform.payoffMatList[2];
            dims=2
        ),
        dims=1
    )
    return [pure_mm1[1], pure_mm2[1]]
end

function brCorrelatMix(normalform::NormalForm,player::Int64,corrMix::Vector)
    payoff=normalform.payoffMatList[player]
    permutation=Tuple(
        unique( [player; 1:normalform.nPlayers] )
    )
    permutedPayoff=permutedims(payoff,permutation)
    permutedPayoff= reshape(permutedPayoff, ( size(permutedPayoff,1) , :))
    return maximum(permutedPayoff*corrMix)
end

function miniMax_i(normalform::NormalForm, player::Int64)
    uBR_i(μ) = brCorrelatMix(normalform, player, μ)
    parametrize_mix(z)=exp.(z)/ sum(exp.(z))
    obj(z) = uBR_i(parametrize_mix(z))

    z0=zeros( prod(normalform.nMoves) ÷ normalform.nMoves[player])
    solution=optimize(obj, z0)
    return solution.minimum
end

function miniMaxProfile(normalform::NormalForm)
    return [miniMax_i(normalform, i) for i=1:normalform.nPlayers]
end

function randomNormalForm(nMovesList::Tuple{Vararg{<:Real}})
    nplayers=length(nMovesList)
    payoffmatlist= Tuple(rand(Float64, nMovesList) for i in 1:nplayers)
    NormalForm(nplayers, nMovesList, payoffmatlist)
end

function greaterThanSet(minPoint::Vector{<:Real})::HPolyhedron
    HPolyhedron(
        -Matrix{Float64}(I,length(minPoint), length(minPoint)),
        -minPoint
    )
end

#PLOTTING

function plotFeasible2!(p,normalform::NormalForm)
    payoff1=normalform.payoffMatList[1]
    payoff2=normalform.payoffMatList[2]

    payoff_points=[
        Vector{Float64}( [payoff1[i], payoff2[i]] )
        for i in eachindex(payoff1)
    ]
    feasible = VPolygon(payoff_points)
    Plots.plot!(p,
        feasible,
        xlabel="Payoff player 1",
        ylabel="Payoff player 2",
        label="Feasible set"
    );

    Plots.scatter!(p,
        [payoff_points[i][1] for i in eachindex(payoff_points)],
        [payoff_points[i][2] for i in eachindex(payoff_points)],
        label="Pure payoffs"
    );
end

function plotFeasible2(normalform::NormalForm)
    p=Plots.plot(title="Normal Form Payoffs")
    plotFeasible2!(p,normalform)
    display(p)
    return p
end

function plotIR_Set2!(p::Plots.Plot, normalform::NormalForm)
    plot!(
        p,
        greaterThanSet( miniMaxProfile(normalform) ),
        label="IR Set",
        legend=:bottomleft
    )
end

function plotMinimax2!(p::Plots.Plot, normalform::NormalForm)
    mm=miniMaxProfile(normalform)
    Plots.scatter!(p,
        [mm[1]],
        [mm[2]],
        label="Minimax"
    );
end

function plotall(normalform::NormalForm)
    p=plot(
        title="Normal form payoffs"
        )
    plotFeasible2!(p,normalform)
    plotIR_Set2!(p,normalform)
    plotMinimax2!(p,normalform)
    display(p)
    return p
end

#= TO DO:
- generate functions to plot 3D
=#

#= CODE TO PLOT 3D FEASIBLE SET
Using LazySets, GLMakie, Polyhedra

test=CH(
    Singleton([1.,1.,3.]),
    Singleton([2.,1.,4.])
)

test=CH(test,
    Singleton([2., 1.5, 3.])
)

figureplot3d(
    overapproximate(test)
) =#

max_stars_repo_name: vfarinhaluz/JupyterJuliaTemplate

max_issues_repo_name: vfarinhaluz/JupyterJuliaTemplate

max_forks_repo_name: vfarinhaluz/JupyterJuliaTemplate


text:
const a = [984,981,976,950,899,890,887,885,880,800,798,790,777,767,750,701,
    697,688,680,678,650,599,589,567,550,501,9,8,7,6,5,4,3,2,1]
const b = [1090,1080,1074,1065,1057,1056,1047,1041,1041,1038,1025,1013,1008,
    992,991,991,991,978,977,959,945,935,925,925,923,915,908,904,901,901,
    900,897,894,882,880,876,866,854,849,849,833,818,818,812,811,809,798,
    794,793,788,772,763,747,746,743,737,736,734,732,730,728,718,714,713,
    706,701,699,691,690,689,681,672,663,656,654,653,652,651,646,644,640,
    637,637,635,634,633,630,625,621]

struct DigitSet
    x::UInt32
end

function DigitSet(x::Integer)
    set = DigitSet()
    absolute = unsigned(abs(x))
    while absolute > 0
        set = push!(set, rem(absolute, 10))
        absolute = oftype(absolute, div(absolute, 10))
    end
    return set
end

DigitSet() = DigitSet(zero(UInt32))
Base.push!(x::DigitSet, y::Integer) = DigitSet(x.x | (UInt32(1) << y))
Base.union(x::DigitSet, y::DigitSet) = DigitSet(x.x | y.x)
Base.hash(x::DigitSet, h::UInt) = hash(x.x, h)
isdisjoint(x::DigitSet, y::DigitSet) = x.x & y.x == zero(UInt32)
Base.issubset(x::DigitSet, y::DigitSet) = (x.x & ~y.x) == zero(UInt32)

function linear_deduplicate(numbers)
    bydigits = Dict{DigitSet, eltype(numbers)}()
    for number in numbers
        set = DigitSet(number)
        existing = get(bydigits, set, typemin(Int))
        if number > existing
            bydigits[set] = number
        end
    end
    return collect(bydigits)
end

function square_deduplicate(numberpairs)
    accepted = Int[]
    for (set, number) in numberpairs
        add = true
        for (set_, number_) in numberpairs
            if issubset(set_, set) && number_ > number
                add = false
                break
            end
        end
        if add
            push!(accepted, number)
        end
    end
    # We sort to prune the tree of recusive possibilities earlier
    sort!(accepted, rev=true)
    return accepted
end

abstract type Solution end

struct SmallSolution <: Solution
    sum::Int64
    x::UInt64
end

mutable struct BigSolution <: Solution
    sum::Int64
    x::Vector{Int}
end

SmallSolution(x) = SmallSolution(x, x)
SmallSolution() = SmallSolution(0)
BigSolution() = BigSolution(0, Int[])

Base.copy(x::SmallSolution) = x
Base.copy(x::BigSolution) = BigSolution(x.sum, copy(x.x))

function pushindex!(x::SmallSolution, numbers, index)
    @inbounds n = numbers[index]
    return SmallSolution(n+x.sum, x.x | UInt64(1) << (index - 1))
end

function pushindex!(x::BigSolution, numbers, index)
    n = numbers[index]
    x.sum += n
    push!(x.x, n)
    return x
end

Base.collect(x::BigSolution, nums::Vector{T}) where T<:Integer = x.x

function Base.collect(x::SmallSolution, nums::Vector{T}) where T<:Integer
    result = Int[]
    @inbounds for i in 1:min(64, unsigned(length(nums)))
        if x.x >>> (i-1) & UInt64(1) == UInt64(1)
            push!(result, nums[i])
        end
    end
    return result
end

function _recurse(current::Solution, numbers, sets, index, set, solution)
    if index > length(numbers)
        if current.sum > solution[].sum
            solution[] = current
        end
        return nothing
    end
    @inbounds nextset = sets[index]
    _recurse(current, numbers, sets, index+1, set, solution)
    if isdisjoint(set, nextset)
        nextcurrent = pushindex!(copy(current), numbers, index)
        _recurse(nextcurrent, numbers, sets, index+1, union(set, nextset), solution)
    end
    return nothing
end

function recurse(numbers)
    numbersarr = collect(numbers)
    T = length(numbersarr) ≤ 64 ? SmallSolution : BigSolution
    solution = Ref(T())
    _recurse(T(), numbersarr, map(DigitSet, numbersarr), 1, DigitSet(), solution)
    return sort!(collect(solution[], numbers))
end

function solve(numbers::Vector{Int})
    filtered = linear_deduplicate(numbers)
    filtered2 = square_deduplicate(filtered)
    return recurse(filtered2)
end

ext: jl

max_stars_repo_name: jakobnissen/play

max_issues_repo_name: jakobnissen/play

max_forks_repo_name: jakobnissen/play


text:
using SpecialFunctions: gamma, erf
using ForwardDiff: derivative
using QuadGK: quadgk
"""
    StudentizedRange(k,ν)

The *Studentized range distribution* has probability density function

```math
f(q; k, \\nu) = \\frac{\\sqrt{2\\pi}k(k-1)\\nu^{\\nu/2}}{\\Gamma(\nu/2)2^{\\nu/2 - 1}}
\\int_{0}^{\\infty}\\varphi(\\sqrt{\\nu}x)\\Big[\\int_{-\\infty}^{\\infty}\\varphi(u)
\\varphi(u-qx)(\\Phi(u)\\Phi(u-qx)^{k-2})du\\Big]dx, \\quad x > 0, \\quad k > 1
```

```julia
StudentizedRange(n, ν)    # Studentized range distribution with number of samples n and degrees of freedom ν

params(d)          # Get the parameters, i.e. (k, ν)
dof(d)             # Get the degrees of freedom, i.e. ν
```

External links

* [Studentized range distribution on Wikipedia](https://en.wikipedia.org/wiki/Studentized_range_distribution)

"""
struct StudentizedRange{T<:Real} <: ContinuousUnivariateDistribution
    k::T
    ν::T
    coeff_pdf::T
    coeff_cdf::T
end

function StudentizedRange{T}(k, ν) where T
    (k, ν) = (Float64(k), Float64(ν))
    coeff_pdf = (√(2π) * k * (k-1) * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
    coeff_cdf = (k * ν^(ν/2)) / (gamma(ν/2) * 2^(ν/2 - 1))
    return StudentizedRange(k, ν, coeff_pdf, coeff_cdf)
end

StudentizedRange(k::T, ν::T) where {T<:Real} = StudentizedRange{T}(k, ν)
StudentizedRange(k::Real, ν::Real) = StudentizedRange(promote(k, ν)...)
StudentizedRange(k::Integer, ν::Integer) = StudentizedRange(Float64(k), Float64(ν))
StudentizedRange(a) = StudentizedRange(a,a)
StudentizedRange() = StudentizedRange(2,2)

@distr_support StudentizedRange 0.0 Inf

### Conversions
function convert(::Type{StudentizedRange{T}}, k::Real, ν::Real) where T<:Real
    StudentizedRange(T(k), T(ν))
end
function convert(::Type{StudentizedRange{T}}, d::StudentizedRange{S}) where {T <: Real, S <: Real}
    StudentizedRange(T(d.k), T(d.ν))
end

### Parameters

params(d::StudentizedRange) = (d.k, d.ν)
dof(d::StudentizedRange) = d.ν

### Statistics

mean(d::StudentizedRange{T}) = T(NaN)
mode(d::StudentizedRange{T}) = T(NaN)
var(d::StudentizedRange{T}) = T(NaN)
skewness(d::StudentizedRange{T}) = T(NaN)
entropy(d::StudentizedRange{T}) = T(NaN)

### Evaluation

# Helper functions for cdf and pdf of standard normal.
function 𝚽(x)
    return (1+erf(x / √2)) / 2
end

function ϕ(x)
    return derivative(𝚽, x)
end

function cdf(d::StudentizedRange, q)
    function outer(x)

        function inner(u)
            return ϕ(u) * (𝚽(u) - 𝚽(u - q*x))^(d.k-1)
        end
        inner_part = quadgk(inner, -Inf, Inf)[1]
        return inner_part * x^(d.ν-1) * exp(-x^2*d.ν / 2)
    end
    integral = quadgk(outer, 0.0, Inf)[1]
    return integral * d.coeff_cdf
end

function pdf(d::StudentizedRange, q)

    function outer(x)

        function inner(u)
            return ϕ(u) * ϕ(u - q*x) * (𝚽(u) - 𝚽(u - q*x))^(d.k-2)
        end
        inner_part = quadgk(inner, -Inf, Inf)[1]
        return inner_part * x^d.ν * ϕ(x*√(d.ν))
    end
    integral = quadgk(outer, 0.0, Inf)[1]
    return integral * d.coeff_pdf
end

logpdf(d::StudentizedRange, q) = log(pdf(d, q))

# To get quantile to work correctly I had to implement my quick naive version of
# the bisection method. I'm not sure why, but trying to do this with Roots.jl
# was WAY too slow. Like 30+ seconds.
function simple_bisection(f::Function, brackets, abstol=10.0^-6, maxeval=1e3)
    if brackets[1] > brackets[2]
        xmax, xmin = brackets
    else
        xmin, xmax = brackets
    end
    @assert f(xmin) * f(xmax) < 0

    a = xmin
    b = xmax
    error = 1
    numeval = 0
    while error > abstol

        numeval += 1
        if numeval > maxeval break end

        fnew = f((a+b)/2)
        error = abs(fnew)
        if fnew * f(a) < 0
            b = (a+b)/2
        elseif fnew * f(b) < 0
            a = (a+b)/2
        else
            throw(BoundsError("Algorithm failed to converge. Sorry :("))
        end

    end
    return (a+b)/2
end

function quantile(d::StudentizedRange, x)
    @assert 0.0 <= x < 1.0
    if x == 0.0 return 0.0 end

    return simple_bisection(y -> cdf(d, y) - x, [0.0, 1000.0])
end

median(d::StudentizedRange) = quantile(d, 0.5)


text:
function glue_piece!(cur_player::AbstractPlayer)
  cur_piece = cur_player.piece

  isnull(cur_piece) && return

  cur_player.clock.lock =
    Nullable{Base.Random.UUID}()

  cur_player.clock.drag =
    Nullable{Base.Random.UUID}()

  cur_player.state.can_hold |= true

  cur_js = """
    \$(".js-shadow-piece, .js-active-piece").removeClass();

    var tmp_cell;
  """

  cur_grid = cur_player.grid

  cur_rows = SortedSet()

  for cur_block in cur_piece.blocks
    (cur_row, cur_col) = calc_block_coords(cur_block)

    cur_js *= """
      tmp_cell = \$(".cs-main-area .cs-row-$(cur_row) td:nth-child($(cur_col))");

      tmp_cell.addClass("cs-color cs-$(cur_piece.color)");
    """

    ( cur_grid.height < cur_row ) &&
      ( cur_grid.height = cur_row )

    push!(cur_rows, cur_row)
    cur_player.grid.table[cur_row, cur_col] = cur_piece.color
  end

  cleared_rows = Int[]

  for cur_row in cur_rows
    any(
      cur_cell -> ( cur_cell == "" ),
      cur_player.grid.table[cur_row, :]
    ) && continue

    push!(cleared_rows, cur_row)
  end

  cleared_count = length(cleared_rows)

  score!(cur_player, "glue", cleared_count)

  cur_player.glues += 1

  if iszero(cleared_count)
    cur_player.combo = 0

    is_valid = false

    for cur_block in cur_piece.blocks
      (cur_row, cur_col) = calc_block_coords(cur_block)
      ( cur_row > cur_grid.rows ) && continue

      is_valid = true
      break
    end

    is_valid || return end_round(cur_player)
  else
    @assert cleared_count < 5

    ( cur_player.combo > 0 ) &&
      score!(cur_player, "combo", cleared_count)

    cur_player.combo += 1

    score!(cur_player, "clear", cleared_count)

    if cleared_count == 4
      cur_player.chain = "tetris"
    else
      cur_player.chain = ""
    end

    cur_js *= """
      \$(".js-score-text").text("$(lpad(cur_player.score, 8, "0"))");
    """

    cur_js *= clear_rows!(cur_player, cleared_rows)
  end

  if is_repl
    cur_string = []
    for cur_row in 1:cur_player.grid.rows
      for cur_col in 1:cur_player.grid.cols
        cur_color = cur_player.grid.table[cur_row, cur_col]
        ( cur_color == "" ) && ( cur_color = "invisible" )

        push!(cur_string, "\x1b[$(cur_player.grid.rows-cur_row+4);$(29+2*cur_col)H")
        push!(cur_string, crayon_dict[cur_color],)
        push!(cur_string, "  ")
        push!(cur_string, inv(crayon_dict[cur_color]))
      end
    end

    push!(cur_string, "\x1b[u")
    print(cur_string...)
  end

  tetris_js(cur_player.game.scope, JSString(cur_js))

  cur_piece.owner = Nullable{AbstractContainer}()
  cur_player.piece = Nullable{AbstractPiece}()

  if cur_player.level <= max_level
    drop_clock(cur_player)
  else
    end_round(cur_player)
  end

end


text:
@testset "tsne()" begin
    @testset "API tests" begin
        iris = dataset("datasets", "iris")
        X = convert(Matrix{Float64}, iris[:, 1:4])
        Y = tsne(X, 2, -1, 10, 15, verbose=true)
        @test size(Y) == (150, 2)
        Y = tsne(X, 3, -1, 10, 15, verbose=false)
        @test size(Y) == (150, 3)
        tsne(X, 2, -1, 10, 15, verbose=true, progress=false)
        tsne(X, 2, -1, 10, 15, verbose=false, progress=false)
        Y = tsne(X, 2, -1, 10, 15, pca_init=true, cheat_scale=1.0, progress=false)
        @test size(Y) == (150, 2)
        Y, beta, kldiv = tsne(X, 2, -1, 10, 15, pca_init=true, cheat_scale=1.0, progress=false, extended_output=true)
        @test size(Y) == (150, 2)
        @test beta isa AbstractVector
        @test length(beta) == 150
        @test isfinite(kldiv)

        @testset "distance = true" begin
            @test_throws ArgumentError tsne(X, 3, -1, 10, 15, distance=true, verbose=false)
            XX = pairwise(CosineDist(), X')
            Y = tsne(XX, 3, -1, 10, 15, distance=true, verbose=false)
            @test size(Y) == (150, 3)
        end

        @testset "distance isa Function" begin
            Y = tsne(X, 3, -1, 10, 15, distance=cityblock, verbose=false)
            @test size(Y) == (150, 3)
        end

        @testset "X isa Vector, distance isa Function" begin
            Y = tsne([view(X, i, :) for i in 1:size(X, 1)],
                     3, -1, 10, 15, distance=cityblock, verbose=false)
            @test size(Y) == (150, 3)
        end

        @testset "distance isa Distances.Metric" begin
            Y = tsne(X, 3, -1, 10, 15, distance=Minkowski(0.5), verbose=false)
            @test size(Y) == (150, 3)
        end
    end

    @testset "Iris dataset" begin
        iris = dataset("datasets", "iris")
        X = convert(Matrix{Float64}, iris[:, 1:4])
        # embed in 3D
        Y3d = tsne(X, 3, -1, 100, 15, progress=false)
        @test size(Y3d) == (150, 3)
        # embed in 2D
        Y2d = tsne(X, 2, 50, 50, 20, progress=false)
        @test size(Y2d) == (150, 2)
    end

    @testset "MNIST.traindata() dataset" begin
        Random.seed!(345678)
        train_data, labels = MLDatasets.with_accept(true) do
                MNIST.traindata(Float64)
        end
        X_labels = labels[1:2500] .+ 1
        X = reshape(permutedims(train_data[:, :, 1:2500], (3, 1, 2)),
                    2500, size(train_data, 1)*size(train_data, 2))
        X .-= mean(X, dims=1)
        X ./= std(X, dims=1)

        Y = tsne(X, 2, 50, 2000, 20, progress=true)
        @test size(Y) == (2500, 2)
    end
end


text:
export latin, check_latin

function funky_mod(i::Int, n::Int)
    v = mod(i,n)
    return v==0 ? n : v
end

"""
`latin(n)` returns a simple `n`-by-`n` Latin square. More generally,
`latin(n,a,b)` generates a Latin square whose `i,j`-entry is
`a*(i-1) + b*(j-1) + 1` (wrapping around `n`, of course).

*Note*: If parameters `n,a,b` do not generate a legitimate Latin square
an error is thrown.
"""
function latin(n::Int, a::Int=1, b::Int=1)::Matrix{Int}
    A = zeros(Int,n,n)
    for i=1:n
        for j=1:n
            A[i,j] = funky_mod(a*(i-1) + b*(j-1) + 1, n)
        end
    end
    @assert check_latin(A) "Parameters n=$n, a=$a, and b=$b do not generate a Latin square"
    return A
end

"""
`check_latin(A)` checks if `A` is a Latin square.
(The entries must be chosen from `1:n` when `A` is an `n`-by-`n`
matrix).
"""
function check_latin(A::Matrix{Int})::Bool
    r,c = size(A)
    if r != c
        return false
    end

    vals = collect(1:r)

    # check rows
    for i=1:r
        if sort(A[:,i]) != vals
            return false
        end
    end

    # check cols
    for i=1:r
        if sort(A[i,:]) != vals
            return false
        end
    end

    return true
end


text:
using GeometricAlgebra

@testset "scalar multiplication" begin
	a1 = Blade{(1,1)}(1, 0b01)
	a6 = Blade{(1,1)}(6, 0b01)
	@test 0a1 == 0a6
	@test a1*6 == a6
	@test 6*a1 == a6
	@test (3*a1)*2 == a6
	@test a6*2 == 3(2a1)*2
	@test 3.0a1 == 0.5a6
	@test 6im*a1 == im*a6
	@test a6/2 == 3a1
	@test a1 == 6.0\a6

	u = Multivector{(1,1,1),1}([1, 2, 4])
	v = Multivector{(1,1,1),1}([12, 24, 48])
	@test 12u == v
	@test 4(u*3) == v
	@test 120f0*u == BigInt(10)v
	@test u/2 == v/24
	@test 6\v == 2u

	A = MixedMultivector{(1,1)}([3, 0, 10, 7])
	B = MixedMultivector{(1,1)}([21, 0, 70, 49])
	@test 7A == B
	@test (7A)*2 == 2B
	@test (7 + 0im)A == (1//1)B
	@test B/7 == A
	@test 1\A == 7\B

end

@testset "addition, subtraction" begin
	sig = (1, 1, 1)
	x, y, z = basis(sig)
	@test x + y == Multivector{sig,1}([1, 1, 0])
	@test x + y == y + x
	@test (x + y) + z == y + x + z
	@test x + 2 == MixedMultivector{sig}([2, 1, 0, 0, 0, 0, 0, 0])
	@test x + 0 == x
	@test x + 3.0y == (3//1)y + x

	u = x + 2y
	@test -x == (-1)x
	@test -u == (-1)u
	@test +u == -(-u)
	@test x - y == x + (-y)
	@test u - 1 == -1 + u
	@test u - x == u - (0 + x)
end

@testset "geometric product" begin
	t, x, y, z = basis((-1, 1, 1, 1))
	@testset "Blade" begin
		@test x*x == y*y == z*z == 1
		@test t*t == -1
		@test x*y == -y*x
		@test (x*y)*(x*y) == -1
		@test (y*z)*(t*x) == t*x*y*z
		@test x*x*x == x
	end

	@testset "Multivector" begin
		u = 3x + 4y
		@test u*x == 3 + 4y*x
		@test 25 == u*u
		@test (t + x)*(t + x) == 0
	end

	@testset "MixedMultivector" begin
		@test (1 + x)*(1 + y) == 1 + x + y + x*y
		@test (1 + x)*(1 - x) == 0
		@test (2 + x*y)*(y + x*y) == 2y + 2x*y + x - 1
	end
end

@testset "multiplicative inverses" begin
	t, x, y, z = basis((-1, 1, 1, 1))
	@test inv(x) == 1/x == x\1 == x
	@test inv(t) == -t
	@test inv(x*y) == y*x
	@test x/x == y\y == 1
	@test (t*x*y*z)/(t*z) == t*x*y*z*t*z == x*y

	@test inv(3x + 4y) == (3x + 4y)/25
	cases = [
		x + 2y + 3z,
		t*x + y*z, # simplest multivector with non-scalar square
		1 + 2x,
		3 + x*y*z,
		1 + t + x*y
	]
	@testset "inv($a)" for a ∈ cases
		@test inv(a)*a ≈ 1
		@test a*inv(a) ≈ 1
	end
end

@testset "integer powers" begin
	t, x, y, z = basis((-1, 1, 1, 1))

	@testset "power_with_scalar_square()" begin
		using GeometricAlgebra: power_with_scalar_square
		pwss = power_with_scalar_square

		@testset "a = $a" for a ∈ float.([x, 2x*y, x*t])
			@testset "a^p, p = $p" for p ∈ [-4:4; 13; 14; 171; -333; 10042]
				a² = a*a
				@assert isscalar(a²)
				@test pwss(a, scalar(a²), p) == prod(fill(p > 0 ? a : inv(a), abs(p)))
			end
		end

	end

	@test x^2 == 1
	@test t^2 == -1
end

@testset "duality operations" begin
	t, x, y, z = basis((-1, 1, 1, 1))
	@test reversion(x*y) == y*x
	@test reversion(42) == 42
	@test ~x == x
	@test ~(x*y + y*z) == y*x + z*y
	@test ~(1 + t*x*y) == 1 - t*x*y
	for a ∈ [x*z + y*z, 1 + x + x*y + x*y*z + t*x*y*z]
		@test ~(~a) == a
	end

	@test involute(4) == 4
	@test involute(4x) == -4x
	@test involute(4x*y) == 4x*y
	@test involute(1 + 4x*y + t) == 1 + 4x*y - t
end

@testset "grade operations" begin
	x, y, z = basis((1, 1, 1))
	@test grade(42) == 0
	@test grade(x) == 1
	@test grade(x*y) == 2
	@test grade(x + y) == 1

	@test grade(42, 1) == 0
	@test grade(x, 1) == x
	@test grade(x, 10) == 0
	@test grade(x + y, 1) == x + y
	@test grade(x + y, 0) == 0
	@test grade(1 + x + y*z, 2) == y*z

	@test x[1] == 1
	@test x[2] == 0
	@test (1 + 3x*y)[2,1] == -3

	@test isscalar(x^2)
	@test isscalar(x^2 + 0)
end

@testset "other products" begin
	t, x, y, z = basis((-1, 1, 1, 1))
	@test x∧x == 0
	@test x∧(y + x) == x*y
	@test (1 + x)∧(y + x) == x + y + x*y
	@test 4∧x == 4x

	@test x⋅y == 0
	@test x⋅x == -t⋅t == 1
	@test z⋅2 == 2z

	u = x - y
	v = t*x + 2y*z
	@test u⋅u == 2
	@test u∧u == 0
	@test u⋅v == -t - 2z
	@test u∧v == 2x*y*z - t*x*y

	@test u∗v == scalar(u*v)
	@test u∗v == scalar(u*v)
end


text:
function matrixfunction(fn::Function, m::Matrix{DoubleFloat{T}}) where {T<:IEEEFloat}
    issquare(m) || throw(ErrorException("matrix must be square"))
    evecs = eigvecs(m)
    rank(evecs) == size(m)[1] || throw(ErrorException("matrix is not diagonalizable"))
    evals = eigenvals(m)
    fnevals = fn.(evals)
    invevecs = inv(evecs)
    diagevals = Diagonal(fnevals)
    result = evecs * diagevals * invevecs
    return result
end

function matrixfunction(fn::Function, m::Matrix{Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
    issquare(m) || throw(ErrorException("matrix must be square"))
    evecs = eigvecs(m)
    rank(evecs) == size(m)[1] || throw(ErrorException("matrix is not diagonalizable"))
    evals = eigenvals(m)
    fnevals = fn.(evals)
    invevecs = inv(evecs)
    diagevals = Diagonal(fnevals)
    result = evecs * diagevals * invevecs
    return result
end

for F in (:sqrt, :cbrt, :log, :exp,
          :sin, :cos, :tan, :csc, :sec, :cot,
          :asin, :acos, :atan, :acsc, :asec, :acot,
          :sinh, :cosh, :tanh, :csch, :sech, :coth,
          :asinh, :acosh, :atanh, :acsch, :asech, :acoth)
  @eval begin
    function $F(m::Matrix{DoubleFloat{T}}) where {T<:IEEEFloat}
        return matrixfunction($F, m)
    end
    function $F(m::Matrix{Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
        return matrixfunction($F, m)
    end
  end
end

function Base.:(^)(m::Matrix{DoubleFloat{T}}, p::Union{IEEEFloat, DoubleFloat{T}}) where {T<:IEEEFloat}
    pw = DoubleFloat{T}(p)
    res = pw * log(m)
    result = exp(res)
    return result
end
function Base.:(^)(m::Matrix{Complex{DoubleFloat{T}}}, p::Union{IEEEFloat, DoubleFloat{T}}) where {T<:IEEEFloat}
    pw = DoubleFloat{T}(p)
    res = pw * log(m)
    result = exp(res)
    return result
end

function Base.:(^)(m::Matrix{DoubleFloat{T}}, p::Union{Complex{T}, Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
    pw = Complex{DoubleFloat{T}}(p)
    res = pw * log(m)
    result = exp(res)
    return result
end
function Base.:(^)(m::Matrix{Complex{DoubleFloat{T}}}, p::Union{Complex{T}, Complex{DoubleFloat{T}}}) where {T<:IEEEFloat}
    pw = Complex{DoubleFloat{T}}(p)
    res = pw * log(m)
    result = exp(res)
    return result
end

# sort to match usual output for `eigvals`

@inline lt_for_sort(a::DoubleFloat{T}, b::DoubleFloat{T}) where {T<:IEEEFloat} = a < b
@inline lt_for_sort(a::Complex{DoubleFloat{T}},b::Complex{DoubleFloat{T}}) where {T<:IEEEFloat} =
    (real(a) < real(b)) || (real(a)==real(b) && imag(a)<imag(b))

function eigenvals(m::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    evals = eigvals(m)
    return sort(evals, lt=lt_for_sort)
end

function eigenvals(m::Array{Complex{DoubleFloat{T}},2}) where {T<:IEEEFloat}
    evals = eigvals(m)
    return sort(evals, lt=lt_for_sort)
end


text:
using Libxc
using Printf

length(ARGS) == 1 || error("Usage: xc-info.jl identifier")

func = Functional(Symbol(ARGS[1]))

for s in (:identifier, :name, :family, :kind)
    @printf "%-10s: %-20s\n" string(s) getproperty(func, s)
end

println()
if is_global_hybrid(func)
    println("This is a global hybrid functional with $(func.exx_coefficient)% " *
            "exact exchange")
elseif is_range_separated(func)
    print("This is a range-separated hybrid functional with range-separation " *
          "constant $(func.omega) and $(100func.alpha + 100func.beta)% short-range " *
          "and $(100func.alpha) long-range exact exchange")
    if :hyb_lcy in func.flags || :hym_camy in func.flags
        println("using the Yukawa kernel.")
    elseif :hyb_lc in func.flags || :hym_cam in func.flags
        println("using the Coulomb kernel.")
    end
else
    println("This is a pure functional with no exact exchange.")
end

println()
println("References:")
for ref in func.references
    println("  *) $(ref.reference)")
    if !isempty(ref.doi)
        println("     doi: $(ref.doi)")
    end
end

println()
println("Implementation has support for:")
:exc in func.flags && println("  *) energy")
:vxc in func.flags && println("  *) 1st derivative (potential)")
:fxc in func.flags && println("  *) 2nd derivative (kernel)")
:kxc in func.flags && println("  *) 3rd derivative")
:lxc in func.flags && println("  *) 4th derivative")


text:
struct AnimationTrack{T}
    name::String
    jstype::String
    frames::Vector{Int}
    values::Vector{T}
end

@with_kw struct AnimationClip
    tracks::Dict{String, AnimationTrack} = Dict{String, AnimationTrack}()
    fps::Int = 30
    name::String = "default"
end

struct Animation
    clips::Dict{Path, AnimationClip}
    default_framerate::Int
end

Animation(fps::Int=30) = Animation(Dict{Path, AnimationClip}(), fps)

function convert_frames_to_video(tar_file_path::AbstractString, output_path::AbstractString="output.mp4"; framerate=60, overwrite=false)
    output_path = abspath(output_path)
    if !isfile(tar_file_path)
        error("Could not find the input file $tar_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        run(unpack_cmd(tar_file_path, tmpdir, ".tar", nothing))
        cmd = `ffmpeg -r $framerate -i %07d.png -vcodec libx264 -preset slow -crf 18`
        if overwrite
            cmd = `$cmd -y`
        end
        cmd = `$cmd $output_path`

        cd(tmpdir) do
            try
                run(cmd)
            catch e
                println("""
Could not call `ffmpeg` to convert your frames into a video.
If you want to convert the frames manually, you can extract the
.tar archive into a directory, cd to that directory, and run:
ffmpeg -r 60 -i %07d.png \\\n\t -vcodec libx264 \\\n\t -preset slow \\\n\t -crf 18 \\\n\t output.mp4""")
                rethrow(e)
            end
        end
    end
    @info("Saved output as $output_path")
    return output_path
end


text:
using MacroTools

# From 1/16/2020 meeting
#
# c1 = Component(A)
# Component(B) # equiv B = Component(B)
#
# x3 = Parameter(a.p1, a.p2, b.p3, default=3, description="asflijasef", visibility=:private)
#
# This creates external param x3, and connects b.p3 and ANY parameter in any child named p1 to it
# AND now no p1 in any child can be connected to anything else. Use Not from the next if you want
# an exception for that
# x3 = Parameter(p1, b.p3, default=3, description="asflijasef", visibility=:private)
#
# x3 = Parameter(p1, p2, Not(c3.p1), b.p3, default=3, description="asflijasef", visibility=:private)
#
# connect(B.p2, c1.v4)
# connect(B.p3, c1.v4)
#
# x2 = Parameter(c2.x2, default=35)
#
# BUBBLE UP PHASE
#
# for p in unique(unbound_parameters)
#   x1 = Parameter(c1.x1)
# end
#
# if any(unbound_parameter) then error("THIS IS WRONG")
#
#
# Expressions to parse in @defcomposite:
#
# 1. name_ = Component(compname_)
# 2. Component(compname_) => (compname = Component(compname_))
# 3. pname_ = Parameter(args__) # args can be: pname, comp.pname, or keyword=value
# 4. connect(a.param, b.var)
#
#

 # splitarg produces a tuple for each arg of the form (arg_name, arg_type, slurp, default)
_arg_name(arg_tup) = arg_tup[1]
_arg_type(arg_tup) = arg_tup[2]
_arg_slurp(arg_tup) = arg_tup[3]
_arg_default(arg_tup) = arg_tup[4]

function _typecheck(obj, expected_type, msg)
    obj isa expected_type || error("$msg must be a $expected_type; got $(typeof(obj)): $obj")
end

"""
    parse_dotted_symbols(expr)

Parse and expression like `a.b.c.d` and return the tuple `(ComponentPath(:a, :b, :c), :d)`,
or `nothing` if the expression is not a series of dotted symbols.
"""
function parse_dotted_symbols(expr)
    global Args = expr
    syms = Symbol[]

    ex = expr
    while @capture(ex, left_.right_) && right isa Symbol
        push!(syms, right)
        ex = left
    end

    if ex isa Symbol
        push!(syms, ex)
    else
        return nothing
    end

    syms = reverse(syms)
    datum_name = pop!(syms)
    return ComponentPath(syms), datum_name
end

#
# Convert @defcomposite "shorthand" statements into Mimi API calls
#
function _parse(expr)
    valid_keys = (:default, :description, :unit)
    result = nothing

    if @capture(expr, newname_ = Component(compname_)) ||
       @capture(expr, Component(compname_))
        # check newname is nothing or Symbol, compname is Symbol
        _typecheck(compname, Symbol, "Referenced component name")

        if newname !== nothing
            _typecheck(newname, Symbol, "Local name for component name")
        end

        newname = (newname === nothing ? compname : newname)
        result = :(Mimi.add_comp!(obj, $compname, $(QuoteNode(newname))))

    elseif @capture(expr, localparname_ = Parameter(args__))
        regargs = []
        keyargs = []

        for arg in args
            if @capture(arg, keywd_ = value_)
                if keywd in valid_keys
                    push!(keyargs, arg)
                else
                    error("Unrecognized Parameter keyword '$keywd'; must be one of $valid_keys")
                end

            elseif @capture(arg, (cname_.pname_ | pname_))
                cname = (cname === nothing ? :(:*) : cname) # wildcard
                push!(regargs, :(obj[$(QuoteNode(cname))] => $(QuoteNode(pname))))

            end
        end
        result = :(Mimi.import_param!(obj, $(QuoteNode(localparname)), $(regargs...);
                                       $(keyargs...)))

    elseif @capture(expr, localvarname_ = Variable(datum_expr_))
        if ((tup = parse_dotted_symbols(datum_expr)) === nothing)
            error("In @defcomposite's Variable(x), x must a Symbol or ",
                  "a dotted series of Symbols. Got :($datum_expr)")
        end
        comppath, varname = tup
        # @info "Variable: $comppath, :$varname"
        _typecheck(localvarname, Symbol, "Local variable name")
        _typecheck(comppath, ComponentPath, "The referenced component")
        _typecheck(varname, Symbol, "Name of referenced variable")

        # import from the added copy of the component, not the template -- thus
        # the lookup of obj[varcomp].
        result = :(Mimi._import_var!(obj, $(QuoteNode(localvarname)), $comppath,
                                     $(QuoteNode(varname))))

    elseif @capture(expr, connect(parcomp_.parname_, varcomp_.varname_))
        # raise error if parameter is already bound
        result = :(Mimi.connect_param!(obj,
                    $(QuoteNode(parcomp)), $(QuoteNode(parname)),
                    $(QuoteNode(varcomp)), $(QuoteNode(varname));
                    # allow_overwrite=false # new keyword to implement
        ))
    else
        error("Unrecognized composite statement: $expr")
    end
    return result
end

# TBD: finish documenting this!
"""
    defcomposite(cc_name::Symbol, ex::Expr)

Define a Mimi CompositeComponentDef `cc_name` with the expressions in `ex`. Expressions
are all shorthand for longer-winded API calls, and include the following:

    p = Parameter(...)
    v = Variable(varname)
    local_name = Component(name)
    Component(name)  # equivalent to `name = Component(name)`
    connect(...)

Variable names are expressed as the component id (which may be prefixed by a module,
e.g., `Mimi.adder`) followed by a `.` and the variable name in that component. So the
form is either `modname.compname.varname` or `compname.varname`, which must be known
in the current module.

Unlike leaf components, composite components do not have user-defined `init` or
`run_timestep` functions; these are defined internally to iterate over constituent
components and call the associated method on each.
"""
macro defcomposite(cc_name, ex)
    @capture(ex, exprs__)

    calling_module = __module__

    # @info "defcomposite calling module: $calling_module"

    stmts = [_parse(expr) for expr in exprs]

    result = :(
        let cc_id = Mimi.ComponentId($calling_module, $(QuoteNode(cc_name))),
            obj = Mimi.CompositeComponentDef(cc_id)

            global $cc_name = obj
            $(stmts...)
            Mimi.import_params!(obj)
            $cc_name
        end
    )
    return esc(result)
end

"""
    import_params!(obj::AbstractCompositeComponentDef)

Imports all unconnected parameters below the given composite `obj` by adding references
to these parameters in `obj`.

N.B. This is also called at the end of code emitted by @defcomposite.
"""
function import_params!(obj::AbstractCompositeComponentDef)

    unconn = unconnected_params(obj)

    # Check for unresolved parameter name collisions.
    # Users must explicitly define any parameters that come from multiple subcomponents.
    all_names = [ref.datum_name for ref in unconn]
    unique_names = unique(all_names)
    _map = Dict([name => count(isequal(name), all_names) for name in unique_names])
    non_unique = [name for (name, val) in _map if val>1]
    isempty(non_unique) || error("Cannot build composite :$(obj.name). There are unresolved parameter name collisions from subcomponents for the following parameter names: $(join(non_unique, ", ")).")

    for param_ref in unconn
        name = param_ref.datum_name
        haskey(obj, name) && error("Cannot build composite :$(obj.name). Failed to auto-import parameter :$name from component :$(param_ref.comp_name), this name has already been defined in the composite component's namespace.")
        obj[name] = CompositeParameterDef(obj, param_ref)
    end
end

# Helper function for finding any field collisions for parameters that want to be joined
function _find_collisions(fields, pairs::Vector{Pair{T, Symbol}}) where T
    collisions = Symbol[]

    pardefs = [comp.namespace[param_name] for (comp, param_name) in pairs]
    for f in fields
        subcomponent_set = Set([getproperty(pardef, f) for pardef in pardefs])
        length(subcomponent_set) > 1 && push!(collisions, f)
    end

    return collisions
end

# `kwargs` contains the keywords specified by the user when defining the composite parameter in @defcomposite.
# If the user does not provide a value for one or any of the possible fields, this function looks at the fields
# of the subcomponents' parameters to use, but errors if any of them are in conflict.
# Note that :dim_names and :datatype can't be specified at the composite level, but must match from the subcomponents.
function _resolve_composite_parameter_kwargs(obj::AbstractCompositeComponentDef, kwargs::Dict{Symbol, Any}, pairs::Vector{Pair{T, Symbol}}, parname::Symbol)  where T <: AbstractComponentDef

    fields = (:default, :description, :unit, :dim_names, :datatype)
    collisions = _find_collisions(fields, pairs)

    # Create a new dictionary of resolved values to return
    new_kwargs = Dict{Symbol, Any}()

    for f in fields
        try
            new_kwargs[f] = kwargs[f] # Get the user specified value for this field if there is one
        catch e
            # If the composite definition does not specify a value, then need to look to subcomponents and resolve or error
            if f in collisions
                error("Cannot build composite parameter :$parname, subcomponents have conflicting values for the \"$f\" field.")
            else
                compdef, curr_parname = pairs[1]
                pardef = compdef[curr_parname]
                new_kwargs[f] = getproperty(pardef, f)
            end
        end
    end

    return new_kwargs
end

# Helper function for detecting whether a specified datum has already been imported or connected
function _is_connected(obj::AbstractCompositeComponentDef, comp_name::Symbol, datum_name::Symbol)
    for (k, item) in obj.namespace
        if isa(item, AbstractCompositeParameterDef)
            for ref in item.refs
                if ref.comp_name == comp_name && ref.datum_name == datum_name
                    return true
                end
            end
        elseif isa(item, AbstractCompositeVariableDef)
            ref = item.ref
            if ref.comp_name == comp_name && ref.datum_name == datum_name
                return true
            end
        end
    end
    return false

    # cannot use the following, because all parameters haven't bubbled up yet
    # return UnnamedReference(comp_name, datum_name) in unconnected_params(obj)
end

# This function creates a CompositeParameterDef in the CompositeComponentDef obj
function import_param!(obj::AbstractCompositeComponentDef, localname::Symbol,
                        pairs::Pair...; kwargs...)

    print_pairs = [(comp.comp_id, name) for (comp, name) in pairs]
    # @info "import_param!($(obj.comp_id), :$localname, $print_pairs)"

    for (comp, pname) in pairs

        if comp == :*       # wild card
            error("Got wildcard component specification (*) for param $pname (Not yet implemented)")
        else
            compname = nameof(comp)
            has_comp(obj, compname) ||
                error("_import_param!: $(obj.comp_id) has no element named $compname")

            _is_connected(obj, compname, pname) &&
                error("Duplicate import of $(comp.name).$pname")
        end
    end
    new_kwargs = _resolve_composite_parameter_kwargs(obj, Dict{Symbol, Any}(kwargs), collect(pairs), localname)

    obj[localname] = CompositeParameterDef(localname, pathof(obj), collect(pairs), new_kwargs)
end

"""
Import a variable from the given subcomponent
"""
function _import_var!(obj::AbstractCompositeComponentDef, localname::Symbol,
                      path::ComponentPath, vname::Symbol)
    if haskey(obj, localname)
        error("Can't import variable; :$localname already exists in component $(obj.comp_id)")
    end

    comp = @or(find_comp(obj, path), error("$path not found from component $(obj.comp_id)"))
    obj[localname] = CompositeVariableDef(localname, pathof(obj), comp, vname)
end

nothing


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LibTorchCAPI"
version = v"1.7.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/data-efficient-ml/ThArrays.jl.git",
              "bb6c6a63ffd0d11e0df42d5faab2cdf9426e1189"),
]

# Bash recipe for building across all platforms
script_file = joinpath(dirname(@__FILE__), "build_dylib.sh")
if !isfile(script_file)
    script_file = joinpath(dirname(@__FILE__), "L/libtorch_capi/build_dylib.sh")
end
script = read(script_file, String)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi = "cxx11")
    # MacOS(:x86_64), # can't build it on MacOS SDK
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libtorch_capi", :libtorch_capi, dont_dlopen = true),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version = v"8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using NativeSVG

export custom_anchor

```
    Mark the location of a custom anchor for Shaper Origin.
    ```
function custom_anchor(io, x, y; xdirection=1, ydirection=-1, size=0.25u"inch")
    path(io;
         style="fill:red; stroke=none;",
	 d=pathd(
	     ["M", x, y],
	     ["h", xdirection * size],
	     ["L", x, y + ydirection * 2 * size],
	     ["z"]
	 ))
end

max_stars_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_issues_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_forks_repo_name: MarkNahabedian/ShaperOriginDesignLib


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## load packages
using Trixi, OrdinaryDiffEq
using Trixi: TimerOutputs
using DelimitedFiles

## set up benchmark code
function compute_relative_performance(to, semi)
    pid(to, semi) = 1e-9 * TimerOutputs.time(to) / (TimerOutputs.ncalls(to) * Trixi.ndofs(semi))

    time_total = pid(to["rhs!"], semi)
    time_volume = pid(to["rhs!"]["volume integral"], semi)
    time_entropy_projection = pid(to["rhs!"]["entropy_projection!"], semi)
    time_surface = pid(to["rhs!"]["surface integral"], semi)
    time_interface = pid(to["rhs!"]["interface flux"], semi)
    time_other = time_total - (time_volume + time_entropy_projection +
                               time_interface + time_surface)

    open(joinpath(@__DIR__, "Gauss_$(ndims(semi))D_relative_timings.dat"), "a") do io
        writedlm(io, hcat(Trixi.polydeg(semi.solver), # polydeg
                          time_volume,
                          time_entropy_projection,
                          time_surface,
                          time_interface,
                          time_other)
                )
    end
end

function generate_timings(elixir_path, polydeg_range)
    # run for a very short time to avoid timing precompilation
    redirect_stdout(devnull) do
        trixi_include(elixir_path, polydeg=first(polydeg_range), tspan=(0.0, 1e-5))
    end

    for polydeg in polydeg_range
        # run for a very short time to avoid timing precompilation
        redirect_stdout(devnull) do
            trixi_include(elixir_path, polydeg=polydeg, tspan=(0.0, 1e-5))
        end

        t = @elapsed begin
            redirect_stdout(devnull) do
                trixi_include(elixir_path, polydeg=polydeg)
            end
        end
        println("polydeg=$polydeg completed in $t seconds")

        if polydeg==first(polydeg_range)
            open(joinpath(@__DIR__, "Gauss_$(ndims(semi))D_relative_timings.dat"), "w") do io
                header = "Polydeg\ttime_volume\ttime_entropy_projection\ttime_surface\ttime_interface\tother"
                println(io, header)
            end
        end
        to = Trixi.timer()
        compute_relative_performance(to, semi)
    end
end

generate_timings(joinpath(@__DIR__, "elixir_khi_2D.jl"), 3:15)
generate_timings(joinpath(@__DIR__, "elixir_tgv_3D.jl"), 3:15)


text:
# Day 9: Smoke Basin
load(::Day{9}, path) = parse.(Int, hcat(split.(readlines(path), "")...)) # transposed!

# Part 1 – find discrete minima in 4-neighborhood
function is_minimum_in_row(data)
    A = diff(data; dims=2)
    return (A[:, 1:(end - 1)] .< 0) .& (A[:, 2:end] .> 0)
end
function is_minimum_in_col(data)
    A = diff(data; dims=1)
    return (A[1:(end - 1), :] .< 0) .& (A[2:end, :] .> 0)
end
function is_minimum(data)
    h, w = size(data)
    padcol = fill(9, h, 1)
    padrow = fill(9, 1, w)

    r = is_minimum_in_row(hcat(padcol, data, padcol))
    c = is_minimum_in_col(vcat(padrow, data, padrow))
    return r .& c
end

solve1(::Day{9}, data) = sum(data[is_minimum(data)] .+ 1)

# Part 2 – naive flood fill
function neighborhood(I::CartesianIndex) # Von Neumann neighborhood
    r, c = Tuple(I)
    return [
        CartesianIndex(r - 1, c),
        CartesianIndex(r + 1, c),
        CartesianIndex(r, c - 1),
        CartesianIndex(r, c + 1),
    ]
end

function solve2(::Day{9}, data)
    h, w = size(data)
    function is_in_basin(I::CartesianIndex)
        r, c = Tuple(I)
        return 1 <= r <= h && 1 <= c <= w && data[I] < 9
    end

    # We make use of the fact that basins never overlap
    function basin_size(I0::CartesianIndex)
        in_basin = Set([I0])
        to_check = Set([I0])
        while !isempty(to_check)
            for I in to_check
                new_in_basin = setdiff(filter(is_in_basin, neighborhood(I)), in_basin)
                union!(in_basin, new_in_basin)
                union!(to_check, new_in_basin)
                pop!(to_check, I)
            end
        end
        return length(in_basin)
    end

    sizes = basin_size.(findall(is_minimum(data)))
    return prod(sort(sizes; rev=true)[1:3])
end

testresult1(::Day{9}) = 15
testresult2(::Day{9}) = 1134


text:
module RocketSwitchMapOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: switch_map()" begin

    println("Testing: operator switch_map()")

    run_proxyshowcheck("SwitchMap", switch_map(Any), args = (check_subscription = true, ))

    run_testset([
        (
            source      = from(1:5) |> switch_map(Int, d -> of(d ^ 2)),
            values      = @ts([ 1, 4, 9, 16, 25, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> switch_map(Float64, d -> of(convert(Float64, d))),
            values      = @ts([ 1.0, 2.0, 3.0, 4.0, 5.0, c ]),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map(Float64, _ -> faulted(Float64, "err")),
            values      = @ts(e("err")),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map(Float64, _ -> completed(Float64)),
            values      = @ts(c),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> switch_map(Float64, _ -> never(Float64)),
            values      = @ts(),
            source_type = Float64
        ),
        (
            source      = completed() |> switch_map(Int, d -> of(1)),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(Int, "e") |> switch_map(String, d -> string(d)),
            values      = @ts(e("e")),
            source_type = String
        ),
        (
            source      = never() |> switch_map(Int, d -> of(1)),
            values      = @ts(),
            source_type = Int
        ),
        (
            source      = from(1:5) |> switch_map(Int, d -> of(1.0)), # Invalid output: Float64 instead of Int
            values      = @ts(),
            source_type = Int,
            throws      = Exception
        ),
        (
            source      = from(1:5) |> safe() |> switch_map(Int, d -> of(1.0)), # Invalid output: Float64 instead of Int
            values      = @ts(e),
            source_type = Int
        ),
        (
            source      = from(1:5) |> async(0) |> switch_map(Int, d -> of(d ^ 2)),
            values      = @ts([ 1 ] ~ [ 4 ] ~ [ 9 ] ~ [ 16 ] ~ [ 25 ] ~ c),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> switch_map(Int),
            values      = @ts([ 1, 2, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> async(0) |> switch_map(Int),
            values      = @ts([ 1 ] ~ [ 2 ] ~ c ),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> switch_map(Int),
            values      = @ts([ 1, e("err") ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> async(0) |> switch_map(Int),
            values      = @ts([ 1 ] ~ e("err")),
            source_type = Int
        ),
        (
            source      = from([ 0, 0 ]) |> switch_map(Int, d -> from([ 1, 2 ]) |> async(0)),
            values      = @ts([ 1 ] ~ [ 2 ] ~ c),
            source_type = Int
        )
    ])

    customsource1 = make(Int) do actor
        subject1 = Subject(Int)
        subject2 = Subject(Int)
        ssubject = Subject(Any)
        source   = ssubject |> switch_map(Int)

        subscribe!(source, lambda(
            on_next     = (d) -> next!(actor, d),
            on_error    = (e) -> error!(actor, e),
            on_complete = () -> complete!(actor)
        ))

        @async begin
            next!(ssubject, subject1)
            next!(subject1, 1)
            next!(subject2, 2)
            @async begin
                next!(ssubject, subject2)
                next!(subject1, 3)
                next!(subject2, 4)
                @async begin
                    next!(ssubject, subject2)
                    next!(subject1, 5)
                    next!(subject2, 6)
                end
            end
        end
    end

    run_testset([ ( source = customsource1, values = @ts([ 1 ] ~ [ 4 ] ~ [ 6 ]) ) ])

    customsource2 = make(Int) do actor
        ssubject = Subject(Any)
        source   = ssubject |> switch_map(Int)

        subscribe!(source, lambda(
            on_next     = (d) -> next!(actor, d),
            on_error    = (e) -> error!(actor, e),
            on_complete = () -> complete!(actor)
        ))

        @async begin
            next!(ssubject, from([ 1, 2, 3 ]))
            complete!(ssubject)
            next!(ssubject, from([ 1, 2, 3 ])) # should be skipped
        end
    end

    run_testset([ ( source = customsource2, values = @ts([ 1, 2, 3, c ]) ) ])

    customsource3 = make(Int) do actor
        ssubject = Subject(Any)
        source   = ssubject |> switch_map(Int)

        subscribe!(source, lambda(
            on_next     = (d) -> next!(actor, d),
            on_error    = (e) -> error!(actor, e),
            on_complete = ()  -> complete!(actor)
        ))

        @async begin
            complete!(ssubject)
            next!(ssubject, from([ 1, 2, 3 ]))
            next!(ssubject, from([ 1, 2, 3 ])) # should be skipped
        end
    end

    run_testset([ ( source = customsource3, values = @ts(c) ) ])

    customsource4 = make(Int) do actor
        subject1 = Subject(Int)
        ssubject = Subject(Any)
        source   = ssubject |> switch_map(Int)

        subscribe!(source, lambda(
            on_next     = (d) -> next!(actor, d),
            on_error    = (e) -> error!(actor, e),
            on_complete = ()  -> complete!(actor)
        ))

        @async begin
            next!(ssubject, of(0))
            next!(ssubject, subject1)
            complete!(ssubject)
            @async begin
                next!(subject1, 1)
                complete!(subject1)
                next!(ssubject, from([ 1, 2, 3 ])) # should be skipped
            end
        end
    end

    run_testset([ ( source = customsource4, values = @ts([ 0 ] ~ [ 1, c ]) ) ])

end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Taskflow_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Taskflow")
JLLWrappers.@generate_main_file("Taskflow", UUID("fadd631c-8d5c-5dbb-9ac9-de413eafe702"))
end  # module Taskflow_jll


text:
include("utils.jl")

include("test_kmeans.jl")
include("test_distances.jl")

max_stars_repo_name: JuliaManifolds/ManifoldML

max_issues_repo_name: JuliaManifolds/ManifoldML

max_forks_repo_name: JuliaManifolds/ManifoldML


text:
# push!(LOAD_PATH,"../src/")
using Documenter
using CUDA
using LIKWID
using Literate

const src = "https://github.com/JuliaPerf/LIKWID.jl"
const ci = get(ENV, "CI", "") == "true"

@info "Building Literate.jl documentation"
cd(@__DIR__) do
    Literate.markdown("src/examples/dynamic_pinning.jl", "src/examples/";
                        repo_root_url="$src/blob/main/docs") #, codefence = "```@repl 1" => "```")
    Literate.markdown("src/examples/perfmon.jl", "src/examples/";
                        repo_root_url="$src/blob/main/docs") #, codefence = "```@repl 1" => "```")
end

@info "Generating Documenter.jl site"
DocMeta.setdocmeta!(LIKWID, :DocTestSetup, :(using LIKWID, CUDA); recursive=true)
makedocs(
    sitename = "LIKWID.jl",
    authors = "Carsten Bauer",
    modules = [LIKWID],
    doctest = ci,
    pages = [
        "LIKWID" => "index.md",
        "Marker API" => [
            "CPU" => "marker.md",
            "GPU" => "marker_gpu.md",
        ],
        "Library" => [
            "CPU topology" => "topo.md",
            "Performance monitoring" => "perfmon.md",
            "GPU topology" => "topo_gpu.md",
            "NVIDIA monitoring" => "nvmon.md",
            "CPU clock timer" => "timer.md",
            "CPU temperature" => "temperature.md",
            "Power / Energy" => "power.md",
            "Affinity" => "affinity.md",
            "HPM / Access" => "access.md",
            "Miscellaneous" => "misc.md",
        ],
        "CLI Tools" => [
            "likwid-pin" => "likwid-pin.md",
        ],
        "Examples" => [
            "SAXPY CPU+GPU" => "examples/saxpy.md",
            "Pinning Julia threads" => "examples/dynamic_pinning.md",
            "Monitoring performance" => "examples/perfmon.md",
        ],
    ],
    # assets = ["assets/custom.css", "assets/custom.js"]
)

if ci
    @info "Deploying documentation to GitHub"
    deploydocs(
        repo = "github.com/JuliaPerf/LIKWID.jl.git",
        devbranch = "main",
        push_preview = true,
        # target = "site",
    )
end


text:
module TimeMarching

export  stegerwarmingrk1!,upwindrk1!

using ..Tools
using ..FiniteDifference

"""
    this function's inputs are t and uu

        t is the current time

        uu has 3 rows
        the 1st row is u (velocity)
        the 2nd row is m (mass flow rate)
        the 3rd row is ̂E (total energy per volume)

    this function's outputs are new t and uunew

        t is the current time

        uunew has the same structure as uu

    this function uses flux splitting scheme (Steger and Warming, 1981)
    and 1st order Euler time marching

    The current function is only good for u<=c (speed of sound)
"""

function stegerwarmingrk1!(t,uu,righthand,shocktubesystem,Δt)

    gamma=shocktubesystem.gamma
    Δx=shocktubesystem.Δx

    #u<=c
    urhoc = uutourhoc(uu,gamma)

    u = urhoc[1,:]
    rho = urhoc[2,:]
    c = urhoc[3,:]

    #initialization of uunew
    uunew =  deepcopy( zeros(size(uu)))

    #initialization of f
    fplus = deepcopy( zeros(size(uu)))
    fminus = deepcopy( zeros(size(uu)))

    #get fplus and fminus
    fplus[1,:] = (rho./(2gamma)).*(2gamma.*u+c-u)
    fplus[2,:] = (rho./(2gamma)).*(2(gamma-1).*u.^2+(u+c).^2)
    fplus[3,:] = (rho./(2gamma)).*((gamma-1).*u.^2+(((u+c).^3)/2)+((3-gamma).*(u+c).*c.^2)./(2(gamma-1)))

    fminus[1,:] = (rho./(2gamma)).*(u-c)
    fminus[2,:] = (rho./(2gamma)).*(u-c).^2
    fminus[3,:] = (rho./(2gamma)).*((((u-c).^3)/2)+((3-gamma).*(u-c).*c.^2)./(2(gamma-1)))

    for i in 2:1:length(uu[1,:])-1 # shrink a grid
        uunew[:,i] = (-(fplus[:,i]-fplus[:,i-1])/Δx)*Δt + (-(fminus[:,i+1]-fminus[:,i])/Δx)*Δt .+ righthand*Δt + uu[:,i]
    end

    # deal with the side
    uunew[:,length(uunew[1,:])] = uu[:,length(uunew[1,:])]
    uunew[:,1] = uu[:,1]

    t = t + Δt
    return t,uunew
end

"""
    this function's inputs are t and w and impulsesystem

        t is the current time

        w has 1 row
        the 1st row is h (enthalpy)

        impulsesystem is an ImpulseSystem struct

    this function's outputs are new t and out

        t is the current time

        out has the same structure as w

    this function uses upwind scheme and 1st order Euler time marching

    The current function is only good for information tavelling to the right
"""

function upwindrk1!(t::Any,w::Any,impulsesystem::Any,Δt::Any)

    ρ = impulsesystem.ρ
    G = impulsesystem.G
    P = impulsesystem.P
    Ac = impulsesystem.Ac
    qw = impulsesystem.qw
    Δz = impulsesystem.Δz

    out = deepcopy(w) .+ (-G/ρ.*upwind(w,Δz) .+ P*qw/Ac/ρ).*Δt

    t = deepcopy(t) .+ Δt
    return t,out
end

end


text:
################################################################################
#
#  Lattice enumeration for finding all vectors
#
################################################################################

################################################################################
#
#  Pseudo cholesky decomposition
#
################################################################################

function _pseudo_cholesky(G::fmpz_mat)
  return _pseudo_cholesky(G, fmpq_mat)
end

function _pseudo_cholesky(G::fmpz_mat, ::Type{fmpq_mat})
  n = ncols(G)
  C = zero_matrix(FlintQQ, n, n)
  return __pseudo_cholesky!(C, G)
end

function _pseudo_cholesky(G::fmpz_mat, ::Type{Matrix{fmpq}})
  n = ncols(G)
  C = Matrix{fmpq}(undef, n, n)
  return __pseudo_cholesky!(C, G)
end

function _pseudo_cholesky(G::fmpz_mat, ::Type{Matrix{S}}) where {T, S <: Union{Rational{T}, UnsafeRational{T}}}
  n = ncols(G)
  z = Matrix{S}(undef, n, n)
  C = _pseudo_cholesky(G, Matrix{fmpq})
  for i in 1:n
    for j in 1:n
      z[i, j] = S(T(numerator(C[i, j])), T(denominator(C[i, j])))
    end
  end
  return z
end

function __pseudo_cholesky!(C::fmpq_mat, G)
  n = ncols(G)
  limit = n

  for i=1:limit
    for j=1:limit
      C[i, j] = G[i, j]
    end
  end
  for i = 1:limit-1
    for j = i+1:limit
      C[j, i] = deepcopy(C[i, j])
      C[i, j] = divexact(C[i, j], C[i, i])
    end
    for k = i+1:limit
      for l = i+1:limit
        C[k,l] = C[k,l] - C[k,i] * C[i,l]
      end
    end
  end
  for j = 1:limit-1
    if C[j,j] <= 0
      throw(LowPrecisionCholesky())
    end
    for i=j+1:limit
      C[i, j] = zero!(C[i, j])
    end
  end
  if C[limit, limit] <= 0
    throw(LowPrecisionCholesky())
  end
  return C
end

function __pseudo_cholesky!(C::Matrix{fmpq}, G)
  n = ncols(G)
  limit = n
  t = fmpq()

  for i=1:limit
    for j=1:limit
      C[i, j] = G[i, j]
    end
  end
  for i = 1:limit-1
    for j = i+1:limit
      C[j, i] = _deepcopy_cheap(C[i, j])
      #C[i, j] = divexact(C[i, j], C[i, i])
      C[i, j] = divexact!(C[i, j], C[i, j], C[i, i])
    end
    for k = i+1:limit
      for l = i+1:limit
        t = mul!(t, C[k, i], C[i, l])
        C[k, l] = sub!(C[k, l], C[k, l], t)
        #C[k,l] = C[k,l] - C[k,i] * C[i,l]
      end
    end
  end
  for j = 1:limit-1
    if C[j,j] <= 0
      throw(LowPrecisionCholesky())
    end
    for i=j+1:limit
      C[i, j] = zero!(C[i, j])
    end
  end
  if C[limit, limit] <= 0
    throw(LowPrecisionCholesky())
  end
  return C
end

################################################################################
#
#  Lattice enumeration
#
################################################################################

function __enumerate_gram(G::fmpz_mat, c::Int)
  # This is the "high" level function
  Q = _pseudo_cholesky(G, Matrix{fmpq})
  n = nrows(G)
  d = lcm([denominator(Q[i]) for i in 1:length(G)])
  @vprint :Lattice 1 "Denominator of pseudo-Cholesky of bit size $(nbits(d))\n"
  k = nbits(d) + nbits(n) + 2
  isbig = Int == Int64
  if 2 * k < (isbig ? 32 : 64)
    @vprint :Lattice 1 "Enumerating using Int\n"
    Qint = Matrix{UnsafeRational{Int}}([Int(numerator(q))//Int(denominator(q)) for q in Q])
    res = __enumerate_cholesky(Qint, c)
    @hassert :Lattice 1 length(__enumerate_gram(G, c, Rational{Int})) == length(res)
  elseif 2 * k < 64
    @vprint :Lattice 1 "Enumerating using Int64\n"
    Qint64 = Matrix{UnsafeRational{Int64}}([Int64(numerator(q))//Int64(denominator(q)) for q in Q])
    res = __enumerate_cholesky(Qint64, c)
    @hassert :Lattice 1 length(__enumerate_gram(G, c, Rational{Int64})) == length(res)
  elseif 2 * k < 128
    Qint128 = Matrix{UnsafeRational{Int128}}([Int128(numerator(q))//Int128(denominator(q)) for q in Q])
    @vprint :Lattice 1 "Enumerating using Int128\n"
    res = __enumerate_cholesky(Qint128, c)
    @hassert :Lattice 1 length(__enumerate_gram(G, c, Rational{Int128})) == length(res)
  else
    @vprint :Lattice 1 "Enumerating using fmpq\n"
    res = __enumerate_cholesky(Q, c)
  end
  @hassert :Lattice 1 length(__enumerate_gram(G, c, fmpq)) == length(res)
  return res
end

function __enumerate_gram(G::fmpz_mat, c::Union{Int, fmpz}, ::Type{T}) where {T}
  Q = _pseudo_cholesky(G, Matrix{T})
  v = __enumerate_cholesky(Q, c)
  return v
end

function __enumerate_cholesky(Q::fmpq_mat, c::Int)
  return __enumerate(Matrix(Q), c)
end

function __enumerate_cholesky(Q::Matrix{fmpq}, c::S) where {S <: Union{Int, fmpz}}
  res = Tuple{Vector{S}, S}[]
  n = nrows(Q)
  i = n
  T = Vector{fmpq}(undef, n)
  U = Vector{fmpq}(undef, n)

  for i in 1:n
    T[i] = fmpq()
    U[i] = fmpq()
  end
  x = Vector{S}(undef, n)
  L = Vector{S}(undef, n)

  T[i] = c
  U[i] = 0

  t1 = fmpq()
  t2 = fmpz()
  t3 = fmpz()
  t4 = fmpz()
  t5 = fmpz()

  Qd = Vector{fmpq}(undef, n) # diagonal of Q
  @inbounds for i in 1:n
    Qd[i] = Q[i, i]
  end

  @label compute_bounds

  # For debugging purposes. This gives the tightest bounds
  #xc = -U[i]
  #while (Q[i, i] * (xc + U[i])^2 <= T[i])
  #  xc -= 1
  #end
  #x[i] = Int(FlintZZ(floor(xc)))

  #upperbound_can = -U[i]
  #while (Q[i, i] * (upperbound_can + U[i])^2 <= T[i])
  #  upperbound_can += 1
  #end
  #L[i] = Int(FlintZZ(ceil(upperbound_can)))

  # For debugging purposes: Fast, but allocating
  #
  #_t = isqrt(FlintZZ(floor(divexact(T[i], Q[i, i]))))
  #_new_upp = Int(FlintZZ(ceil(_t + 2 - U[i])))
  #_new_low = Int(FlintZZ(floor(-(_t + 2) - U[i]))) - 1

  @inbounds _new_upp, _new_low = _compute_bounds(T[i], Qd[i], U[i], t1, t2, t3, t4, t5)

  @inbounds x[i] = _new_low
  @inbounds L[i] = _new_upp

  @label main_loop

  @inbounds x[i] = x[i] + 1

  @inbounds if x[i] > L[i]
    i = i + 1
    @goto main_loop
  else
    if i > 1
      #T[i - 1] = T[i] - Q[i, i] * (x[i] + U[i])^2
      @inbounds update_T!(T, i, Qd, x, U, t1)
      @inbounds if isnegative(T[i - 1])
        @goto main_loop
      end
      i = i - 1
      #U[i] = sum(Q[i, j] * x[j] for j in (i + 1):n)
      update_U!(U, Q, i, n, x, t1, t2)
      @goto compute_bounds
    end
  end

  if iszero(x)
    return res
  else
    #len = c - T[1] + Q[1, 1]*(x[1] + U[1])^2
    t1 = compute_len!(t1, c, T, Q, x, U)
    t2 = numerator!(t2, t1)
    _short_enough = false
    if S === Int
      _len = ccall((:fmpz_get_si, libflint), Int, (Ref{fmpz}, ), t2)
      _short_enough = _len <= c
      len = _len
    else
      _short_enough = t2 <= c
      len = deepcopy(t2)
    end

    if _short_enough
      # t1 must be a UInt
      #z = fmpq()
      #ccall((:fmpq_set, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}), z, t1)
      # Todo
      y = Vector{S}(undef, n)
      if S === fmpz
        for i in 1:n
          @inbounds y[i] = deepcopy(x[i])
        end
      else
        for i in 1:n
          @inbounds y[i] = x[i]
        end
      end
      push!(res, (y, len))
    end
    @goto main_loop
  end
end

function __enumerate_cholesky(Q::Matrix{S}, c::Int) where {S <: Union{Rational, UnsafeRational}}
  res = Tuple{Vector{Int}, Int}[]
  n = nrows(Q)
  i = n
  T = Vector{S}(undef, n)
  U = Vector{S}(undef, n)
  L = Vector{S}(undef, n)
  #for i in 1:n
  #  T[i] = 0
  #  U[i] = 0
  #end
  x = Vector{Int}(undef, n)

  T[i] = c
  U[i] = 0

  Qd = Vector{S}(undef, n)
  for i in 1:n
    Qd[i] = Q[i, i]
  end

  #@show d = lcm(Int[denominator(a) for a in Q])
  #@show nbits(d)

  @label compute_bounds

  # For debugging purposes. This gives the tightest bounds
  #xc = -U[i]
  #while (Q[i, i] * (xc + U[i])^2 <= T[i])
  #  xc -= 1
  #end
  #x[i] = Int(FlintZZ(floor(xc)))

  #upperbound_can = -U[i]
  #while (Q[i, i] * (upperbound_can + U[i])^2 <= T[i])
  #  upperbound_can += 1
  #end
  #L[i] = Int(FlintZZ(ceil(upperbound_can)))

  # For debugging purposes: Fast, but allocating
  #
  #_t = isqrt(FlintZZ(floor(divexact(T[i], Q[i, i]))))
  #_new_upp = Int(FlintZZ(ceil(_t + 2 - U[i])))
  #_new_low = Int(FlintZZ(floor(-(_t + 2) - U[i]))) - 1

  _new_upp, _new_low = @inbounds _compute_bounds(T[i], Qd[i], U[i])

  @inbounds x[i] = _new_low
  @inbounds L[i] = _new_upp

  @label main_loop

  @inbounds x[i] = x[i] + 1

  @inbounds if x[i] > L[i]
    i = i + 1
    @goto main_loop
  else
    if i > 1
      #T[i - 1] = T[i] - Q[i, i] * (x[i] + U[i])^2
      @inbounds update_T!(T, i, Qd, x, U)
      @inbounds if isnegative(T[i - 1])
        @goto main_loop
      end
      i = i - 1
      #U[i] = sum(Q[i, j] * x[j] for j in (i + 1):n)
      update_U!(U, Q, i, n, x)

      @goto compute_bounds
    end
  end

  if iszero(x)
    return res
  else
    #len = c - T[1] + Q[1, 1]*(x[1] + U[1])^2
    len = Int(compute_len!(c, T, Q, x, U))
    if len <= c
      y = Vector{Int}(undef, n)
      @inbounds for i in 1:n
        y[i] = x[i]
      end
      push!(res, (y, len))
    end
    @goto main_loop
  end
end

@inline function _compute_bounds(Ti, Qi, Ui, t1 = fmpq(),
                                             t2 = fmpz(),
                                             t3 = fmpz(),
                                             t4 = fmpz(),
                                             t5 = fmpz())
  t1 = divexact!(t1, Ti, Qi)
  t2 = floor!(t2, t1, t3, t4)
  t2 = isqrt!(t2, t2)
  # Now t2 = ceil(sqrt(Ti/Qi))
  t2 = add_two!(t2, t2)
  # Now t2 = ceil(sqrt(Ti/Qi)) + 2

  t1 = sub!(t1, Ui, t2)
  t1 = neg!(t1, t1)
  t5 = ceil!(t5, t1, t3, t4)
  ub = ccall((:fmpz_get_si, libflint), Int, (Ref{fmpz}, ), t5)

  t1 = add!(t1, Ui, t2)
  t1 = neg!(t1, t1)
  t5 = floor!(t5, t1, t3, t4)
  lb = ccall((:fmpz_get_si, libflint), Int, (Ref{fmpz}, ), t5) - 1

  return ub, lb
end

@inline function _compute_bounds(Ti::Union{UnsafeRational, Rational{<:Integer}}, Qi, Ui)
  t1 = Ti//Qi
  t2 = Int(floor(t1))
  t2 = isqrt(t2)
  # Now t2 = ceil(sqrt(Ti/Qi))
  t2 += 2
  # Now t2 = ceil(sqrt(Ti/Qi)) + 2

  t1 = Ui - t2
  t1 = -t1
  t5 = Int(ceil(t1))
  ub = t5

  t1 = Ui + t2
  t1 = -t1
  t5 = Int(floor(t1))
  lb = t5

  return ub, lb - 1
end

@inline function update_T!(T, i, Qd, x, U, t = fmpq())
  @inbounds t = add!(t, U[i], x[i])
  t = mul!(t, t, t)
  t = mul!(t, t, Qd[i])
  sub!(T[i - 1], T[i], t)
end

@inline function update_T!(T::Vector{S}, i, Qd, x, U) where S <: Union{Rational{<:Integer}, UnsafeRational}
  # T[i - 1] = T[i] - QQ[i, i] * (U[i] + x[i])^2
  @inbounds t = U[i] + x[i]
  t = t * t
  @inbounds t = t * Qd[i]
  @inbounds T[i - 1] = T[i] - t
  #@assert Rational(T[i - 1]) == Rational(T[i]) - Rational(QQ[i, i]) * (Rational(U[i]) + Rational(x[i]))^2
end

@inline function update_U!(U, QQ, i, n, x, t = fmpq(), t2 = fmpz())
  # U[i] = sum(QQ[i, j] * x[j] for j in (i + 1):n)
  s = @inbounds U[i]
  @inbounds ccall((:fmpz_set_si, libflint), Cvoid, (Ref{fmpz}, Int), t2, x[i + 1])
  @inbounds s = mul!(s, QQ[i, i + 1], t2)
  for j in (i + 2):n
    @inbounds ccall((:fmpz_set_si, libflint), Cvoid, (Ref{fmpz}, Int), t2, x[j])
    @inbounds mul!(t, QQ[i, j], t2)
    s = add!(s, s, t)
    #addmul!(s, QQ[i, j], t2, t)
  end
  return s
end

@inline function update_U!(U, QQ::Matrix{S}, i, n, x) where S <: Union{Rational{<:Integer}, UnsafeRational}
  # U[i] = sum(QQ[i, j] * x[j] for j in (i + 1):n)
  @inbounds s = QQ[i, i + 1] * x[i + 1]
  @inbounds for j in (i + 2):n
    t = QQ[i, j] * x[j]
    s = s + t
    #addmul!(s, QQ[i, j], t2, t)
  end
  U[i] = s
  #@assert Rational(U[i]) == sum(Rational(QQ[i, j]) * Rational(x[j]) for j in (i + 1):n)
end

@inline function compute_len!(t, c, T, QQ, x, U)
#len = c - T[1] + QQ[1, 1]*(x[1] + U[1])^2
  t = @inbounds add!(t, U[1], x[1])
  t = mul!(t, t, t)
  t = @inbounds mul!(t, t, QQ[1, 1])
  t = add!(t, t, c)
  t = @inbounds sub!(t, t, T[1])
  return t
end

@inline function compute_len!(c, T::Vector{S}, QQ, x, U) where S <: Union{Rational{<:Integer}, UnsafeRational}
#len = c - T[1] + QQ[1, 1]*(x[1] + U[1])^2
  s = @inbounds begin
  t = U[1] + x[1]
  t = t * t
  t = t * QQ[1, 1]
  t = t + c
  t = t - T[1]
  end
  return s
end

function _deepcopy_cheap(x::fmpq)
  z = fmpq()
  ccall((:fmpq_set, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}), z, x)
  return z
end

function isnegative(x::fmpq)
  c = ccall((:fmpq_sgn, libflint), Cint, (Ref{fmpq}, ), x)
  return c < 0
end

function islessorequal(x::fmpq, y::UInt)
  c = ccall((:fmpq_cmp_ui, libflint), Cint, (Ref{fmpq}, UInt), x, y)
  return c <= 0
end

################################################################################
#
#  Enumeration functions on gram matrices
#
################################################################################

# This is the work horse
# Assumes that G is a gram matrix, finds all vectors v with length l
# lb <= v < ub. Also applies the transformation transform
#
# If transform === nothing, no transform
#
# Return value is Vector{Tuple{Vector{Int}, fmpz}}
function _short_vectors_gram_nolll_integral(G, lb, ub, transform)
  n = nrows(G)
  max = _round_down(fmpz, ub)
  # G is integral, so q(x) <= ub is equivalent to q(x) <= floor(ub)

  if ub isa fmpz && fits(Int, ub)
    V = __enumerate_gram(G, Int(ub), fmpq)
  else
    V = __enumerate_gram(G, ub, fmpq)
  end

  if ub isa fmpz
    W = Vector{Tuple{Vector{fmpz}, fmpz}}()
  else
    W = Vector{Tuple{Vector{Int}, fmpz}}()
  end

  for (v, l) in V
    if l < lb
      continue
    end

    m = _transform(v, transform)

    positive = false
    for k in 1:n
      if iszero(m[k])
        continue
      end
      if m[k] > 0
        positive = true
      end
      break
    end

    if positive
      push!(W, (m, FlintZZ(l)))
    else
      push!(W, (-m, FlintZZ(l)))
    end
  end
  return W
end

# No assumption on _G, algorithm applies LLL
function _short_vectors_gram(_G, lb, ub; hard::Bool = false)
  d = denominator(_G)
  G = change_base_ring(FlintZZ, d * _G)
  if hard
    Glll, T = lll_gram_with_transform(G, lll_ctx(0.99999999999999, 0.500000000001, :gram))
  else
    Glll, T = lll_gram_with_transform(G, lll_ctx(0.9999, 0.5001, :gram))
  end
  if isone(T)
    V = _short_vectors_gram_nolll_integral(Glll, d * lb, d * ub, nothing)
  else
    V = _short_vectors_gram_nolll_integral(Glll, d * lb, d * ub, T)
  end
  W = Vector{Tuple{Vector{Int}, elem_type(base_ring(_G))}}(undef, length(V))
  for i in 1:length(V)
    W[i] = (V[i][1], V[i][2]//d)
  end
  return W
end

# No assumption on _G, algorithm applies LLL
function _short_vectors_gram(_G, ub; hard::Bool = false)
  return _short_vectors_gram(_G, fmpz(0), ub, hard = hard)
end

# No assumption on _G, algorithm applies LLL
function _shortest_vectors_gram(_G)
  d = denominator(_G)
  G = change_base_ring(FlintZZ, d * _G)
  Glll, T = lll_gram_with_transform(G)
  max = maximum([Glll[i, i] for i in 1:nrows(G)])
  @assert max > 0
  if isone(T)
    V = _short_vectors_gram_nolll_integral(Glll, 0, max, nothing)
  else
    V = _short_vectors_gram_nolll_integral(Glll, 0, max, T)
  end
  min = minimum(v[2] for v in V)
  return min//d, Vector{fmpz}[ v[1] for v in V if v[2] == min]
end

function _short_vectors_gram_integral(_G, ub; hard = false)
  if hard
    Glll, T = lll_gram_with_transform(_G, lll_ctx(0.99999999999999, 0.500000000001, :gram))
  else
    Glll, T = lll_gram_with_transform(_G)
  end
  if isone(T)
    V = _short_vectors_gram_nolll_integral(Glll, 0, ub, nothing)
  else
    V = _short_vectors_gram_nolll_integral(Glll, 0, ub, T)
  end

  return V
end

function _shortest_vectors_gram_integral(_G)
  Glll, T = lll_gram_with_transform(_G)
  max = maximum([Glll[i, i] for i in 1:nrows(Glll)])
  @assert max > 0
  if isone(T)
    V = _short_vectors_gram_nolll_integral(Glll, 0, max, nothing)
  else
    V = _short_vectors_gram_nolll_integral(Glll, 0, max, T)
  end
  min = minimum(v[2] for v in V)
  return min, [ v for v in V if v[2] == min]
end

_transform(m::fmpz_mat, T) = m * T

_transform(m, ::Nothing) = m

_transform(m::Vector{Int}, ::Nothing) = m

function _transform(m::Vector{Int}, T)
  n = matrix(FlintZZ, 1, nrows(T), m) * T
  return fmpz[n[1, i] for i in 1:nrows(T)]
end

function _transform(m::Vector{fmpz}, T)
  n = matrix(FlintZZ, 1, nrows(T), m) * T
  return fmpz[n[1, i] for i in 1:nrows(T)]
end

_round_up(::Type{fmpz}, x::fmpz) = x

_round_up(::Type{fmpz}, x::fmpq) = FlintZZ(ceil(x))

_round_down(::Type{fmpz}, x::fmpz) = x

_round_down(::Type{fmpz}, x::fmpq) = FlintZZ(floor(x))

################################################################################
#
#  Additional inplace operations
#
################################################################################

@inline function numerator!(z::fmpz, y::fmpq)
  ccall((:fmpq_numerator, libflint), Cvoid, (Ref{fmpz}, Ref{fmpq}), z, y)
  return z
end

@inline function denominator!(z::fmpz, y::fmpq)
  ccall((:fmpq_denominator, libflint), Cvoid, (Ref{fmpz}, Ref{fmpq}), z, y)
  return z
end

@inline function floor!(z::fmpz, y::fmpq, t0::fmpz, t1::fmpz)
  numerator!(t0, y)
  denominator!(t1, y)
  ccall((:fmpz_fdiv_q, libflint), Cvoid, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, t0, t1)
  return z
end

@inline function ceil!(z::fmpz, y::fmpq, t0::fmpz, t1::fmpz)
  numerator!(t0, y)
  denominator!(t1, y)
  ccall((:fmpz_cdiv_q, libflint), Cvoid, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, t0, t1)
  return z
end

@inline function divexact!(z::fmpq, a::fmpq, b::fmpq)
  ccall((:fmpq_div, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
  return z
end

@inline function add_two!(z::fmpz, x::fmpz)
  ccall((:fmpz_add_ui, libflint), Cvoid, (Ref{fmpz}, Ref{fmpz}, Int), z, x, 2)
  return z
end

@inline function sub!(z::fmpq, a::fmpq, b::fmpq)
  ccall((:fmpq_sub, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
  return z
end

@inline function sub!(z::fmpq, a::fmpq, b::fmpz)
  ccall((:fmpq_sub_fmpz, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
  return z
end

@inline function add!(z::fmpq, a::fmpq, b::fmpz)
  ccall((:fmpq_add_fmpz, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
  return z
end

@inline function mul!(z::fmpq, a::fmpq, b::fmpz)
  ccall((:fmpq_mul_fmpz, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
  return z
end

@inline function neg!(z::fmpq, a::fmpq)
  ccall((:fmpq_neg, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}), z, a)
  return z
end

@inline function isqrt!(z::fmpz, a::fmpz)
  ccall((:fmpz_sqrt, libflint), Cvoid, (Ref{fmpz}, Ref{fmpz}), z, a)
  return z
end

@inline function add!(z::fmpq, a::fmpq, b::Int)
  ccall((:fmpq_add_si, libflint), Cvoid, (Ref{fmpq}, Ref{fmpq}, Int), z, a, b)
  return z
end

divexact!(z::Rational{Int}, x::Rational{Int}, y::Rational{Int}) = divexact(x, y)

floor!(z::Int, x::Rational{Int}, y::Int, w::Int) = Int(floor(x))

isqrt!(z::Int, x::Int) = isqrt(x)

add_two!(z::Int, x::Int) = x + 2

sub!(z::Rational{Int}, x::Rational{Int}, y::Int) = x - y

neg!(z::Rational{Int}, x::Rational{Int}) = -x

ceil!(z::Int, x::Rational{Int}, y::Int, w::Int) = Int(ceil(x))

add!(z::Rational{Int}, x::Rational{Int}, y::Int) = x + y

mul!(z::Rational{Int}, x::Rational{Int}, y::Int) = x * y

numerator!(z::Int, x::Rational{Int}) = numerator(x)

isnegative(x::Rational) = x.num < 0


text:
# This doesn't really belong here, but it's convenient.

adapt_(T, x) = x

adapt(T, x) = adapt_(T, x)

adapt(T, x::RowVector) = RowVector(adapt(T, x.vec))


text:
function dl_msldzl()
    url = "http://www.dushu369.com/zhongguomingzhu/HTML"
    pf = open("msl.md", "w")
    for page in 61706:-1:61685
        download("$url/$page.html", "tmp.txt")
        run(`iconv -cf gbk -t utf8 tmp.txt -o page.txt`)
        chapter = parsehtml(read("page.txt", String))
        for t in eachmatch(sel"title", chapter.root)
            println(pf, "## ", text(t), "\n\n")
        end

        for c in eachmatch(Selector(".content"), chapter.root)
            println(pf, text(c), "\n")
        end
    end
    close(pf)
end


text:
@testset "514.freedom-trail.jl" begin
    @test find_rotate_steps("godding", "gd") == 4
    @test find_rotate_steps("godding", "godding") == 13
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Nemo
matrix(T, A) = MatrixSpace(T, size(A)...)(A)
pl(x...) = println(x...)
macro s(expr)
    :(println($(string(expr)), " ="); show(stdout, "text/plain", $expr); println("\n"))
end

F, a = FiniteField(7, 2, "a")
P, x = F["x"]
X = matrix(F, (0:6)' .* a .+ (0:6))
A = F[a+1 a+2; a+3 3a+4]
v = F[5a+2; 3a+6]

@show F P; pl()
@show [F(3)^k for k in 1:6]; pl()
@show [a^k for k in 1:48]; pl()
@s X
@s rref(X)[2]
@show det(X) tr(X) rank(X) charpoly(P, X) minpoly(P, X); pl()
@s A
@s A'
@s A'A
@show det(A) tr(A) rank(A) charpoly(P, A) lu(A) inv(A); pl()
@show v A*v; pl()

# %%
methodswith(typeof(A))

# %%

max_stars_repo_path: 0011/FiniteField with linear algebra.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0011/FiniteField with linear algebra.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0011/FiniteField with linear algebra.jl

max_forks_repo_name: genkuroki/public


text:
module TestExamples3DPart2

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "3D-Part2" begin

# Run basic tests
@testset "Examples 3D" begin
  # MHD
  include("test_examples_3d_mhd.jl")

  # Lattice-Boltzmann
  include("test_examples_3d_lbm.jl")

  # Curved mesh
  include("test_examples_3d_structured.jl")

  # P4estMesh
  include("test_examples_3d_p4est.jl")
end

@testset "Additional tests in 3D" begin
  @testset "ideal GLM MHD" begin
    eqn = IdealGlmMhdEquations3D(1.4)
    u = [1.0, 2.0, 3.0, 4.0, 20.0, 0.1, 0.2, 0.3, 1.5]

    @test isapprox(density(u, eqn), 1.0)
    @test isapprox(pressure(u, eqn), 1.7219999999999995)
    @test isapprox(density_pressure(u, eqn), 1.7219999999999995)

    @test isapprox(Trixi.entropy_thermodynamic(u, eqn), 0.5434864060055388)
    @test isapprox(Trixi.entropy_math(u, eqn), -1.3587160150138473)
    @test isapprox(Trixi.entropy(u, eqn), -1.3587160150138473)

    @test isapprox(energy_total(u, eqn), 20.0)
    @test isapprox(energy_kinetic(u, eqn), 14.5)
    @test isapprox(energy_magnetic(u, eqn), 0.07)
    @test isapprox(energy_internal(u, eqn), 4.305)

    @test isapprox(cross_helicity(u, eqn), 2.0)
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # 3D-Part2

end #module


text:
include("load_data.jl")

using Knet

#=
Network consists of 3 layers:

input layer consisting of linear units (taking the input matrix)
hidden layer including user defined number of logistic units
Output layer with 10 units all of which are a softmax group
=#

@knet function forward_net(train_data; hidden=100)
  input_to_hidden_weights = par(init=Gaussian(0, 0.1), dims = (hidden,256))
  input_to_hidden = input_to_hidden_weights * train_data

  hidden_activation = sigm(input_to_hidden)
  hidden_to_output_weights = par(init=Gaussian(0, 0.1), dims= (10,hidden))

  input_to_output = hidden_to_output_weights * hidden_activation
  return soft(input_to_output)
end

function train(f, data, loss)
    for (x,y) in data
        forw(f, x)
        back(f, y, loss)
        update!(f)
    end
end

function test(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in data
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    sumloss / numloss
end

function main(numofepochs = 100)
   # Load the data
   (train_input,
    train_target,
    valid_input,
    valid_target,
    test_input,
    test_target) = load_data()

    # Split data into minibatches
    batchsize = 100
    train_data = minibatch(train_input, train_target, batchsize)
    test_data = minibatch(test_input, test_target, batchsize)

   net = compile(:forward_net)
   setp(net; lr=0.1)

   println("Test error before training: ", test(net, test_data, zeroone))

   for i=1:numofepochs
     train(net, train_data, softloss)
     println("Test error after epoch ", i, ": ",test(net, test_data, zeroone))
   end

end
main()

max_stars_repo_name: kirnap/Machine-Learning-in-Julia

max_issues_repo_name: kirnap/Machine-Learning-in-Julia

max_forks_repo_name: kirnap/Machine-Learning-in-Julia


text:
abstract type Prolog end

mutable struct Swipl <: Prolog
    initiated::Bool
end

Swipl() = Swipl(false)