text:
using Rimu: SplittablesThreading
using Rimu.StochasticStyles: ThresholdCompression, NoCompression

"""
    WorkingMemory

This structure is used when spawning or performing matrix-vector multiplications to
ensure different threads access different regions in memory.

The memory is a large rectangular array of dictionaries ([`Storage`](@ref)), optinally
distributed over MPI. The local part of the memory is segmented into blocks, where each
block belongs to a rank. Each block has `n×n` entries, where `n` is [`num_threads`](@ref).

The purpose of this segmentation into entries is that each entry belongs to a unique
combination of thread and MPI rank.

## Workflow:

```julia
empty!(w)
spawns!(operation, w, src)
merge_rows!(w)
synchronize!(w)
collect!(operation, dst, w)
```
"""
struct WorkingMemory{W,S,NR,ID}
    entries::Matrix{W} # num_segments * num_ranks × num_segments
    style::S
end

rank_id(w::WorkingMemory{<:Any,<:Any,<:Any,ID}) where {ID} = ID

"""
    height(::WorkingMemory)

The number of blocks of a working memory is the same as the number of MPI ranks over which
the vector is distributed.
"""
num_ranks(w::WorkingMemory{<:Any,<:Any,NR}) where {NR} = NR

"""
    height(::WorkingMemory)

The height of a working memory is the number of entries in each of its columns. This is
equivalent to the number of segments (accross all ranks) in the vector.
"""
height(w::WorkingMemory) = size(w.entries, 1)

"""
    num_threads(::WorkingMemory)

The num_threads of a working memory is the number of entries in each of its rows. This is
equivalent to the number of segments in the local part of the vector.
"""
num_threads(w::WorkingMemory) = size(w.entries, 2)

Base.length(w::WorkingMemory) = sum(length, w.entries)

function WorkingMemory(t::TVec; style=t.style)
    nsegs = num_segments(t)
    nranks = num_ranks(t)
    blocks = [empty(t.segments[1]) for i in 1:nsegs * nranks, j in 1:nsegs]

    return WorkingMemory{eltype(blocks),typeof(style),num_ranks(t),rank_id(t)}(
        blocks, style
    )
end
function Rimu.deposit!(w::WorkingMemory, index, key, value, parent=nothing)
    h = hash(key)
    seg = mod1(h, height(w))
    deposit!(w.entries[seg, index], key, value, h, parent)
end
#function Base.getindex(w::WorkingMemory, index, key)
#    h = hash(key)
#    seg = mod1(h, height(w))
#    return w.entries[seg, index][key]
#end

"""
    get_diagonal(w::WorkingMemory, i)

Get the `i`-th diagonal entry in working memory. This is the entry that is transferred back
to the vector in the [`merge_and_compress!`](@ref) step.

NOTE: diagonal is diagonal in the block.
"""
function get_diagonal(w::WorkingMemory, index)
    return w.entries[index + rank_id(w) * num_threads(w), index]
end
function Base.empty!(w::WorkingMemory)
    foreach(empty!, w.entries)
    return w
end

Base.getindex(w::WorkingMemory, r::Int, c::Int) = w.entries[r, c]
Base.getindex(w::WorkingMemory, ::Colon, c) = WorkingMemoryColumn(w, c)

"""
    WorkingMemoryColumn

A column in a [`WorkingMemory`](@ref). Used to allow spawning from a vector in a thread-safe
manner. Supports enough of the `AbstractDVec` interface to be usable as a target for
[`fciqmc_col!`](@ref).
"""
struct WorkingMemoryColumn{W}
    working_memory::W
    index::Int
end
function Rimu.deposit!(s::WorkingMemoryColumn, args...)
    return deposit!(s.working_memory, s.index, args...)
end
function Rimu.getindex(s::WorkingMemoryColumn, key)
    return s.working_memory[s.index, key]
end
function Base.empty!(s::WorkingMemoryColumn)
    return foreach(empty!, view(s.working_memory.entries, :, s.index))
end
Base.keytype(s::WorkingMemoryColumn) = keytype(eltype(s.working_memory.entries))
Base.valtype(s::WorkingMemoryColumn) = valtype(eltype(s.working_memory.entries))

###
### Operations
###
function merge_rows!(w::WorkingMemory)
    nlocal = num_threads(w)

    Folds.foreach(1:height(w)) do i
        diag_index = mod1(i, nlocal)
        for j in 1:nlocal
            j == diag_index && continue
            add!(w.entries[i, diag_index], w.entries[i, j])
        end
    end
end

function exchange!(w::WorkingMemory, rank_id, thread_id)
    @assert num_threads(w) ≥ 2 # or we can't find the recieving buffer

    # Maybe this would be more efficient if all ranks (asynchronously) sent first, then all
    # rank recieved, then everything got moved around.

    row = rank_id * num_threads(w) + thread_id
    send_arr = w[row, thread_id].pairs
    recv_arr = w[row, mod1(thread_id + 1, num_threads(w))].pairs

    MPI.Isend(MPI.Buffer(send_arr), rank_id, thread_id, MPI.COMM_WORLD)

    recv_len = MPI.Get_count(
        MPI.Probe(rank_id, thread_id, MPI.COMM_WORLD), eltype(recv_arr)
    )
    resize!(recv_arr, recv_len)
    MPI.Recv!(recv_arr, rank_id, thread_id, MPI.COMM_WORLD)

    # Move to diagonal
    diag = get_diagonal(w, thread_id)
    Threads.@spawn for (k, v) in recv_arr
        deposit!(diag, k, v)
    end
end
function mpi_send!(w::WorkingMemory, rank_id, thread_id)
    @assert num_threads(w) ≥ 2 # or we can't find the recieving buffer

    # Maybe this would be more efficient if all ranks (asynchronously) sent first, then all
    # rank recieved, then everything got moved around.
    row = rank_id * num_threads(w) + thread_id
    send_arr = w[row, thread_id].pairs

    MPI.Isend(MPI.Buffer(send_arr), rank_id, thread_id, MPI.COMM_WORLD)
end
function mpi_recv!(w::WorkingMemory, rank_id, thread_id)
    row = rank_id * num_threads(w) + thread_id
    recv_arr = w[row, mod1(thread_id + 1, num_threads(w))].pairs

    recv_len = MPI.Get_count(
        MPI.Probe(rank_id, thread_id, MPI.COMM_WORLD), eltype(recv_arr)
    )
    resize!(recv_arr, recv_len)
    MPI.Recv!(recv_arr, rank_id, thread_id, MPI.COMM_WORLD)
end
function mpi_collect!(w::WorkingMemory, rank_id, thread_id)
    # Move to diagonal
    row = rank_id * num_threads(w) + thread_id
    recv_arr = w[row, mod1(thread_id + 1, num_threads(w))].pairs

    diag = get_diagonal(w, thread_id)
    for (k, v) in recv_arr
        deposit!(diag, k, v)
    end
end

function synchronize!(w::WorkingMemory)
    if num_ranks(w) > 1
        foreach(1:num_threads(w)) do thread_id
            for rank in 0:(num_ranks(w) - 1)
                rank == rank_id(w) && continue
                mpi_send!(w, rank, thread_id)
            end
        end
        foreach(1:num_threads(w)) do thread_id
            for rank in 0:(num_ranks(w) - 1)
                rank == rank_id(w) && continue
                mpi_recv!(w, rank, thread_id)
            end
        end
        Folds.foreach(1:num_threads(w)) do thread_id
            for rank in 0:(num_ranks(w) - 1)
                rank == rank_id(w) && continue
                mpi_collect!(w, rank, thread_id)
            end
        end
    end
    return w
end

function move_and_compress!(t::ThresholdCompression, dst::TVec, src::WorkingMemory)
    Folds.foreach(1:num_segments(dst)) do i
        dst_seg = dst.segments[i]
        src_seg = get_diagonal(src, i)
        empty!(dst_seg)
        for (add, val) in pairs(src_seg)
            prob = abs(val) / t.threshold
            if prob < 1 && prob > cRand()
                dst_seg[add] = t.threshold * sign(val)
            elseif prob ≥ 1
                dst_seg[add] = val
            end
        end
    end
    return dst
end
function move_and_compress!(::NoCompression, dst::TVec, src::WorkingMemory)
    Folds.foreach(1:num_segments(dst)) do i
        dst_seg = dst.segments[i]
        src_seg = get_diagonal(src, i)
        copy!(dst_seg, src_seg)
    end
    return dst
end

function LinearAlgebra.mul!(dst, op, src, w, style=StochasticStyle(src))
    T = valtype(dst)
    # Perform spawns. Note that setting shift to 1 and dτ to -1 turns this into regular
    # matrix-vector multiply.
    Folds.foreach(1:num_segments(src)) do i
        strip = w[:, i]
        empty!(strip)
        foreach(pairs(src.segments[i])) do (add, val)
            fciqmc_col!(style, strip, op, add, val, one(T), -one(T))
        end
    end
    merge_rows!(w)
    synchronize!(w)
    move_and_compress!(CompressionStrategy(style), dst, w)
    return dst
end

function Base.:*(op::AbstractHamiltonian, tv::TVec)
    wm = WorkingMemory(tv)
    dst = similar(tv, promote_type(eltype(op), valtype(tv)))
    mul!(dst, op, tv, wm)
end

###
### Rimu compat.
###
function Rimu.working_memory(::SplittablesThreading, t::TVec)
    return WorkingMemory(t)
end
function Rimu.fciqmc_step!(
    ::SplittablesThreading, w::WorkingMemory, ham, src::TVec, shift, dτ
)
    stat_names, stats = step_stats(src, Val(1))
    style = StochasticStyle(src)
    result = Folds.sum(1:num_segments(src); init=stats) do i
        strip = w[:, i]
        empty!(strip)
        sum(pairs(src.segments[i]); init=stats) do (add, val)
            fciqmc_col!(style, strip, ham, add, val, shift, dτ)
        end
    end
    return stat_names, result
end
function Rimu.apply_memory_noise!(w::WorkingMemory, t::TVec, args...)
    return 0.0
end
function Rimu.sort_into_targets!(dst::TVec, w::WorkingMemory, stats)
    merge_rows!(w)
    synchronize!(w)
    move_and_compress!(CompressionStrategy(StochasticStyle(dst)), dst, w)
    return dst, w, stats
end
function Rimu.StochasticStyles.compress!(::ThresholdCompression, t::TVec)
    return t
end


text:
## Looking at the results
using Plots

ppes = map(v -> getPPE(getVariable(newfg, v)).suggested, ls(newfg))
x = map(ppe -> ppe[1], ppes); y = map(ppe -> ppe[2], ppes); h = map(ppe -> ppe[3], ppes)
Plots.plot(x, y, title="Path Plot", lw=3)

## Stuff
using Optim

cost(tf, im1, im2) = evaluateTransform(im1,im2, tf )

# Plotting
xrange = -100.0:1.0:100.0
hrange = -pi:0.1:pi
val = reshape(
  [sweepx(sweeps[10],sweeps[11],xrange); sweepx(sweep_original[10],sweep_original[11],xrange)],
  length(xrange), 2)
Plots.plot(xrange,val)
# Heading
val = reshape(
  [sweeph(sweeps[10],sweeps[11],hrange); sweeph(sweep_original[10],sweep_original[11],hrange)],
  length(hrange), 2)
Plots.plot(hrange,val)

corr_func = (a,b)->sqrt(sum((a .- 0.5).*(b .- 0.5)))
val = reshape(
  [sweepx(sweeps[10],sweeps[11],xrange,diff_func=corr_func);
  sweepx(sweep_original[10],sweep_original[11],xrange,diff_func=corr_func)],
  length(xrange), 2)
Plots.plot(xrange,val)

## Sweep plotting
# sanity check: identity transform should yield zero cost
# @assert evaluateTransform(sweeps[11],sweeps[11],0.,0.,0.) == 0 "There's error with no transform!"

# let's try small displacements:
# sweepx(im1, im2, xrange) = (x->@show evaluateTransform(im1,im2,x,0.,0.)).(xrange)
# sweepy(im1, im2, yrange) = (y->@show evaluateTransform(im1,im2,0.,y,0.)).(yrange)
# sweeph(im1, im2, hrange) = (h->@show evaluateTransform(im1,im2,0.,0.,h)).(hrange)

# using Plots
# xrange = -10:0.1:10
# hrange = -pi:0.1:pi
# Plots.plot(xrange,sweepx(sweeps[10],sweeps[11],xrange))
# Plots.plot(xrange,sweepy(sweeps[10],sweeps[11],xrange))
# Plots.plot(hrange,sweeph(sweeps[10],sweeps[11],hrange))

# fs10 = imfilter(sweeps[10],Kernel.gaussian(3))
# fs11 = imfilter(sweeps[11],Kernel.gaussian(3))
# ffs10 = imfilter(fs10,Kernel.gaussian(3))
# ffs11 = imfilter(fs11,Kernel.gaussian(3))
#
# Plots.plot(xrange,sweepx(ffs10,ffs11,xrange))
# Plots.plot(xrange,sweepy(fs10,fs11,xrange))

ext: jl


