text:
@testset "Problems" begin
    F = equations(katsura(5))
    G = equations(cyclic(6))

    P1 = StartTargetInput(G, F)
    (PP1, start1) = problem_startsolutions(P1, [rand(ComplexF64, 6), rand(ComplexF64, 6)])
    @test PP1 isa Problem{AffineTracking}
    @test length(start1) == 2

    P2 = TargetSystemInput(G)
    (PP2, start2) = problem_startsolutions(P2)
    @test PP2 isa Problem{AffineTracking}
    @test length(start2) == 720

    @test homotopy(PP2) isa AbstractHomotopy
    @test homvars(PP2) == nothing

    @polyvar x y z
    F = [x^2+y^2+2z^2, x+y+3z]
    P, _ = problem_startsolutions(TargetSystemInput(F), nothing, homvar=y)
    @test homvars(P) == (2,)

    P, _ = problem_startsolutions(TargetSystemInput(F), nothing)
    @test homvars(P) === nothing

    @test_throws ArgumentError problem_startsolutions(StartTargetInput(
        [x^2+y^2+z^2, x^4+y^4+z^3], [x^3+z^3,y^3-z^3]), [rand(ComplexF64, 3)])

    P, _ = problem_startsolutions(StartTargetInput(
        [x^2+y^2+z^2, x^4+y^4+z^4], [x^3+z^3,y^3-z^3]), [rand(ComplexF64, 3)])
    @test P isa Problem{ProjectiveTracking}
    @test homvars(P) === nothing

    P, _ = problem_startsolutions(StartTargetInput(
        [x^2+y^2+z^2, x^4+y^4+z^2], [x^3+z^3,y^3-z^2]), [rand(ComplexF64, 3)])
    @test P isa Problem{AffineTracking}

    # not supported right now
    @test_throws MethodError problem_startsolutions(
        StartTargetInput(SPSystem([x^2+y^2+z^2, x^4+y^4+z^2]), SPSystem([x^3+z^3,y^3-z^2])),
        [rand(ComplexF64, 3)])

    P, starts = problem_startsolutions(HC.input_startsolutions(
        [x^2+y^2+1, x^4+y^4+1], [x^3+1,y^3-1], rand(ComplexF64, 3))...)
    @test P isa Problem{AffineTracking}
    @test length(starts) == 1

    P, _ = problem_startsolutions(StartTargetInput(
        [x^2+y^2+z^2, x^4+y^4+z^4], [x^3+z^3,y^3-z^3]), [rand(ComplexF64, 3)], homvar=z)
    @test homvars(P) === (3,)

    F = FPSystem(homogenize(equations(cyclic(6))))
    P, startvals = problem_startsolutions(TargetSystemInput(F))
    @test P isa Problem{ProjectiveTracking}
    @test homvars(P) == nothing

    F = homogenize(equations(cyclic(6)))
    P, _ = problem_startsolutions(TargetSystemInput(F))
    @test P isa Problem{ProjectiveTracking}

    F = FPSystem(homogenize(equations(cyclic(6))))
    P, startvals = problem_startsolutions(TargetSystemInput(F), homvar=5)
    @test P isa Problem{ProjectiveTracking}
    @test homvars(P) == (5,)
    @test length(startvals) == 720
end


text:
using ClusterStructureAnalysis
using Test

@testset "ClusterStructureAnalysis.jl" begin
    # Write your tests here.
end

ext: jl


text:
using Documenter
using TerminalPager

makedocs(
    modules = [TerminalPager],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://ronisbr.github.io/TerminalPager.jl/stable/",
        edit_link = "main"
    ),
    sitename = "Terminal Pager",
    authors = "Ronan Arraes Jardim Chagas",
    pages = [
        "Introduction" => "index.md",
        "Usage"        => "man/usage.md",
    ]
)

deploydocs(
    repo = "github.com/ronisbr/TerminalPager.jl.git",
    devbranch = "main",
    target = "build",
)


text:
"""
Operation - Google Code Jam 2017 Finals, Problem B
(https://code.google.com/codejam/contest/6314486/dashboard#s=p1&a=0)

Run the script:

    julia --project src/Operation.jl <in.txt >out.txt
"""

module Operation

using Combinatorics: permutations

p(s::AbstractString) = parse(Int, s)

OP_MAP = Dict(
    "+" => +,
    "-" => -,
    "*" => *,
    "/" => //,
)

function solve!(start::BigInt, card_map::Dict{String, Vector{BigInt}})::Rational{BigInt}
    cards = Tuple{Function, BigInt}[]

    for op in ["+", "-"]
        push!(cards, (OP_MAP[op], sum(card_map[op])))
    end

    z_ind = findall(iszero, card_map["*"])
    if !isempty(z_ind)
        push!(cards, (OP_MAP["*"], 0))
        deleteat!(card_map["*"], z_ind)
    end

    for op in ["*", "/"]
        pd = prod(card_map[op])
        neg_list = filter(x -> x < 0, card_map[op])

        if length(neg_list) == 1
            m = neg_list[1]
            pd //= m
            push!(cards, (OP_MAP[op], m))
        elseif length(neg_list) >= 2
            mi, ma = minimum(neg_list), maximum(neg_list)
            pd //= mi * ma
            push!(cards, (OP_MAP[op], mi))
            push!(cards, (OP_MAP[op], ma))
        end

        push!(cards, (OP_MAP[op], pd))
    end

    reduce(max, (reduce((acc, (op, x)) -> op(acc, x), cs; init = start // big(1))
        for cs in permutations(cards)))
end

function main()
    n_tests = p(readline())
    for t in 1:n_tests
        (start, n_cards) = p.(split(readline()))

        card_map = Dict(op => BigInt[] for (op, _) in OP_MAP)
        for _ in 1:n_cards
            (op, x) = begin
                (op, x) = split(readline())
                (String(op), p(x))
            end
            if op == "+" && x < 0
                (op, x) = ("-", -x)
            elseif op == "-" && x < 0
                (op, x) = ("+", -x)
            end
            push!(card_map[op], x)
        end

        res = solve!(big(start), card_map)
        println("Case #$(t): $(res.num) $(res.den)")
        flush(stdout)
    end
end

main()

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
struct HoldOut{T <: Persa.AbstractDataset} <: ResamplingStrategy{T}
    dataset::T
    index::Array{Int}
    k::Float64
end

function HoldOut(dataset::Persa.AbstractDataset; is_shuffle::Bool = true, margin::Float64 = 0.9)
    if is_shuffle == true
        return HoldOut(dataset, shuffle(collect(1:length(dataset))), margin)
    else
        return HoldOut(dataset, collect(1:length(dataset)), margin)
    end
end

Base.length(ho::HoldOut) = 1

function getTrainData(holdout::HoldOut{T}, _::Int) :: T where T
    index = findall(r -> r <= length(holdout.dataset) * holdout.k, holdout.index)
    return Persa.sample(holdout.dataset, index)
end

function getTestData(holdout::HoldOut{T}, _::Int) :: T where T
    index = findall(r -> r > length(holdout.dataset) * holdout.k, holdout.index)
    return Persa.sample(holdout.dataset, index)
end


text:
using MultipleViewGeometry, Test, LinearAlgebra
using MultipleViewGeometry.ModuleTypes
using StaticArrays, Calculus, GeometryTypes
using MultipleViewGeometry.ModuleSyntheticData

# Fix random seed.
Random.seed!(1234)
# Construct two camera matrices and parametrise two planar surfaces.
f = 50
image_width = 640 / 10
image_height = 480 / 10
𝐊₁ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]
𝐑₁ = SMatrix{3,3,Float64,9}(rotxyz(0, 25*(pi/180), 0))
𝐭₁ = [-30.0, 0.0, -5.0]

𝐊₂ = @SMatrix [f 0 0 ;
               0 f 0 ;
               0 0 1 ]

𝐑₂ = SMatrix{3,3,Float64,9}(rotxyz(0, -25*(pi/180), 0))
𝐭₂ = [30.0, 0.0, 5.0]

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

# Normals and distance from origin
𝐧₁ = [0.0, 0.0, 1.0]
d₁ = 55.0

𝐧₂ = [0.5, -0.2, 2.0]
d₂ = 145.0

𝒳₁ = generate_planar_points(𝐧₁,d₁, 20, 50)
𝒳₂ = generate_planar_points(𝐧₂,d₂, 20, 50)

world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
picture_basis = (Point(0.0, 0.0), Vec(-1.0, 0.0), Vec(0.0, -1.0))

camera₁ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
camera₂ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
relocate!(camera₁, 𝐑₁, 𝐭₁)
relocate!(camera₂, 𝐑₂, 𝐭₂)

𝐑₁′, 𝐭₁′ = ascertain_pose(camera₁, world_basis... )
𝐊₁′ = obtain_intrinsics(camera₁, CartesianSystem())
𝐑₂′, 𝐭₂′ = ascertain_pose(camera₂, world_basis... )
𝐊₂′ = obtain_intrinsics(camera₂, CartesianSystem())

𝐏₁ = construct(ProjectionMatrix(),𝐊₁′,𝐑₁′,𝐭₁′)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂′,𝐑₂′,𝐭₂′)

# Set of corresponding points.
ℳ₁ = project(camera₁,𝐏₁,𝒳₁)
ℳ₁ʹ= project(camera₂,𝐏₂,𝒳₁)
ℳ₂ = project(camera₁,𝐏₁,𝒳₂)
ℳ₂ʹ= project(camera₂,𝐏₂,𝒳₂)

𝐇₁ = estimate(HomographyMatrix(), DirectLinearTransform(), (ℳ₁, ℳ₁ʹ))
𝐇₂ = estimate(HomographyMatrix(), DirectLinearTransform(), (ℳ₂, ℳ₂ʹ))

for i = zip(ℳ₁,ℳ₁ʹ)
    m, mʹ =  i
    𝐦  = hom(m)
    𝐦ʹ = hom(mʹ)
    residual = vec2antisym(𝐦ʹ)*𝐇₁*𝐦
    @test isapprox(sum(residual), 0.0; atol = 1e-7)
end

for i = zip(ℳ₂,ℳ₂ʹ)
    m, mʹ =  i
    𝐦  = hom(m)
    𝐦ʹ = hom(mʹ)
    residual = vec2antisym(𝐦ʹ)*𝐇₂*𝐦
    @test isapprox(sum(residual), 0.0; atol = 1e-7)
end

# Tests for homography matrix estimation
#
# # Normal to the planar surface.
# 𝐧 = [0.0, 0.0, -1.0]
# # Distance of the plane from the origin.
# d = 100
# # Sample points on the planar surface.
# 𝒳 = [Point3D(x,y,d) for x = -100:5:100 for y = -100:5:100]
# 𝒳 = 𝒳[1:50:end]
# # Intrinsic and extrinsic parameters of camera one.
# 𝐊₁ = Matrix{Float64}(I,3,3)
# 𝐑₁ = Matrix{Float64}(I,3,3)
# 𝐭₁ = [0.0, 0.0, 0.0]
#
# # Intrinsic and extrinsic parameters of camera two.
# 𝐊₂ = Matrix{Float64}(I,3,3)
# 𝐑₂ = Matrix{Float64}(I,3,3)
# 𝐭₂ = [100.0, 2.0, -100.0]
#
# # Camera projection matrices.
# 𝐏₁ = construct(ProjectionMatrix(),𝐊₁,𝐑₁,𝐭₁)
# 𝐏₂ = construct(ProjectionMatrix(),𝐊₂,𝐑₂,𝐭₂)
#
# # Set of corresponding points.
# ℳ = project(Pinhole(),𝐏₁,𝒳)
# ℳʹ = project(Pinhole(),𝐏₂,𝒳)
#
# # Estimate of the homography matrix and the true homography matrix.
# 𝐇 = estimate(HomographyMatrix(), DirectLinearTransform(), (ℳ, ℳʹ))
# 𝐇₀ = construct(HomographyMatrix(),𝐊₁,𝐑₁,𝐭₁,𝐊₂,𝐑₂,𝐭₂,𝐧,d)
#
# 𝐇₀ = 𝐇₀ / norm(𝐇₀)
# 𝐇₀ = 𝐇₀ / sign(𝐇₀[3,3])
#
# 𝐇 = 𝐇 / norm(𝐇)
# 𝐇 = 𝐇 / sign(𝐇[3,3])
#
# for i = zip(ℳ,ℳʹ)
#     m, mʹ =  i
#     𝐦  = hom(m)
#     𝐦ʹ = hom(mʹ)
#     residual = vec2antisym(𝐦ʹ)*𝐇₀*𝐦
#     @test isapprox(sum(residual), 0.0; atol = 1e-7)
# end
#
# for i = zip(ℳ,ℳʹ)
#     m, mʹ =  i
#     𝐦  = hom(m)
#     𝐦ʹ = hom(mʹ)
#     residual = vec2antisym(𝐦ʹ)*𝐇*𝐦
#     @show residual
#     @test isapprox(sum(residual), 0.0; atol = 1e-7)
# end
#

#@test isapprox(sum(residual), 0.0; atol = 1e-7)
#dot(𝒳[10],𝐧) + d


text:
using JLD2, FileIO
using Test

function better_success(cmd)
    fn1, _ = mktemp()
    fn2, _ = mktemp()
    try
       run(pipeline(cmd, stdout=fn1, stderr=fn2))
    catch
        println(String(read(fn1)))
        println(String(read(fn2)))
        return false
    end
    return true
end

include("bufferedreader.jl")
include("lookup3.jl")
include("internal.jl")
include("rw.jl")
include("append.jl")
include("groups.jl")
include("consistency.jl")
include("loadsave.jl")
include("modules.jl")
include("modules-nested.jl")
include("isreconstructed.jl")
include("backwards_compatibility.jl")
include("inlineunion.jl")
include("customserialization.jl")
include("compression.jl")


text:
function sendCmd(cl::TCPSocket, cmd::ASCIIString)
  println(cl, cmd)
  cmd == "QUIT" ?  close(cl) : readline(cl)