text:
@testsection "Bases" begin
    for T ∈ (Sym, Float64)
        @testsection "type $T" begin
            sq2 = √T(2)
            v = T[1 0 0; 0 1 0; 0 1 1]
            bv = Basis(v)
            @test eltype(bv) == T
            @test vecbasis(bv) == T[1 0 0; 0 1 0; 0 1 1]
            @test vecbasis(bv, :cont) == T[1 0 0; 0 1 -1; 0 0 1]
            @test metric(bv) == T[1 0 0; 0 2 1; 0 1 1]
            @test metric(bv, :cont) == bv.gⁱʲ == T[1 0 0; 0 1 -1; 0 -1 2]
            @test vecbasis(bv, :cont)' * vecbasis(bv, :cov) == I
            @test metric(bv, :cont)' * metric(bv, :cov) == I
            @test !isorthogonal(bv)

            bw = Basis(vecbasis(bv, :cont), :cont)
            @test eltype(bw) == T
            @test vecbasis(bw) == T[1 0 0; 0 1 0; 0 1 1]
            @test vecbasis(bw, :cont) == T[1 0 0; 0 1 -1; 0 0 1]
            @test metric(bw) == T[1 0 0; 0 2 1; 0 1 1]
            @test metric(bw, :cont) == T[1 0 0; 0 1 -1; 0 -1 2]
            @test vecbasis(bw, :cont)' * vecbasis(bw, :cov) == I
            @test metric(bw, :cont)' * metric(bw, :cov) == I
            @test !isorthogonal(bw)

            nb = normalize(bv)
            @test opequal(vecbasis(nb), vecbasis(nb, :cov))
            @test opequal(vecbasis(nb), T[1 0 0; 0 sq2/2 0; 0 sq2/2 1])
            θ, ϕ, ψ = T == Sym ? symbols("θ, ϕ, ψ", real = true) : rand(T, 3)
            br = Basis(θ, ϕ, ψ)
            @test br isa RotatedBasis
            @test opequal(metric(br), I)
            @test opequal(vecbasis(br, :cont), vecbasis(br, :cov))
            if T == Sym
                @test angles(br) == (θ = θ, ϕ = ϕ, ψ = ψ)
            else
                angbr = angles(br)
                @test angbr.θ ≈ θ && angbr.ϕ ≈ ϕ && angbr.ψ ≈ ψ
                angbr = angles(Array(br))
                @test angbr.θ ≈ θ && angbr.ϕ ≈ ϕ && angbr.ψ ≈ ψ
            end
            λ = T == Sym ? [symbols("λ$i", real = true) for i in 1:getdim(br)] : rand(getdim(br))
            e = vecbasis(br, :cov) .* λ'
            bo = Basis(e)
            @test bo isa OrthogonalBasis
            @test bo isa AllOrthogonalBasis

            br = Basis(θ)
            @test opequal(metric(br), I)
            @test opequal(vecbasis(br, :cont), vecbasis(br, :cov))
            if T == Sym
                @test angles(br) == (θ = θ,)
            else
                @test angles(br).θ ≈ θ
                @test angles(Array(br)).θ ≈ θ
            end
            λ = T == Sym ? [symbols("λ$i", real = true) for i in 1:getdim(br)] : rand(getdim(br))
            e = vecbasis(br, :cov) .* λ'
            bo = Basis(e)
            @test bo isa OrthogonalBasis
            @test bo isa AllOrthogonalBasis

            b2 = Basis(zero(T))
            @test b2 isa CanonicalBasis
            @test getdim(b2) == 2
            b3 = Basis(zero(T), zero(T), zero(T))
            @test b3 isa CanonicalBasis
            @test getdim(b3) == 3

            if T ≠ Sym
                v = rand(T, 3, 3)
                while det(v) ≈ 0
                    v = rand(T, 3, 3)
                end
                @test eltype(bv) == T
                @test vecbasis(bv, :cont)' * vecbasis(bv, :cov) == I
                @test metric(bv, :cont)' * metric(bv, :cov) == I
            end

        end
    end
end


text:
function matrixfma(output, mtx, input::AbstractVector, bias, isize, osize)
  for idx = 1:osize
    @inbounds output[idx] = bias[idx]
    @inbounds for jdx = 1:isize
      output[idx] += mtx[idx, jdx] * input[jdx]
    end
  end
end

function matrixfma(output, mtx, input::AbstractMatrix, bias, isize, osize)
  for mdx = 1:size(input, 2)
    for idx = 1:osize
      @inbounds output[idx, mdx] = bias[idx]
      @inbounds for jdx = 1:isize
        output[idx, mdx] += mtx[idx, jdx] * input[jdx, mdx]
      end
    end
  end
end

function evaluate!{F, i, o, tf}(output::AbstractArray{F}, fcl::FullyConnectedLayer{F, i, o, tf}, input::AbstractArray)
  matrixfma(output, fcl.transition, input, fcl.bias, i, o)
  for idx = 1:o
    output[idx] = tf(output[idx])
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
module Woody

export @runtimes,
       @runduring,
       Collector,
       Reporter,
       Controller,
       report,
       finalize

include("core.jl")

end


text:
using Random

@testset "SpatialFields" begin

σ = rand()
g = Gaussian(σ,0,1)
@test g(σ) ≈ exp(-1)/sqrt(π*σ^2)

@test radius(g) == σ
@test center(g) == 0
@test strength(g) == 1

g = EmptySpatialField()
@test g(rand(2)...) == 0.0

g = SpatialGaussian(σ,σ,0.0,0.5,1)
@test g(0,0.5+σ) ≈ g(σ,0.5) ≈ g(-σ,0.5) ≈ g(0,0.5-σ) ≈ exp(-1)/(π*σ^2)

g = EmptySpatialField()
for x in [-0.5,0,0.5], y in [-0.5,0,0.5]
  g += SpatialGaussian(0.2,0.5,x,y,1)
end
@test g(0.5,0) ≈ 5.5357580598289235

g2 = -g
@test g2(0.5,0) ≈ -5.5357580598289235

gr = PhysicalGrid((-3.0,3.0),(-2.0,2.0),0.03)
w = Nodes(Dual,size(gr))
xg, yg = coordinates(w,gr)

gfield = GeneratedField(w,g,gr)

@test gfield()[104,70] == g(xg[104],yg[70])
@test datatype(gfield) == typeof(w)

wfield = PulseField(gfield,0.5,0.1)

@test maximum(wfield(0.5)) ≈ maximum(gfield()) ≈ 5.532960088678624
@test maximum(wfield(2)) ≈ 0.0
@test datatype(wfield) == datatype(gfield) <: ScalarGridData

q = Edges(Primal,size(gr))
gauss = SpatialGaussian(0.5,0,0,1)
gaussfield = GeneratedField(q,gauss,EmptySpatialField(),gr)
q .= gaussfield()

ffield = PulseField(gaussfield,0.5,0.1)

@test maximum(ffield(0.5).u) ≈ maximum(q.u) ≈ 1.272094144652253
@test maximum(ffield(2).u) ≈ 0.0
@test maximum(ffield(0.5).v) ≈ maximum(q.v) ≈ 0.0
@test maximum(ffield(2).v) ≈ maximum(q.v) ≈ 0.0

@test datatype(ffield) == datatype(gaussfield) <: VectorGridData

gq = EdgeGradient(Dual,size(gr))
tfield = GeneratedField(gq,[g,gauss,EmptySpatialField(),EmptySpatialField()],gr)
gq .= tfield()
@test maximum(gq.dudx) ≈ 5.532960088678624
@test maximum(gq.dudy) ≈ 1.2732395447351625
@test maximum(gq.dvdx) == maximum(gq.dvdy) == 0.0

ffield = PulseField(tfield,0.5,0.1)
@test datatype(ffield) <: TensorGridData

end


text:
"""
This module contains all the constants used in the package, as well as the macro
that they were created by.
"""
module Constants

using ..Exceptions
using SyntaxTree

"""
    @constant(name::Symbol, docstring::String, value::Union{Expr, String}) -> Expr

Create a constant with a documentation string and export it.

# Arguments
- `name::Symbol`: name of the constant
- `docstring::String`: documentation string
- `value::Union{Expr, String}`: value of the constant

# Returns
- `Expr`: constant definition

# Example
```jldoctest; output = false
using PDFHighlights
using Suppressor
using SyntaxTree

name = :CONSTANT_NAME
docstring = "Docstring"
value = "Value"

d1 = @capture_out @macroexpand(
    PDFHighlights.Internal.Constants.@constant(
        CONSTANT_NAME,
        "Docstring",
        "Value"
    )
) |> linefilter! |> dump

macro oof()
    return esc(
        quote
            @doc \$(docstring)
            const \$(name) = \$(value)
            export \$(name)
        end
    )
end

d2 = @capture_out @macroexpand(@oof) |> linefilter! |> dump

d1 == d2

# output

true
```
"""
macro constant(name::Symbol, docstring::String, value::Union{Expr, String})
    name_string = "$(name)"
    if endswith(name_string, "_EXCEPTION")
        docstring = """
        A documentation constant for the [`$(docstring)`](@ref) exception.
        """
    end
    if typeof(value) == String
        docstring = """
            $(name) = "$(value)"

        $(docstring)
        """
    end
    return esc(
        quote
            @doc $(docstring)
            const $(name) = $(value)
            export $(name)
        end
    )
end

# For the `CSV` module

@constant(
    HEADER,
    "Default (expected) header for CSV files.",
    "Highlight,Title,Author,Note,Location",
)

# Documentation constants

## For exceptions

@constant(
    INCORRECT_HEADER_EXCEPTION,
    "IncorrectHeader",
    "the specified file has an incorrect header",
)
@constant(
    LAST_ELEMENT_IS_NOT_AN_INTEGER_EXCEPTION,
    "LastElementIsNotAnInteger",
    "the last element in the line is not an integer",
)
@constant(
    NOT_FIVE_COLUMNS_EXCEPTION,
    "NotFiveColumns",
    "the row does not represent elements for 5 columns",
)
@constant(
    INTEGRITY_CHECK_FAILED_EXCEPTION,
    "IntegrityCheckFailed",
    "another exception was thrown while checking the integrity of the table",
)

@constant(
    NOT_CSV_EXCEPTION,
    "NotCSV",
    "the specified path does not end in `.csv`",
)
@constant(
    NOT_CSV_OR_DIR_EXCEPTION,
    "NotCSVorDir",
    "the specified path does not end in `.csv` and is not a directory",
)
@constant(
    NOT_CSV_OR_PDF_OR_DIR_EXCEPTION,
    "NotCSVorPDForDir",
    "the specified path does not end in `.csv` or `.pdf` and is not a directory",
)

@constant(
    FILE_DOES_NOT_EXIST_EXCEPTION,
    "FileDoesNotExist",
    "the specified file doesn't exist",
)
@constant(
    NOT_PDF_EXCEPTION,
    "NotPDF",
    "the specified path does not end in `.pdf`",
)
@constant(
    DIRECTORY_DOES_NOT_EXIST_EXCEPTION,
    "DirectoryDoesNotExist",
    "the specified directory doesn't exist",
)
@constant(
    DOES_NOT_EXIST_EXCEPTION,
    "DoesNotExist",
    "the specified file or directory doesn't exist",
)

## For arguments

@constant(
    CSV_ARGUMENT,
    "A documentation constant for the `csv` argument.",
    "absolute or relative path to the CSV file",
)
@constant(
    PDF_ARGUMENT,
    "A documentation constant for the `pdf` argument.",
    "absolute or relative path to the PDF file",
)
@constant(
    TARGET_CSV_ARGUMENT,
    "A documentation constant for the `target` argument that mentions CSV.",
    "a CSV file or a directory with PDF files",
)
@constant(
    TARGET_PDF_ARGUMENT,
    "A documentation constant for the `target` argument that mentions PDF.",
    "a PDF file or a directory with PDF files",
)
@constant(
    HIGHLIGHTS_ARGUMENT,
    "A documentation constant for the `highlights` argument.",
    "the highlights vector",
)
@constant(
    COMMENTS_ARGUMENT,
    "A documentation constant for the `comments` argument.",
    "the comments vector",
)
@constant(
    PAGES_ARGUMENT,
    "A documentation constant for the `pages` argument.",
    "the pages vector",
)
@constant(
    DIR_ARGUMENT,
    "A documentation constant for the `dir` argument.",
    "a directory with PDF files",
)

## For keywords

@constant(
    CONCATENATE_KEYWORD,
    "A documentation constant for the `concatenate` keyword.",
    "if `true`, concatenate the highlights",
)
@constant(
    CONCATENATE_BOTH_KEYWORD,
    "A documentation constant for the `concatenate` keyword (for the `Both` module).",
    "if `true`, concatenate the highlights (only for PDFs)",
)

end


text:
module ConstrainedDynamics

using LinearAlgebra
using StaticArrays
using StaticArrays: SUnitRange
using Rotations

using CoordinateTransformations
using GeometryTypes: GeometryTypes, Vec, Point, GLUVMesh
using Blink
using Colors: RGBA, RGB
using FileIO, MeshIO
using MeshCat
using LightXML

using Plots