end

function getParticles(cl::TCPSocket, lbl::ASCIIString)
  res = sendCmd(cl, "GETPARTICLES $(lbl)")
  rows = split(res[1:(end-1)],';')
  V = readdlm(IOBuffer(rows[1]),',')
  for i in 2:length(rows)
    v = readdlm(IOBuffer(rows[i]),',')
    V = [V;v]
  end
  return V
end

cl = connect(60001)
println("connected")
sendCmd(cl, "INIT")
sendCmd(cl, "ODOMETRY 1 2 10.0 0.0 0.0 0.1 0 0 0.05 0 0.05")
sendCmd(cl, "ODOMETRY 2 3 10.0 0.0 0.0 0.1 0 0 0.05 0 0.05")
Vb = getParticles(cl, "x3");
sendCmd(cl, "SOLVETREE")
Va = getParticles(cl, "x3");
@show size(Va)
sendCmd(cl, "LANDMBR 3 4 $(pi/2) 10.0 0.01 0 0.5")
sendCmd(cl, "SOLVETREE")
l1 = getParticles(cl, "l1");
@show size(l1)

sendCmd(cl, "QUIT")
close(cl)


text:
using Printf: @printf

print("Sum of first 100,000 primes:  ")
println(Iterators.sum(Iterators.take(PrimesPaged(), 100000)))
print("First 20 primes:  ( ")
foreach((p->@printf("%d ", p)), Iterators.take(PrimesPaged(), 20))
println(")")
print("Primes between 100 and 150:  ( ")
for p in Iterators.filter((p->p>=100), PrimesPaged()) p > 150 && break; @printf("%d ", p)) end
println(")")
let cnt = 0
    for p in PrimesPaged()
        p > 8000 && break; if p > 7700 cnt += 1 end
    end; println("Number of primes between 7700 and 8000:  ", cnt)
end
@printf("The 10,000th prime:  %d\n", Iterators.first(Iterators.drop(PrimesPaged(), 9999)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test, NiLangCore
using NiLangCore: get_argname, get_ftype, match_function, MyOrderedDict

@testset "match function" begin
    ex = match_function(:(function f(x) x end))
    @test ex[1] == nothing
    @test ex[2] == :f
    @test ex[3] == [:x]
    @test ex[4] == []
    @test length(filter(x->!(x isa LineNumberNode), ex[5])) == 1
    ex = match_function(:(@inline function f(x; y) x end))
    @test ex[1][1] == Symbol("@inline")
    @test ex[1][2] isa LineNumberNode
    @test ex[2] == :f
    @test ex[3] == [Expr(:parameters, :y), :x]
    @test length(filter(x->!(x isa LineNumberNode), ex[5])) == 1
    @test ex[4] == []
    ex = match_function(:(function f(x::T) where T x end))
    @test ex[2] == :f
    @test ex[3] == [:(x::T)]
    @test length(filter(x->!(x isa LineNumberNode), ex[5])) == 1
    @test ex[4] == [:T]
    ex = match_function(:(f(x)=x))
    @test ex[2] == :f
    @test ex[3] == [:x]
    @test length(ex[5]) == 2
    @test ex[4] == []
end

@testset "argname and type" begin
    @test get_argname(:(y=3)) == :y
    @test get_argname(:(y::Int)) == :y
    @test get_argname(:(y::Int=3)) == :y
    @test get_argname(:(f(; k::Int=4)).args[2]) == :(f(; k::Int=4)).args[2]
end

@testset "my ordered dict" begin
    od = MyOrderedDict{Any, Any}()
    od[:a] = 2
    od[:b] = 4
    od[:c] = 7
    @test length(od) == 3
    @test od[:b] == 4
    od[:b] = 1
    @test od[:b] == 1
    delete!(od, :b)
    @test_throws KeyError od[:b]
    @test pop!(od) == (:c, 7)
    @test length(od) == 1
end

@testset "unzipped broadcast" begin
    x = [1, 2, 3.0]
    res = NiLangCore.unzipped_broadcast(exp, x)
    @test res === x
    @test res ≈ exp.([1, 2, 3.0])
end


text:
using Compat

const hippath = Base.@__DIR__;
const lib_hip_jl = Libdl.dlopen(hippath * "/hip_jl.so")

const hipTranspose = Libdl.dlsym(lib_hip_jl, :tranpose)
const hipAlloc = Libdl.dlsym(lib_hip_jl, :hipArray)
const hipAllocUninitialized = Libdl.dlsym(lib_hip_jl, :hipArrayUninit)
const hipFree = Libdl.dlsym(lib_hip_jl, :hipDelete)
const hipDeviceToHost = Libdl.dlsym(lib_hip_jl, :hipMemDevToHost)
const sync = Libdl.dlsym(lib_hip_jl, :synchronize)
const setDevice = Libdl.dlsym(lib_hip_jl, :setDevice)
const sgemmNaive = Libdl.dlsym(lib_hip_jl, :naive_sgemm)
const sgemmTile = Libdl.dlsym(lib_hip_jl, :tiled_sgemm)
const sgemmHip = Libdl.dlsym(lib_hip_jl, :hipSgemm)

struct HipArray{N} <: AbstractArray{Float32,N}
    ptr::Ptr{Void}
    size::NTuple{N,Int}
    length::Int
end
function HipArray(x::AbstractArray{T,N}) where {T,N}
    HipArray(Array(Float32.(x)))
end
function HipArray(x::Array{Float32,N}) where {N}
    y = ccall(hipAlloc, Ptr{Cvoid}, (Ref{Cfloat}, Cint), x, Cint(length(x)))
    HipArray{N}(y, size(x), length(x))
end
function HipArray(x::NTuple{N,T}) where {N, T<: Integer}
    l = prod(x)
    y = ccall(hipAllocUninitialized, Ptr{Cvoid}, (Cint,), Cint(l))
    HipArray{N}(y, x, l)
end

function transpose!(x::HipArray{2}, y::HipArray{2})
    @assert y.size[1] == x.size[2]
    @assert y.size[2] == x.size[1]
    ccall(hipTranspose, Cvoid, (Ref{Cvoid}, Ref{Cvoid}, Cint, Cint), x.ptr, y.ptr, Cint(y.size[1]), Cint(y.size[2]))
end

function sgemm!(C::HipArray{2}, alpha::Cfloat, A::HipArray{2}, B::HipArray{2},
                    beta::Cfloat, transA = false, transB = false)
    @assert A.size[2] == B.size[1]
    @assert C.size[1] == A.size[1]
    @assert C.size[2] == B.size[2]
    ccall(sgemmHip, Cvoid, (Bool, Bool, Cint, Cint, Cint, Cfloat, Ref{Cvoid}, Ref{Cvoid}, Cfloat, Ref{Cvoid}), transA, transB, Cint(size(A, 1)), Cint(size(B,2)), Cint(size(A,2)), alpha, A.ptr, B.ptr, beta, C.ptr)
end

function naiveMul!(C::HipArray{2}, A::HipArray{2}, B::HipArray{2})
    @assert A.size[2] == B.size[1]
    @assert C.size[1] == A.size[1]
    @assert C.size[2] == B.size[2]
    ccall(sgemmNaive, Cvoid, (Cint, Cint, Cint, Cfloat, Ref{Cvoid}, Ref{Cvoid}, Cfloat, Ref{Cvoid}),
        Cint(size(A, 1)), Cint(size(B,2)), Cint(size(A,2)), 1f0, A.ptr, B.ptr, 0f0, C.ptr)
end
function naivesgemm!(C::HipArray{2}, alpha::Cfloat, A::HipArray{2}, B::HipArray{2}, beta::Cfloat)
    @assert A.size[2] == B.size[1]
    @assert C.size[1] == A.size[1]
    @assert C.size[2] == B.size[2]
    ccall(sgemmNaive, Cvoid, (Cint, Cint, Cint, Cfloat, Ref{Cvoid}, Ref{Cvoid}, Cfloat, Ref{Cvoid}),
        Cint(size(A, 1)), Cint(size(B,2)), Cint(size(A,2)), alpha, A.ptr, B.ptr, beta, C.ptr)
end
function tilesgemm!(C::HipArray{2}, alpha::Cfloat, A::HipArray{2}, B::HipArray{2}, beta::Cfloat)
    @assert A.size[2] == B.size[1]
    @assert C.size[1] == A.size[1]
    @assert C.size[2] == B.size[2]
    ccall(sgemmTile, Cvoid, (Cint, Cint, Cint, Cfloat, Ref{Cvoid}, Ref{Cvoid}, Cfloat, Ref{Cvoid}),
        Cint(size(A, 1)), Cint(size(B,2)), Cint(size(A,2)), alpha, A.ptr, B.ptr, beta, C.ptr)
end

Base.size(y::HipArray) = y.size
Base.length(y::HipArray) = y.length
Base.show(io::IO, y::HipArray) = println(io, "HipArray of size: $(y.size).")
Base.getindex(y::HipArray, i) = throw("Get index not defined for HipArray.")
Base.setindex!(y::HipArray, i, v) = throw("Set index not defined for HipArray.")

synchronize() = ccall(sync, Cvoid, ())
function set_device(i)
    ccall(setDevice, Cvoid, (Cint,), Cint(i))
end

function free!(y::HipArray)
    ccall(hipFree, Cvoid, (Ref{Cvoid},), y.ptr)
end
function Base.copy!(x::Array{Float32}, y::HipArray)#dims do not have to match.
    @assert length(x) == length(y)#but lengths do.
    ccall(hipDeviceToHost, Cvoid, (Ref{Cfloat}, Ref{Cvoid}, Cint), x, y.ptr, Cint(length(y)))
end
function Base.Array(y::HipArray{N}) where N
    out = Array{Float32,N}(uninitialized, y.size...)
    ccall(hipDeviceToHost, Cvoid, (Ref{Cfloat}, Ref{Cvoid}, Cint), out, y.ptr, Cint(length(y)))
    out
end


text:
export Identified, LeftIdentified, RightIdentitified, Unidentified,
       identification

"""
Trait to indicate that a binary operation • over set S indeed has an identity
element. Only methods of • with the signature •(x::S, y::S) are to be
considered.

That is to say, ``\exists I \in S. \forall x ∈ S. I \bullet x = x \bullet I =
x``. Note that an identity exists if and only if both left- and right-identities
exist (and if so, they are necessarily equal).
"""
abstract Identified

"""
Trait to indicate that a binary operation • over set S has a left-identity
element. Only methods of • with the signature •(x::S, y::S) are to be
considered.

That is to say, ``\exists I \in S. \forall x ∈ S. I \bullet x = x``. Note that
this trait does not imply that the left-identity could not be a right-identity
(and hence an identity element) also.
"""
abstract LeftIdentitified

"""
Trait to indicate that a binary operation • over set S has a right-identity
element. Only methods of • with the signature •(x::S, y::S) are to be
considered.

That is to say, ``\exists I \in S. \forall x ∈ S. x \bullet I = x``. Note that
this trait does not imply that the right-identity could not be a left-identity
(and hence an identity element) also.
"""
abstract RightIdentitified

"""
Trait to indicate that a binary operation • over set S is not known to have a
(possibly one-directional) identity element. Note that this trait does not imply
there there could not possibly be an identity element.
"""
abstract Unidentified

"""
Return the trait describing the existence of a full identity element of the set
S under the binary operation •.
"""
identification(::Type, ::Function) = WithoutIdentity

# Assume all number types have 0 and 1 (possibly bad assumption)
identification{T<:Number}(::Type{T}, ::typeof(+)) = Identified
identification{T<:Number}(::Type{T}, ::typeof(*)) = Identified

# Empty string is an identity element
identification(::Type{String}, ::typeof(*)) = Identified


text:
using Test

@testset "WrappedRNG" begin
    using RandomNumbers.Xorshifts

    @info "Testing WrappedRNG"
    seed = RandomNumbers.gen_seed(UInt64, 2)
    r = Xoroshiro128(seed)

    r1 = WrappedRNG(r, UInt32)
    r2 = WrappedRNG(Xoroshiro128, UInt32, seed)

    @test rand(r1, UInt64) == rand(r2, UInt64)

    r3 = WrappedRNG(r2, UInt128)

    rand(r2, UInt64)

    @test rand(r3, UInt128) == rand(r2, UInt128)

    @test copyto!(copy(r1), r1) == r1
end


text:
recA(::Type{T},::Chebyshev,k) where {T} = k == 0 ? one(T) : 2one(T)
recB(::Type{T},::Chebyshev,_) where {T} = zero(T)
recC(::Type{T},::Chebyshev,k) where {T} = one(T)   # one(T) ensures we get correct type

recα(::Type{T},::Chebyshev,_) where {T} = zero(T)
recβ(::Type{T},::Chebyshev,k) where {T} = ifelse(k==1,one(T),one(T)/2)   # one(T) ensures we get correct type,ifelse ensures inlining
recγ(::Type{T},::Chebyshev,k) where {T} = one(T)/2   # one(T) ensures we get correct type

## Evaluation

Evaluation(S::Chebyshev,x::typeof(leftendpoint),o::Integer) =
    ConcreteEvaluation(S,x,o)
Evaluation(S::Chebyshev,x::typeof(rightendpoint),o::Integer) =
    ConcreteEvaluation(S,x,o)

Evaluation(S::Chebyshev,x::Number,o::Integer) =
    o==0 ? ConcreteEvaluation(S,x,o) : EvaluationWrapper(S,x,o,Evaluation(x)*Derivative(S,o))

function evaluatechebyshev(n::Integer,x::T) where T<:Number
    if n == 1
        [one(T)]
    else
        p = zeros(T,n)
        p[1] = one(T)
        p[2] = x

        for j=2:n-1
            p[j+1] = 2x*p[j] - p[j-1]
        end

        p
    end
end

function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
    T=eltype(op)
    if op.order == 0
        ifelse(isodd(j),  # right rule
            one(T),
            -one(T))
    else
        #TODO: Fast version
        op[j:j][1]
    end
end

function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
    T=eltype(op)
    if op.order == 0
        one(T)
    else
        #TODO: Fast version
        op[j:j][1]
    end
end

function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
    T=eltype(op)
    x = op.x
    d = domain(op)
    p = op.order
    cst = convert(T,(2/complexlength(d))^p)
    n=length(k)

    ret = Array{T}(undef, n)
    k1=1-first(k)
    @simd for j=k
        @inbounds ret[j+k1]=(-1)^(p+1)*(-one(T))^j
    end

    for m=0:p-1
        k1=1-first(k)
        @simd for j=k
            @inbounds ret[j+k1] *= (j-1)^2-m^2
        end
        scal!(convert(T,1/(2m+1)), ret)
    end

    scal!(cst,ret)
end

function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractRange) where {DD<:IntervalOrSegment,RR}
    T=eltype(op)
    x = op.x
    d = domain(op)
    p = op.order
    cst = convert(T,(2/complexlength(d))^p)
    n=length(k)

    ret = fill(one(T),n)

    for m=0:p-1
        k1=1-first(k)
        @simd for j=k
            @inbounds ret[j+k1] *= (j-1)^2-m^2
        end
        scal!(convert(T,1/(2m+1)), ret)
    end

    scal!(cst,ret)
end

function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
                j::Integer) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
    if op.order == 0
        convert(T,evaluatechebyshev(j,tocanonical(domain(op),op.x))[end])
    else
        error("Only zero–second order implemented")
    end
end

function getindex(op::ConcreteEvaluation{Chebyshev{DD,RR},M,OT,T},
                                             k::AbstractRange) where {DD<:IntervalOrSegment,RR,M<:Real,OT,T}
    if op.order == 0
        Array{T}(evaluatechebyshev(k[end],tocanonical(domain(op),op.x))[k])
    else
        error("Only zero–second order implemented")
    end
end

function Dirichlet(S::Chebyshev,order)
    order == 0 && return ConcreteDirichlet(S,ArraySpace([ConstantSpace.(Point.(endpoints(domain(S))))...]),0)
    default_Dirichlet(S,order)
end

function getindex(op::ConcreteDirichlet{<:Chebyshev},
                                             k::Integer,j::Integer)
    if op.order == 0
        k == 1 && iseven(j) && return -one(eltype(op))
        return one(eltype(op))
    else
        error("Only zero Dirichlet conditions implemented")
    end
end

function Matrix(S::SubOperator{T,ConcreteDirichlet{C,V,T},
                                Tuple{UnitRange{Int},UnitRange{Int}}}) where {C<:Chebyshev,V,T}
    ret = Array{T}(undef, size(S)...)
    kr,jr = parentindices(S)
    isempty(kr) && return ret
    isempty(jr) && return ret
    if first(kr) == 1
        if isodd(jr[1])
            ret[1,1:2:end] .= one(T)
            ret[1,2:2:end] .= -one(T)
        else
            ret[1,1:2:end] .= -one(T)
            ret[1,2:2:end] .= one(T)
        end
    end
    if last(kr) == 2
        ret[end,:] .= one(T)
    end
    return ret
end
#

# Multiplication

Base.stride(M::ConcreteMultiplication{U,V}) where {U<:Chebyshev,V<:Chebyshev} =
    stride(M.f)

getindex(M::ConcreteMultiplication{C,C,T},k::Integer,j::Integer) where {T,C<:Chebyshev} =
    chebmult_getindex(coefficients(M.f),k,j)

getindex(M::ConcreteMultiplication{C,PS,T},k::Integer,j::Integer) where {PS<:PolynomialSpace,T,C<:Chebyshev} =
    M[k:k,j:j][1,1]

function BandedMatrix(S::SubOperator{T,ConcreteMultiplication{C,C,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {C<:Chebyshev,T}
    ret = BandedMatrix(Zeros, S)

    kr,jr=parentindices(S)
    cfs=parent(S).f.coefficients

    isempty(cfs) && return ret

    # Toeplitz part
    sym_toeplitz_axpy!(1.0,0.5,cfs,kr,jr,ret)

    #Hankel part
    hankel_axpy!(0.5,cfs,kr,jr,ret)

    # divide first row by half
    if first(kr)==1
        if first(jr)==1
            ret[1,1]+=0.5cfs[1]
        end

        for j=1:min(1+ret.u,size(ret,2))
            ret[1,j]/=2
        end
    end

    ret
end

## Derivative

Derivative(sp::Chebyshev{DD},order::Integer) where {DD<:IntervalOrSegment} =
    ConcreteDerivative(sp,order)

rangespace(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} =
    Ultraspherical(D.order,domain(D))

bandwidths(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} = -D.order,D.order
Base.stride(D::ConcreteDerivative{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} = D.order

function getindex(D::ConcreteDerivative{Chebyshev{DD,RR},K,T},k::Integer,j::Integer) where {DD<:IntervalOrSegment,RR,K,T}
    m=D.order
    d=domain(D)

    if j==k+m
        C=convert(T,pochhammer(one(T),m-1)/2*(4/complexlength(d))^m)
        convert(T,C*(m+k-one(T)))
    else
        zero(T)
    end
end

linesum(f::Fun{Chebyshev{DD,RR}}) where {DD<:IntervalOrSegment,RR} =
    sum(setcanonicaldomain(f))*arclength(domain(f))/2

## Clenshaw-Curtis functional

for (Func,Len) in ((:DefiniteIntegral,:complexlength),(:DefiniteLineIntegral,:arclength))
    ConcFunc = Meta.parse("Concrete"*string(Func))
    @eval begin
        $Func(S::Chebyshev{D}) where {D<:IntervalOrSegment} = $ConcFunc(S)
        function getindex(Σ::$ConcFunc{Chebyshev{D,R},T},k::Integer) where {D<:IntervalOrSegment,R,T}
            d = domain(Σ)
            C = $Len(d)/2

            isodd(k) ? convert(T,2C/(k*(2-k))) : zero(T)
        end
    end
end

ReverseOrientation(S::Chebyshev) = ReverseOrientationWrapper(NegateEven(S,reverseorientation(S)))
Reverse(S::Chebyshev) = ReverseWrapper(NegateEven(S,S))


text:
# IntervalStreams
# ===============
#
# Often we'd rather avoid reading a large interval dataset into an
# IntervalCollection. We can still do efficient intersections if we can assume
# the data is sorted

type IntervalStreamIntersectIterator{S, T, SS, TS}
    a::SS
    b::TS
    seqname_isless::Function

    # buffered intervals from a
    a_buffer::StreamBuffer{Interval{S}}
    b_interval::Interval{T}

    function IntervalStreamIntersectIterator(a::SS,
                                             b::TS,
                                             seqname_isless::Function)
        return new(a, b, seqname_isless, StreamBuffer{Interval{S}}())
    end
end

Base.iteratorsize(::IntervalStreamIntersectIterator) = Base.SizeUnknown()

immutable IntervalStreamIntersectIteratorState{U, V}
    a_buffer_pos::Int
    a_state::U
    b_state::V
end

"""
Intersect two `IntervalStreams` returning an iterator over all pairs of
intersecting intervals.
"""
function Base.intersect{SS<:IntervalStreamOrArray,TS<:IntervalStreamOrArray}(
        a::SS, b::TS, seqname_isless::Function=isless)
    S = metadatatype(a)
    T = metadatatype(b)
    return IntervalStreamIntersectIterator{S, T, SS, TS}(a, b, seqname_isless)
end

function first_intersection{S, T, SS, TS, U, V}(
                                        a::SS,
                                        b::TS,
                                        a_state::U,
                                        b_state::V,
                                        a_interval::Interval{S},
                                        b_interval::Interval{T},
                                        seqname_isless::Function)

    while !done(a, a_state) || !done(b, b_state)
        if precedes(a_interval, b_interval, seqname_isless) && !done(a, a_state)
            a_interval_next, a_state = next(a, a_state)
            if isless(a_interval_next, a_interval, seqname_isless)
                error("Interval streams must be sorted to perform intersection.")
            end
            a_interval = a_interval_next
        elseif precedes(b_interval, a_interval, seqname_isless) && !done(b, b_state)
            b_interval_next, b_state = next(b, b_state)
            if isless(b_interval_next, b_interval, seqname_isless)
                error("Interval streams must be sorted to perform intersection.")
            end
            b_interval = b_interval_next
        else # intersection
            break
        end
    end

    return a_interval, b_interval, a_state, b_state
end

function Base.start{S,T,SS,TS}(it::IntervalStreamIntersectIterator{S,T,SS,TS})
    a_state = start(it.a)
    b_state = start(it.b)

    if done(it.a, a_state) || done(it.b, b_state)
        return IntervalStreamIntersectIteratorState(a_state, b_state)
    end

    a_interval, a_state = next(it.a, a_state)
    b_interval, b_state = next(it.b, b_state)
    a_interval, b_interval, a_state, b_state = first_intersection(
        it.a, it.b, a_state, b_state, a_interval, b_interval, it.seqname_isless)

    if isoverlapping(a_interval, b_interval)
        push!(it.a_buffer, a_interval)
        it.b_interval = b_interval
    end

    return IntervalStreamIntersectIteratorState(1, a_state, b_state)
end

# Do we need a state if this is going to be stateful. State can maybe just
# be the current interval in b
function Base.next{S,T,SS,TS,U,V}(it::IntervalStreamIntersectIterator{S,T,SS,TS},
                                  state::IntervalStreamIntersectIteratorState{U,V})
    value = (it.a_buffer[state.a_buffer_pos], it.b_interval)
    a_state = state.a_state
    b_state = state.b_state

    # another intersection in a_buffer?
    a_buffer_pos = state.a_buffer_pos + 1
    while a_buffer_pos <= length(it.a_buffer)
        if precedes(it.b_interval, it.a_buffer[a_buffer_pos], it.seqname_isless)
            a_buffer_pos = length(it.a_buffer) + 1
        elseif isoverlapping(it.b_interval, it.a_buffer[a_buffer_pos])
            break
        else
            a_buffer_pos += 1
        end
    end

    # look for an intersection for b_interval by reading more entries from a
    b_interval = it.b_interval
    if a_buffer_pos > length(it.a_buffer) && !done(it.a, a_state)
        while true
            a_interval, a_state = next(it.a, a_state)
            if !isordered(it.a_buffer[end], a_interval, it.seqname_isless)
                error("Interval streams must be sorted to perform intersection.")
            end
            push!(it.a_buffer, a_interval)

            if isoverlapping(b_interval, it.a_buffer[end]) ||
               precedes(b_interval, it.a_buffer[end], it.seqname_isless) ||
               done(it.a, a_state)
                break
            end
            a_buffer_pos += 1
        end

        if !isoverlapping(it.a_buffer[a_buffer_pos], b_interval)
            a_buffer_pos += 1
        end
    end

    # no more intersections against b_interval, try next interval from b
    while (a_buffer_pos > length(it.a_buffer) ||
           !isoverlapping(it.a_buffer[a_buffer_pos], b_interval)) && !done(it.b, b_state)

        b_interval_next, b_state = next(it.b, b_state)
        if !isordered(b_interval, b_interval_next, it.seqname_isless)
            error("Interval streams must be sorted to perform intersection.")
        end
        b_interval = b_interval_next

        # evict as much as possible from a_buffer
        a_interval = it.a_buffer[1]
        while !isempty(it.a_buffer) && precedes(it.a_buffer[1], b_interval, it.seqname_isless)
            a_interval = shift!(it.a_buffer)
        end

        if isempty(it.a_buffer)
            #a_interval, b_interval, a_state, b_state = first_intersection(
                #it.a, it.b, a_state, b_state, a_interval, b_interval, it.seqname_isless)
            # See JuliaLang/julia#11181
            first_intersection_value = first_intersection(
                it.a, it.b, a_state, b_state, a_interval, b_interval, it.seqname_isless)
            a_interval = first_intersection_value[1]
            b_interval = first_intersection_value[2]
            a_state = first_intersection_value[3]
            b_state = first_intersection_value[4]

            it.b_interval = b_interval
            if isoverlapping(a_interval, b_interval)
                push!(it.a_buffer, a_interval)
                a_buffer_pos = 1
            end
            break
        elseif isoverlapping(it.a_buffer[1], b_interval)
            a_buffer_pos = 1
            it.b_interval = b_interval
        end
    end

    return value, IntervalStreamIntersectIteratorState{U, V}(
                        a_buffer_pos, a_state, b_state)
end

function Base.done{S,T,SS,TS,U,V}(it::IntervalStreamIntersectIterator{S,T,SS,TS},
                                  state::IntervalStreamIntersectIteratorState{U,V})
    return state.a_buffer_pos > length(it.a_buffer)
end

# TODO: IntervalCollection(stream::IntervalStream) constructor.

# Helper function for coverage. Process remaining interval end points after
# all intervals have been read.
function coverage_process_lasts_heap!(cov::IntervalCollection{UInt32},
                                      current_coverage, coverage_seqname,
                                      coverage_first, lasts)
    while !isempty(lasts)
        pos = heappop!(lasts)
        if pos == coverage_first - 1
            current_coverage -= 1
        else
            @assert pos >= coverage_first
            push!(cov, Interval{UInt32}(coverage_seqname, coverage_first,
                                        pos, STRAND_BOTH, current_coverage))
            current_coverage -= 1
            coverage_first = pos + 1
        end
    end
    @assert current_coverage == 0
end

"""
    coverage(intervals)

Compute the coverage of a collection of intervals and return an
`IntervalCollection` that contains run-length encoded coverage data.

For example, given intervals like:

    [------]     [------------]
       [---------------]

This function would return a new set of disjoint intervals with annotated
coverage like:

    [1][-2-][-1-][--2--][--1--]
"""
function coverage(stream::Union{IntervalStreamOrArray,IntervalTrees.IntervalTree},
                  seqname_isless::Function=isless)
    cov = IntervalCollection{UInt32}()
    lasts = Int64[]

    stream_state = start(stream)
    if done(stream, stream_state)
        return cov
    end

    current_coverage = 0
    coverage_seqname = ""
    coverage_first = 0
    last_interval_first = 0
    interval, stream_state = next(stream, stream_state)
    while true
        if interval.seqname != coverage_seqname
            coverage_process_lasts_heap!(cov, current_coverage, coverage_seqname,
                                         coverage_first, lasts)
            if !(isempty(coverage_seqname) || seqname_isless(coverage_seqname, interval.seqname))
                error("Intervals must be sorted to compute coverage.")
            end

            coverage_seqname = interval.seqname
            current_coverage = 0
            coverage_first = 0
            last_interval_first = 0
        end

        if interval.first < last_interval_first
            error("Intervals must be sorted to compute coverage.")
        end

        if !isempty(lasts) && lasts[1] < first(interval)
            pos = heappop!(lasts)
            if first(interval) == pos + 1
                heappush!(lasts, last(interval))
                if done(stream, stream_state)
                    break
                end
                last_interval_first = first(interval)
                interval, stream_state = next(stream, stream_state)
            elseif pos == coverage_first - 1
                current_coverage -= 1
            else
                @assert pos >= coverage_first
                push!(cov, Interval{UInt32}(coverage_seqname, coverage_first,
                                            pos, STRAND_BOTH, current_coverage))
                current_coverage -= 1
                coverage_first = pos + 1
            end
        else
            if coverage_first == 0
                coverage_first = first(interval)
                current_coverage = 1
            elseif coverage_first == first(interval)
                current_coverage += 1
            else
                if current_coverage > 0
                    push!(cov, Interval{UInt32}(coverage_seqname, coverage_first,
                                                first(interval) - 1, STRAND_BOTH,
                                                current_coverage))
                end
                current_coverage += 1
                coverage_first = first(interval)
            end

            heappush!(lasts, last(interval))
            if done(stream, stream_state)
                break
            end
            last_interval_first = first(interval)
            interval, stream_state = next(stream, stream_state)
        end
    end

    coverage_process_lasts_heap!(cov, current_coverage, coverage_seqname,
                                 coverage_first, lasts)

    return cov
end

function coverage(ic::IntervalCollection)
    return coverage(ic, isless)
end


text:
#glut.h
@getf libglut mgl_create_graph_glut(draw::DrawCallback, title::Ptr{Cchar}, par::Ptr{Void}, load::LoadCallback)::HMGL
@getf libglut mgl_glut_toggle_alpha(gr::HMGL)::Void
@getf libglut mgl_glut_toggle_light(gr::HMGL)::Void
@getf libglut mgl_glut_toggle_no(gr::HMGL)::Void
@getf libglut mgl_glut_update(gr::HMGL)::Void
@getf libglut mgl_glut_reload(gr::HMGL)::Void
@getf libglut mgl_glut_next_frame(gr::HMGL)::Void
@getf libglut mgl_glut_prev_frame(gr::HMGL)::Void
@getf libglut mgl_glut_animation(gr::HMGL)::Void


text:
using AdventOfCode2019
using Test

@testset "Day 1" begin
    @test AdventOfCode2019.Day01.day01() == [3442987, 5161601]
end

@testset "Day 2" begin
    @test AdventOfCode2019.Day02.day02() == [3101844, 8478]
end

@testset "Day 3" begin
    inp1 = "R75,D30,R83,U83,L12,D49,R71,U7,L72\nU62,R66,U55,R34,D71,R55,D58,R83"
    @test AdventOfCode2019.Day03.day03(inp1) == [159, 610]
    inp2 = "R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51\nU98,R91,D20,R16,D67,R40,U7,R15,U6,R7"
    @test AdventOfCode2019.Day03.day03(inp2) == [135, 410]
    @test AdventOfCode2019.Day03.day03() == [1983, 107754]
end

@testset "Day 4" begin
    @test AdventOfCode2019.Day04.day04() == [1675, 1142]
end

@testset "Day 5" begin
    inp = parse.(Int, split("3,9,8,9,10,9,4,9,99,-1,8", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [8])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [9])[end] == 0

    inp = parse.(Int, split("3,9,7,9,10,9,4,9,99,-1,8", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [7])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [8])[end] == 0
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [9])[end] == 0

    inp = parse.(Int, split("3,3,1108,-1,8,3,4,3,99", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [8])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [9])[end] == 0

    inp = parse.(Int, split("3,3,1107,-1,8,3,4,3,99", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [7])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [8])[end] == 0
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [9])[end] == 0

    inp = parse.(Int, split("3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [0])[end] == 0
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [-1])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [1])[end] == 1

    inp = parse.(Int, split("3,3,1105,-1,9,1101,0,0,12,4,12,99,1", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [0])[end] == 0
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [-1])[end] == 1
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [1])[end] == 1

    inp = parse.(Int, split("3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99", ","))
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [7])[end] == 999
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [8])[end] == 1000
    @test AdventOfCode2019.IntCode.run_program!(copy(inp), [9])[end] == 1001

    @test AdventOfCode2019.Day05.day05() == [7566643, 9265694]