export Box,
    Cylinder,
    Sphere,
    Mesh,

    Quaternion,
    Origin,
    Body,
    EqualityConstraint,
    InequalityConstraint,
    Mechanism,
    PID,

    OriginConnection,
    Prismatic,
    Spherical,
    Cylindrical,
    Revolute,
    Planar,
    Fixed,
    FixedOrientation,
    CylindricalFree,

    Impact,
    Friction,

    setPosition!,
    setVelocity!,
    setForce!,
    simulate!,
    plotθ,
    plotλ,
    visualize!,

    disassemble,
    getcomponent,
    getbody,
    geteqconstraint,
    getineqconstraint,
    minimalCoordinates,

    RotX,
    RotY,
    RotZ,
    RGBA

include(joinpath("util", "util.jl"))
include(joinpath("util", "customdict.jl"))
include(joinpath("util", "quaternion.jl"))
include(joinpath("util", "shapes.jl"))
include(joinpath("components", "component.jl"))

include(joinpath("components", "body.jl"))
include(joinpath("joints", "joint.jl"))
include(joinpath("bounds", "bound.jl"))

include(joinpath("joints", "translational.jl"))
include(joinpath("joints", "translational0.jl"))
include(joinpath("joints", "translational1.jl"))
include(joinpath("joints", "translational2.jl"))
include(joinpath("joints", "translational3.jl"))
include(joinpath("joints", "rotational.jl"))
include(joinpath("joints", "rotational0.jl"))
include(joinpath("joints", "rotational1.jl"))
include(joinpath("joints", "rotational2.jl"))
include(joinpath("joints", "rotational3.jl"))
include(joinpath("joints", "prototypes.jl"))

include(joinpath("bounds", "impact.jl"))
include(joinpath("bounds", "friction.jl"))

include(joinpath("components", "equalityconstraint.jl"))
include(joinpath("components", "inequalityconstraint.jl"))

include(joinpath("util", "graph.jl"))
include(joinpath("util", "storage.jl"))

include(joinpath("control", "controller.jl"))

include(joinpath("solver", "sparseldu.jl"))
include(joinpath("components", "mechanism.jl"))
include(joinpath("components", "mechanism_functions.jl"))
include(joinpath("components", "initialize.jl"))
include(joinpath("solver", "solverfunctions.jl"))

include(joinpath("util", "urdf.jl"))

include(joinpath("solver", "newton.jl"))

include(joinpath("util", "visualize.jl"))

end


text:
using Test
using ModelingToolkit, OrdinaryDiffEq