end

@testset "Day 6" begin
    @test AdventOfCode2019.Day06.day06("COM)B\nB)C\nC)D\nD)E\nE)F\nB)G\nG)H\nD)I\nE)J\nJ)K\nK)L\nK)YOU\nI)SAN") == [54, 4]
    @test AdventOfCode2019.Day06.day06() == [194721, 316]
end

@testset "Day 7" begin
    inp = parse.(Int, split("3,15,3,16,1002,16,10,16,1,16,15,15,4,15,99,0,0", ","))
    @test AdventOfCode2019.Day07.run_amplifiers(inp, [4, 3, 2, 1, 0], 0) == 43210

    inp = parse.(Int, split("3,23,3,24,1002,24,10,24,1002,23,-1,23,101,5,23,23,1,24,23,23,4,23,99,0,0", ","))
    @test AdventOfCode2019.Day07.run_amplifiers(inp, [0, 1, 2, 3, 4], 0) == 54321

    inp = parse.(Int, split("3,31,3,32,1002,32,10,32,1001,31,-2,31,1007,31,0,33,1002,33,7,33,1,33,31,31,1,32,31,31,4,31,99,0,0,0", ","))
    @test AdventOfCode2019.Day07.run_amplifiers(inp, [1, 0, 4, 3, 2], 0) == 65210

    inp = parse.(Int, split("3,26,1001,26,-4,26,3,27,1002,27,2,27,1,27,26,27,4,27,1001,28,-1,28,1005,28,6,99,0,0,5", ","))
    @test AdventOfCode2019.Day07.run_feedback_loop(inp, [9, 8, 7, 6, 5], 0) == 139629729

    inp = parse.(Int, split("3,52,1001,52,-5,52,3,53,1,52,56,54,1007,54,5,55,1005,55,26,1001,54,-5,54,1105,1,12,1,53,54,53,1008,54,0,55,1001,55,1,55,2,53,55,53,4,53,1001,56,-1,56,1005,56,6,99,0,0,0,0,10", ","))
    @test AdventOfCode2019.Day07.run_feedback_loop(inp, [9, 7, 8, 5, 6], 0) == 18216

    @test AdventOfCode2019.Day07.day07() == [844468, 4215746]
end

@testset "Day 8" begin
    part1, imagecode = AdventOfCode2019.Day08.day08()
    @test part1 == 2375
    image = AdventOfCode2019.Day08.generate_image(imagecode)
    expected = "███  █  █ █  █ ███  █   █\n" *
               "█  █ █ █  █  █ █  █ █   █\n" *
               "█  █ ██   ████ █  █  █ █ \n" *
               "███  █ █  █  █ ███    █  \n" *
               "█ █  █ █  █  █ █ █    █  \n" *
               "█  █ █  █ █  █ █  █   █  \n"
    @test image == expected
end

@testset "Day 9" begin
    for T in (Int, BigInt)
        inp = parse.(T, split("109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99", ","))
        @test AdventOfCode2019.IntCode.run_program!(inp, nothing, nothing, nothing) == [109,1,204,-1,1001,100,1,100,1008,100,16,101,1006,101,0,99]

        inp = parse.(T, split("1102,34915192,34915192,7,4,7,99,0", ","))
        @test AdventOfCode2019.IntCode.run_program!(inp, nothing, nothing, nothing)[end] == 1219070632396864

        inp = parse.(T, split("104,1125899906842624,99", ","))
        @test AdventOfCode2019.IntCode.run_program!(inp, nothing, nothing, nothing)[end] == 1125899906842624
    end
    @test AdventOfCode2019.Day09.day09() == [3601950151, 64236]
end

@testset "Day 10" begin
    inp = ".#..#\n" *
          ".....\n" *
          "#####\n" *
          "....#\n" *
          "...##\n"
    map = AdventOfCode2019.Day10._generate_map(inp)
    ast, count = AdventOfCode2019.Day10._evaluate_best_asteroid(map)
    @test AdventOfCode2019.Day10._to_puzzle_coordinates(ast) == CartesianIndex(3, 4)
    @test count == 8

    inp = "......#.#.\n" *
          "#..#.#....\n" *
          "..#######.\n" *
          ".#.#.###..\n" *
          ".#..#.....\n" *
          "..#....#.#\n" *
          "#..#....#.\n" *
          ".##.#..###\n" *
          "##...#..#.\n" *
          ".#....####\n"
    map = AdventOfCode2019.Day10._generate_map(inp)
    ast, count = AdventOfCode2019.Day10._evaluate_best_asteroid(map)
    @test AdventOfCode2019.Day10._to_puzzle_coordinates(ast) == CartesianIndex(5, 8)
    @test count == 33

    inp = "#.#...#.#.\n" *
          ".###....#.\n" *
          ".#....#...\n" *
          "##.#.#.#.#\n" *
          "....#.#.#.\n" *
          ".##..###.#\n" *
          "..#...##..\n" *
          "..##....##\n" *
          "......#...\n" *
          ".####.###.\n"
    map = AdventOfCode2019.Day10._generate_map(inp)
    ast, count = AdventOfCode2019.Day10._evaluate_best_asteroid(map)
    @test AdventOfCode2019.Day10._to_puzzle_coordinates(ast) == CartesianIndex(1, 2)
    @test count == 35

    inp = ".#..#..###\n" *
          "####.###.#\n" *
          "....###.#.\n" *
          "..###.##.#\n" *
          "##.##.#.#.\n" *
          "....###..#\n" *
          "..#.#..#.#\n" *
          "#..#.#.###\n" *
          ".##...##.#\n" *
          ".....#.#..\n"
    map = AdventOfCode2019.Day10._generate_map(inp)
    ast, count = AdventOfCode2019.Day10._evaluate_best_asteroid(map)
    @test AdventOfCode2019.Day10._to_puzzle_coordinates(ast) == CartesianIndex(6, 3)
    @test count == 41

    inp = ".#..##.###...#######\n" *
          "##.############..##.\n" *
          ".#.######.########.#\n" *
          ".###.#######.####.#.\n" *
          "#####.##.#.##.###.##\n" *
          "..#####..#.#########\n" *
          "####################\n" *
          "#.####....###.#.#.##\n" *
          "##.#################\n" *
          "#####.##.###..####..\n" *
          "..######..##.#######\n" *
          "####.##.####...##..#\n" *
          ".#####..#.######.###\n" *
          "##...#.##########...\n" *
          "#.##########.#######\n" *
          ".####.#.###.###.#.##\n" *
          "....##.##.###..#####\n" *
          ".#.#.###########.###\n" *
          "#.#.#.#####.####.###\n" *
          "###.##.####.##.#..##\n"
    map = AdventOfCode2019.Day10._generate_map(inp)
    ast, count = AdventOfCode2019.Day10._evaluate_best_asteroid(map)
    @test AdventOfCode2019.Day10._to_puzzle_coordinates(ast) == CartesianIndex(11, 13)
    @test count == 210
    ns = [1, 2, 3, 10, 20, 50, 100, 199, 200, 201, 299]
    expected = [CartesianIndex(a) for a in [(11,12), (12,1), (12,2), (12,8), (16,0), (16,9), (10,16), (9,6), (8,2), (10,9), (11,1)]]
    for (n, e) in zip(ns, expected)
        res = AdventOfCode2019.Day10._vaporize(map, ast, n)
        @test AdventOfCode2019.Day10._to_puzzle_coordinates(res) == e
    end

    @test AdventOfCode2019.Day10.day10() == [329, 512]
end

@testset "Day 11" begin
    p1, p2 = AdventOfCode2019.Day11.day11()
    @test p1 == 2088
    img = " █  █ ███   ██   ██  ████ █     ██  ███    \n" *
          " █  █ █  █ █  █ █  █ █    █    █  █ █  █   \n" *
          " █  █ █  █ █    █  █ ███  █    █    █  █   \n" *
          " █  █ ███  █    ████ █    █    █    ███    \n" *
          " █  █ █ █  █  █ █  █ █    █    █  █ █      \n" *
          "  ██  █  █  ██  █  █ █    ████  ██  █      \n"
    @test AdventOfCode2019.Day08.generate_image(p2) == img
end

@testset "Day 12" begin
    inp = "<x=-1, y=0, z=2>\n" *
          "<x=2, y=-10, z=-7>\n" *
          "<x=4, y=-8, z=8>\n" *
          "<x=3, y=5, z=-1>"
    @test AdventOfCode2019.Day12.day12(inp) == [183, 2772]

    inp = "<x=-8, y=-10, z=0>\n" *
          "<x=5, y=5, z=10>\n" *
          "<x=2, y=-7, z=3>\n" *
          "<x=9, y=-8, z=-3>"
    @test AdventOfCode2019.Day12.day12(inp) == [14645, 4686774924]

    @test AdventOfCode2019.Day12.day12() == [10635, 583523031727256]
end