@isdefined(t) || @parameters t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Inductor(; name, L = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L=L
    D = Differential(t)
    eqs = [
           D(i) ~ v / L
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

@connector function HeatPort(;name)
    @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
    ODESystem(Equation[], t, [T, Q_flow], [], name=name)
end

function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
    @named p = Pin()
    @named n = Pin()
    @named h = HeatPort()
    @variables v(t) RTherm(t)
    @parameters R=R TAmbient=TAmbient alpha=alpha
    eqs = [
           RTherm ~ R*(1 + alpha*(h.T - TAmbient))
           v ~ p.i * RTherm
           h.Q_flow ~ -v * p.i # -LossPower
           v ~ p.v - n.v
           0 ~ p.i + n.i
          ]
    compose(ODESystem(
        eqs, t, [v, RTherm], [R, TAmbient, alpha],
        name=name,
    ), p, n, h)
end

function HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)
    @parameters rho=rho V=V cp=cp
    C = rho*V*cp
    @named h = HeatPort()
    D = Differential(t)
    eqs = [
           D(h.T) ~ h.Q_flow / C
          ]
    compose(ODESystem(
        eqs, t, [], [rho, V, cp],
        name=name,
    ), h)
end


text:
include("typeDefine.jl")
include("maxThroughputPrediction.jl")
include("lockPrediction.jl")
include("findClosestValue.jl")
include("getConcurrencyLevel.jl")
include("cfFlushRateApprox.jl")
include("probOfBeingChosenAtLeastOnce.jl")
include("recpow.jl")
include("mae.jl")
include("mre.jl")
include("useLockModel.jl")
include("prepareLockModel.jl")
include("Div.jl")
include("bmult.jl")
include("Max.jl")
include("basicThomasian.jl")
include("valueToString.jl")
include("implicitUseLockModel.jl")
include("barzanCurveFit.jl")
include("barzanLinInvoke.jl")
include("barzanLinSolve.jl")
include("barzanRegressTreeLearn.jl")
include("barzanRegressTreeInvoke.jl")
include("DoSmooth.jl")
include("localmax.jl")
include("findMaxThroughput.jl")

tc=PredictionConfig()
tc.io_conf=[]
tc.lock_conf=[]
tc.transactionType=[1 2 3 4 5]
tc.groupingStrategy=[]
tc.initialized=1
tc.configSummary="{0-101"
tc.currentLockWait=[0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0]
tc.lockWaitTime=[0
0.0220000000000000
0.0160000000000000
0.0150000000000000
0.0160000000000000
0.0290000000000000
0.0870000000000000
0.00700000000000000
0.0250000000000000
0.0530000000000000
0.0340000000000000
0.0680000000000000
0.0160000000000000
0
0.0190000000000000
0.00400000000000000
0.0110000000000000
0.0180000000000000
0
0.00200000000000000
0.00700000000000000
0.00300000000000000
0.0110000000000000
0.00500000000000000
0.0270000000000000
0.0240000000000000
0.0120000000000000
0.0290000000000000
0.00900000000000000
0.0200000000000000
0.0110000000000000
0.00200000000000000
0.00900000000000000
0.0320000000000000
0
0.0120000000000000
0.00300000000000000
0.0100000000000000
0
0
0
0
0.0200000000000000
0.00900000000000000
0.0110000000000000
0.00200000000000000
0
0
0.0220000000000000
0.00800000000000000
0.0230000000000000
0.00900000000000000
0
0
0
0
0]
tc.pagesFlushed=[0
0
29.4000000000000
63.2857142857143
82.0000000000000
98.5555555555555
115.111111111111
131.777777777778
147.888888888889
147.888888888889
147.777777777778
148.000000000000
147.888888888889
148.111111111111
147.888888888889
147.888888888889
147.888888888889
148.222222222222
148.111111111111
148.333333333333
148.222222222222
148.333333333333
147.666666666667
147.666666666667
147.666666666667
147.333333333333
147.222222222222
147.444444444444
147.333333333333
147.333333333333
144.000000000000
140.777777777778
137.666666666667
134.111111111111
130.444444444444
126.777777777778
123.444444444444
119.777777777778
116.222222222222
116.222222222222
116.444444444444
116.111111111111
116
116.444444444444
116.555555555556
116.333333333333
116.444444444444
116.222222222222
116.000000000000
115.888888888889
127.111111111111
167.888888888889
208.111111111111
234.571428571429
282.200000000000
393.333333333333
480]
tc.TPS=[99
101
99
101
99
100
101
99
101
99
101
99
100
101
100
99
101
99
101
100
100
100
100
100
100
100
99
101
100
101
98
101
100
100
100
101
98
101
100
99
101
100
100
100
101
99
101
98
101
100
101
100
100
100
26
0
0]
tc.TPSUngrouped=[99
101
99
101
99
100
101
99
101
99
101
99
100
101
100
99
101
99
101
100
100
100
100
100
100
100
99
101
100
101
98
101
100
100
100
101
98
101
100
99
101
100
100
100
101
99
101
98
101
100
101
100
100
100
26
0
0]
tc.transactionCount=[40	53	5	1	0
41	42	10	4	4
38	47	5	6	3
38	52	6	3	2
44	43	6	3	3
52	36	2	6	4
46	43	7	1	4
37	49	6	4	3
49	38	4	5	5
48	45	2	3	1
44	41	6	7	3
44	37	4	4	10
43	40	7	8	2
44	44	7	4	2
54	38	2	4	2
42	45	4	5	3
46	39	2	6	8
50	42	0	3	4
45	47	2	3	4
46	45	2	2	5
43	46	5	3	3
48	41	4	3	4
52	35	4	6	3
42	39	7	6	6
40	45	5	3	7
42	41	8	7	2
42	46	6	2	3
53	39	2	6	1
46	46	4	2	2
43	47	3	3	5
45	35	3	7	8
52	40	3	3	3
42	44	6	3	5
47	38	5	4	6
50	44	2	3	1
45	45	3	5	3
49	41	2	1	5
51	41	3	4	2
43	43	5	4	5
52	38	2	3	4
52	41	2	4	2
48	39	6	3	4
57	30	6	1	6
47	48	1	0	4
42	43	9	4	3
47	36	4	8	4
49	37	6	2	7
34	43	8	7	6
40	49	3	5	4
48	42	3	3	4
41	44	1	7	8
45	47	3	5	0
45	40	5	6	4
36	49	7	2	6
13	8	1	2	2
0	0	0	0	0
0	0	0	0	0]
tc.averageCpuUsage=[4.11956250000000
4.07018750000000
5.16787500000000
4.38475000000000
4.39287500000000
4.41987500000000
4.11575000000000
3.56131250000000
4.31768750000000
4.43718750000000
4.18150000000000
5.19593750000000
3.18137500000000
3.49325000000000
3.94462500000000
9.88331250000000
10.1602500000000
8.01743750000000
4.62093750000000
3.87143750000000
3.25200000000000
3.86825000000000
4.49712500000000
3.87268750000000
4.05700000000000
4.25318750000000
3.86637500000000
3.81187500000000
3.80843750000000
3.63162500000000
4.31850000000000
4.13968750000000
3.85312500000000
4.00525000000000
3.76906250000000
3.22900000000000
3.62962500000000
5.34856250000000
3.38562500000000
4.10531250000000
3.44781250000000
3.55137500000000
4.69587500000000
4.31081250000000
3.87256250000000
3.17537500000000
3.74956250000000
3.12956250000000
3.25768750000000
3.17718750000000
3.63018750000000
2.49456250000000
3.57793750000000
2.98962500000000
4.61468750000000
0.375625000000000
0.500000000000000]

pc=PredictionCenter()

pc.TEST_MODE_MIXTURE_TPS=0
pc.TEST_MODE_DATASET=1
pc.NUM_TPS_SAMPLES=100
pc.testVar=[]
pc.taskName="WhatIfAnalysisLatency"
pc.workloadName=[]
pc.lockType="waitTime"
pc.learnLock=0
pc.whichTransactionToPlot=1
pc.ioConf=[700000 1000 10]
pc.lockConf=[0.200000000000000 10 100 1]
pc.trainConfig = tc
pc.testMode=0
pc.testConfig=[]
pc.testMixture=[0.453657637737240 0.423069836005159 0.043486272341994 0.040353786622443 0.039432467293164]
pc.testMinTPS=0
pc.testMaxTPS=202
pc.testWorkloadRatio=2
pc.testSampleTPS=[0.232513540079662
1.35649349233245
7.22407212036182
13.7345392309421
18.3463037290628
18.4049196646801
21.3370990724624
24.0813372929466
24.1180420550557
30.8496824885228
31.0386569534439
35.5266319680732
38.1097193118812
38.7325416032832
42.2048543300566
45.2560862264922
47.7185765527469
48.6228211669924
49.0556188602702
51.4676116325951
51.8010804302878
53.4853633480773
53.8272411373726
54.3504405105757
56.7630711118419
58.0746309593585
58.3910434782443
64.2509632471739
65.1393050517293
65.6794277277531
68.9071706239201
69.5814070828105
70.7440387151032
78.1277812932389
84.0640351738989
85.5374896304731
85.7185060427058
85.9972260579799
88.8971980783476
90.9278147871110
92.3997218689102
92.6625497170713
93.1051059378507
93.4147301669504
95.2141450495436
95.6441705789947
96.8515714396753
104.646525888943
105.544051766975
106.482833835652
107.329449126266
109.229346099131
109.864958667927
110.032654326406
110.411809147245
116.394294893927
117.614331234500
121.638438491523
122.675396018498
122.692621181206
123.413649066733
123.919068836201
128.428915523845
128.817237810579
129.142026130102
130.177441849130
130.401478744450
130.756919019212
130.818761294882
133.712839884942
133.725938868545
134.902210856771
135.705249413671
136.417077280894
136.576705380478
137.161384326827
138.039375145520
140.418380909451
143.274903947530
145.651417357122
145.795122745003
149.162221675342
154.631003357531
155.572265178939
155.597673990339
157.318052848467
157.664969851734
158.517337541670
168.249197453194
170.069688843644
170.567215618495
175.728088536317
176.825062918045
185.030231314548
185.319717093975
189.845553009674
190.611770158774
190.925170848099
193.454175848000
200.728334072412]
pc.testSampleTransactionCount=[0.105481543334462	0.0983694652704815	0.0101111471271052	0.00938280178320359	0.00916858256440899
0.615383633337480	0.573891479343157	0.0589888454377113	0.0547396489443168	0.0534898852697890
3.27725549296680	3.05628700725092	0.314147967644258	0.291518664890223	0.284862987609624
6.23077862291865	5.81066926004112	0.597263913488547	0.554240665483016	0.541586769010801
8.32294081093653	7.76176770995543	0.797812360430959	0.740342825993136	0.723440021746717
8.34953237782243	7.78656634422433	0.800361348970793	0.742708200951711	0.725751392710804
9.67973796137875	9.02708300541252	0.927870901253203	0.861032743112083	0.841374461305871
10.9246825898719	10.1880874193119	1.04720759188049	0.971773146702655	0.949586545179763
10.9413339855440	10.2036160969980	1.04880374516181	0.973254322840833	0.951033904511134
13.9951940826872	13.0515701112306	1.34153769435994	1.24490150451198	1.21647909573316
14.0809237920359	13.1315195071139	1.34975548940718	1.25252733974649	1.22393082514041
16.1169279353964	15.0302463605484	1.54492079315741	1.43363412585370	1.40090275311732
17.2887652378573	16.1230726994802	1.65724963287340	1.53787148135286	1.50276026031741
17.5713133273048	16.3865700241640	1.68433385265797	1.56300471920380	1.52731967995257
19.1465545164178	17.8556008000387	1.83533178955102	1.70312568606641	1.66424153798270
20.5307691707433	19.1464849780774	1.96801849077799	1.82625444695076	1.78455913994275
21.6478967151028	20.1882903565703	2.07510301574503	1.92562525613628	1.88166120919253
22.0581141907380	20.5708489772277	2.11442524330389	1.96211495035403	1.91731780536878
22.2544561698885	20.7539526263461	2.13324600166275	1.97957997611925	1.93438408625351
23.3486751132208	21.7743940129792	2.23813457624699	2.07691301778852	2.02949491235956
23.4999557802411	21.9154746025319	2.25263589120102	2.09036974649586	2.04264440981787
24.2640435900067	22.6280439003474	2.32587907686498	2.15833693997216	2.10905984088604
24.4191390602932	22.7726820805984	2.34074606751795	2.17213300333234	2.12254092563069
24.6564924520061	22.9940319536174	2.36349805795022	2.19324607919957	2.14317196780232
25.7510007513092	24.0147431864362	2.46841436933752	2.29060485968185	2.23830794507724
26.3459998934846	24.5696245960360	2.52544921805945	2.34353126591110	2.29002598586746
26.4895428493528	24.7034891885109	2.53920881902813	2.35629970918288	2.30250291216957
29.1479402090551	27.1826444841554	2.79403488600204	2.59276966116291	2.53357400679846
29.5509432536130	27.5584751057252	2.83266555964771	2.62861761679173	2.56860351595174
29.7959740309062	27.7869847176932	2.85615348143537	2.65041361200994	2.58990188570834
31.2602642484045	29.1525453754414	2.99651598806802	2.78066526011397	2.71717975189219
31.5661367676211	29.4377944835329	3.02583601834223	2.80787325430910	2.74376655900524
32.0935734874856	29.9296688575414	3.07639453413753	2.85478984311915	2.78961199281962
35.4432647031425	33.0535076191775	3.39748597479351	3.15275181559228	3.08077118053310
38.1362916156512	35.5649575749534	3.65563152773911	3.39230213802909	3.31485231752614
38.8047354837341	36.1883317102573	3.71970656952122	3.45176160476757	3.37295426219297
38.8868549616992	36.2649142940948	3.72757829852194	3.45906630244197	3.38009218594786
39.0132984254186	36.3828323252482	3.73969879301332	3.47031371046575	3.39108280383411
40.3288928816827	37.6097230123247	3.86580776607519	3.58733856258672	3.50543585567835
41.2500976609300	38.4688156902906	3.95411171729467	3.66928163596413	3.58550808263161
41.9178395506278	39.0915351780022	4.01811946951590	3.72867866027111	3.64354901049323
42.0370734113561	39.2027297126213	4.02954887290010	3.73928475917424	3.65391296101958
42.2378424210408	39.3899619003695	4.04879399324356	3.75714357847601	3.67136404472085
42.3783058174004	39.5209545721979	4.06225839679386	3.76963808855024	3.68357329200799
43.1946241223468	40.2822327314819	4.14050824243455	3.84225129276765	3.75452866051268
43.3897084882043	40.4641635617048	4.15920844972227	3.85960445122533	3.77148562813799
43.9374551104626	40.9749784458254	4.21171381237578	3.90832764792499	3.81909642308652
47.4736957321868	44.2727885463448	4.55068732445008	4.22288357650240	4.12647070945897
47.8808652018227	44.6525046723740	4.58971737921614	4.25910214427261	4.16186236928922
48.3067508574488	45.0496750482138	4.63054151192442	4.29698555555698	4.19888086250773
48.6908243502612	45.4078524403736	4.66735765502098	4.33114968834574	4.23226499226478
49.5527271229150	46.2116415411101	4.74997709220469	4.40781772539334	4.30718261750764
49.8410776193913	46.4805500463535	4.77761751347536	4.43346709936909	4.33224638933783
49.9171540356757	46.5514970210850	4.78490997255055	4.44023425418886	4.33885904290601
50.0891605160342	46.7119059889580	4.80139800234934	4.45553458692587	4.35381005297779
52.8031608676708	49.2429152527100	5.06155400681166	4.69695054021929	4.58971422651566
53.3566396718885	49.7590758272364	5.11460883938492	4.74618362637838	4.63782326961175
55.1822066641107	51.4615542245322	5.28960226349721	4.90857159197410	4.79650374740848
55.6526303662323	51.9002596754140	5.33469568092235	4.95041675475421	4.83739354117535
55.6604446928560	51.9075471221761	5.33544473903900	4.95111185529467	4.83807277184045
55.9875445001466	52.2125922714608	5.36679955403517	4.98020806073603	4.86650468035393
56.2168320388293	52.4264201304435	5.38877837577730	5.00060366226792	4.88643462888280
58.2627584337030	54.3344002289935	5.58489479705683	5.18259305320104	5.06426901089051
58.4389238049836	54.4986876751593	5.60178148577422	5.19826332790065	5.07958151676137
58.5862665067829	54.6360958162362	5.61590531909048	5.21136976644413	5.09238872154815
59.0559907559532	55.0741489746826	5.66093168903532	5.25315271143532	5.13321771802355
59.1576268046500	55.1689322272446	5.67067421847985	5.26219344850461	5.14205204557071
59.3188750000552	55.3193082859978	5.68613099106946	5.27653680950938	5.15606793257994
59.3469302207481	55.3454718874239	5.68882028111152	5.27903237950603	5.15850652609265
60.6598510773404	56.5698692419064	5.81467297085798	5.39581940939787	5.27262718543901
60.6657935312987	56.5754110267514	5.81524259682636	5.39634800298718	5.27314371068153
61.1994183028136	57.0730562239074	5.86639428085460	5.44381503180999	5.31952701738510
61.5637228775491	57.4127976144812	5.90131543424110	5.47622067838475	5.35119280901524
61.8866490262688	57.7139505135309	5.93227017473576	5.50494562825056	5.37926193810785
61.9590655328426	57.7814843474437	5.93921180574770	5.51138722652011	5.38555646792376
62.2243096024779	58.0288443733913	5.96463731364126	5.53498123596371	5.40861180135267
62.6226168432412	58.4002958050698	6.00281786149672	5.57041149011772	5.44323314559448
63.7018709782694	59.4067813834714	6.10627195405020	5.66641338108895	5.53704321257094
64.9977544718664	60.6152901167365	6.23049149283529	5.78168490224969	5.64968296384217
66.0758779313127	61.6207212551966	6.33383720218920	5.87758621728574	5.74339475113766
66.1410709781088	61.6815186700804	6.34008641382359	5.88338527384477	5.74906140914512
67.6685811248742	63.1060366623523	6.48650899491077	6.01926046561635	5.88183442758858
70.1495357041165	65.4197132317837	6.72432592452132	6.23994651470411	6.09748198240493
70.5765463185089	65.8179327162049	6.76525789243221	6.27792999340108	6.13459825839192
70.5880732198665	65.8286824178771	6.76636282692465	6.27895533515466	6.13560019051642
71.3685362286576	66.5565228192518	6.84117569048059	6.34837913650529	6.20343897357138
71.5259177768509	66.7032929389316	6.85626181776475	6.36237855123085	6.21711876695617
71.9126008895506	67.0639039977287	6.89332811126480	6.39677481511437	6.25072972801130
76.3275334678022	71.1811603745223	7.31653042177146	6.78949221342351	6.63448097567412
77.1534132915151	71.9513553685291	7.39569680617286	6.86295593454176	6.70626744288555
77.3791201129050	72.1618439395735	7.41733239100146	6.88303302385305	6.72588615116234
79.7203895294663	74.3452535985600	7.64175951622829	7.09129378836439	6.92939210369853
80.2180403361392	74.8093503703394	7.68946284294430	7.13556085849492	6.97264851012746
83.9403776481328	78.2807096182425	8.04627503044653	7.46667047316860	7.29619854455745
84.0717050829864	78.4031823194706	8.05886368788984	7.47835232054184	7.30761368308655
86.1248851132886	80.3179269781116	8.25567542109509	7.66098693737214	7.48607855980657
86.4724853751433	80.6420903417258	8.28899534871398	7.69190670071339	7.51629239247793
86.6146619915275	80.7746807199623	8.30262397644212	7.70455360525773	7.52865055490938
87.7619644256082	81.8446264505266	8.41260097662207	7.80660853339082	7.62837546185221
91.0619418622220	84.9221033776042	8.72892700222761	8.10014836223664	7.91521346812165]
pc.throttleLatencyType=0
pc.throttleTargetLatency=0
pc.throttleTargetTransactionIndex=0
pc.throttlePenalty=[]
pc.throttleIndividualTransactions=0

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
# BEGIN 0.1 deprecation

@deprecate restrict(A::AbstractArray, region::Vector{Int}) restrict(A, (region...,))

@deprecate meanfinite(A, region) meanfinite(A; dims=region)

# END 0.1 deprecation


text:
### CONCRETE TYPE: DIRECT PROX EVALUATION
# prox! is computed using a QR factorization of A'.

using LinearAlgebra: QRCompactWY

struct IndAffineDirect{R <: Real, T <: RealOrComplex{R}, M <: AbstractMatrix{T}, V <: AbstractVector{T}, F <: Factorization} <: IndAffine
  A::M
  b::V
  fact::F
  res::V
  function IndAffineDirect{R, T, M, V, F}(A::M, b::V) where {R <: Real, T <: RealOrComplex{R}, M <: AbstractMatrix{T}, V <: AbstractVector{T}, F <: Factorization}
    if size(A, 1) > size(A, 2)
      error("A must be full row rank")
    end
    normrowsinv = 1 ./ vec(sqrt.(sum(abs2.(A); dims=2)))
    A = normrowsinv.*A # normalize rows of A
    b = normrowsinv.*b # and b accordingly
    fact = qr(M(A'))
    new(A, b, fact, similar(b))
  end
end

is_cone(f::IndAffineDirect) = norm(f.b) == 0.0

IndAffineDirect(A::M, b::V) where {R <: Real, T <: RealOrComplex{R}, M <: DenseMatrix{T}, V <: AbstractVector{T}} = IndAffineDirect{R, T, M, V, QRCompactWY{T, M}}(A, b)

IndAffineDirect(A::M, b::V) where {R <: Real, T <: RealOrComplex{R}, I <: Integer, M <: SparseMatrixCSC{T, I}, V <: AbstractVector{T}} = IndAffineDirect{R, T, M, V, SuiteSparse.SPQR.Factorization{T}}(A, b)

IndAffineDirect(a::V, b::T) where {R <: Real, T <: RealOrComplex{R}, V <: AbstractVector{T}} = IndAffineDirect(reshape(a,1,:), [b])

function (f::IndAffineDirect{R, T, M, V, F})(x::V) where {R, T, M, V, F}
  mul!(f.res, f.A, x)
  f.res .= f.b .- f.res
  # the tolerance in the following line should be customizable
  if norm(f.res, Inf) <= 1e-12
    return zero(R)
  end
  return typemax(R)
end

function prox!(y::V, f::IndAffineDirect{R, T, M, V, F}, x::V, gamma::R=one(R)) where {R, T, M, V, F <: QRCompactWY}
  mul!(f.res, f.A, x)
  f.res .= f.b .- f.res
  Rfact = view(f.fact.factors, 1:length(f.b), 1:length(f.b))
  LinearAlgebra.LAPACK.trtrs!('U', 'C', 'N', Rfact, f.res)
  LinearAlgebra.LAPACK.trtrs!('U', 'N', 'N', Rfact, f.res)
  mul!(y, adjoint(f.A), f.res)
  y .+= x
  return zero(R)
end

function prox!(y::V, f::IndAffineDirect{R, T, M, V, F}, x::V, gamma::R=one(R)) where {R, T, M, V, F <: SuiteSparse.SPQR.Factorization}
  mul!(f.res, f.A, x)
  f.res .= f.b .- f.res
  ##############################################################################
  # We need to solve for z: AA'z = res
  # We have: QR=PA'S so A'=P'QRS' and AA'=SR'Q'PP'QRS'=SR'RS'
  # So: z = S'\R\R'\S\res
  # TODO: the following lines should be made more efficient
  temp = f.res[f.fact.pcol]
  temp = LowerTriangular(adjoint(f.fact.R))\temp
  temp = UpperTriangular(f.fact.R)\temp
  RRres = similar(temp)
  RRres[f.fact.pcol] .= temp
  ##############################################################################
  mul!(y, adjoint(f.A), RRres)
  y .+= x
  return zero(R)
end

function prox_naive(f::IndAffineDirect, x::AbstractArray{T,1}, gamma::R=one(R)) where {R <: Real, T <: RealOrComplex{R}}
  y = x + f.A'*((f.A*f.A')\(f.b - f.A*x))
  return y, zero(R)
end


text:
mutable struct robot_urdf_t <: LCMType
    utime::Int64
    robot_name::String
    urdf_xml_string::String

    left_hand::Int8
    right_hand::Int8
end

@lcmtypesetup(robot_urdf_t)


text:
@testset "Szegedy model test" begin
  @testset "Generalized neighbor function" begin
    g = DiGraph(3)
    @test QuantumWalk.outneighbors_looped(g, 3) == [3]

    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)

    @test QuantumWalk.outneighbors_looped(g, 1) == [2, 3]
    @test QuantumWalk.outneighbors_looped(g, 3) == [3]

    @test QuantumWalk.outneighbors_looped(complete_graph(3), 3) == [1, 2]
    @test QuantumWalk.outneighbors_looped(complete_digraph(3), 3) == [1, 2]

    @test_throws BoundsError QuantumWalk.outneighbors_looped(g, 0)
    @test_throws BoundsError QuantumWalk.outneighbors_looped(g, 4)
  end

  @testset "Graph stochastic matrix checker" begin
  g = smallgraph(:bull)

  @test QuantumWalk.graphstochasticcheck(g, QuantumWalk.default_stochastic(g)) == nothing

  g = DiGraph(3)
  add_edge!(g, 1, 3)

  @test QuantumWalk.graphstochasticcheck(g, QuantumWalk.default_stochastic(g)) == nothing

   # negative cases
   g = smallgraph(:bull)
   @test_throws AssertionError QuantumWalk.graphstochasticcheck(g, adjacency_matrix(g))
   @test_throws AssertionError QuantumWalk.graphstochasticcheck(complete_graph(nv(g)+1), QuantumWalk.default_stochastic(g))
   @test_throws AssertionError QuantumWalk.graphstochasticcheck(DiGraph(nv(g)), adjacency_matrix(g))

   g = DiGraph(3)
   add_edge!(g, 1, 2)
   add_edge!(g, 1, 3)

   @test_throws AssertionError QuantumWalk.graphstochasticcheck(g, sparse([0 0 0; -1 1 0; 2 0 1]))
   #dense not allowed
   @test_throws MethodError QuantumWalk.graphstochasticcheck(g, Matrix(QuantumWalk.default_stochastic(g)))
  end

  @testset "Default stochastic matrix" begin
    g = smallgraph(:bull)
    result = [0.0  1. /3  1. /3  0.0  0.0;
              0.5  0.0   1. /3  1.0  0.0;
              0.5  1. /3  0.0   0.0  1.0;
              0.0  1. /3  0.0   0.0  0.0;
              0.0  0.0   1. /3  0.0  0.0]

    @test isapprox(QuantumWalk.default_stochastic(g), sparse(result), atol=1e-5)

    g = DiGraph(3)
    add_edge!(g, 1, 3)
    @test isapprox(QuantumWalk.default_stochastic(g), [0 0 0; 0 1 0; 1 0 1], atol=1e-5)
  end

  @testset "Szegedy walk operators" begin
    g = smallgraph(:bull)
    sz = Szegedy(g, QuantumWalk.default_stochastic(g), true)

    @test sqrtstochastic(sz) == sz.sqrtstochastic

    r1, r2 = QuantumWalk.szegedy_walk_operators(sz)
    @test isapprox(norm(r1*r1'-I, Inf), 0, atol=1e-8)
    @test isapprox(norm(r2*r2'-I, Inf), 0, atol=1e-8)

    q1, q2 = QuantumWalk.szegedy_oracle_operators(sz, [2, 3])

    @test all(typeof(m) == SparseMatrixCSC{Float64,Int} for m = [r1, r2, q1, q2])

    @test isapprox(norm(q1*q1'-I, Inf), 0, atol=1e-8)
    @test isapprox(norm(q2*q2'-I, Inf), 0, atol=1e-8)
    @test isdiag(q1)
    @test isdiag(q2)

    # is this the right place for this functions
    v = [1/3,1/3,1/3]
    @test QuantumWalk.proj(v) == v*v'

  end

  @testset "Szegedy initial state and measurement" begin
    g = smallgraph(:bull)
    n = nv(g)
    qws = QWSearch(Szegedy(g), [1])

    @test measure(qws, initial_state(qws), collect(1:n)) ≈ fill(1/n, (n))
  end

  @testset "Szegedy conversion" begin
    szegedy = Szegedy(smallgraph(:bull))
    qws = QWSearch(szegedy, [1], 0.5)

    function is_equal_mine(q1::QWSearch, q2::QWSearch)
      model(q1) == model(q2) &&
      all(parameters(q1)[:operators] .≈ parameters(q2)[:operators]) &&
      penalty(q1) == penalty(q2) &&
      marked(q1) == marked(q2)
    end

    qss1 = QWSearch(qws, marked = [1])
    qss2 = QWSearch(qws, penalty = 0.5)
    qss3 = QWSearch(qws, marked = [2])
    qss4 = QWSearch(qws, penalty = 1.)
    qss5 = QWSearch(qws, marked = [2], penalty = 1.)
    qss6 = QWSearch(qws, marked = [2], penalty = 2.)

    qss1_ref = QWSearch(szegedy, [1], 0.5)
    qss2_ref = QWSearch(szegedy, [1], 0.5)
    qss3_ref = QWSearch(szegedy, [2], 0.5)
    qss4_ref = QWSearch(szegedy, [1], 1.)
    qss5_ref = QWSearch(szegedy, [2], 1.)

    @test is_equal_mine(qss1, qss1_ref)
    @test is_equal_mine(qss2, qss2_ref)
    @test is_equal_mine(qss3, qss3_ref)
    @test is_equal_mine(qss4, qss4_ref)
    @test is_equal_mine(qss5, qss5_ref)
    @test !is_equal_mine(qss6, qss5_ref)
  end

end


text:
origin(prim::HyperRectangle) = prim.origin
maximum(prim::HyperRectangle) = origin(prim) + widths(prim)
minimum(prim::HyperRectangle) = origin(prim)
length(prim::HyperRectangle{N, T}) where {T, N} = N
widths(prim::HyperRectangle) = prim.widths

width(prim::HyperRectangle) = prim.widths[1]
height(prim::HyperRectangle) = prim.widths[2]

"""
Splits an HyperRectangle into two along an axis at a given location.
"""
split(b::HyperRectangle, axis, value::Integer) = _split(b, axis, value)
split(b::HyperRectangle, axis, value::Number) = _split(b, axis, value)
function _split(b::H, axis, value) where H<:HyperRectangle
    bmin = minimum(b)
    bmax = maximum(b)
    b1max = setindex(bmax, value, axis)
    b2min = setindex(bmin, value, axis)

    return H(bmin, b1max-bmin),
           H(b2min, bmax-b2min)
end

# empty constructor such that update will always include the first point
Rect() = Rect{2, Float32}()

TRect{T}() where {T} = Rect{2, T}()

Rect{N}() where {N} = Rect{N, Float32}()

function Rect{N, T}() where {T,N}
    return Rect(Vec{N,T}(typemax(T)), Vec{N,T}(typemin(T)))
end

# conversion from other HyperRectangles
function Rect{N,T1}(a::Rect{N, T2}) where {N,T1,T2}
    return Rect(Vec{N, T1}(minimum(a)), Vec{N, T1}(widths(a)))
end

function Rect(v1::Vec{N, T1}, v2::Vec{N, T2}) where {N,T1,T2}
    T = promote_type(T1, T2)
    return Rect{N,T}(Vec{N, T}(v1), Vec{N, T}(v2))
end

function TRect{T}(v1::VecTypes{N, T1}, v2::VecTypes{N, T2}) where {N,T,T1,T2}
    if T <: Integer
         Rect{N, T}(round.(T, v1), round.(T, v2))
    else
        return Rect{N, T}(Vec{N, T}(v1), Vec{N, T}(v2))
    end
end

function Rect{N}(v1::VecTypes{N, T1}, v2::VecTypes{N, T2}) where {N,T1,T2}
    T = promote_type(T1, T2)
    return Rect{N,T}(Vec{N, T}(v1), Vec{N, T}(v2))
end

"""
```
HyperRectangle(vals::Number...)
```
HyperRectangle constructor for indidually specified intervals.
e.g. HyperRectangle(0,0,1,2) has origin == Vec(0,0) and
width == Vec(1,2)
"""
@generated function Rect(vals::Number...)
    # Generated so we get goodish codegen on each signature
    n = length(vals)
    @assert iseven(n)
    mid = div(n,2)
    v1 = Expr(:call, :Vec)
    v2 = Expr(:call, :Vec)
    # TODO this can be inbounds
    append!(v1.args, [:(vals[$i]) for i = 1:mid])
    append!(v2.args, [:(vals[$i]) for i = mid+1:length(vals)])
    return Expr(:call, :HyperRectangle, v1, v2)
end

Rect3D(a::Vararg{Number, 6}) = Rect3D(Vec{3}(a[1], a[2], a[3]), Vec{3}(a[4], a[5], a[6]))
Rect3D(args::Vararg{Number, 4}) = Rect3D(Rect{2}(args...))
#=
From different args
=#
function (Rect)(args::Vararg{Number, 4})
    args_prom = promote(args...)
    return Rect2D{typeof(args_prom[1])}(args_prom...)
end

function (Rect2D)(args::Vararg{Number, 4})
    args_prom = promote(args...)
    return Rect2D{typeof(args_prom[1])}(args_prom...)
end

function (Rect{2, T})(args::Vararg{Number, 4}) where {T}
    x, y, w, h = T <: Integer ? round.(T, args) : args
    return Rect2D{T}(Vec{2, T}(x, y), Vec{2, T}(w, h))
end

function TRect{T}(args::Vararg{Number, 4}) where {T}
    x, y, w, h = T <: Integer ? round.(T, args) : args
    return Rect2D{T}(Vec{2, T}(x, y), Vec{2, T}(w, h))
end

Rect(r::SimpleRectangle{T}) where T = HyperRectangle{2, T}(r)
Rect{N}(r::SimpleRectangle{T}) where {N, T} = Rect{N, T}(r)

function Rect{N, T}(r::SimpleRectangle) where {N, T}
    if N === 2
        return Rect(Vec{N, T}(T(r.x), T(r.y)), Vec{N, T}(T(r.w), T(r.h)))
    else
        return Rect(
            Vec{N, T}(r.x, r.y, Vec{N - 2}(zero(T))...),
            Vec{N, T}(r.w, r.h, Vec{N - 2}(zero(T))...)
        )
    end
end

#=
From other types
=#
function Rect2D(xy::NamedTuple{(:x, :y)}, wh::NamedTuple{(:width, :height)})
    return Rect2D(xy.x, xy.y, wh.width, wh.height)
end

function FRect3D(x::Rect2D{T}) where T
    return Rect{3, T}(Vec{3, T}(minimum(x)..., 0), Vec{3, T}(widths(x)..., 0.0))
end

function Rect2D{T}(a::Rect2D) where {T}
    return Rect2D{T}(minimum(a), widths(a))
end

function TRect{T}(a::Rect2D) where {T}
    return Rect2D{T}(minimum(a), widths(a))
end

function Rect{N, T}(a::GeometryPrimitive) where {N, T}
    return Rect{N, T}(Vec{N, T}(minimum(a)), Vec{N, T}(widths(a)))
end

function Rect2D(xy::VecTypes{2}, w::Number, h::Number)
    return Rect2D(xy..., w, h)
end

function Rect2D(x::Number, y::Number, wh::VecTypes{2})
    return Rect2D(x, y, wh...)
end

function TRect{T}(xy::VecTypes{2}, w::Number, h::Number) where {T}
    return Rect2D{T}(xy..., w, h)
end

function TRect{T}(x::Number, y::Number, wh::VecTypes{2}) where {T}
    return Rect2D{T}(x, y, wh...)
end

#=
From limits
=#
function FRect3D(x::Tuple{Tuple{<: Number, <: Number}, Tuple{<: Number, <: Number}})
    FRect3D(Vec3f0(x[1]..., 0), Vec3f0(x[2]..., 0))
end
function FRect3D(x::Tuple{Tuple{<: Number, <: Number, <: Number}, Tuple{<: Number, <: Number, <: Number}})
    FRect3D(Vec3f0(x[1]...), Vec3f0(x[2]...))
end

"""
Transform a `HyperRectangle` using a matrix. Maintains axis-align properties
so a significantly larger HyperRectangle may be generated.
"""
function *(m::Mat{N1,N1,T1}, h::HyperRectangle{N2, T2}) where {N1,N2,T1,T2}

    # TypeVar constants
    T = promote_type(T1, T2)
    D = N1 - N2

    # get all points on the HyperRectangle
    d = decompose(Point, h)
    # make sure our points are sized for the tranform
    pts = (Vec{N1, T}[vcat(pt, ones(Vec{D, T})) for pt in d]...,)::NTuple{2^N2,Vec{N1,T}}

    vmin = Vec{N1, T}(typemax(T))
    vmax = Vec{N1, T}(typemin(T))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * pt
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    HyperRectangle{N2, T}(vmin, vmax - vmin)
end

function -(h::HyperRectangle{N, T}, move::Vec{N}) where {N,T}
    HyperRectangle{N, T}(minimum(h) .- move, widths(h))
end
function +(h::HyperRectangle{N, T}, move::Vec{N}) where {N,T}
    HyperRectangle{N, T}(minimum(h) .+ move, widths(h))
end

function -(h::SimpleRectangle{T}, move::Vec{2}) where T
    SimpleRectangle{T}((minimum(h) .- move)..., widths(h)...)
end
function +(h::SimpleRectangle{T}, move::Vec{2}) where T
    SimpleRectangle{T}((minimum(h) .+ move)..., widths(h)...)
end

function *(m::Mat{N,N,T1}, h::HyperRectangle{N,T2}) where {N,T1,T2}
    # equal dimension case

    # TypeVar constants
    T = promote_type(T1,T2)

    # get all points on the HyperRectangle
    pts = decompose(Point, h)

    # make sure our points are sized for the tranform
    vmin = Vec{N, T}(typemax(T))
    vmax = Vec{N, T}(typemin(T))

    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * Vec(pt)
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    HyperRectangle{N,T}(vmin, vmax-vmin)
end

# fast path. TODO make other versions fast without code duplications like now
function *(m::Mat{4,4,T}, h::HyperRectangle{3,T}) where T
    # equal dimension case

    # get all points on the HyperRectangle
    pts = (
        Vec{4,T}(0.0,0.0,0.0, 1.0),
        Vec{4,T}(1.0,0.0,0.0, 1.0),
        Vec{4,T}(0.0,1.0,0.0, 1.0),
        Vec{4,T}(1.0,1.0,0.0, 1.0),
        Vec{4,T}(0.0,0.0,1.0, 1.0),
        Vec{4,T}(1.0,0.0,1.0, 1.0),
        Vec{4,T}(0.0,1.0,1.0, 1.0),
        Vec{4,T}(1.0,1.0,1.0, 1.0)
    )

    # make sure our points are sized for the tranform
    vmin = Vec{4, T}(typemax(T))
    vmax = Vec{4, T}(typemin(T))
    o, w = origin(h), widths(h)
    _o = Vec{4, T}(o[1], o[2], o[3], T(0))
    _w = Vec{4, T}(w[1], w[2], w[3], T(1))
    # tranform all points, tracking min and max points
    for pt in pts
        pn = m * (_o + (pt .* _w))
        vmin = min.(pn, vmin)
        vmax = max.(pn, vmax)
    end
    _vmin = Vec{3, T}(vmin[1], vmin[2], vmin[3])
    _vmax = Vec{3, T}(vmax[1], vmax[2], vmax[3])
    HyperRectangle{3,T}(_vmin, _vmax - _vmin)
end

function Rect(geometry::AbstractArray{<: Point{N, T}}) where {N,T}
    return Rect{N,T}(geometry)
end

@inline function minmax(p::StaticVector, vmin, vmax)
    isnan(p) && return (vmin, vmax)
    min.(p, vmin), max.(p, vmax)
end

# Annoying special case for view(Vector{Point}, Vector{Face})
@inline function minmax(tup::Tuple, vmin, vmax)
    for p in tup
        isnan(p) && continue
        vmin = min.(p, vmin)
        vmax = max.(p, vmax)
    end
    vmin, vmax
end

"""
Construct a HyperRectangle enclosing all points.
"""
function (t::Type{Rect{N1, T1}})(
        geometry::AbstractArray{PT}
    ) where {N1, T1, PT <: Point}
    N2, T2 = length(PT), eltype(PT)
    @assert N1 >= N2
    vmin = Point{N2, T2}(typemax(T2))
    vmax = Point{N2, T2}(typemin(T2))
    for p in geometry
        vmin, vmax = minmax(p, vmin, vmax)
    end
    o = vmin
    w = vmax - vmin
    if N1 > N2
        z = zero(Vec{N1-N2, T1})
        return Rect{N1, T1}(vcat(o, z),
                                     vcat(w, z))
    else
        return Rect{N1, T1}(o, w)
   end
end

xwidth(a::SimpleRectangle)  = a.w + a.x
width(a::SimpleRectangle)  = a.w
yheight(a::SimpleRectangle) = a.h + a.y
height(a::SimpleRectangle) = a.h
area(a::SimpleRectangle) = a.w*a.h
widths(a::SimpleRectangle{T}) where {T} = Point{2,T}(a.w, a.h)
maximum(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x + widths(a)[1], a.y +widths(a)[2])
minimum(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x, a.y)
origin(a::SimpleRectangle{T}) where {T} = Point{2, T}(a.x, a.y)

SimpleRectangle(val::Vec{2, T}) where {T} = SimpleRectangle{T}(0, 0, val...)
function SimpleRectangle(position::Vec{2,T}, width::Vec{2,T}) where T
    SimpleRectangle{T}(position..., width...)
end

function Base.to_indices(A::AbstractArray{T, 2}, I::Tuple{<: Union{Rect2D, SimpleRectangle}}) where T
    rect = I[1]
    mini = minimum(rect)
    wh = widths(rect)
    return (mini[1] + 1 : (mini[1] + wh[1]), mini[2] + 1 : (mini[2] + wh[2]))
end

AABB(a) = AABB{Float32}(a)

function Rect{T}(a::Pyramid) where T
    w,h = a.width/T(2), a.length
    m = Vec{3,T}(a.middle)
    Rect{T}(m-Vec{3,T}(w,w,0), m+Vec{3,T}(w, w, h))
end

Rect{T}(a::Cube) where T = Rect{T}(origin(a), widths(a))

Rect{T}(a::AbstractMesh) where T = Rect{T}(vertices(a))
Rect{N, T}(a::AbstractMesh) where {N, T} = Rect{N, T}(vertices(a))

function positive_widths(rect::Rect{N, T}) where {N, T}
    mini, maxi = minimum(rect), maximum(rect)
    realmin = min.(mini, maxi)
    realmax = max.(mini, maxi)
    Rect{N, T}(realmin, realmax .- realmin)
end

# set operations

"""
Perform a union between two HyperRectangles.
"""
function union(h1::Rect{N}, h2::Rect{N}) where N
    m = min.(minimum(h1), minimum(h2))
    mm = max.(maximum(h1), maximum(h2))
    Rect{N}(m, mm - m)
end

"""
Perform a difference between two HyperRectangles.
"""
diff(h1::Rect, h2::Rect) = h1

"""
Perform a intersection between two HyperRectangles.
"""
function intersect(h1::Rect{N}, h2::Rect{N}) where N
    m = max.(minimum(h1), minimum(h2))
    mm = min.(maximum(h1), maximum(h2))
    Rect{N}(m, mm - m)
end

function intersect(a::SimpleRectangle, b::SimpleRectangle)
    min_n = max.(minimum(a), minimum(b))
    max_n = min.(maximum(a), maximum(b))
    w = max_n - min_n
    SimpleRectangle(min_n[1], min_n[2], w[1], w[2])
end

function update(b::HyperRectangle{N, T}, v::Vec{N, T2}) where {N, T, T2}
    update(b, Vec{N, T}(v))
end

function update(b::HyperRectangle{N, T}, v::Vec{N, T}) where {N, T}
    m = min.(minimum(b), v)
    maxi = maximum(b)
    mm = if isnan(maxi)
        v-m
    else
        max.(v, maxi) - m
    end
    HyperRectangle{N, T}(m, mm)
end

# Min maximum distance functions between hrectangle and point for a given dimension
@inline function min_dist_dim(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int) where {N, T}
    max(zero(T), max(minimum(rect)[dim] - p[dim], p[dim] - maximum(rect)[dim]))
end

@inline function max_dist_dim(rect::HyperRectangle{N, T}, p::Vec{N, T}, dim::Int) where {N, T}
    max(maximum(rect)[dim] - p[dim], p[dim] - minimum(rect)[dim])
end

@inline function min_dist_dim(rect1::HyperRectangle{N, T},
                              rect2::HyperRectangle{N, T},
                              dim::Int) where {N, T}
    max(zero(T), max(
        minimum(rect1)[dim] - maximum(rect2)[dim],
        minimum(rect2)[dim] - maximum(rect1)[dim]
    ))
end

@inline function max_dist_dim(rect1::HyperRectangle{N, T},
                              rect2::HyperRectangle{N, T},
                              dim::Int) where {N, T}
    max(
        maximum(rect1)[dim] - minimum(rect2)[dim],
        maximum(rect2)[dim] - minimum(rect1)[dim]
    )
end

# Total minimum maximum distance functions
@inline function min_euclideansq(rect::HyperRectangle{N, T},
                                 p::Union{Vec{N, T},
                                 HyperRectangle{N, T}}) where {N, T}
    minimum_dist = T(0.0)
    for dim in 1:length(p)
        d = min_dist_dim(rect, p, dim)
        minimum_dist += d*d
    end
    return minimum_dist
end

@inline function max_euclideansq(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    maximum_dist = T(0.0)
    for dim in 1:length(p)
        d = max_dist_dim(rect, p, dim)
        maximum_dist += d*d
    end
    return maximum_dist
end
function max_euclidean(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    sqrt(max_euclideansq(rect, p))
end

# Functions that return both minimum and maximum for convenience
@inline function minmax_dist_dim(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}, dim::Int) where {N, T}
    minimum_d = min_dist_dim(rect, p, dim)
    maximum_d = max_dist_dim(rect, p, dim)
    return minimum_d, maximum_d
end

@inline function minmax_euclideansq(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    minimum_dist = min_euclideansq(rect, p)
    maximum_dist = max_euclideansq(rect, p)
    return minimum_dist, maximum_dist
end

function minmax_euclidean(rect::HyperRectangle{N, T}, p::Union{Vec{N, T}, HyperRectangle{N, T}}) where {N, T}
    minimumsq, maximumsq = minmax_euclideansq(rect, p)
    return sqrt(minimumsq), sqrt(maximumsq)
end

#
# http://en.wikipedia.org/wiki/Allen%27s_interval_algebra
#

function before(b1::Rect{N}, b2::Rect{N}) where N
    for i = 1:N
        maximum(b1)[i] < minimum(b2)[i] || return false
    end
    true
end

meets(b1::Rect{N}, b2::Rect{N}) where N = maximum(b1) == minimum(b2)

function overlaps(b1::Rect{N}, b2::Rect{N}) where N
    for i = 1:N
        maximum(b2)[i] > maximum(b1)[i] > minimum(b2)[i] &&
        minimum(b1)[i] < minimum(b2)[i] || return false
    end
    true
end

function starts(b1::Rect{N}, b2::Rect{N}) where N
    if minimum(b1) == minimum(b2)
        for i = 1:N
            maximum(b1)[i] < maximum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

function during(b1::Rect{N}, b2::Rect{N}) where N
    for i = 1:N
        maximum(b1)[i] < maximum(b2)[i] && minimum(b1)[i] > minimum(b2)[i] ||
        return false
    end
    true
end

function finishes(b1::Rect{N}, b2::Rect{N}) where N
    if maximum(b1) == maximum(b2)
        for i = 1:N
            minimum(b1)[i] > minimum(b2)[i] || return false
        end
        return true
    else
        return false
    end
end

#
# Containment
#
function isinside(rect::SimpleRectangle, x::Real, y::Real)
    rect.x <= x && rect.y <= y && rect.x + rect.w >= x && rect.y + rect.h >= y
end

# TODO only have point in c and deprecate above methods
in(x::AbstractPoint{2}, c::Circle) = isinside(c, x...)
in(x::AbstractPoint{2}, c::SimpleRectangle) = isinside(c, x...)

"""
Check if HyperRectangles are contained in each other. This does not use
strict inequality, so HyperRectangles may share faces and this will still
return true.
"""
function contains(b1::Rect{N}, b2::Rect{N}) where N
    for i = 1:N
        maximum(b2)[i] <= maximum(b1)[i] &&
            minimum(b2)[i] >= minimum(b1)[i] ||
            return false
    end
    return true
end

"""
Check if a point is contained in a Rect. This will return true if
the point is on a face of the Rect.
"""
function contains(b1::Rect{N, T}, pt::Union{FixedVector, AbstractVector}) where {T, N}
    for i = 1:N
        pt[i] <= maximum(b1)[i] && pt[i] >= minimum(b1)[i] || return false
    end
    return true
end

"""
Check if HyperRectangles are contained in each other. This does not use
strict inequality, so HyperRectangles may share faces and this will still
return true.
"""
in(b1::Rect, b2::Rect) = contains(b2, b1)

"""
Check if a point is contained in a Rect. This will return true if
the point is on a face of the Rect.
"""
in(pt::Union{FixedVector, AbstractVector}, b1::Rect) = contains(b1, pt)

#
# Equality
#
(==)(b1::Rect, b2::Rect) = minimum(b1) == minimum(b2) && widths(b1) == widths(b2)

isequal(b1::Rect, b2::Rect) = b1 == b2

isless(a::SimpleRectangle, b::SimpleRectangle) = isless(area(a), area(b))

centered(R::Type{HyperRectangle{N,T}}) where {N, T} = R(Vec{N,T}(-0.5), Vec{N,T}(1))
centered(::Type{T}) where {T <: HyperRectangle} = centered(HyperRectangle{ndims_or(T, 3), eltype_or(T, Float32)})


text:
# # KenKen

# ## Introduction

# [KenKen](https://en.wikipedia.org/wiki/KenKen) is a number puzzle similar to the Sodoku where the goal is to place the digits 1 to ``N`` on an ``N \times N`` grid, with some of the digits already filled in.

# Our solution must satisfy the following rules:
#   * The digits 1 to ``N`` must appear once each in each row.
#   * The digits 1 to ``N`` must appear once each in each column.
#   * The digits in each shape outlined by thicker lines (called a cage) must equal the number in the corner of the cage when combined using the cage's operation.

# We'll define some structs to hold information about the puzzle and the cages.

struct Cage
    cells::Vector{Tuple{Integer, Integer}}
    operation::Function
    target::Number
end

struct Puzzle
    size::Integer
    cages::Vector{Cage}
end

# Define an easy ``3 \times 3`` puzzle using only the ``+`` operator.

puzzle = Puzzle(
    3,
    [
        Cage([(1, 1), (2, 1)], +, 5),
        Cage([(1, 2), (1, 3)], +, 3),
        Cage([(2, 2)], +, 3),
        Cage([(2, 3), (3, 3)], +, 4),
        Cage([(3, 1), (3, 2)], +, 3),
    ],
)

# ## Solver

# We will use the integer solver from [Cbc.jl](https://github.com/JuliaOpt/Cbc.jl) for this example.

using JuMP
using Cbc

# This isn't an optimization problem, it's actually a feasibility problem.
# We want to find any solution that satsifies the rules of the puzzle, so we set `maxSolutions` to 1 so that the solver stops as soon as it finds a solution.

## Create a JuMP optimiser
kenken = Model(with_optimizer(Cbc.Optimizer, maxSolutions = 1))

# We will define a binary variable for each possible number in each possible cell, i.e., the variable ``x_{i, j, k}`` will be 1 if the cell at position ``(i, j)`` contains the digit ``k`` and 0 otherwise.

@variable(kenken, x[i = 1:(puzzle.size), j = 1:(puzzle.size), k = 1:(puzzle.size)], Bin)

# Now we can add out first constraint: that each cell must contain exactly one digit.
# For each cell at position ``(i, j)``, the sum across all possible digits must be 1.

for i in 1:(puzzle.size), j in 1:(puzzle.size)
    @constraint(kenken, sum(x[i, j, digit] for digit in 1:(puzzle.size)) == 1)
end

# Next we'll add the constraints that each row and column must contain each digit once.
# We loop over all the columns using the variable `index` and over all possible digits using `digit`.
# In each value of `index`, we constrain that column to only contains the digit `digit` once.
# The constraints for the rows are very similar, so we can add them at the same time.

for index in 1:(puzzle.size), digit in 1:(puzzle.size)
    ## The sum over each column must be 1
    @constraint(kenken, sum(x[index, j, digit] for j in 1:(puzzle.size)) == 1)
    ## The sum over each row must be 1
    @constraint(kenken, sum(x[i, index, digit] for i in 1:(puzzle.size)) == 1)
end

# Finally, we have to make sure that the digits in each cage combine to that cage's target number.
# For each cage we'll use the `reduce` function to combine the digits in the cage's cells.

# TODO: explain this more

# We constrain this to be equal to the cage's target value.

for cage in puzzle.cages
    @constraint(
        kenken,
        reduce(
            cage.operation,
            [sum(x[i, j, k] * k for k in 1:(puzzle.size)) for (i, j) in cage.cells],
        ) == cage.target
    )
end

# Now we have all of the constraints in place we can solve the problem!

optimize!(kenken)

# Extract the values of ``x``.
# Integer programs are solved as a series of linear programs so the values might not be precisely 0 and 1.
# We can just round them to the nearest integer to fix this.

x_solution = round.(Int, value.(x))

# To display the solution as a matrix we need to look for the values of `x[i, j, k]` that are equal to 1.

matrix_solution = zeros(Int, puzzle.size, puzzle.size)
for i in 1:(puzzle.size), j in 1:(puzzle.size), k in 1:(puzzle.size)
    if x_solution[i, j, k] == 1
        matrix_solution[i, j] = k
    end
end

matrix_solution


text:
#################### Geweke Diagnostic ####################
"""
    gewekediag(x::Vector{T}; first::Real=0.1, last::Real=0.5,
                etype=:imse, args...) where {T<:Real}

"""
function gewekediag(
    x::Vector{T};
    first::Real = 0.1,
    last::Real = 0.5,
    etype = :imse,
    args...,
) where {T<:Real}
    if !(0.0 < first < 1.0)
        throw(ArgumentError("first is not in (0, 1)"))
    elseif !(0.0 < last < 1.0)
        throw(ArgumentError("last is not in (0, 1)"))
    elseif first + last > 1.0
        throw(ArgumentError("first and last proportions overlap"))
    end
    n = length(x)
    x1 = x[1:round(Int, first * n)]
    x2 = x[round(Int, n - last * n + 1):n]
    z =
        (mean(x1) - mean(x2)) /
        sqrt(mcse(x1, etype; args...)^2 + mcse(x2, etype; args...)^2)
    [round(z, digits = 3), round(1.0 - erf(abs(z) / sqrt(2.0)), digits = 4)]
end
"""
    gewekediag(c::AbstractChains; first::Real=0.1, last::Real=0.5,
                etype=:imse, args...)

Compute the convergence diagnostic of Geweke [37] for MCMC sampler output. The diagnostic is designed to asses convergence of posterior means estimated with autocorrelated samples. It computes a normal-based test statistic comparing the sample means in two windows containing proportions of the first and last iterations. Users should ensure that there is sufficient separation between the two windows to assume that their samples are independent. A non-significant test p-value indicates convergence. Significant p-values indicate non-convergence and the possible need to discard initial samples as a burn-in sequence or to simulate additional samples.

Returns a `ChainSummary` type object with parameters contained in the rows of the `value` field, and test Z-scores and p-values in the first and second columns. Results are chain-specific.

* `c` : sampler output on which to perform calculations.

* `first` : proportion of iterations to include in the first window.

* `last` : proportion of iterations to include in the last window.

* `etype` : method for computing Monte Carlo standard errors. See `mcse()` for options.

* `args...` : additional arguments to be passed to the `etype` method.

* `x` : vector on which to perform calculations.
"""
function gewekediag(
    c::AbstractChains;
    first::Real = 0.1,
    last::Real = 0.5,
    etype = :imse,
    args...,
)
    _, p, m = size(c.value)
    vals = Array{Float64}(undef, p, 2, m)
    for j = 1:p, k = 1:m
        vals[j, :, k] =
            gewekediag(c.value[:, j, k], first = first, last = last, etype = etype; args...)
    end
    hdr =
        header(c) *
        "\nGeweke Diagnostic:\nFirst Window Fraction = $first\n" *
        "Second Window Fraction = $last\n"
    ChainSummary(vals, c.names, ["Z-score", "p-value"], hdr)
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
include("Tree.jl")

mutable struct NodePool
    population::Array{Tree,1}
    nodes::Array{Tree,1}
    numNodes::Array{UInt,1}
    heights::Array{UInt,1}
    semantics::Array{Array{<:Real,1}}
    hasZerosInSemantics::Array{Bool,1}
    toleranceValue::Real
    NodePool() = new([],[],[],[],[],[], 0.);
end

function checkIntegrity(sol::NodePool)
    @assert (length(sol.nodes)==length(sol.numNodes));
    @assert (length(sol.nodes)==length(sol.heights));
    @assert (length(sol.nodes)==length(sol.semantics));
    for numNode = 1:length(sol.nodes)
        @assert (numNodes(sol.nodes[numNode])==sol.numNodes[numNode]);
        @assert (height(sol.nodes[numNode])==sol.heights[numNode]);
    end;
end;

function removeNodesFromPool!(sol::NodePool, indexNodes::Array{Int,1})
    numberOfNodes = length(sol.nodes);
    indexRemaining = setdiff(1:numberOfNodes, indexNodes);

    sol.nodes = sol.nodes[indexRemaining];
    sol.numNodes = sol.numNodes[indexRemaining];
    sol.heights = sol.heights[indexRemaining];
    sol.semantics = sol.semantics[indexRemaining];
    return indexRemaining;
end;

function addNodeToPool!(sol::NodePool, node::Tree, numNodes::UInt, height::UInt)
    semantic = node.semantic;
    @assert (semantic!=nothing);
    if (sol.toleranceValue!=0.)
        add = all( sqrt.( (x->sum(x.^2)).([semantic] .- sol.semantics)) .> sol.toleranceValue );
    else
        add = all( [semantic] .!= sol.semantics );
    end;

    if (add)
        push!(sol.nodes, node);
        push!(sol.numNodes, numNodes);
        push!(sol.heights, height);
        push!(sol.semantics, semantic);
        # push!(sol.hasZerosInSemantics, any(semantic.==0));
        push!(sol.hasZerosInSemantics, any0(semantic));
        return true;
    else
        return false;
    end;
end;


text:
# Defines the configuration of the Master.
# used in the creation of the master side algorithm.
# Allows the definition of algorithm parameters such as n, m, k, distance function

struct MasterConfiguration
    n::Int64
    k::Int64
    m::Int64

    # Distance function type, used by the Slave to calculate distance betweeen 2 points
    dist::String
    # report::Bool

    # for future implementations
    #port::Int64
    #conntype::String

    function MasterConfiguration(n::T, k::T, m::T, dist::S) where {T<:Int64, S<:String}
        new(n, k, m, dist)
    end
    MasterConfiguration() = new()
end


text:
module EPExperiments

using Reexport
@reexport using DifferentialEquations
using DiffEqSensitivity
@reexport using EPT
using StatsFuns: logsumexp, logistic
using DataFrames
@reexport using Parameters
using Logging
using LoggingExtras
using LaTeXStrings
import LinearAlgebra

export sir_ode!,
    base_reproduction_rate,
    predict,
    neg_bin_params2,
    bayes_sir,
    geomspace,
    SIRParameters,
    plot_predictive,
    process_results,
    COST_FUNCTIONS,
    plot_is_samples,
    MHTransitionConfig,
    MHFixedTuple,
    MHThreeLevel,
    MHFunction,
    AnISConfig,
    ISConfig,
    AnISHMCConfig,
    AnISHMCGranularConfig,
    AnISMHConfig,
    SIRExperimentConfig,
    convergence_plot,
    ISGroundTruthConfig,
    get_anis_alg,
    radon_hierarchical,
    preprocess_radon_df,
    POST_PRED_COUNTY_IDX,
    POST_PRED_FLOOR,
    RADON_COST_FNS,
    radon_target_f,
    plot_intermediate_ess,
    plot_intermediate_weights,
    plot_intermediate_acceptance,
    get_median_and_quantiles,
    check_diagnostics,
    ess_batch,
    make_ess_plot,
    relative_squared_error,
    plot_joint_dist,
    plot_sir_samples,
    post_pred,
    POST_PRED_Y_OBSERVED,
    POST_PRED_GROUND_TRUTH,
    POST_PRED_DIMENSION,
    compute_ess,
    banana,
    banana_f

# Constants for the loss function for radon problem.
const POST_PRED_COUNTY_IDX = 10
const POST_PRED_FLOOR = 0

include("configs.jl")
include("cost_functions.jl")
include("check_diagnostics.jl")
include("post_pred_model.jl")

include("radon_model.jl")

function sir_ode!(du,u,p,t)
    (S,I,R,C) = u
    (β,γ) = p
    N = S+I+R
    infection = β*I/N*S
    recovery = γ*I
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
        du[4] = infection
    end
    nothing
end

function base_reproduction_rate(β, γ)
    return β / γ
end

function predict(y,chain)
    # Length of data
    l = length(y)
    # Length of chain
    m = length(chain)
    # Choose random
    idx = sample(1:m)
    i₀ = Array(chain[:i₀])[idx]
    β = Array(chain[:β])[idx]
    I = i₀*10.0
    u0=[1000.0-I,I,0.0,0.0]
    p=[β,0.25]
    tspan = (0.0,float(l))
    prob = ODEProblem(sir_ode!,
            u0,
            tspan,
            p)
    sol = solve(prob,
                Tsit5(),
                saveat = 1.0)
    out = Array(sol)
    sol_X = [0.0; out[4,2:end] - out[4,1:(end-1)]]
    return hcat(sol.t, out', sol_X)
end

neg_bin_r(mean, var) = mean^2 / (var - mean)
neg_bin_p(r, mean) = r / (r + mean)

function neg_bin_params(mean, var)
    r = neg_bin_r(mean, var)
    return r, neg_bin_p(r, mean)
end

neg_bin_params2(mean, phi) = neg_bin_params(mean, mean + mean^2 / phi)

@expectation function bayes_sir(y, total_population, γ, cost_fn)
    # Calculate number of timepoints
    l = length(y)
    i₀ ~ truncated(Normal(10, 10), 0, total_population/10)
    β ~ truncated(Normal(2, 1.5), 0, Inf)

    I = i₀ * 10
    u0 = [total_population-I, I, 0.0, 0.0]
    p = [β, γ]

    tspan = (0.0, float(l))
    prob = ODEProblem(sir_ode!, u0, tspan, p)
    sol = solve(prob, Tsit5(), saveat = 1.0)
    sol_C = Array(sol)[4,:] # Cumulative cases
    sol_X = sol_C[2:end] - sol_C[1:(end-1)]

    l = length(y)
    if any(sol_X .< 0) || length(sol_X) != l
        # Check if we have negative cumulative cases
        @warn "Negative cumulative cases" β i₀
        Turing.acclogp!(_varinfo, -Inf)
        return l
    end

    phi = 0.5
    variance = sol_X .+ sol_X.^2 ./ phi
    rs = neg_bin_r.(sol_X, variance)
    ps = neg_bin_p.(rs, sol_X)
    if !all(rs .> 0) || !all(0 .< ps .<= 1)
        Turing.acclogp!(_varinfo, -Inf)
        return l
    end

    y ~ arraydist(NegativeBinomial.(rs, ps))
    #y ~ arraydist([Poisson(x) for x in sol_X])
    return cost_fn(β, γ)
end

function geomspace(start, stop, length)
    logstart = log10(start)
    logstop = log10(stop)
    points = 10 .^ range(logstart, logstop; length=length)
    points[1] = start
    points[end] = stop
    return points
end

function relative_squared_error(mu_true, mu_hat)
    return (mu_true - mu_hat)^2 / mu_true^2
end

function compute_ess(log_weights)
    denominator = logsumexp(2 * log_weights)
    numerator = 2 * logsumexp(log_weights)
    return exp(numerator - denominator)
end

include("plots.jl")

function est_exp(chn, f)
    log_weights = Array(chn[:log_weight])
    i₀ = Array(chn[:i₀])
    β = Array(chn[:β])

    normalisation = exp(logsumexp(log_weights))
    weighted_sum = sum(zip(log_weights, i₀, β)) do (lw, i, b)
        exp(lw) * f(i, b)
    end

    return weighted_sum / normalisation
end

function process_results(
    result_folder,
    logger;
    ode_tabi,
    R0_estimate,
    tabi,
    X,
    Y,
    obstimes,
    i_0_true,
    β_true,
    cost_fn
)
    betas = tabi.Z2_alg.inference_algorithm.betas

    savefig(plot_intermediate_weights(
        ode_tabi[:Z2_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "intermediate_weights.png")
    )
    savefig(plot_intermediate_weights(
        ode_tabi[:Z2_info].info[:intermediate_log_weights], betas; ixs=1:14),
        joinpath(result_folder, "intermediate_weights_first40.png")
    )
    savefig(plot_intermediate_ess(
        ode_tabi[:Z2_info].info[:intermediate_log_weights], betas),
        joinpath(result_folder, "intermediate_ess.png")
    )
    savefig(plot_intermediate_ess(
        ode_tabi[:Z2_info].info[:intermediate_log_weights], betas; ixs=1:14),
        joinpath(result_folder, "intermediate_ess_first40.png")
    )
    savefig(
        plot_predictive(obstimes, X, Y, ode_tabi[:Z2_info]),
        joinpath(result_folder, "anis_post_pred.png")
    )
    savefig(
        plot_joint_dist(ode_tabi[:Z2_info], β_true, i_0_true),
        joinpath(result_folder, "joint_samples.png")
    )

    anis_estimate = est_exp(ode_tabi[:Z2_info], cost_fn)

    with_logger(logger) do
        @info "Expectation estimate TABI: $(R0_estimate)"
        @info "Expectation estimate AnIS: $(anis_estimate)"

        @info "Estimation algorithm: $(tabi)"

        for (k, est_results) in pairs(ode_tabi)
            if isnothing(est_results)
                continue
            end
            msg = string([
                "$(string(k)):\n",
                "ESS: $(est_results.info[:ess])\n",
                "Log evidence: $(est_results.logevidence)\n",
                "Log weights: $(get(est_results, :log_weight)[:log_weight])\n"
            ]...)

            @info "$(msg)"
        end
    end
end

end


text:
module day15

# Part 1
function playgame(turns)
    numbers = [0,8,15,2,12,1,4]
    numturns = Dict(num => idx for (idx, num) in enumerate(numbers) if idx < lastindex(numbers))
    lastnum = numbers[end]
    for i in (lastindex(numbers)+1):turns
        if lastnum ∈ keys(numturns)
            num = (i - 1) - numturns[lastnum]
        else
            num = 0
        end
        numturns[lastnum] = i - 1
        lastnum = num
    end
    lastnum
end
println("Part 1: 2020th number = $(playgame(2020))")

# Part 2
println("Part 2: 30,000,000th number = $(playgame(30000000))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
"""
$(SIGNATURES)
Computes correlation matrix (R) * n on GPU with cuBLAS library. n is not removed because of performance choice.

# Arguments
- `a` : standardized phenotype matrix
- `b` : standardized genotype matrix

# Output:
returns correlation matrix R on GPU device (of type CuArray).

"""
function calculate_nr(a::CuArray,b::CuArray)
    return CUDA.CUBLAS.gemm('T', 'N', a,b);
end

"""
$(SIGNATURES)
Computes the number of blocks for chuncking large phenotype matrix, to fit into GPU memory.

# Arguments
- `n` : number of individuals
- `m` : number of phenotypes
- `p` : number of genotype markers
- datatype: datatype of phenotype. (Float32, or Float64)

# Output:
returns
- `num_block` : number of blocks to disect phenotype
- `block_size` : size of each block.
"""
function get_pheno_block_size(n::Int, m::Int, p::Int, datatype::DataType)
    total_data_size = (n*m + n*p + m*p) * sizeof(datatype) # get the number of bytes in total
    gpu_mem = CUDA.available_memory() * 0.89 # can not use all of gpu memory, need to leave some for intermediate result.
    #if m is too big for gpu memory, I need to seperate m into several blocks to process
    block_size = Int(ceil((gpu_mem - (n*p))/((n+p) * sizeof(datatype))))
    num_block = Int(ceil(m/block_size))
    return (num_block, block_size)
end

"""
$(SIGNATURES)
Sets up GPU environment and calls GPU kernels.

# Arguments
- `d_nr` : Output matrix (correlation matrix) in CuArray format
- `n` : number of individuals
- `m` : number of phenotypes
- `p` : number of genotype markers
- `export_matrix` : boolean value to indicate user desire output

# Output:
returns correlation matrix R, or maximum of R if `export_matrix` is set to false

"""
function gpu_square_lod(d_nr::CuArray{<:Real,2},n,m,p, export_matrix::Bool)
    #Get total number of threads
    ndrange = prod(size(d_nr))
    #Get maximum number of threads per block
    dev = CUDA.device()
    threads = CUDA.warpsize(dev)
    blocks = min(Int(ceil(ndrange/threads)), attribute(dev, CUDA.CU_DEVICE_ATTRIBUTE_MAX_GRID_DIM_X))
    @cuda blocks=blocks threads=threads lod_kernel(d_nr, ndrange,n)
    if !export_matrix
        @cuda blocks=blocks threads=threads reduce_kernel(d_nr,m,p)
    end
    return
end

################
## GPU kernel ##
################
"""
$(SIGNATURES)
GPU kernel that computes LOD scores.

# Arguments
- `nr` : correlation matrix. will be updated with output results. (Type CuArray)
- `MAX` : number of total GPU threads, also called ndrange.
- `n` : number of individuals

# Output:
Results written back to input matrix `nr`. Returns the LOD scores.

"""
function lod_kernel(nr, MAX,n)
    # TODO: n, 1, 2, do they need to be converted to correct data type for better performance?
    tid = (blockIdx().x-1) * blockDim().x + threadIdx().x
    if(tid < MAX+1)
        r_square = (nr[tid]/n)^2
        nr[tid] = (-n/2.0) * CUDA.log10(1.0-r_square)
    end
    return
end

"""
$(SIGNATURES)
Computes the maximum LOD score of each row. Results written back to the first two columns.

# Arguments
- `input` : LOD score
- `rows` : number of rows for `input`
- `cols` : number of cols for `input`

# Output:
Output is written to the `input` matrix, first column of `input` matrix is the index of maximum, and second column is the maximum value of that row.

"""
function reduce_kernel(input, rows, cols)
    tid = (blockIdx().x-1) * blockDim().x + threadIdx().x

    # Trying for simplest kernel
    if(tid <= rows)
        temp_max = input[tid, 1]
        max_idx = 1
        for i in 1:cols
            if temp_max < input[tid,i]
                temp_max = input[tid,i]
                max_idx = i
            end
        end
        input[tid,1] = max_idx
        input[tid,2] = temp_max
    end

    return
end

"""
$(SIGNATURES)

# Arguments:
- `Y` : a matrix of phenotypes.
- `G` : a matrix of genotypes.
- `X` : a matrix of covariates. Default is `nothing`. If `nothing`, scan is run without covariates.
- `maf_threshold`: a floating point number to indicate the maf_threshold. Default is 0.05. Set to 0 if no maf filtering should be done.
- `export_matrix` : a boolean value that determines whether the result should be the maximum value of LOD score of each phenotype and its corresponding index, or the whole LOD score matrix.
- `lod_or_pval`: a string value of either `lod` or `pval` to indicate the desired output.
- `timing_file`: a string that indicates the file location for the timing outputs. Default is nothing.
# Output:
returns LOD score, in vector or matrix format depending on value of `export_matrix`. P value for GPU is not currently implementated.
"""
function gpurun(pheno::Array{<:Real,2}, geno::Array{<:Real,2}, X::Union{Array{<:Real, 2}, Nothing}=nothing; maf_threshold=0.05, export_matrix=false, timing_file="")
    start = time_ns()
    Y = pheno
    G = geno
    if maf_threshold > 0
        G = filter_maf(geno, maf_threshold=maf_threshold)
    end
    m = size(Y,2)
    p = size(G,2)
    n = size(G,1)
    (num_block, block_size) = get_pheno_block_size(n,m,p, typeof(Y[1,1]))
    println("Seperated into $num_block blocks")

    data_transfer_time = 0.0
    result_reorg_time =0.0
    compute_time = 0.0
    pval_time = 0.0

    # Output array:
    if export_matrix
        result_dim = p
    else
        result_dim = 2
    end
    result_reorg_time += @elapsed lod = convert(Array{typeof(Y[1,1]), 2},zeros(0,result_dim))
    # Transfer Genotype array to GPU
    data_transfer_time += CUDA.@elapsed d_g = CuArray(G);

    if !isnothing(X)
        # Calculate px matrix for covar (X) on GPU
        data_transfer_time += CUDA.@elapsed d_px = CuArray(calculate_px(X))
        # Calculate g hat on GPU
        compute_time += CUDA.@elapsed d_g_hat = CUDA.CUBLAS.gemm('N', 'N', d_px, d_g)
        # Calculate tilda on GPU
        compute_time += CUDA.@elapsed d_g = d_g .- d_g_hat
    end
    # Standardizing matrix on GPU
    compute_time += CUDA.@elapsed d_g_std = get_standardized_matrix_gpu(d_g)

    for i = 1:num_block
        # i = 1
        begining = block_size * (i-1) +1
        ending = i * block_size
        if (i == num_block)
            ending = size(Y)[2]
        end
        # println("processing $begining to $ending...")

        data_transfer_time += @elapsed y_block = Y[:, begining : ending]
        # Transfer phenotype array to GPU
        data_transfer_time += CUDA.@elapsed d_y = CuArray(y_block)

        compute_time += CUDA.@elapsed begin
            if !isnothing(X)
                # Calulate y hat on GPU
                d_y_hat = CUDA.CUBLAS.gemm('N', 'N', d_px, d_y)
                # Calculate y tilda on GPU
                d_y = d_y .- d_y_hat
            end
            d_y_std = get_standardized_matrix_gpu(d_y)

            # calculate correlation matrix
            d_nr = calculate_nr(d_y_std,d_g_std);
            actual_block_size = ending - begining + 1 #it is only different from block size at the last loop since we are calculating the left over block not a whole block.
            gpu_square_lod(d_nr,n,actual_block_size,p, export_matrix)
        end

        if export_matrix
            data_transfer_time += CUDA.@elapsed res = collect(d_nr)
        else
            data_transfer_time += CUDA.@elapsed res = collect(d_nr[:, 1:2])
        end
        result_reorg_time += @elapsed lod = vcat(lod, res)
    end
    stop = time_ns()
    elapsed_total = (stop - start) * 1e-9

    if timing_file != ""
        open(timing_file, "a") do io
            write(io, "$(now()),GPU,$data_transfer_time,$(compute_time),$(result_reorg_time),$(pval_time),$(elapsed_total)\n")
        end
    end
    return lod

end