@testset "Day 13" begin
    @test AdventOfCode2019.Day13.day13() == [363, 17159]
end

@testset "Day 14" begin
    inp = "10 ORE => 10 A\n" *
          "1 ORE => 1 B\n" *
          "7 A, 1 B => 1 C\n" *
          "7 A, 1 C => 1 D\n" *
          "7 A, 1 D => 1 E\n" *
          "7 A, 1 E => 1 FUEL"
    @test day14(inp) == [31, 34482758620]

    inp = "9 ORE => 2 A\n" *
          "8 ORE => 3 B\n" *
          "7 ORE => 5 C\n" *
          "3 A, 4 B => 1 AB\n" *
          "5 B, 7 C => 1 BC\n" *
          "4 C, 1 A => 1 CA\n" *
          "2 AB, 3 BC, 4 CA => 1 FUEL"
    @test day14(inp) == [165, 6323777403]

    inp = "157 ORE => 5 NZVS\n" *
          "165 ORE => 6 DCFZ\n" *
          "44 XJWVT, 5 KHKGT, 1 QDVJ, 29 NZVS, 9 GPVTF, 48 HKGWZ => 1 FUEL\n" *
          "12 HKGWZ, 1 GPVTF, 8 PSHF => 9 QDVJ\n" *
          "179 ORE => 7 PSHF\n" *
          "177 ORE => 5 HKGWZ\n" *
          "7 DCFZ, 7 PSHF => 2 XJWVT\n" *
          "165 ORE => 2 GPVTF\n" *
          "3 DCFZ, 7 NZVS, 5 HKGWZ, 10 PSHF => 8 KHKGT"
    @test day14(inp) == [13312, 82892753]

    inp = "2 VPVL, 7 FWMGM, 2 CXFTF, 11 MNCFX => 1 STKFG\n" *
          "17 NVRVD, 3 JNWZP => 8 VPVL\n" *
          "53 STKFG, 6 MNCFX, 46 VJHF, 81 HVMC, 68 CXFTF, 25 GNMV => 1 FUEL\n" *
          "22 VJHF, 37 MNCFX => 5 FWMGM\n" *
          "139 ORE => 4 NVRVD\n" *
          "144 ORE => 7 JNWZP\n" *
          "5 MNCFX, 7 RFSQX, 2 FWMGM, 2 VPVL, 19 CXFTF => 3 HVMC\n" *
          "5 VJHF, 7 MNCFX, 9 VPVL, 37 CXFTF => 6 GNMV\n" *
          "145 ORE => 6 MNCFX\n" *
          "1 NVRVD => 8 CXFTF\n" *
          "1 VJHF, 6 MNCFX => 4 RFSQX\n" *
          "176 ORE => 6 VJHF"
    @test day14(inp) == [180697, 5586022]

    inp = "171 ORE => 8 CNZTR\n" *
          "7 ZLQW, 3 BMBT, 9 XCVML, 26 XMNCP, 1 WPTQ, 2 MZWV, 1 RJRHP => 4 PLWSL\n" *
          "114 ORE => 4 BHXH\n" *
          "14 VRPVC => 6 BMBT\n" *
          "6 BHXH, 18 KTJDG, 12 WPTQ, 7 PLWSL, 31 FHTLT, 37 ZDVW => 1 FUEL\n" *
          "6 WPTQ, 2 BMBT, 8 ZLQW, 18 KTJDG, 1 XMNCP, 6 MZWV, 1 RJRHP => 6 FHTLT\n" *
          "15 XDBXC, 2 LTCX, 1 VRPVC => 6 ZLQW\n" *
          "13 WPTQ, 10 LTCX, 3 RJRHP, 14 XMNCP, 2 MZWV, 1 ZLQW => 1 ZDVW\n" *
          "5 BMBT => 4 WPTQ\n" *
          "189 ORE => 9 KTJDG\n" *
          "1 MZWV, 17 XDBXC, 3 XCVML => 2 XMNCP\n" *
          "12 VRPVC, 27 CNZTR => 2 XDBXC\n" *
          "15 KTJDG, 12 BHXH => 5 XCVML\n" *
          "3 BHXH, 2 VRPVC => 7 MZWV\n" *
          "121 ORE => 7 VRPVC\n" *
          "7 XCVML => 6 RJRHP\n" *
          "5 BHXH, 4 VRPVC => 5 LTCX"
    @test day14(inp) == [2210736, 460664]

    @test day14() == [371695, 4052920]
end

@testset "Day 15" begin
    @test day15() == [374, 482]
end

@testset "Day 16" begin
    for (inp, out) in zip(("80871224585914546619083218645595",
                           "19617804207202209144916044189917",
                           "69317163492948606335995924319873"), (
                           24176176, 73745418, 52432133))
        digits = parse.(Int8, split(strip(inp), ""))
        @test AdventOfCode2019.Day16.part1!(digits) == out
    end
    for (inp, out) in zip(("03036732577212944063491565474664",
                           "02935109699940807407585447034323",
                           "03081770884921959731165446850517"), (
                           84462026, 78725270, 53553731))
        digits = parse.(Int8, split(strip(inp), ""))
        @test AdventOfCode2019.Day16.part2(digits) == out
    end
    @test day16() == [25131128, 53201602]
end

@testset "Day 17" begin
    @test day17() == [3936, 785733]
end

@testset "Day 18" begin
    inp = "#########\n" *
          "#b.A.@.a#\n" *
          "#########\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part1(maze) == 8

    inp = "########################\n" *
          "#f.D.E.e.C.b.A.@.a.B.c.#\n" *
          "######################.#\n" *
          "#d.....................#\n" *
          "########################\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part1(maze) == 86

    inp = "########################\n" *
          "#...............b.C.D.f#\n" *
          "#.######################\n" *
          "#.....@.a.B.c.d.A.e.F.g#\n" *
          "########################\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part1(maze) == 132

    inp = "#################\n" *
          "#i.G..c...e..H.p#\n" *
          "########.########\n" *
          "#j.A..b...f..D.o#\n" *
          "########@########\n" *
          "#k.E..a...g..B.n#\n" *
          "########.########\n" *
          "#l.F..d...h..C.m#\n" *
          "#################\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part1(maze) == 136

    inp = "########################\n" *
          "#@..............ac.GI.b#\n" *
          "###d#e#f################\n" *
          "###A#B#C################\n" *
          "###g#h#i################\n" *
          "########################\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part1(maze) == 81

    inp = "#######\n" *
          "#a.#Cd#\n" *
          "##1#2##\n" *
          "#######\n" *
          "##4#3##\n" *
          "#cB#Ab#\n" *
          "#######\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part2(maze) == 8

    inp = "###############\n" *
          "#d.ABC.#.....a#\n" *
          "######1#2######\n" *
          "###############\n" *
          "######3#4######\n" *
          "#b.....#.....c#\n" *
          "###############\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part2(maze) == 24

    inp = "#############\n" *
          "#DcBa.#.GhKl#\n" *
          "#.###1#2#I###\n" *
          "#e#d#####j#k#\n" *
          "###C#3#4###J#\n" *
          "#fEbA.#.FgHi#\n" *
          "#############\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part2(maze) == 32

    inp = "#############\n" *
          "#g#f.D#..h#l#\n" *
          "#F###e#E###.#\n" *
          "#dCba1#2BcIJ#\n" *
          "#############\n" *
          "#nK.L3#4G...#\n" *
          "#M###N#H###.#\n" *
          "#o#m..#i#jk.#\n" *
          "#############\n"
    maze = AdventOfCode2019.Day18.parse_input(inp)
    @test AdventOfCode2019.Day18.part2(maze) == 72

    @test day18() == [3962, 1844]
end

@testset "Day 19" begin
    @test day19() == [166, 3790981]
end

@testset "Day 20" begin
    inp = "         A          \n" *
          "         A          \n" *
          "  #######.######### \n" *
          "  #######.........# \n" *
          "  #######.#######.# \n" *
          "  #######.#######.# \n" *
          "  #######.#######.# \n" *
          "  #####  B    ###.# \n" *
          "BC...##  C    ###.# \n" *
          "  ##.##       ###.# \n" *
          "  ##...DE  F  ###.# \n" *
          "  #####    G  ###.# \n" *
          "  #########.#####.# \n" *
          "DE..#######...###.# \n" *
          "  #.#########.###.# \n" *
          "FG..#########.....# \n" *
          "  ###########.##### \n" *
          "             Z      \n" *
          "             Z      \n"
    @test day20(inp) == [23, 26]

    inp = "                   A               \n" *
          "                   A               \n" *
          "  #################.#############  \n" *
          "  #.#...#...................#.#.#  \n" *
          "  #.#.#.###.###.###.#########.#.#  \n" *
          "  #.#.#.......#...#.....#.#.#...#  \n" *
          "  #.#########.###.#####.#.#.###.#  \n" *
          "  #.............#.#.....#.......#  \n" *
          "  ###.###########.###.#####.#.#.#  \n" *
          "  #.....#        A   C    #.#.#.#  \n" *
          "  #######        S   P    #####.#  \n" *
          "  #.#...#                 #......VT\n" *
          "  #.#.#.#                 #.#####  \n" *
          "  #...#.#               YN....#.#  \n" *
          "  #.###.#                 #####.#  \n" *
          "DI....#.#                 #.....#  \n" *
          "  #####.#                 #.###.#  \n" *
          "ZZ......#               QG....#..AS\n" *
          "  ###.###                 #######  \n" *
          "JO..#.#.#                 #.....#  \n" *
          "  #.#.#.#                 ###.#.#  \n" *
          "  #...#..DI             BU....#..LF\n" *
          "  #####.#                 #.#####  \n" *
          "YN......#               VT..#....QG\n" *
          "  #.###.#                 #.###.#  \n" *
          "  #.#...#                 #.....#  \n" *
          "  ###.###    J L     J    #.#.###  \n" *
          "  #.....#    O F     P    #.#...#  \n" *
          "  #.###.#####.#.#####.#####.###.#  \n" *
          "  #...#.#.#...#.....#.....#.#...#  \n" *
          "  #.#####.###.###.#.#.#########.#  \n" *
          "  #...#.#.....#...#.#.#.#.....#.#  \n" *
          "  #.###.#####.###.###.#.#.#######  \n" *
          "  #.#.........#...#.............#  \n" *
          "  #########.###.###.#############  \n" *
          "           B   J   C               \n" *
          "           U   P   P               \n"
    @test AdventOfCode2019.Day20.solve(inp, (true, false)) == [58]

    inp = "             Z L X W       C                 \n" *
          "             Z P Q B       K                 \n" *
          "  ###########.#.#.#.#######.###############  \n" *
          "  #...#.......#.#.......#.#.......#.#.#...#  \n" *
          "  ###.#.#.#.#.#.#.#.###.#.#.#######.#.#.###  \n" *
          "  #.#...#.#.#...#.#.#...#...#...#.#.......#  \n" *
          "  #.###.#######.###.###.#.###.###.#.#######  \n" *
          "  #...#.......#.#...#...#.............#...#  \n" *
          "  #.#########.#######.#.#######.#######.###  \n" *
          "  #...#.#    F       R I       Z    #.#.#.#  \n" *
          "  #.###.#    D       E C       H    #.#.#.#  \n" *
          "  #.#...#                           #...#.#  \n" *
          "  #.###.#                           #.###.#  \n" *
          "  #.#....OA                       WB..#.#..ZH\n" *
          "  #.###.#                           #.#.#.#  \n" *
          "CJ......#                           #.....#  \n" *
          "  #######                           #######  \n" *
          "  #.#....CK                         #......IC\n" *
          "  #.###.#                           #.###.#  \n" *
          "  #.....#                           #...#.#  \n" *
          "  ###.###                           #.#.#.#  \n" *
          "XF....#.#                         RF..#.#.#  \n" *
          "  #####.#                           #######  \n" *
          "  #......CJ                       NM..#...#  \n" *
          "  ###.#.#                           #.###.#  \n" *
          "RE....#.#                           #......RF\n" *
          "  ###.###        X   X       L      #.#.#.#  \n" *
          "  #.....#        F   Q       P      #.#.#.#  \n" *
          "  ###.###########.###.#######.#########.###  \n" *
          "  #.....#...#.....#.......#...#.....#.#...#  \n" *
          "  #####.#.###.#######.#######.###.###.#.#.#  \n" *
          "  #.......#.......#.#.#.#.#...#...#...#.#.#  \n" *
          "  #####.###.#####.#.#.#.#.###.###.#.###.###  \n" *
          "  #.......#.....#.#...#...............#...#  \n" *
          "  #############.#.#.###.###################  \n" *
          "               A O F   N                     \n" *
          "               A A D   M                     \n"
    @test day20(inp) == [77, 396]

    @test day20() == [422, 5040]
end

@testset "Day 21" begin
    @test day21() == [19353074, 1147582556]
end

@testset "Day 22" begin
    @test day22() == [5540, 6821410630991]
end

@testset "Day 23" begin
    @test day23() == [17286, 11249]
end

@testset "Day 24" begin
    input = "....#\n" *
            "#..#.\n" *
            "#..##\n" *
            "..#..\n" *
            "#....\n"
    layout = AdventOfCode2019.Day24.parse_input(input)
    @test AdventOfCode2019.Day24.part1(layout) == 2129920
    @test AdventOfCode2019.Day24.part2(layout; niter = 10) == 99

    @test day24() == [18350099, 2037]
end

@testset "Day 25" begin
    @test day25() == 262848
end


text:
cd("dummy") do

  Julz.generate("function", "roo")

end

module Woof
  using Julz
  Julz.include_all_files("$(pwd())/dummy/src/functions", is_revised=false)
  Julz.@export_all_files "$(pwd())/dummy/src/functions"
end

@testset "Export All Files Macro Tests" begin

  @test isdefined(Julz, Symbol("@export_all_files")) == true

  @test isdefined(Woof, :roo) == true

  cd("dummy") do

    Julz.destroy("function", "roo")

  end

end


text:
using DocStringExtensions
using KernelAbstractions
using CLIMA.Mesh.Geometry

import CLIMA.Mesh.Grids:
    _ξ1x1, _ξ2x1, _ξ3x1, _ξ1x2, _ξ2x2, _ξ3x2, _ξ1x3, _ξ2x3, _ξ3x3

using CLIMA.VariableTemplates
import CLIMA.VariableTemplates.varsindex

include("diagnostic_fields_kernels.jl")
"""
    VecGrad{FT <: AbstractFloat,
         FTA2D <: AbstractArray{FT, 2},
         FTA3D <: AbstractArray{FT, 3},}

This data structure stores the gradients of a velocity field.

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    VecGrad(Npl, Nel, ::Type{FT}) where {FT <: AbstractFloat}

# Arguments for the innter constructor
 - 'Npl': Number of local dof in a spectral element
 - 'Nel": Number of spectral elements
 - 'FT': Floating point precision
"""
struct VecGrad{
    FT <: AbstractFloat,
    FTA2D <: AbstractArray{FT, 2},
    FTA3D <: AbstractArray{FT, 3},
}
    "Device array storing the gradient data"
    data::FTA3D
    "∂u₁∂x₁"
    ∂₁u₁::FTA2D
    "∂u₁∂x₂"
    ∂₂u₁::FTA2D
    "∂u₁∂x₃"
    ∂₃u₁::FTA2D
    "∂u₂∂x₁"
    ∂₁u₂::FTA2D
    "∂u₂∂x₂"
    ∂₂u₂::FTA2D
    "∂u₂∂x₃"
    ∂₃u₂::FTA2D
    "∂u₃∂x₁"
    ∂₁u₃::FTA2D
    "∂u₃∂x₂"
    ∂₂u₃::FTA2D
    "∂u₃∂x₃"
    ∂₃u₃::FTA2D

    function VecGrad(Npl, Nel, ::Type{FT}) where {FT <: AbstractFloat}
        DA = CLIMA.array_type()
        data = DA{FT}(undef, Npl, 9, Nel)

        ∂₁u₁, ∂₂u₁, ∂₃u₁ =
            view(data, :, 1, :), view(data, :, 2, :), view(data, :, 3, :)
        ∂₁u₂, ∂₂u₂, ∂₃u₂ =
            view(data, :, 4, :), view(data, :, 5, :), view(data, :, 6, :)
        ∂₁u₃, ∂₂u₃, ∂₃u₃ =
            view(data, :, 7, :), view(data, :, 8, :), view(data, :, 9, :)

        return new{FT, typeof(∂₁u₁), typeof(data)}(
            data,
            ∂₁u₁,
            ∂₂u₁,
            ∂₃u₁,
            ∂₁u₂,
            ∂₂u₂,
            ∂₃u₂,
            ∂₁u₃,
            ∂₂u₃,
            ∂₃u₃,
        )
    end
end
#--------------------------------------------------------------------------------------------------
"""
    Vorticity{FT <: AbstractFloat, FTA2D <: AbstractArray{FT, 2}, FTA3D <: AbstractArray{FT, 3}}

This data structure stores the vorticity of a velocity field.

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    Vorticity(Npl, Nel, ::Type{FT}) where {FT <: AbstractFloat}

# Arguments for the innter constructor
 - 'Npl': Number of local dof in a spectral element
 - 'Nel": Number of spectral elements
 - 'FT': Floating point precision
"""
struct Vorticity{
    FT <: AbstractFloat,
    FTA2D <: AbstractArray{FT, 2},
    FTA3D <: AbstractArray{FT, 3},
}
    "Device array storing the vorticity data"
    data::FTA3D
    "x1 component of vorticity"
    Ω₁::FTA2D
    "x2 component of vorticity"
    Ω₂::FTA2D
    "x3 component of vorticity"
    Ω₃::FTA2D
    function Vorticity(Npl, Nel, ::Type{FT}) where {FT <: AbstractFloat}
        DA = CLIMA.array_type()
        data = DA{FT}(undef, Npl, 3, Nel)
        Ω₁ = view(data, :, 1, :)
        Ω₂ = view(data, :, 2, :)
        Ω₃ = view(data, :, 3, :)
        return new{FT, typeof(Ω₁), typeof(data)}(data, Ω₁, Ω₂, Ω₃)
    end
end
#--------------------------------------------------------------------------------------------------
"""
    compute_vec_grad(model::BalanceLaw,
                         Q::AbstractArray{FT},
                        dg::DGModel,
                     vgrad::VecGrad{FT},) where {FT <: AbstractFloat}

This function computes the spatial gradients of the velocity field.

# Arguments
 - `model`: BalanceLaw
 - `Q`: State array
 - `dg`: DGmodel
 - `vgrad`: Vector gradients
"""
function compute_vec_grad(
    model::BalanceLaw,
    Q::AbstractArray{FT},
    dg::DGModel,
) where {FT <: AbstractFloat}
    Nel = size(Q.realdata, 3) # # of spectral elements
    Npl = size(Q.realdata, 1) # # of dof per element
    qm1 = size(dg.grid.D, 2)  # poly order + 1
    DA = CLIMA.array_type()

    vgrad = VecGrad(Npl, Nel, FT)
    ind = [
        varsindex(vars_state(model, FT), :ρ)
        varsindex(vars_state(model, FT), :ρu)
    ]
    _ρ, _ρu, _ρv, _ρw = ind[1], ind[2], ind[3], ind[4]

    vgrad_data = vgrad.data
    vgeo = dg.grid.vgeo
    D = dg.grid.D
    sv = Q.data

    device = DA <: Array ? CPU() : CUDA()
    g = DA{FT}(undef, Npl, Nel, 3, 3)

    workgroup = (qm1, qm1)
    ndrange = (Nel * qm1, qm1)

    kernel = compute_vec_grad_kernel!(device, workgroup)
    event = kernel(
        sv,
        D,
        vgeo,
        g,
        vgrad_data,
        _ρ,
        _ρu,
        _ρv,
        _ρw,
        Val(qm1),
        ndrange = ndrange,
    )
    wait(event)

    return vgrad
end
#--------------------------------------------------------------------------------------------------
"""
    compute_vorticity(dg::DGModel,
                   vgrad::VecGrad{FT},
                    vort::Vorticity{FT},) where {FT <: AbstractFloat}

This function computes the spatial gradients of the velocity field.

# Arguments
 - `dg`: DGmodel
 - `vgrad`: Velocity gradients
 - `vort`: Vorticity
"""
function compute_vorticity(
    dg::DGModel,
    vgrad::VecGrad{FT},
) where {FT <: AbstractFloat}
    Npl = size(vgrad.∂₁u₁, 1)
    Nel = size(vgrad.∂₁u₁, 2)
    qm1 = size(dg.grid.D, 2)  # poly order + 1
    DA = CLIMA.array_type()
    device = DA <: Array ? CPU() : CUDA()

    vort = Vorticity(Npl, Nel, FT)

    vgrad_data = vgrad.data
    vort_data = vort.data

    Ω₁, Ω₂, Ω₃ = vort.Ω₁, vort.Ω₂, vort.Ω₃

    workgroup = (qm1, qm1)
    ndrange = (Nel * qm1, qm1)

    kernel = compute_vorticity_kernel!(device, workgroup)
    event = kernel(vgrad_data, vort_data, Val(qm1), ndrange = ndrange)
    wait(event)

    return vort
end
#--------------------------------------------------------------------------------------------------

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
module Images

using FileType

TypeJpeg     = FileType.Types("jpg", MIME("image/jpeg"))
TypeJpeg2000 = FileType.Types("jp2", MIME("image/jp2"))
TypePng      = FileType.Types("png", MIME("image/png"))
TypeGif      = FileType.Types("gif", MIME("image/gif"))
TypeWebp     = FileType.Types("webp", MIME("image/webp"))
TypeCR2      = FileType.Types("cr2", MIME("image/x-canon-cr2"))
TypeTiff     = FileType.Types("tif", MIME("image/tiff"))
TypeBmp      = FileType.Types("bmp", MIME("image/bmp"))
TypeJxr      = FileType.Types("jxr", MIME("image/vnd.ms-photo"))
TypePsd      = FileType.Types("psd", MIME("image/vnd.adobe.photoshop"))
TypeIco      = FileType.Types("ico", MIME("image/vnd.microsoft.icon"))
TypeDwg      = FileType.Types("dwg", MIME("image/vnd.dwg"))

function Jpeg(buf)::Bool
	return length(buf) > 2 &&
		buf[1] == 0xFF &&
		buf[2] == 0xD8 &&
		buf[3] == 0xFF
end

function Jpeg2000(buf)::Bool
	return length(buf) > 12 &&
		buf[1] == 0x0 &&
		buf[2] == 0x0 &&
		buf[3] == 0x0 &&
		buf[4] == 0xC &&
		buf[5] == 0x6A &&
		buf[6] == 0x50 &&
		buf[7] == 0x20 &&
		buf[8] == 0x20 &&
		buf[9] == 0xD &&
		buf[10] == 0xA &&
		buf[11] == 0x87 &&
		buf[12] == 0xA &&
		buf[13] == 0x0
end

function Png(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x89 && buf[2] == 0x50 &&
		buf[3] == 0x4E && buf[4] == 0x47
end

function Gif(buf)::Bool
	return length(buf) > 2 &&
		buf[1] == 0x47 && buf[2] == 0x49 && buf[3] == 0x46
end

function Webp(buf)::Bool
	return length(buf) > 11 &&
		buf[9] == 0x57 && buf[10] == 0x45 &&
		buf[11] == 0x42 && buf[12] == 0x50
end

function CR2(buf)::Bool
	return length(buf) > 10 &&
		((buf[1] == 0x49 && buf[2] == 0x49 && buf[3] == 0x2A && buf[4] == 0x0) || # Little Endian
			(buf[1] == 0x4D && buf[2] == 0x4D && buf[3] == 0x0 && buf[4] == 0x2A)) && # Big Endian
		buf[9] == 0x43 && buf[10] == 0x52 && # CR2 magic word
		buf[11] == 0x02 # CR2 major version
end

function Tiff(buf)::Bool
	return length(buf) > 10 &&
		((buf[1] == 0x49 && buf[2] == 0x49 && buf[3] == 0x2A && buf[4] == 0x0) || # Little Endian
			(buf[1] == 0x4D && buf[2] == 0x4D && buf[3] == 0x0 && buf[4] == 0x2A)) && # Big Endian
		!CR2(buf) # To avoid conflicts differentiate Tiff from CR2
end

function Bmp(buf)::Bool
	return length(buf) > 1 &&
		buf[1] == 0x42 &&
		buf[2] == 0x4D
end

function Jxr(buf)::Bool
	return length(buf) > 2 &&
		buf[1] == 0x49 &&
		buf[2] == 0x49 &&
		buf[3] == 0xBC
end

function Psd(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x38 && buf[2] == 0x42 &&
		buf[3] == 0x50 && buf[4] == 0x53
end

function Ico(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x00 && buf[2] == 0x00 &&
		buf[3] == 0x01 && buf[4] == 0x00
end

function Dwg(buf)::Bool
	return length(buf) > 3 &&
		buf[1] == 0x41 && buf[2] == 0x43 &&
		buf[3] == 0x31 && buf[4] == 0x30
end

Image = Dict(
	TypeJpeg     => Jpeg,
	TypeJpeg2000 => Jpeg2000,
	TypePng      => Png,
	TypeGif      => Gif,
	TypeWebp     => Webp,
	TypeCR2      => CR2,
	TypeTiff     => Tiff,
	TypeBmp      => Bmp,
	TypeJxr      => Jxr,
	TypePsd      => Psd,
	TypeIco      => Ico,
	TypeDwg      => Dwg,
)

end


text:
# coding: utf-8

# In[1]:

#fractal is one of the interesting topics in geometry
#it is usually described by a recursive function
#voila,here we are!
using Plots

# In[2]:

#create rectangle shape
rectangle(w,h,x,y)=Shape(x.+[0,w,w,0],y.-[0,0,h,h])

#compute euclidean distance
function euclidean_distance(point1,point2)
    return √((point1[1]-point2[1])^2+(point1[2]-point2[2])^2)
end

# In[3]:

#recursively plot jerusalem cross
#it kinda looks like flag of georgia (n=2)
#i mean the eurasian country not a yankee state
#i call it jerusalem cross but it is aka cross menger square,jerusalem square
#it is a 2d version of jerusalem cube
#a good reference to jerusalem cube
# https://robertdickau.com/jerusalemcube.html
#a good understanding of sierpiński carpet is helpful as well
# https://github.com/je-suis-tm/recursion-and-dynamic-programming/blob/master/sierpi%C5%84ski%20carpet.py
#do not confuse it with quadratic cross,which creates new crosses from the tips
# https://onlinemathtools.com/generate-quadratic-cross-fractal
#or fibonacci snowflakes,which is more like koch snowflake
# http://www.slabbe.org/Publications/2011-fibo-snowflakes.pdf
#or vicsek fractal,which is more similar to crosslet cross
# https://en.wikipedia.org/wiki/Vicsek_fractal
function jerusalem_cross(top_left,top_right,bottom_left,bottom_right,n)

    if n<=0
        return
    else

        #compute the width
        width=euclidean_distance(top_left,top_right)

        #create the cross
        plot!(rectangle(width*(√(2)-1)^2,
                width*(1-2*((√(2)-1)^2)),
                top_left[1]+width*(√(2)-1),
                top_left[2]-width*(√(2)-1)^2),
                color="black")
        plot!(rectangle(width*(1-2*((√(2)-1)^2)),
                width*(√(2)-1)^2,
                top_left[1]+width*(√(2)-1)^2,
                top_left[2]-width*(√(2)-1)),
                color="black")

        #top left corner recursion
        jerusalem_cross(top_left,(top_left[1]+width*(√(2)-1),top_left[2]),
                        (top_left[1],top_left[2]-width*(√(2)-1)),
                        (top_left[1]+width*(√(2)-1),
                        top_left[2]-width*(√(2)-1)),n-1)

        #top right corner recursion
        jerusalem_cross((top_right[1]-width*(√(2)-1),top_left[2]),top_right,
                        (top_right[1]-width*(√(2)-1),
                         top_left[2]-width*(√(2)-1)),
                        (top_right[1],top_left[2]-width*(√(2)-1)),n-1)

        #bottom left corner recursion
        jerusalem_cross((bottom_left[1],bottom_left[2]+width*(√(2)-1)),
                        (bottom_left[1]+width*(√(2)-1),
                         bottom_left[2]+width*(√(2)-1)),
                        bottom_left,
                        (bottom_left[1]+width*(√(2)-1),bottom_left[2]),n-1)

        #bottom right corner recursion
        jerusalem_cross((bottom_right[1]-width*(√(2)-1),
                         bottom_right[2]+width*(√(2)-1)),
                        (bottom_right[1],
                         bottom_right[2]+width*(√(2)-1)),
                        (bottom_right[1]-width*(√(2)-1),
                         bottom_right[2]),
                        bottom_right,n-1)

        #top mid corner recursion
        jerusalem_cross((top_left[1]+width*(√(2)-1),top_left[2]),
                        (top_right[1]-width*(√(2)-1),top_left[2]),
                        (top_left[1]+width*(√(2)-1),
                         top_left[2]-width*(√(2)-1)^2),
                        (top_right[1]-width*(√(2)-1),
                         top_left[2]-width*(√(2)-1)^2),n-2)

        #bottom mid corner recursion
        jerusalem_cross((bottom_left[1]+width*(√(2)-1),
                         bottom_left[2]+width*(√(2)-1)^2),
                        (bottom_right[1]-width*(√(2)-1),
                         bottom_left[2]+width*(√(2)-1)^2),
                        (bottom_left[1]+width*(√(2)-1),
                         bottom_left[2]),
                        (bottom_right[1]-width*(√(2)-1),
                         bottom_left[2]),n-2)

        #left mid corner recursion
        jerusalem_cross((bottom_left[1],
                         top_left[2]-width*(√(2)-1)),
                        (bottom_left[1]+width*(√(2)-1)^2,
                         top_left[2]-width*(√(2)-1)),
                        (bottom_left[1],bottom_left[2]+width*(√(2)-1)),
                        (bottom_left[1]+width*(√(2)-1)^2,
                         bottom_left[2]+width*(√(2)-1)),n-2)

        #right mid corner recursion
        jerusalem_cross((bottom_right[1]-width*(√(2)-1)^2,
                         top_right[2]-width*(√(2)-1)),
                        (bottom_right[1],
                         top_right[2]-width*(√(2)-1)),
                        (bottom_right[1]-width*(√(2)-1)^2,
                         bottom_right[2]+width*(√(2)-1)),
                        (bottom_right[1],bottom_right[2]+width*(√(2)-1)),
                        n-2)

    end

end

# In[4]:

#initialize
top_left=(0,0)
top_right=(1,0)
bottom_left=(0,-1)
bottom_right=(1,-1)
n=5;

# In[5]:

#viz
gr(size=(500,500))
fig=plot(legend=false,grid=false,axis=false,ticks=false,

    xlim=(top_left[1],top_right[1]),
    ylim=(bottom_right[2],top_right[2]))
jerusalem_cross(top_left,top_right,bottom_left,bottom_right,n)
fig

max_stars_repo_path: jerusalem cross.jl

max_issues_repo_path: jerusalem cross.jl

max_forks_repo_path: jerusalem cross.jl


text:
#=
max_heap:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-25
=#

module MyMaxHeap

    # heap list
    heap_list = [0]
    current_size = 0
    input_list = [9, 21, 5, 14, 33, 19, 17, 11, 27, 18]

    function min_child_idx(idx)
        if (idx * 2 + 1) > current_size
            return idx * 2
        else
            if heap_list[idx * 2] < heap_list[idx * 2 + 1]
                return idx * 2 + 1
            else
                return idx * 2
            end
        end
    end

    # method to swap the root with its biggest
    # child bigger than the root. For percolating
    # a node down the tree.
    function perc_down(idx)
        while (idx * 2) <= current_size
            mc_idx = min_child_idx(idx)
            if heap_list[idx] < heap_list[mc_idx]
                tmp = heap_list[idx]
                heap_list[idx] = heap_list[mc_idx]
                heap_list[mc_idx] = tmp
            end
            idx = mc_idx
        end
    end

    # method to make as it moves the node
    # in an input list into their proper position
    function build_heap(list)
        println("Building a Heap from List $(list)")
        i = length(list) ÷ 2
        global current_size = length(list)
        global heap_list = list
        while i > 0
            perc_down(i)
            println("i = $(i): $(heap_list)")
            i -= 1
        end
    end

    # method to percolate a new item
    # as far up in the tree as it needs to
    # go to maintain the heap property
    function perc_up(idx)
        while (idx ÷ 2) > 0
            if heap_list[idx] > heap_list[idx ÷ 2]
                tmp = heap_list[idx ÷ 2]
                heap_list[idx ÷ 2] = heap_list[idx]
                heap_list[idx] = tmp
            end
            println("i = $(idx): $(heap_list)")
            idx = idx ÷ 2
        end
    end

    # method to append an item to the end of list
    # regain the heap structure property by
    # comparing the newly added item with its parent
    function insert(element)
        println("Insert $(element) to heap")
        push!(heap_list, element)
        println(heap_list)
        global current_size += 1
        perc_up(current_size)
    end

    # method to remove maximum item from root of heap
    # after removed, restore full compliance with
    # heap structure and heap order policy
    function remove_max()
        max_item = heap_list[1]
        println("Remove max item $(max_item) from heap")
        heap_list[1] = heap_list[current_size]
        global current_size -= 1
        pop!(heap_list)
        perc_down(1)
        println(heap_list)
    end

    function main()
        build_heap(input_list)
        println("")
        insert(50)
        println("")
        remove_max()
    end

end

if abspath(PROGRAM_FILE) == @__FILE__
    using .MyMaxHeap
    MyMaxHeap.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Agenie
using Test

@testset "Agenie.jl" begin
    # Write your tests here.
end

max_stars_repo_name: aminya/Agenie

max_issues_repo_name: aminya/Agenie

max_forks_repo_name: aminya/Agenie


text:
export percept_lengths

function findlengths(x)
  indices = [1; findall(diff(x) .!= 0) .+ 1; length(x)+1]
  filter(!iszero,diff(indices)), x[indices[1:end-1]]
end

mergelengths(threshold) = lens -> mergelengths(lens,threshold)
function mergelengths(lens,vals,threshold)
  mergedlen = similar(lens)
  mergedval = similar(vals)

  i = 0
  for (len,val) in zip(lens,vals)
    if i > 0 && (len < threshold || mergedval[i] == val)
      mergedlen[i] += len
    else
      i += 1
      mergedlen[i] = len
      mergedval[i] = val
    end
  end

  mergedlen[1:i], mergedval[1:i]
end

function percept_lengths(counts; threshold = 0.45,
                         min_length=750ms)
  lens,vals = findlengths(Array(counts .< threshold))
  slens = lens * ustrip(uconvert(s,Δt(counts)))

  mergelengths(slens,vals,ustrip(asseconds(min_length)))
end

function percept_lengths(spmask::AbstractMatrix, sp::AbstractMatrix,
                         settings)
  settings = read_settings(settings)
  startHz, stopHz = asHz.(settings.rates.freqbound.freq_limits)
  ratio, = bandwidth_ratio(spmask, sp[:,startHz .. stopHz], settings)
  percept_lengths(ratio, settings)
end

function percept_lengths(mask,params,settings)
  stim = audiospect_stimulus(params,settings,cache=true)
  percept_lengths(mask,stim,settings)
end

function percept_lengths(result::NamedTuple,settings)
  settings = read_settings(settings)
  startHz, stopHz = asHz.(settings.rates.freqbound.freq_limits)
  ratio, = bandwidth_ratio(result.primary_source,
                           result.input[:,startHz .. stopHz],
                           settings)
  percept_lengths(ratio,settings)
end

percept_lengths(ratio::AbstractVector,settings) =
  percept_lengths(ratio;settings.percept_lengths...)

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DensityWithShape{D<:AbstractDensity,S<:AbstractValueShape} <: AbstractDensity
    density::D
    shape::S
end

Base.parent(density::DensityWithShape) = density.density

var_bounds(density::DensityWithShape) = var_bounds(parent(density))

ValueShapes.varshape(density::DensityWithShape) = density.shape

function eval_logval_unchecked(density::DensityWithShape, v::Any)
    eval_logval_unchecked(parent(density), v)
end

struct ReshapedDensity{D<:AbstractDensity,S<:AbstractValueShape} <: AbstractDensity
    density::D
    shape::S
end

Base.parent(density::ReshapedDensity) = density.density

var_bounds(density::ReshapedDensity) = var_bounds(density.density)

ValueShapes.varshape(density::ReshapedDensity) = density.shape

function eval_logval(density::ReshapedDensity, v::Any, T::Type{<:Real})
    v_shaped = fixup_variate(varshape(density), v)
    orig_density = parent(density)
    orig_varshape = varshape(orig_density)
    v_reshaped = reshape_variate(orig_varshape, v_shaped)
    eval_logval(orig_density, v_reshaped, T)
end

ValueShapes.unshaped(density::AbstractDensity) = _unshaped_density(density, varshape(density))

_unshaped_density(density::AbstractDensity, ::ArrayShape{<:Real,1}) = density
_unshaped_density(density::AbstractDensity, ::AbstractValueShape) = ReshapedDensity(density, ArrayShape{Real}(totalndof(density)))

(shape::AbstractValueShape)(density::AbstractDensity) = _reshaped_density(density, shape, varshape(density))

function _reshaped_density(density::AbstractDensity, new_shape::VS, orig_shape::VS) where {VS<:AbstractValueShape}
    if orig_shape == new_shape
        density
    else
        throw(ArgumentError("Value shapes are incompatible"))
    end
end

function _reshaped_density(density::AbstractDensity, new_shape::AbstractValueShape, orig_shape::AbstractValueShape)
    if totalndof(orig_shape) == totalndof(new_shape)
        ReshapedDensity(density, new_shape)
    else
        throw(ArgumentError("Value shapes are incompatible"))
    end
end


text:
type TensorDesc
    ptr::Ptr{Void}

    function TensorDesc{T,N}(x::CuArray{T,N})
        csize = Cint[size(x,i) for i=N:-1:1]
        cstrides = Cint[stride(x,i) for i=N:-1:1]
        p = Ptr{Void}[0]
        cudnnCreateTensorDescriptor(p)
        cudnnSetTensorNdDescriptor(p[1], datatype(T), N, csize, cstrides)
        desc = new(p[1])
        finalizer(desc, cudnnDestroyTensorDescriptor)
        desc
    end
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::TensorDesc) = desc.ptr


text:
#####
##### 3D modes on tetrahedra
#####

"""
    simplex_3D(a, b, c, i, j, k)

Evaluate 3D "Legendre" basis phi_ijk at (a,b,c) coordinates on the [-1,1] cube
"""
function simplex_3D(a, b, c, i, j, k)
    h1 = jacobiP(a,0,0,i)
    h2 = jacobiP(b,2*i+1,0,j)
    h3 = jacobiP(c,2*(i+j)+2,0,k)
    return @. 2*sqrt(2)*h1*h2*((1-b)^i)*h3*((1-c)^(i+j))
end

"""
    grad_simplex_3D(a, b, c, id, jd, kd)

Evalute the partial derivative w.r.t. (r,s,t) of 3D "Legendre" polynomial with
index, or order, (id,jd,kd) at (a,b,c)
"""

function grad_simplex_3D(a, b, c, id, jd, kd)
    fa = jacobiP(a,0,0,id)
    gb = jacobiP(b,2*id+1,0,jd)
    hc = jacobiP(c,2*(id+jd)+2,0,kd)

    dfa = grad_jacobiP(a,0,0,id)
    dgb = grad_jacobiP(b,2*id+1,0,jd)
    dhc = grad_jacobiP(c,2*(id+jd)+2,0,kd)

    # r-derivative
    V3Dr = @. dfa.*(gb.*hc)
    if id>0
        V3Dr = @. V3Dr*((0.5*(1-b))^(id-1))
    end
    if id+jd>0
        V3Dr = @. V3Dr*((0.5*(1-c))^(id+jd-1))
    end

    # s-derivative
    V3Ds = @. 0.5*(1+a)*V3Dr
    tmp = @. dgb*((0.5*(1-b))^id)
    if id>0
        tmp = @. tmp+(-0.5*id)*(gb*(0.5*(1-b))^(id-1))
    end
    if id+jd>0
        tmp = @. tmp*((0.5*(1-c))^(id+jd-1))
    end
    tmp = @. fa*(tmp*hc)
    V3Ds = @. V3Ds+tmp

    # t-derivative
    V3Dt = @. 0.5*(1+a).*V3Dr+0.5*(1+b)*tmp
    tmp = @. dhc*((0.5*(1-c))^(id+jd))
    if id+jd>0
      tmp = @. tmp-0.5*(id+jd)*(hc*((0.5*(1-c))^(id+jd-1)))
    end
    tmp = @. fa*(gb*tmp)
    tmp = @. tmp*((0.5*(1-b))^id)
    V3Dt = @. V3Dt+tmp

    # normalize
    V3Dr = @. V3Dr*(2^(2*id+jd+1.5))
    V3Ds = @. V3Ds*(2^(2*id+jd+1.5))
    V3Dt = @. V3Dt*(2^(2*id+jd+1.5))
    return V3Dr,V3Ds,V3Dt
end

function equi_nodes(elem::Tet,N)
    Np = (N+1)*(N+2)*(N+3) ÷ 6
    r1D = LinRange(-1,1,N+1)
    r,s,t = ntuple(x->zeros(Np),3)
    sk = 1
    #for i = 0:N, j = 0:N-i, k = 0:N-i-j
    for k = 0:N, j = 0:N-k, i = 0:N-j-k
        r[sk] = r1D[i+1]
        s[sk] = r1D[j+1]
        t[sk] = r1D[k+1]
        sk += 1
    end

    return r,s,t
end

"""
    quad_nodes_tet(N)

Returns quadrature nodes and weights which exactly integrate degree N polynomials
"""
function quad_nodes_tet(N)

    if N<16
        rstw::Matrix{Float64} = readdlm(string(@__DIR__,"/QuadratureData/quad_nodes_tet_N", N, ".txt"),' ', Float64, '\n')
        r = rstw[:,1]
        s = rstw[:,2]
        t = rstw[:,3]
        w = rstw[:,4]
    else
        cubN = convert(Int,ceil((N+1)/2))
        r,s,t,w = stroud_quad_nodes(Tet(),cubN)
    end

    return r,s,t,w
end

function stroud_quad_nodes(elem::Tet,N)
    cubA,cubWA = gauss_quad(0,0,N)
    cubB,cubWB = gauss_quad(1,0,N)
    cubC,cubWC = gauss_quad(2,0,N)

    a,b,c = vec.(meshgrid(cubA,cubB,cubC))
    wa,wb,wc = vec.(meshgrid(cubWA,cubWB,cubWC))

    r = @. .5*(1+a)*.5*(1-b)*(1-c)-1
    s = @. .5*(1+b)*(1-c)-1
    t = c
    w = @. wa*wb*wc
    w = (4/3)*w./sum(w)
    return r,s,t,w
end

quad_nodes(elem::Tet,N) = quad_nodes_tet(2*N)

function basis(elem::Tet,N,r,s,t)
    Np = (N+1)*(N+2)*(N+3)÷6

    V,Vr,Vs,Vt = ntuple(x->zeros(length(r),Np),4)

    a = @. 2*(1+r)/(-s-t)-1
    b = @. 2*(1+s)/(1-t)-1
    c = t

    tol = 1e-14
    ida = @. abs(s+t) < tol
    idb = @. abs(1-t) < tol
    a[ida] .= -1
    b[idb] .= -1

    sk = 1
    for i = 0:N
        for j = 0:N-i
            for k = 0:N-i-j
                V[:, sk] .= simplex_3D(a, b, c, i, j, k)
                # should have fewer allocations than previous implementation:
                # `Vr[:,sk], Vs[:,sk], Vt[:,sk] = grad_simplex(...)`
                map((out, x) -> out .= x, view.((Vr, Vs, Vt), :, sk), grad_simplex_3D(a, b, c, i, j, k))
                sk += 1
            end
        end
    end
    return V, Vr, Vs, Vt
end


text:
using Compat
import BinDeps: download_cmd, unpack_cmd, splittarpath
if is_apple()
    using Homebrew
end

function download_and_unpack(baseurl, filename)
    downloadurl = baseurl * filename
    info("Downloading $filename from $downloadurl\nTo avoid this " *
        "download, install git manually and add it to your path " *
        "before\nrunning Pkg.add(\"Git\") or Pkg.build(\"Git\")")
    dest = joinpath("usr", string(Sys.MACHINE))
    for dir in ("downloads", "usr", dest)
        isdir(dir) || mkdir(dir)
    end
    filename = joinpath("downloads", filename)
    isfile(filename) || run(download_cmd(downloadurl, filename))
    # TODO: checksum validation
    (b, ext, sec_ext) = splittarpath(filename)
    run(unpack_cmd(filename, dest, is_windows() ? ".7z" : ext, sec_ext))
    # TODO: make this less noisy on windows, see how WinRPM does it
end

gitcmd = `git`
gitver = "notfound"
try
    gitver = readchomp(`$gitcmd --version`)
end
if gitver == "notfound"
    if is_apple()
        # we could allow other options, but lots of other packages already
        # depend on Homebrew.jl on mac and it needs a working git to function
        error("Working git not found on path, try running\nPkg.build(\"Homebrew\")")
    end
    baseurl = ""
    if is_linux() && (Sys.ARCH in (:x86_64, :i686, :i586, :i486, :i386))
        # use conda for a non-root option on x86/amd64 linux
        # TODO? use conda-forge when we no longer build julia on centos 5
        gitver = "2.6.4"
        plat = "download/linux-$(Sys.WORD_SIZE)/"
        baseurl = "http://anaconda.org/anaconda/git/$gitver/$plat"
        download_and_unpack(baseurl, "git-$gitver-0.tar.bz2")
        # dependencies for linux
        sslver = "1.0.2h"
        sslbase = "http://anaconda.org/anaconda/openssl/$sslver/$plat"
        download_and_unpack(sslbase, "openssl-$sslver-1.tar.bz2")
        zlibver = "1.2.8"
        zlibbase = "http://anaconda.org/anaconda/zlib/$zlibver/$plat"
        download_and_unpack(zlibbase, "zlib-$zlibver-3.tar.bz2")
    elseif is_windows()
        # download and extract portablegit
        gitver = "2.9.0"
        baseurl = "https://github.com/git-for-windows/git/releases/download/"
        download_and_unpack(baseurl * "v$gitver.windows.1/",
            "PortableGit-$gitver-$(Sys.WORD_SIZE)-bit.7z.exe")
    end
    if !isempty(baseurl)
        gitpath = joinpath(dirname(@__FILE__), "usr", string(Sys.MACHINE), "bin", "git")
        gitcmd = `$gitpath`
    end
    try
        gitver = readchomp(`$gitcmd --version`)
        info("Successfully installed $gitver to $gitcmd")
        # TODO: fix a warning about missing /templates here on linux
        # by setting an environment variable in deps.jl
    catch err
        error("Could not automatically install git, error was: $err\n" *
            (is_windows() ? "Report an issue at https://github.com/JuliaPackaging/Git.jl/issues/new" :
            "Try installing git via your system package manager then running\nPkg.build(\"Git\")"))
    end
else
    try
        # this is in a try because some environments like centos 7
        # docker containers don't have `which` installed by default
        gitpath = chomp(readlines(is_windows() ? `where git` : `which git`)[1])
        gitcmd = `$gitpath`
    end
    info("Using $gitver found on path" * (gitcmd == `git` ?
        "" : " at $gitcmd"))
end
open("deps.jl", "w") do f
    println(f, "gitcmd = $gitcmd")
end


text:
###
### Conversion & Promotion
###
###
### Conversion methods for biological sequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

function Base.convert(::Type{S}, seq::BioSequence) where {S<:AbstractString}
    return S([Char(x) for x in seq])
end
Base.String(seq::BioSequence) = convert(String, seq)

Base.convert(::Type{Vector{DNA}}, seq::BioSequence{<:DNAAlphabet}) = collect(seq)
Base.convert(::Type{Vector{RNA}}, seq::BioSequence{<:RNAAlphabet}) = collect(seq)
Base.convert(::Type{Vector}, seq::BioSequence) = collect(seq)
Base.convert(::Type{Vector{AminoAcid}}, seq::AminoAcidSeq) = collect(seq)


text:
using Base.Test
using MT1DGeneticInversion

include("./utils.jl")
include("./calc_response.jl")


text:
using Test
using WaveProp.Utils
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Mesh
using WaveProp.Integration
using QuadGK
using LinearAlgebra
using Plots

## test
pde   = Laplace(;dim=2)
h     = 0.1
τ     = LagrangeLine((0,0),(h,0))
plot(τ)
Γᶜ    = Nystrom.auxiliary_domain_ldim(τ,0.1)
[plot!(Γ) for Γ in Γᶜ]
plot!()
# plot!(size = (2000,2000))

k     = SingleLayerKernel(pde)
n     = 10
qrule  = GaussLegendre(n)
qinner = GaussLegendre(n)
s     = 0.53
x     = τ(s)

# interpolation surface
# xi,wi,νi = push_forward_quad_with_normal(τ,qrule)
# # auxiliary integration surface
# xq,wq,νq = _auxiliary_quadrature(τ,qinner)
# scatter(xi)
# scatter!(xq)
# scatter!([x],m=:x)

g = (y) -> cos(y[1])*exp(y[2])

# quadgk reference
f = (u) -> begin
    y = τ(u)
    jac = jacobian(τ,u)
    μ = transpose(jac)*jac |> det |> sqrt
    k(x,y)*g(y)*μ
end
I = quadgk(f,0,1)[1]

# naive intgration
x̂,ŵ   = qrule()
xq,wq   = push_forward_quad(τ,qrule)
I1 = integrate(y->k(x,y)*g(y),xq,wq)
@show (I - I1)/I |> abs

# modified integration
w̃ = singular_weights(k,τ,x,qrule,qinner)
I2 = integrate(y->g(y),xq,w̃)
@show (I - I2)/I |> abs

# cov = IMT{10,1}()
# x,w,ν = Nystrom._auxiliary_quadrature(τ,qrule,cov)
# sum(w)

# d⃗    = -20*h/n*normal(τ,0.5)
# Ω    = extrude(τ,d⃗)
# Γ    = boundary(Ω)

# Σ    = Geometry.circle(;radius=10*h,center=SVector(h/2,0))

# x̂, ŵ = qrule()
# mesh = GenericQuadrature{2,Float64}()
# x̃,w̃ = Integration._push_forward_quad(cov,x̂,ŵ)
# for i in 1:4
#     if i == 2
#         x,w,n⃗ = Integration._push_forward_quad_with_normal(Γ[i],x̃,w̃)
#     elseif i == 4
#         xr = -x̃ + 1 |> reverse
#         wr = reverse(w̃)
#         x,w,n⃗ = Integration._push_forward_quad_with_normal(Γ[i],xr,wr)
#     else
#         x,w,n⃗ = Integration._push_forward_quad_with_normal(Γ[i],x̂,ŵ)
#     end
#     append!(mesh.nodes,x)
#     append!(mesh.weights,w)
#     append!(mesh.normals,n⃗)
# end

# plot(mesh.nodes,m=:x,label="")

# # mesh = GenericQuadrature{2,Float64}(Γ,qrule)

# ##
# # plot(mesh.nodes)
# # side = Γ[2]
# # x,w = quadgen(side,qrule)
# # cov = IMT()
# # x̂,ŵ = qrule()
# # x,w = push_forward_map(cov,x̂,ŵ)
# # sum(w)
# # scatter(x,zero(x))
# # scatter!(x)
# ##

# # tmp  = Geometry.circle(;radius=h,center=SVector(0,0))
# # tmp2  = Geometry.circle(;radius=h,center=SVector(3*h,0))
# # mesh = GenericQuadrature{2,Float64}((tmp,tmp2),qrule)

# xs = SVector(10,10)

# G = SingleLayerKernel(pde)
# dG = DoubleLayerKernel(pde)
# γ₀Φ = [G(xs,y) for y in qnodes(mesh)]
# γ₁Φ = [dG(xs,y,ny) for (y,ny) in zip(qnodes(mesh),qnormals(mesh))]

# S = SingleLayerOperator(pde,mesh)
# D = DoubleLayerOperator(pde,mesh)

# ee = D*γ₀Φ - S*γ₁Φ + γ₀Φ/2
# @info norm(ee[1:n],Inf)

# xx = [pt[1] for pt in mesh.nodes[1:n]]
# plot(xx,ee[1:n],m=:x)

# ee[1:n]

# # ee = Utils.error_interior_green_identity(S,D,γ₀Φ,γ₁Φ) / norm(γ₀Φ,Inf)
# # norm(ee,Inf)

# # @testset "Kernels" begin
# #     pde  = Laplace(;dim=2)
# #     h    = 0.01
# #     k    = SingleLayerKernel(pde)
# #     Γ    = ParametricLine() do u
# #         SVector(h*u,0.)
# #     end
# #     N = 1
# #     qrule    = GaussLegendre{2*N}()
# #     x̂, ŵ     = qrule()
# #     xi,wi,ni    = quadgen_with_normal(Γ,qrule)
# #     l        = sum(wi) # approximate length of Γ
# #     d⃗        = 5*l*normal(Γ,0.5)
# #     bnd     = extrude(Γ,d⃗)
# #     els     = bnd[2:end]
# #     xq,wq,nq = quadgen_with_normal(els,qrule)
# #     xs = xi[N]
# #     w  = singular_quadrature(k,Γ,xs,qrule)
# #     f  = u-> k(xs,Γ(u))*sqrt(det(transpose(jacobian(Γ,u))*jacobian(Γ,u)))
# #     ff = (x) -> k(xs,x)
# #     val = quadgk(f,0,1)[1]
# #     val1 = ff.(xi) .* wi |> sum
# #     val2 = w  |> sum
# #     e1 = val - val1 |> abs
# #     e1 = e1 / abs(val)
# #     e2 = val - val2 |> abs
# #     e2 = e2 / abs(val)
# #     @show val, val1, val2
# #     @show e1, e2
# # end

# # xx = [x[1] for x in xi]; yy = [x[2] for x in xi]; scatter(xx,yy,m=:x);p
# # xx = [x[1] for x in xq]; yy = [x[2] for x in xq]; scatter!(xx,yy,m=:o);
# # xx = [x[1] for x in xs]; yy = [x[2] for x in xs]; scatter!(xx,yy,m=:+)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# TODO: doc is missing

OLSbeta(y,x) = inv(x' * x) * (x' * y)

# TODO: make it a doctest string
# a = [10, 14,17,21,20, 18, 42, 51, 77, 11, 91]
# b = [0.11,0.7,0.2,0.19,0.09,0.8, 0.71,0.1, 0.6,0.3,0.81]
# OLSbeta(a,b)

# df = DataFrame(a=a, b=b)
# ols = GLM.lm(@formula(a~b), df)
# coef(ols) # print coefficients

# # Predict
# predict(ols) # predict the model

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia