text:
__precompile__()

"""
    printsummary()

    printsummary(principal, rate, amortization, frequency, compounding, startdate, termString, paymentSize, firstPaymentDate, numberOfPayments, accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, table)

Print out a seven-line summary of the given and calculated mortgage parameters, eg:

julia> printsummary(mortgage(startdate=Date(2016,6,29))...)

Principal: 100000; Annual Interest Rate: 10.0%; Payment frequency: monthly
Compounding: monthly (American) compounding
Your mortgage starts on 2016-06-29 and is amortized over (ie would be fully paid off in) 25.0 years.
Your first payment of 908.70 will be on 2016-07-30.
During the term of 5 years, you will make 60 payments, with a final payment of 908.70 on 2021-06-29.
At the end of the term, the balance remaining will be 94163.77. You will have paid a total of 54522.04
of which 48685.81, or 89.3%, will be interest. This represents 48.7% of the principal amount.

# usage
`printsummary(mortgage()...)`

where mortgage() takes arguments as defined for that function.
"""
function printsummary(principal, rate, amortization, frequency, compounding, startdate, termString, paymentSize, firstPaymentDate, numberOfPayments, accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, table)
	# print out basic mortgage data
	println()
	println("Principal: $(principal); Annual Interest Rate: $(rate)%; Payment frequency: $(frequencyDescriptors[frequency])")
	println("Compounding: $(compoundingDescriptors[compounding])")
	println(@sprintf "Your mortgage starts on %s and is amortized over (ie would be fully paid off in) %0.1f years." startdate amortization)
	println("Your first payment of $(@sprintf "%0.2f" paymentSize) will be on $(firstPaymentDate).")
	println(@sprintf "During the term of %s, you will make %d payments, with a final payment of %0.2f on %s." termString numberOfPayments finalPayment finalPaymentDate)
	println(@sprintf "At the end of the term, the balance remaining will be %0.2f. You will have paid a total of %0.2f" balance accumulatedTotal)
	println(@sprintf "of which %0.2f, or %0.1f%%, will be interest. This represents %0.1f%% of the principal amount." accumulatedInterest accumulatedInterest / accumulatedTotal * 100 accumulatedInterest / principal * 100)
end

"""
    function printtable

    printtable(principal, rate, amortization, frequency, compounding, startdate, termString, paymentSize, firstPaymentDate, numberOfPayments, accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, table)

Print out a header, then a table with a row for each payment, and then a summary footer, eg:

julia> printtable(mortgage(startdate=Date(2016,6,29), term="1/4")...)
                       Mortgage Schedule of Payments

Mortgage Principal: 100000.00 at 10.00% annual interest, amortized over 25.0 years,
using monthly (American) compounding, to be repaid with 3 payments of 908.70 each,
for a term of 3 months. The mortgage start date is 2016-06-29.

                                                   Accumulated   Accumulated   Accumulated
Payment         Date    Principal     Interest       Principal      Interest         Total         Balance

      1   2016-07-30        75.37       833.33           75.37        833.33        908.70        99924.63
      2   2016-08-30        76.00       832.71          151.36       1666.04       1817.40        99848.64
      3   2016-09-29        76.63       832.07          227.99       2498.11       2726.10        99772.01

At the end of the term, the outstanding balance will be 99772.01.
The last payment of the term will be on 2016-09-29 in the amount of 908.70.
Over the term, you will have paid a total of 2726.10 of which 2498.11, or 91.64%, is interest. This
represents 2.50% of the principal amount.

# usage
`printtable(mortgage()...)`

where mortgage() takes arguments as defined for that function.
"""
function printtable(principal, rate, amortization, frequency, compounding, startdate, termString, paymentSize, firstPaymentDate, numberOfPayments, accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, table)
	# used when a tabular listing for each payment is desired
	printtableheader(principal, rate, amortization, compounding, numberOfPayments, paymentSize, termString, startdate)
	printtablecolumnheadings()
	printtablerows(numberOfPayments, table)
	printtablefooter(accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, principal) # use "splat ..." to allow the tuple returned by one function to be used a multiple arguments in another function
end

function printtableheader(principal, rate, amortization, compounding, numberOfPayments, paymentSize, termString, startdate)
	# println()
	print_with_color(:red, "                       Mortgage Schedule of Payments")
	println()
	println()
	s1 = @sprintf "Mortgage Principal: %0.2f at %0.2f%% annual interest, amortized over %0.1f years,\nusing %s, to be repaid with %d payments of %0.2f each,\nfor a term of %s. The mortgage start date is %s." principal rate amortization compoundingDescriptors[compounding] numberOfPayments paymentSize termString startdate
	println(s1)
end

function printtablecolumnheadings()
	println()
	println("                                                   Accumulated   Accumulated   Accumulated")
	println("Payment         Date    Principal     Interest       Principal      Interest         Total         Balance")
	println()
end

function printtablerows(numberOfPayments, table)
	for i = 0:numberOfPayments - 1
		j = i * 8 + 1
		s1 = @sprintf "%7d %12s %12.2f %12.2f" table[j] table[j + 1] table[j + 2] table[j + 3]
		s2 = @sprintf "%16.2f %13.2f %13.2f %15.2f" table[j + 4] table[j + 5] table[j + 6] table[j + 7]
		println(s1, s2)
	end
end

function printtablefooter(accumulatedPrincipal, accumulatedInterest, accumulatedTotal, balance, finalPayment, finalPaymentDate, principal)
	s1 = @sprintf "At the end of the term, the outstanding balance will be %0.2f.\n" balance
	s2 = @sprintf "The last payment of the term will be on %s in the amount of %0.2f.\n" finalPaymentDate finalPayment
	s3 = @sprintf "Over the term, you will have paid a total of %0.2f of which %0.2f, or %0.2f%%, is interest. This\nrepresents %0.2f%% of the principal amount." accumulatedTotal accumulatedInterest accumulatedInterest / accumulatedTotal * 100 accumulatedInterest / principal * 100
	println("\n", s1, s2, s3, "\n\n\n")
end


text:
# This file is a part of SimilaritySearch.jl
using Intersections

export maxlength, maxlength, getpair, getdist, getid, initialstate, idview, distview, reuse!

export KnnResult

"""
    KnnResult(ksearch::Integer)

Creates a priority queue with fixed capacity (`ksearch`) representing a knn result set.
It starts with zero items and grows with [`push!(res, id, dist)`](@ref) calls until `ksearch`
size is reached. After this only the smallest items based on distance are preserved.
"""
struct KnnResult # <: AbstractVector{Tuple{IdType,DistType}}
    id::Vector{Int32}
    dist::Vector{Float32}
    k::Int  # number of neighbors
end

function KnnResult(k::Integer)
    @assert k > 0
    res = KnnResult(Vector{Int32}(undef, 0), Vector{Float32}(undef, 0), k)
    sizehint!(res.id, k)
    sizehint!(res.dist, k)
    res
end

"""
    _shifted_fixorder!(res, shift=0)

Sorts the result in place; the possible element out of order is on the last entry always.
It implements a kind of insertion sort that it is efficient due to the expected
distribution of the items being inserted (it is expected just a few elements smaller than the current ones)
"""
function _shifted_fixorder!(res, shift=0)
    sp = shift + 1
    pos = N = lastindex(res.id)
    id = res.id
    dist = res.dist
    id_, dist_ = res.id[end], res.dist[end]

    #pos = doublingsearch(dist, dist_, sp, N)
    #pos = binarysearch(dist, dist_, sp, N)
    #if N > 16
    #    pos = doublingsearchrev(dist, dist_, sp, N)::Int
    #else
        @inbounds while pos > sp && dist_ < dist[pos-1]
            pos -= 1
        end
    #end

    @inbounds if pos < N
        while N > pos
            id[N] = id[N-1]
            dist[N] = dist[N-1]
            N -= 1
        end

        dist[N] = dist_
        id[N] = id_
    end

    nothing
end

"""
    push!(res::KnnResult, item::Pair)
    push!(res::KnnResult, id::Integer, dist::Real)

Appends an item into the result set
"""
@inline function Base.push!(res::KnnResult, id::Integer, dist::Real)
    if length(res) < maxlength(res)
        k = res.k
        push!(res.id, id)
        push!(res.dist, dist)

        _shifted_fixorder!(res)
        return true
    end

    dist >= last(res.dist) && return false

    @inbounds res.id[end], res.dist[end] = id, dist
    _shifted_fixorder!(res)
    #_shifted_fixorder!(res.shift, res.id, res.dist)
    true
end

#@inline Base.push!(res::KnnResult, id::Integer, dist::Real) = push!(res, convert(Int32, id), convert(Float32, dist))
@inline Base.push!(res::KnnResult, p::Pair) = push!(res, p.first, p.second)

"""
    popfirst!(p::KnnResult)

Removes and returns the nearest neeighboor pair from the pool, an O(length(p.pool)) operation
"""
@inline function Base.popfirst!(res::KnnResult)
    popfirst!(res.id) => popfirst!(res.dist)
end

"""
    pop!(res::KnnResult)

Removes and returns the last item in the pool, it is an O(1) operation
"""
@inline function Base.pop!(res::KnnResult)
    pop!(res.id) => pop!(res.dist)
end

"""
    maxlength(res::KnnResult)

The maximum allowed cardinality (the k of knn)
"""
@inline maxlength(res::KnnResult) = res.k
@inline Base.length(res::KnnResult) = length(res.id)

"""
    reuse!(res::KnnResult)
    reuse!(res::KnnResult, k::Integer)

Returns a result set and a new initial state; reuse the memory buffers
"""
@inline function reuse!(res::KnnResult, k::Integer=res.k)
    @assert k > 0
    empty!(res.id)
    empty!(res.dist)
    if k > res.k
        sizehint!(res.id, k)
        sizehint!(res.dist, k)
    end
    KnnResult(res.id, res.dist, k)
end

"""
    getindex(res::KnnResult, i)

Access the i-th item in `res`
"""
@inline function getpair(res::KnnResult, i)
    @inbounds res.id[i] => res.dist[i]
end

@inline getid(res::KnnResult, i) = @inbounds res.id[i]
@inline getdist(res::KnnResult, i) = @inbounds res.dist[i]

@inline Base.last(res::KnnResult) = last(res.id) => last(res.dist)
@inline Base.first(res::KnnResult) = @inbounds res.id[1] => res.dist[1]
@inline Base.maximum(res::KnnResult) = last(res.dist)
@inline Base.minimum(res::KnnResult) = @inbounds res.dist[1]
@inline Base.argmax(res::KnnResult) = last(res.id)
@inline Base.argmin(res::KnnResult) = @inbounds res.id[1]

@inline idview(res::KnnResult) = res.id
@inline distview(res::KnnResult) = res.dist

@inline Base.eachindex(res::KnnResult) = 1:length(res)
Base.eltype(res::KnnResult) = Pair{Int32,Float32}

##### iterator interface
### KnnResult
"""
    Base.iterate(res::KnnResult, state::Int=1)

Support for iteration
"""
function Base.iterate(res::KnnResult, i::Int=1)
    n = length(res)
    (n == 0 || i > n) && return nothing
    @inbounds res.id[i] => res.dist[i], i+1
end

ext: jl

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
"""
    read(input;
         delim=',',
         quotes=missing,
         escape=missing,
         comment=missing,
         encodings=Dict{String, Any}(),
         header=0,
         skiprows=Vector{Int}(),
         types=Dict{Int,DataType}(),
         allowmissing=Dict{Int,Bool}(),
         coltypes=Vector,
         colparsers=Dict{Int,Function}(),
         typeparsers=Dict{DataType, Function}(),
         typedetectrows=1,
         skipmalformed=false,
         trimwhitespace=false)

Take an input file or IO source and user-defined parsing rules and return:
1. a `Vector{Any}` containing the parsed columns
2. a `Vector{String}` containing the header (column names)

# Arguments
- `input`
    - the path to a local file, or an open IO source from which to read data
- `delim`
    - a `Char` or `String` that separates fields in the dataset
    - default: `delim=','`
        - for CSV files
    - frequently used:
        - `delim='\\t'`
        - `delim=' '`
        - `delim='|'`
- `quotes`
    - a `Char` used for quoting fields in the dataset
    - default: `quotes=missing`
        - by default, the parser does not check for quotes
    - frequently used:
        - `quotes='"'`
- `escape`
    - a `Char` used for escaping other reserved parsing characters
    - default: `escape=missing`
        - by default, the parser does not check for escapes
    - frequently used:
        - `escape='"'`
            - double-quotes within quotes, e.g. `"firstname ""nickname"" lastname"`
        - `escape='\\\\'`
            - note that the first backslash is just to escape the second backslash
            - e.g. `"firstname \\\"nickname\\\" lastname"`
- `comment`
    - a `Char` or `String` at the beginning of lines that should be skipped as comments
        - note that skipped comment lines do not contribute to the line count for the header
          (if the user requests parsing a header on a specific row) or for skiprows
    - default: `comment=missing`
        - by default, the parser does not check for comments
    - frequently used:
        - `comment='#'`
        - `comment='!'`
        - `comment="#!"`
- `encodings`
    - a `Dict{String, Any}` mapping parsed fields to Julia values
        - if your dataset has booleans that are not represented as `"true"` and `"false"` or missing values that you'd like to read as `missing`s, you'll need to use this!
    - default: `encodings=Dict{String, Any}()`
        - by default, the parser does not check for any reserved fields
    - frequently used:
        - `encodings=Dict("" => missing)`
        - `encodings=Dict("NA" => missing)`
        - `encodings=Dict("N/A" => missing)`
        - `encodings=Dict("NULL" => missing)`
        - `encodings=Dict("TRUE" => true, "FALSE" => false)`
        - `encodings=Dict("True" => true, "False" => false)`
        - `encodings=Dict("T" => true, "F" => false)`
        - `encodings=Dict("yes" => true, "no" => false)`
        - ... your encodings here ...
            - can include any number of `String` => value mappings
            - note that if the user requests `quotes`, `escapes`, or `trimwhitespace`, these requests
              will be applied (removed) the raw string *BEFORE* checking whether the field matches
              any strings in in the `encodings` argument
- `header`
    - an `Int` indicating which line of the dataset contains column names or a `Vector{String}` of column names
        - note that commented lines and blank lines do not contribute to this value
          e.g. if the first 3 lines of your dataset are comments, you'll still need to
          set `header=1` to interpret the first line of parsed data as the header
    - default: `header=0`
        - no header is checked for by default
    - frequently used:
        - `header=1`
- `skiprows`
    - a `Range` or `Vector` of `Int`s indicating which rows to skip in the dataset
        - note that this is 1-based in reference to the first row *AFTER* the header.
          if `header=0` or is provided by the user, this will be the first non-empty
          line in the dataset. otherwise `skiprows=1:1` will skip the `header+1`-nth line
          in the file
    - default: `skiprows=Vector{Int}()`
        - no rows are skipped
- `types`
    - declare the types of the columns
        - scalar, e.g. `types=Bool`
            - scalars will be broadcast to apply to every column of the dataset
        - vector, e.g. `types=[Bool, Int, Float64, String, Symbol, Date, DateTime]`
            - the vector length must match the number of parsed columns
        - dictionary, e.g. `types=("column1" => Bool)` or `types=(1 => Union{Int, Missing})`
            - users can refer to the columns by name (only if a header is provided or
              parsed!) or by index
    - default:
        - `types=Dict{Int,DataType}()`
            - column-types will be interpreted from the dataset
    - built-in support for parsing the following:
        - `Int`
        - `Float64`
        - `String`
        - `Symbol`
        - `Date` -- only the default date format will work
        - `DateTime` -- only the default datetime format will work
        - for other types or unsupported formats, see `colparsers` and `typeparsers`
- `allowmissing`
    - declare whether columns should have element-type `Union{T, Missing} where T`
        - boolean scalar, e.g. `allowmissing=true`
            - scalars will be broadcast to apply to every column of the dataset
        - vector, e.g. `allowmissing=[true, false, true, true]`
            - the vector length must match the number of parsed columns
        - dictionary, e.g. `allowmissing=("column1" => true)` or `allowmissing=(17 => true)`
            - users can refer to the columns by name (only if a header is provided or
              parsed!) or by index
    - default: `allowmissing=Dict{Int,Bool}()`
        - Allowing missing values is determined by type detection in rows `1:typedetectrows`
- `coltypes`
    - declare the type of vector that should be used for columns
    - should work for any AbstractVector that allows `push!`ing values
        - scalar, e.g. `coltypes=CategoricalVector`
            - scalars will be broadcast to apply to every column of the dataset
        - vector, e.g. `coltypes=[CategoricalVector, Vector, CategoricalVector]`
            - the vector length must match the number of parsed columns
        - dictionary, e.g. `coltypes=("column1" => CategoricalVector)` or `coltypes=(17 => CategoricalVector)`
            - users can refer to the columns by name (only if a header is provided or
              parsed!) or by index
    - default: `coltypes=Vector`
        - all columns are returned as standard julia `Vector`s
- `colparsers`
    - provide custom functions for converting parsed strings to values by column
        - scalar, e.g. `colparsers=(x -> parse(Float64, replace(x, ',', '.')))`
            - scalars will be broadcast to apply to every column of the dataset
        - vector, e.g. `colparsers=[x -> mydateparser(x), x -> mytimeparser(x)]`
            - the vector length must match the number of parsed columns
        - dictionary, e.g. `colparsers=("column1" => x -> mydateparser(x))`
            - users can refer to the columns by name (only if a header is provided or
              parsed!) or by index
    - default: `colparsers=Dict{Int,Function}()`
        - column parsers are determined based on user-specified types and those
          detected from the data
- `typeparsers`
    - provide custom functions for converting parsed strings to values by column type
        - *NOTE* must be used with `coltypes`. If you supply a custom Int parser you'd like to
          use to parse column 6, you'll need to set `coltypes=dict(6 => Int)` for it to work
    - default: `colparsers=Dict{DataType, Function}()`
        - column parsers are determined based on user-specified types and those
          detected from the data
    - frequently used:
        - `typeparsers=Dict(Float64 => x -> parse(Float64, replace(x, ',' => '.')))` # decimal-comma floats!
- `typedetectrows`
    - specify how many rows of data to read before interpretting the values that each
      column should take on
    - default: `typedetectrows=1`
        - must be >= 1
        - commented, skipped, and empty lines are not counted when determining
          which rows are used for type detection, e.g. setting `typedetectrows=10` and
          `skiprows=1:5` means type detection will occur on rows `6:15`
- `skipmalformed`
    - specify whether the parser should skip a line or fail with an error if a line is
      parsed but does not contain the expected number of rows
    - default: `skipmalformed=false`
        - malformed lines result in an error
- `trimwhitespace`
    - specify whether should extra whitespace be removed from the beginning and ends of fields.
        - e.g `...., myfield ,...`
            - `trimwhitespace=false` -> `" myfield "`
            - `trimwhitespace=true`  -> `"myfield"`
    - leading and trailing whitespace *OUTSIDE* of quoted fields is trimmed by default.
        - e.g. `...., " myfield " ,...` -> `" myfield "` when `quotes='"'`
    - `trimwhitespace=true` will also trim leading and trailing whitespace *WITHIN* quotes
    - default: `trimwhitespace=false`
"""
function read(source::IO;
              delim::Union{Char,String}=',',
              quotes::Union{Char,Missing}=missing,
              escape::Union{Char,Missing}=missing,
              comment::Union{Char,String,Missing}=missing,
              encodings::Dict{String,T} where T=Dict{String,Any}(),
              header::Union{Integer,Vector{String}}=0,
              skiprows::AbstractVector{Int}=Vector{Int}(),
              types::Union{T1,COLMAP{T1},Vector{T1}} where {T1<:Type}=Dict{Int,DataType}(),
              allowmissing::Union{Bool,COLMAP{Bool},Vector{Bool}}=Dict{Int,Bool}(),
              coltypes::Union{Type{<:AbstractVector},COLMAP{UnionAll},Vector{UnionAll}}=Vector,
              colparsers::Union{F1, COLMAP{F1}, Vector{F1}} where F1=Dict{Int,Function}(),
              typeparsers::Dict{T2, F2} where {T2<:Type, F2}=Dict{DataType, Function}(),
              typedetectrows::Int=1,
              skipmalformed::Bool=false,
              trimwhitespace::Bool=false)

        reserved = [x for x in (delim, quotes, escape, comment) if !ismissing(x)]
        if !ismissing(quotes) && isequal(quotes, escape)
            @assert length(unique(string.(reserved))) == length(reserved) - 1
        else
            @assert length(unique(string.(reserved))) == length(reserved)
        end
        if isa(coltypes, Vector{UnionAll})
            @assert all(x -> x <: AbstractVector, coltypes)
        elseif isa(coltypes, COLMAP{UnionAll})
            @assert all(x -> x <: AbstractVector, collect(values(coltypes)))
        end
        @assert typedetectrows >= 1
        if typedetectrows > 100
            @warn """
                  Large values for `typedetectrows` will reduce performance. Consider manually declaring the types of columns using the `types` argument instead.
                  """
        end
        return parsesource(source, delim, quotes, escape, comment, encodings, header,
                           skiprows, types, allowmissing, coltypes, colparsers, typeparsers,
                           typedetectrows, skipmalformed, trimwhitespace)
end

function read(fullpath::String; kwargs...)
    read(open(fullpath); kwargs...)
end


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

JLLWrappers.@generate_main_file_header("HarfBuzz")
JLLWrappers.@generate_main_file("HarfBuzz", UUID("2e76f6c2-a576-52d4-95c1-20adfe4de566"))
end  # module HarfBuzz_jll


text:
using IncrementalInference
using Test

## a new factor that is broken
struct BrokenFactor{T<: SamplableBelief} <: AbstractRelativeRoots
    Z::T
end

IncrementalInference.getSample(cf::CalcFactor{<:BrokenFactor}, N::Int=1) = (reshape(rand(cf.factor.Z, N),:,N), )

function (s::CalcFactor{<:BrokenFactor})(z,
                                         wxi,
                                         wxj)
    #
    error("User factor has a bug.")
end

# # FIXME consolidate with CalcFactor according to #467
# function (s::BrokenFactor{<:IIF.ParametricTypes})(X1::AbstractArray{<:Real},
#                                                   X2::AbstractArray{<:Real};
#                                                   userdata::Union{Nothing,FactorMetadata}=nothing )
#     error("User factor has a bug -- USE NEW CalcFactor API INSTEAD, v0.21.")
# end

@testset "Test CSM monitor/watchdog on errors" begin

# create a factor graph
fg = generateCanonicalFG_lineStep(10;
                                  poseEvery=1,
                                  landmarkEvery=5,
                                  posePriorsAt=[0],
                                  sightDistance=4,
                                  solverParams=SolverParams(algorithms=[:default, :parametric]))

ensureAllInitialized!(fg)

#TODO test FSM watchdog
# add a broken factor - mid
addFactor!(fg, [:x9, :lm10], BrokenFactor(Normal()); graphinit=false)
smtasks = Task[]
@test_throws CompositeException tree, smt, hist = IIF.solveTree!(fg; smtasks=smtasks);
sleep(0.1)

# Test parametric solve also
addFactor!(fg, [:x9, :lm10], BrokenFactor(Normal()); graphinit=false)
@test_throws CompositeException tree2, smt, hist = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)
sleep(0.1)

deleteFactor!(fg, :x9lm10f2)
# add a broken factor - leave
addFactor!(fg, [:x10, :lm10], BrokenFactor(Normal()); graphinit=false)
@test_throws CompositeException tree2, smt, hist = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)
sleep(0.1)

deleteFactor!(fg, :x10lm10f2)
# add a broken factor - root
addFactor!(fg, [:x7, :lm10], BrokenFactor(Normal()); graphinit=false)
@test_throws CompositeException tree2, smt, hist = IIF.solveTree!(fg; smtasks=smtasks, algorithm = :parametric)

end

@testset "test CSM debug options work" begin

## create a factor graph

fg = generateCanonicalFG_lineStep(10;
                                  poseEvery=1,
                                  landmarkEvery=5,
                                  posePriorsAt=[0],
                                  sightDistance=4 )
#

smtasks = Task[];
solveTree!(fg, smtasks=smtasks, recordcliqs=ls(fg));

## make sure we fetch the results

hists = fetchCliqHistoryAll!(smtasks);

printCSMHistoryLogical(hists)

printCSMHistorySequential(hists)

## test CSM resolve steps

@info "test repeatCSMStep"
csmc_ = repeatCSMStep!(hists, 1, 1)

##

end


text:
__precompile__()

module Whippet

import BioSymbols
import BioSequences
import BioAlignments
import FASTX
import XAM

import DataStructures: SortedSet

using StatsBase: midpoints
using BioSymbols
using BioSequences
using BioAlignments
using GenomicFeatures
using FASTX
using XAM
using BufferedStreams
using Distributions
using FMIndexes
using IntArrays
using IntervalTrees
using Libz
using LinearAlgebra
using Nullables
using Printf
using Random

include("types.jl")
include("timer.jl")
include("sgkmer.jl")
include("fmindex_patch.jl")
include("bam.jl")
include("refset.jl")
include("graph.jl")
include("bias.jl")
include("edges.jl")
include("index.jl")
include("record.jl")
include("align.jl")
include("quant.jl")
include("reads.jl")
include("paired.jl")
include("events.jl")
include("io.jl")
include("diff.jl")

export SpliceGraph,
       SpliceGraphQuant,
       GraphLib,
       GraphLibQuant,
       GeneInfo,
       TxInfo,
       SGNode,
       Edges,
       EdgeType,
       EDGETYPE_LS,
       EDGETYPE_SL,
       EDGETYPE_RS,
       EDGETYPE_SR,
       EDGETYPE_LR,
       EDGETYPE_LL,
       EDGETYPE_RR,
       Kmer,
       SGKmer,
       SGNucleotide,
       SGNucleotideSequence,
       SGSequence,
       SGAlignSingle,
       SGAlignPaired,
       @sg_str,
       reverse_complement,
       kmer,
       NucleotideSequence,
       RefSet, RefTx, RefGene,
       MultiMapping,
       DefaultBiasMod,
       DefaultCounter,
       PrimerBiasMod,
       PrimerBiasCounter,
       GCBiasMod,
       GCBiasCounter,
       JointBiasMod,
       JointBiasCounter,
       AlignParam,
       ungapped_align,
       make_fqparser,
       ident_to_fastq_url,
       hasextension,
       isgzipped,
       load_refflat,
       load_gtf,
       fasta_to_index,
       process_reads!,
       process_paired_reads!,
       build_equivalence_classes!,
       set_global_readcount,
       total_multi,
       primer_adjust!,
       gc_adjust!,
       primer_normalize!,
       gc_normalize!,
       calculate_tpm!,
       gene_em!,
       set_gene_tpm!,
       output_tpm,
       output_stats,
       output_junctions,
       output_exons,
       assign_ambig!,
       effective_lengths!,
       global_bias,
       process_events,
       fixpath,
       @timer,
       timer_print,
       open_streams,
       process_psi_files
end


text:
using Symbolics
using SymbolicUtils, Test
using Symbolics: symtype, shape, wrap, unwrap, Unknown, Arr, arrterm, jacobian, @variables, value, get_variables
using Base: Slice
using SymbolicUtils: Sym, term, operation

@testset "arrays" begin
    @variables X[1:5, 1:5] Y[1:5, 1:5]
    @test_throws BoundsError X[1000]
    @test typeof(X) <: Arr
    @test shape(X) == Slice.((1:5, 1:5))
    @test shape(Y) == Slice.((1:5, 1:5))

    A = Y[2, :]
    @test typeof(A) <: Arr{Num, 1}
    @test axes(A) == (1:5,)

    B = A[3:5]
    @test axes(B) == (Slice(1:3),)

    i = Sym{Int}(:i)
    j = Sym{Int}(:j)
    @test symtype(X[i,j]) == Real
    @test symtype(X[1,j]) == Real

    @variables t x[1:2](t)
    @test isequal(get_variables(0 ~ x[1]), [x[1]])
    @test Set(get_variables(2x)) == Set(collect(x)) # both array elements are present
    @test isequal(get_variables(2x[1]), [x[1]])
end

@testset "getindex" begin
    @variables X[1:5, 1:5] Y[1:5, 1:5]

    @test isequal(X[1,1], wrap(term(getindex, unwrap(X), 1,1)))

    XX = unwrap(X)
    @test isequal(unwrap(X[1, :]), Symbolics.@arrayop(XX[1,:], (j,), XX[1, j]))
    @test isequal(unwrap(X[:, 2]), Symbolics.@arrayop(XX[:,2], (i,), XX[i, 2]))
    @test isequal(unwrap(X[:, 2:3]), Symbolics.@arrayop(XX[:,2:3], (i,j), XX[i, j], (+), (j in 2:3)))

    @variables t x[1:4](t)
    @syms i::Int
    @test isequal(x[i], operation(unwrap(x[i]))(t))
end

getdef(v) = getmetadata(v, Symbolics.VariableDefaultValue)
@testset "broadcast & scalarize" begin
    @variables A[1:5,1:3]=42 b[1:3]=[2, 3, 5] t x[1:4](t)
    AA = Symbolics.scalarize(A)
    bb = Symbolics.scalarize(b)
    @test all(isequal(42), getdef.(AA))
    @test getdef.(bb) == [2, 3, 5]
    @test isequal(Symbolics.scalarize([b.*1; b.*1]), [bb; bb])
    @test isequal(Symbolics.scalarize(b.^1), bb)
    c = A*b

    @test isequal(collect(sin.(x)),
                  sin.([x[i] for i in 1:4]))

    @test isequal(Symbolics.scalarize(sin.(A.*c)[1,1]),
                  sin(A[1, 1]*(b[1]*A[1, 1] +
                               b[2]*A[1, 2] +
                               b[3]*A[1, 3])))

    D = Differential(t)
    @test isequal(collect(D.(x) ~ x), map(i->D(x[i]) ~ x[i], eachindex(x)))
    @test_throws ArgumentError A ~ t
end

@testset "Parent" begin
    @variables t x[1:4](t)
    x = unwrap(x)
    @test Symbolics.getparent(collect(x)[1]).metadata === x.metadata
end

n = 2
A = randn(n,n)
foo(x) = A*x # a function to represent symbolically, note, if this function is defined inside the testset, it's not found by the function fun_eval = eval(fun_ex)
function Symbolics.propagate_ndims(::typeof(foo), x)
    ndims(x)
end
function Symbolics.propagate_shape(::typeof(foo), x)
    shape(x)
end
@wrapped function foo(x::AbstractVector)
    t = arrterm(foo, x)
    setmetadata(t, Symbolics.ScalarizeCache, Ref{Any}(nothing))
end

#=
The following two testsets test jacobians for symbolic functions of symbolic arrays. Getting it to work currently requires the user to manually specify propagate_ndims and propagate_shape. Making use of `@register` or `@syms` to specify symbolic functions fail to propagate array shape
=#

@testset "Functions and Jacobians using @syms" begin
    @variables x[1:n]

    function symbolic_call(x)
        @syms foo(x::Symbolics.Arr{Num, 1})::Symbolics.Arr{Num, 1} # symbolic foo can not be created in global scope due to conflict with function foo
        foo(x) # return a symbolic call to foo
    end

    x0 = randn(n)
    @test foo(x0) == A*x0
    ex = symbolic_call(x)

    fun_genf = build_function(ex, x, expression=Val{false})
    @test_broken fun_genf(x0) == A*x0# UndefVarError: foo not defined

    # Generate an expression instead and eval it manually
    fun_ex = build_function(ex, x, expression=Val{true})
    fun_eval = eval(fun_ex)
    @test fun_eval(x0) == foo(x0)

    # Try to provide the hidden argument `expression_module` to solve the scoping issue
    @test_skip begin
        fun_genf = build_function(ex, x, expression=Val{false}, expression_module=Main) # UndefVarError: #_RGF_ModTag not defined
        fun_genf(x0) == A*x0
    end

    ## Jacobians
    @test Symbolics.value.(Symbolics.jacobian(foo(x), x)) == A
    @test_skip Symbolics.value.(Symbolics.jacobian(ex , x)) == A #ERROR: axes of foo(x[1:2]) not known
end

@testset "Functions and Jacobians using manual @wrapped" begin
    @variables x[1:n]

    x0 = randn(n)
    @test foo(x0) == A*x0
    ex = foo(x)

    @test shape(ex) == shape(x)

    fun_genf = build_function(ex, x, expression=Val{false})
    @test fun_genf(x0) == A*x0

    # Generate an expression instead and eval it manually
    fun_ex = build_function(ex, x, expression=Val{true})
    fun_eval = eval(fun_ex)
    @test fun_eval(x0) == foo(x0)

    ## Jacobians
    @test value.(jacobian(foo(x), x)) == A
    @test value.(jacobian(ex , x)) == A
end


text:
module CollectionBenchmarks

using BenchmarkTools
using Random
using StableRNGs

const SUITE = BenchmarkGroup()

const RNG = StableRNG(1)

const iterlen = 1000
const ints = rand(RNG, 1:iterlen, iterlen)
const strings = [randstring(RNG, rand(RNG, 1:30)) for i in 1:iterlen]
const anys = (shuffle!(RNG, [ints; strings;])[1:iterlen])::Vector{Any}
const sortedints = sort(ints)

# random element, in the collections iff 2nd parameter is true
let anyshuffled = shuffle(RNG, anys),
    filter = Iterators.filter
    d = Dict((Int,    true ) => first(filter(i->isa(i, Int),    anyshuffled))::Int,
             (String, true ) => first(filter(i->isa(i, String), anyshuffled))::String,
             (Int,    false) => first(filter(i->!in(i, ints), rand(RNG, 1:iterlen) for _ in 1:typemax(Int)))::Int,
             (String, false) => randstring(RNG, 40)::String)
    global randelt
    randelt(C, T, bool) = (T == Any && (T = Int); C === Dict || C === IdDict ? (d[T, bool] => d[T, bool]) : d[T, bool])
end

const coll = Dict((Vector, Int)    => ints,
                  (Vector, String) => strings,
                  (Vector, Any)    => anys,
                  (Pair, Int)      => map(Pair, ints, ints),
                  (Pair, String)   => map(Pair, strings, strings),
                  (Pair, Any)      => map(Pair, anys, anys),
                  (BitSet, Int)    => BitSet(ints))

# return a collection that can serve to initialize a container C of type T
initcoll(C, T) = let initmap = Dict(Vector => Vector,
                                    Dict   => Pair,
                                    IdDict => Pair,
                                    Set    => Vector,
                                    BitSet => Vector)
    coll[initmap[C], T]
end

for C in (Dict, IdDict, Set), T in (Int, String, Any)
    coll[C, T] = C(initcoll(C, T))
end

@inline newcoll(::Type{C}, ::Type{T}) where {C,T} = C{T}()
@inline newcoll(::Type{Dict}, ::Type{T}) where {T} = Dict{T,T}()
@inline newcoll(::Type{IdDict}, ::Type{T}) where {T} = IdDict{T,T}()
@inline newcoll(::Type{BitSet}, ::Type{Int}) = BitSet()

askey(::Type, elt) = elt
askey(::Type{Dict}, elt) = first(elt)
askey(::Type{IdDict}, elt) = first(elt)

objstr(T) = T === Int ? "Int" : T === Vector ? "Vector" : string(T)

function foreach_container(bench; T = (Int, String, Any), C = (Vector, Dict, IdDict, Set, BitSet))
    for _C in C
        cstr = objstr(_C)
        for _T in T
            _C === BitSet && _T !== Int && continue
            c = coll[_C, _T]
            bench(_C, cstr, _T, objstr(_T), c)
        end
    end
end

function set_tolerance!(g, tol=0.25)
    for b in values(g)
        b.params.time_tolerance = tol
    end
end

##################
# Initialization #
##################

g = addgroup!(SUITE, "initialization", ["AbstractVector", "AbstractSet", "Associative"])

function perf_push!(c, elts)
    for e in elts
        push!(c, e)
    end
end

foreach_container(C = (Vector, Dict, IdDict, Set)) do C, cstr, T, tstr, c
    g[cstr, tstr, "iterator"] = @benchmarkable $C($(initcoll(C, T)))
    g[cstr, tstr, "loop"]     = @benchmarkable perf_push!(newcoll($C, $T), $(initcoll(C, T)))
    g[cstr, tstr, "loop", "sizehint!"] =
        @benchmarkable perf_push!(sizehint!(newcoll($C, $T), iterlen), $(initcoll(C, T)))
end

g["BitSet", "Int", "unsorted", "iterator"]          = @benchmarkable BitSet($ints)
g["BitSet", "Int", "sorted",   "iterator"]          = @benchmarkable BitSet($sortedints)
g["BitSet", "Int", "unsorted", "loop"]              = @benchmarkable perf_push!(BitSet(), $ints)
g["BitSet", "Int", "sorted",   "loop"]              = @benchmarkable perf_push!(BitSet(), $sortedints)
g["BitSet", "Int", "unsorted", "loop", "sizehint!"] = @benchmarkable perf_push!(sizehint!(BitSet(), iterlen), $ints)
g["BitSet", "Int", "sorted",   "loop", "sizehint!"] = @benchmarkable perf_push!(sizehint!(BitSet(), iterlen), $sortedints)

set_tolerance!(g)

#############
# Iteration #
#############

g = addgroup!(SUITE, "iteration", ["AbstractVector", "AbstractSet", "Associative"])

foreach_container() do C, cstr, T, tstr, c
    g[cstr, tstr, "iterate"] = @benchmarkable iterate($c)
    g[cstr, tstr, "iterate second"]  = @benchmarkable iterate($c, $(iterate(c)[2]))
end

set_tolerance!(g)

############
# Deletion #
############

g = addgroup!(SUITE, "deletion", ["AbstractVector", "AbstractSet", "Associative"])

function perf_pop!(c)
    while !isempty(c)
        pop!(c)
    end
end

pred(::Type{C}, ::Type{Any})    where {C} = x -> isa(askey(C, x), Int)
pred(::Type{C}, ::Type{String}) where {C} = x -> Int(askey(C, x)[1]) < 90
pred(::Type{C}, ::Type{Int})    where {C} = x -> iseven(askey(C, x))

foreach_container() do C, cstr, T, tstr, c
    if C !== IdDict
        g[cstr, tstr, "pop!"] = @benchmarkable perf_pop!(d) setup=(d=copy($c)) evals=1
    end
    C === BitSet && return
    if C !== Dict
        g[cstr, tstr, "filter!"] = @benchmarkable filter!($(pred(C, T)), d) setup=(d=copy($c)) evals=1
        g[cstr, tstr, "filter"] =  @benchmarkable filter( $(pred(C, T)), $c)
    end
end

set_tolerance!(g)

#####################
# Queries & Updates #
#####################

g = addgroup!(SUITE, "queries & updates", ["AbstractVector", "AbstractSet", "Associative"])

foreach_container() do C, cstr, T, tstr, c
    if T === Int # seems unnecessary to run those with all types
        g[cstr, tstr, "length"] = @benchmarkable length($c) # probably useful only for BitSet
        g[cstr, tstr, "first"]  = @benchmarkable first($c)
        if C in (Vector, BitSet)
            g[cstr, tstr, "last"] = @benchmarkable last($c)
        end
    end
    eltin, eltout = randelt(C, T, true), randelt(C, T, false)
    g[cstr, tstr, "in", "true"]  = @benchmarkable in($eltin, $c)
    g[cstr, tstr, "in", "false"] = @benchmarkable in($eltout, $c)
    if C === Vector
        g[cstr, tstr, "getindex"]  = @benchmarkable $c[$(iterlen÷2)]
        g[cstr, tstr, "setindex!"] = @benchmarkable d[$(iterlen÷2)] = $eltout setup=(d=copy($c))
    elseif C === Dict || C === IdDict
        keyin  = askey(C, eltin)
        keyout = askey(C, eltout)
        g[cstr, tstr, "getindex"]               = @benchmarkable $c[$keyin]
        g[cstr, tstr, "setindex!", "overwrite"] = @benchmarkable d[$keyin]  = $keyout setup=(d=copy($c))
        g[cstr, tstr, "setindex!", "new"]       = @benchmarkable d[$keyout] = $keyout setup=(d=copy($c)) evals=1
    end
    if C === Vector
        g[cstr, tstr, "push!"] = @benchmarkable push!(d, $eltout) setup=(d=copy($c))
    else
        g[cstr, tstr, "push!", "overwrite"] = @benchmarkable push!(d, $eltin)             setup=(d=copy($c))
        g[cstr, tstr, "push!", "new"]       = @benchmarkable push!(d, $eltout)            setup=(d=copy($c)) evals=1
        # despite the `evals=1` parameter, this gets evaluated more than once, leading to a KeyError
        # g[cstr, tstr, "pop!",  "specified"] = @benchmarkable  pop!(d, $(askey(C, eltin))) setup=(d=copy($c)) evals=1
        # so we add a push!
        g[cstr, tstr, "pop!",  "specified"] = @benchmarkable  pop!(push!(d, $eltin), $(askey(C, eltin))) setup=(d=copy($c))

    end
    if C !== IdDict
        # same remark as above, we need to push! a value to avoid ending up with an empty container
        g[cstr, tstr, "pop!", "unspecified"] = @benchmarkable push!(d, pop!(d)) setup=(d=copy($c))
    end
end

set_tolerance!(g)

##################
# Set operations #
##################

g = addgroup!(SUITE, "set operations", ["AbstractSet", "Array"])
const newints = [rand(RNG, ints, 10); rand(RNG, 1:iterlen, 10); rand(RNG, iterlen:2iterlen, 10);]

foreach_container(C = (BitSet, Set, Vector), T = (Int,)) do C, cstr, T, tstr, c
    g[cstr, tstr, "union"]     = @benchmarkable union($c)
    g[cstr, tstr, "intersect"] = @benchmarkable intersect($c)
    g[cstr, tstr, "symdiff"]   = @benchmarkable symdiff($c)
    for C2 in (BitSet, Set, Vector)
        c2 = C2(newints)
        c2str = objstr(C2)
        g[cstr, tstr, "union",     c2str] = @benchmarkable union($c, $c2)
        g[cstr, tstr, "intersect", c2str] = @benchmarkable intersect($c, $c2)
        g[cstr, tstr, "symdiff",   c2str] = @benchmarkable symdiff($c, $c2)
        g[cstr, tstr, "setdiff",   c2str] = @benchmarkable setdiff($c, $c2)

        g[cstr, tstr, "union", c2str, c2str]     = @benchmarkable union($c, $c2, $c2)
        g[cstr, tstr, "intersect", c2str, c2str] = @benchmarkable intersect($c, $c2, $c2)
        g[cstr, tstr, "symdiff", c2str, c2str]   = @benchmarkable symdiff($c, $c2, $c2)

        if C === BitSet && C2 === BitSet
            g[cstr, tstr, "intersect!", c2str] = @benchmarkable intersect!(d, $c2) setup=(d=copy($c)) evals=1
        end
        if C === BitSet
            g[cstr, tstr, "symdiff!",   c2str] = @benchmarkable symdiff!(d, $c2)   setup=(d=copy($c)) evals=1
        end
        if C in (BitSet, Set)
            g[cstr, tstr, "setdiff!",   c2str] = @benchmarkable setdiff!(d, $c2)   setup=(d=copy($c)) evals=1
            g[cstr, tstr, "union!",     c2str] = @benchmarkable union!(d, $c2)     setup=(d=copy($c)) evals=1
        end
        g[cstr, tstr, "⊆", c2str] = @benchmarkable ⊆($c, $c2)
    end
    C === Vector && return
    empty = newcoll(C, Int)
    c2  = C(newints)
    g[cstr,    tstr, "<",  cstr]   = @benchmarkable <($c, $c2)
    g["empty", tstr, "<",  cstr]   = @benchmarkable <($empty, $c)
    g["empty", tstr, "⊆",  cstr]   = @benchmarkable ⊆($empty, $c)
    g[cstr,    tstr, "⊆",  "self"] = @benchmarkable ⊆($c, $c)
    g[cstr,    tstr, "==", cstr]   = @benchmarkable ==($c, $c2)
    g[cstr,    tstr, "==", "self"] = @benchmarkable ==($c, $c)
end

# test BitSet with very large values
const small, big = BitSet(2), BitSet(2^18)
g["BitSet", "Int", "union!",     "big"]   = @benchmarkable union!(b, $small)     setup=(b=BitSet(2^18)) evals=1
g["BitSet", "Int", "union!",     "small"] = @benchmarkable union!(s, $big)       setup=(s=BitSet(2))    evals=1
g["BitSet", "Int", "intersect!", "big"]   = @benchmarkable intersect!(b, $small) setup=(b=BitSet(2^18)) evals=1
g["BitSet", "Int", "intersect!", "small"] = @benchmarkable intersect!(s, $big)   setup=(s=BitSet(2))    evals=1
g["BitSet", "Int", "setdiff!",   "big"]   = @benchmarkable setdiff!(b, $small)   setup=(b=BitSet(2^18)) evals=1
g["BitSet", "Int", "setdiff!",   "small"] = @benchmarkable setdiff!(s, $big)     setup=(s=BitSet(2))    evals=1
g["BitSet", "Int", "symdiff!",   "big"]   = @benchmarkable symdiff!(b, $small)   setup=(b=BitSet(2^18)) evals=1
g["BitSet", "Int", "symdiff!",   "small"] = @benchmarkable symdiff!(s, $big)     setup=(s=BitSet(2))    evals=1

set_tolerance!(g)

#############################################
# Optimizations for types with "few" values #
#############################################

# cf. issue #20903 and PR #21964
g = addgroup!(SUITE, "optimizations", ["Dict", "IdDict", "Set", "BitSet", "Vector"])

for T in (Nothing, Bool, Int8, UInt16)
    local v
    v::Vector{T} = T === Nothing ? Vector{Nothing}(undef, 100000) :
                                   rand(RNG, one(T):typemax(T), 100000)
    tstr = string(T)
    g["Dict", "abstract", tstr] = @benchmarkable Dict($(map(Pair, v, v)))
    g["Dict", "concrete", tstr] = @benchmarkable Dict{$T,$T}($(map(Pair, v, v)))
    g["IdDict", "abstract", tstr] = @benchmarkable IdDict($(map(Pair, v, v)))
    g["IdDict", "concrete", tstr] = @benchmarkable IdDict{$T,$T}($(map(Pair, v, v)))
    g["Set",  "abstract", tstr] = @benchmarkable Set($v)
    g["Set",  "concrete", tstr] = @benchmarkable Set{$T}($v)
    if T === Nothing
        g["Vector", "abstract", tstr] = @benchmarkable Vector($v)
        g["Vector", "concrete", tstr] = @benchmarkable Vector{$T}($v)
    elseif T !== Bool && (!(T <: Unsigned) || VERSION >= v"0.7.0-")
        # there is a bug on 0.6, will probably get fixed on v6.2, cf. #24365
        g["BitSet", tstr] = @benchmarkable BitSet($v)
    end
end

set_tolerance!(g)

end # module


text:
# TODO The mutability of the coefficient type `T` is currently not exploited yet.
#      We first need to make sure that it is copied with `MA.copy_if_mutable` when it
#      is passed from one function to a mutable one.

MA.mutability(::Type{<:TypedLike}) = MA.IsMutable()

function MA.mutable_copy(func::MOI.ScalarAffineFunction)
    terms = [
        MOI.ScalarAffineTerm(
            MA.copy_if_mutable(t.coefficient),
            t.variable_index,
        ) for t in func.terms
    ]
    return MOI.ScalarAffineFunction(terms, MA.copy_if_mutable(func.constant))
end
function MA.mutable_copy(func::MOI.ScalarQuadraticFunction)
    affine_terms = [
        MOI.ScalarAffineTerm(
            MA.copy_if_mutable(t.coefficient),
            t.variable_index,
        ) for t in func.affine_terms
    ]
    quadratic_terms = [
        MOI.ScalarQuadraticTerm(
            MA.copy_if_mutable(t.coefficient),
            t.variable_index_1,
            t.variable_index_2,
        ) for t in func.quadratic_terms
    ]
    return MOI.ScalarQuadraticFunction(
        affine_terms,
        quadratic_terms,
        MA.copy_if_mutable(func.constant),
    )
end

function MA.isequal_canonical(
    f::F,
    g::F,
) where {F<:Union{
    MOI.ScalarAffineFunction,
    MOI.ScalarQuadraticFunction,
    MOI.VectorAffineFunction,
    MOI.VectorQuadraticFunction,
}}
    return MA.isequal_canonical(MOI.constant(f), MOI.constant(g)) &&
           all(MOI.dict_compare.(
        MOI._dicts(f),
        MOI._dicts(g),
        MA.isequal_canonical,
    ))
end

function MA.iszero!(f::TypedScalarLike)
    return iszero(MOI.constant(f)) && _is_constant(canonicalize!(f))
end

function MA.scaling(f::TypedScalarLike{T}) where {T}
    g = canonical(f)
    if !_is_constant(g)
        throw(InexactError(:convert, T, f))
    end
    return MA.scaling(MOI.constant(g))
end

function MA.promote_operation(
    ::Union{typeof(zero),typeof(one)},
    F::Type{<:TypedScalarLike},
)
    return F
end

# To avoid type piracy, we add at least one `ScalarLike` outside of the `...`.
const PROMOTE_IMPLEMENTED_OP = Union{typeof(+),typeof(-),typeof(*),typeof(/)}
function MA.promote_operation(
    op::PROMOTE_IMPLEMENTED_OP,
    F::Type{<:ScalarLike{T}},
    G::Type{<:ScalarLike{T}},
) where {T,N}
    return promote_operation(op, T, F, G)
end
function MA.promote_operation(
    op::PROMOTE_IMPLEMENTED_OP,
    F::Type{T},
    G::Type{<:TypedLike{T}},
) where {T,N}
    return promote_operation(op, T, F, G)
end
function MA.promote_operation(
    op::PROMOTE_IMPLEMENTED_OP,
    F::Type{<:TypedLike{T}},
    G::Type{T},
) where {T}
    return promote_operation(op, T, F, G)
end
function MA.promote_operation(
    op::PROMOTE_IMPLEMENTED_OP,
    F::Type{<:Number},
    G::Type{<:Union{MOI.SingleVariable,MOI.VectorOfVariables}},
)
    return promote_operation(op, F, F, G)
end
function MA.promote_operation(
    op::PROMOTE_IMPLEMENTED_OP,
    F::Type{<:Union{MOI.SingleVariable,MOI.VectorOfVariables}},
    G::Type{<:Number},
)
    return promote_operation(op, G, F, G)
end

function MA.mutable_operate!(
    op::Union{typeof(zero),typeof(one)},
    f::MOI.ScalarAffineFunction,
)
    empty!(f.terms)
    f.constant = op(f.constant)
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(zero),typeof(one)},
    f::MOI.ScalarQuadraticFunction,
)
    empty!(f.affine_terms)
    empty!(f.quadratic_terms)
    f.constant = op(f.constant)
    return f
end

function MA.mutable_operate!(::typeof(-), f::MOI.ScalarQuadraticFunction)
    operate_terms!(-, f.quadratic_terms)
    operate_terms!(-, f.affine_terms)
    f.constant = -f.constant
    return f
end
function MA.mutable_operate!(::typeof(-), f::MOI.ScalarAffineFunction)
    operate_terms!(-, f.terms)
    f.constant = -f.constant
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarAffineFunction{T},
    g::T,
) where {T}
    f.constant = op(f.constant, g)
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarAffineFunction{T},
    g::MOI.SingleVariable,
) where {T}
    push!(f.terms, MOI.ScalarAffineTerm(op(one(T)), g.variable))
    return f
end
function MA.mutable_operate_to!(
    output::MOI.ScalarAffineFunction{T},
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarAffineFunction{T},
    g::MOI.ScalarAffineFunction{T},
) where {T}
    empty!(output.terms)
    append!(output.terms, f.terms)
    append!(output.terms, operate_terms(op, g.terms))
    output.constant = op(f.constant, g.constant)
    return output
end
function MA.mutable_operate_to!(
    output::MOI.ScalarQuadraticFunction{T},
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarQuadraticFunction{T},
    g::MOI.ScalarQuadraticFunction{T},
) where {T}
    empty!(output.affine_terms)
    append!(output.affine_terms, f.affine_terms)
    append!(output.affine_terms, operate_terms(op, g.affine_terms))
    empty!(output.quadratic_terms)
    append!(output.quadratic_terms, f.quadratic_terms)
    append!(output.quadratic_terms, operate_terms(op, g.quadratic_terms))
    output.constant = op(f.constant, g.constant)
    return output
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarAffineFunction{T},
    g::MOI.ScalarAffineFunction{T},
) where {T}
    append!(f.terms, operate_terms(op, g.terms))
    f.constant = op(f.constant, g.constant)
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarQuadraticFunction{T},
    g::T,
) where {T}
    f.constant = op(f.constant, g)
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarQuadraticFunction{T},
    g::MOI.SingleVariable,
) where {T}
    push!(f.affine_terms, MOI.ScalarAffineTerm(op(one(T)), g.variable))
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarQuadraticFunction{T},
    g::MOI.ScalarAffineFunction{T},
) where {T}
    append!(f.affine_terms, operate_terms(op, g.terms))
    f.constant = op(f.constant, g.constant)
    return f
end
function MA.mutable_operate!(
    op::Union{typeof(+),typeof(-)},
    f::MOI.ScalarQuadraticFunction{T},
    g::MOI.ScalarQuadraticFunction{T},
) where {T}
    append!(f.affine_terms, operate_terms(op, g.affine_terms))
    append!(f.quadratic_terms, operate_terms(op, g.quadratic_terms))
    f.constant = op(f.constant, g.constant)
    return f
end

_constant(::Type{T}, α::T) where {T} = α
_constant(::Type{T}, ::MOI.SingleVariable) where {T} = zero(T)
_constant(::Type{T}, func::TypedScalarLike{T}) where {T} = MOI.constant(func)

_affine_terms(f::MOI.ScalarAffineFunction) = f.terms
_affine_terms(f::MOI.ScalarQuadraticFunction) = f.affine_terms

function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    α::T,
    f::MOI.SingleVariable,
    β::T,
) where {T}
    push!(terms, MOI.ScalarAffineTerm(op(α * β), f.variable))
    return
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    f::MOI.SingleVariable,
    β::T,
) where {T}
    push!(terms, MOI.ScalarAffineTerm(op(β), f.variable))
    return
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    α::T,
    f::MOI.SingleVariable,
) where {T}
    push!(terms, MOI.ScalarAffineTerm(op(α), f.variable))
    return
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    f::MOI.SingleVariable,
) where {T}
    push!(terms, MOI.ScalarAffineTerm(op(one(T)), f.variable))
    return
end

function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    α::T,
    f::TypedScalarLike{T},
    β::T,
) where {T}
    for t in _affine_terms(f)
        push!(terms, operate_term(op, operate_term(*, α, t, β)))
    end
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    f::TypedScalarLike{T},
    β::T,
) where {T}
    for t in _affine_terms(f)
        push!(terms, operate_term(op, operate_term(*, t, β)))
    end
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    α::T,
    f::TypedScalarLike{T},
) where {T}
    for t in _affine_terms(f)
        push!(terms, operate_term(op, operate_term(*, α, t)))
    end
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    f::TypedScalarLike{T},
) where {T}
    append!(terms, operate_terms(op, _affine_terms(f)))
    return
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    args::Vararg{T,N},
) where {T,N}
    return
end
function _add_sub_affine_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarAffineTerm{T}},
    α::T,
    β::T,
    args::Vararg{ScalarQuadraticLike,N},
) where {T,N}
    return _add_sub_affine_terms(op, terms, α * β, args...)
end

function MA.mutable_operate!(
    op::MA.AddSubMul,
    f::MOI.ScalarAffineFunction{T},
    args::Vararg{ScalarAffineLike{T},N},
) where {T,N}
    f.constant = op(f.constant, _constant.(T, args)...)
    _add_sub_affine_terms(MA.add_sub_op(op), f.terms, args...)
    return f
end

function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    α::ScalarAffineLike{T},
    f::MOI.ScalarQuadraticFunction{T},
    β::ScalarAffineLike{T},
) where {T}
    for t in f.quadratic_terms
        push!(
            terms,
            operate_term(
                op,
                operate_term(*, _constant(T, α), t, _constant(T, β)),
            ),
        )
    end
end
function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    f::MOI.ScalarQuadraticFunction{T},
    β::ScalarAffineLike{T},
) where {T}
    for t in f.quadratic_terms
        push!(terms, operate_term(op, operate_term(*, t, _constant(T, β))))
    end
end
function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    α::ScalarAffineLike{T},
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    for t in f.quadratic_terms
        push!(terms, operate_term(op, operate_term(*, _constant(T, α), t)))
    end
end
function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    append!(terms, operate_terms(op, f.quadratic_terms))
    return
end
function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    # Compiler fails in StackOverflowError on Julia v1.1
    #args::Vararg{ScalarAffineLike{T}, N}) where {T, N}
    args::ScalarAffineLike{T},
) where {T}
    return
end
function _merge_constants(
    ::Type{T},
    α::ScalarAffineLike{T},
    β::ScalarAffineLike{T},
    args::Vararg{Any,N},
) where {T,N}
    return (_constant(T, α) * _constant(T, β), args...)
end
function _add_quadratic_terms(
    op::Union{typeof(+),typeof(-)},
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
    args::Vararg{Any,N},
) where {T,N}
    return _add_quadratic_terms(op, terms, _merge_constants(T, args...)...)
end

_num_function_with_terms(::Type{T}, ::T) where {T} = 0
_num_function_with_terms(::Type{T}, ::ScalarLike{T}) where {T} = 1
function _num_function_with_terms(
    ::Type{T},
    f::ScalarQuadraticLike{T},
    args::Vararg{ScalarQuadraticLike{T},N},
) where {T,N}
    return _num_function_with_terms(T, f) + _num_function_with_terms(T, args...)
end
function MA.mutable_operate!(
    op::MA.AddSubMul,
    f::MOI.ScalarQuadraticFunction{T},
    args::Vararg{ScalarQuadraticLike{T},N},
) where {T,N}
    if isone(_num_function_with_terms(T, args...))
        f.constant = op(f.constant, _constant.(T, args)...)
        _add_sub_affine_terms(MA.add_sub_op(op), f.affine_terms, args...)
        _add_quadratic_terms(MA.add_sub_op(op), f.quadratic_terms, args...)
        return f
    else
        return MA.mutable_operate!(MA.add_sub_op(op), f, *(args...))
    end
end
# `args` could be `(x', a)` where `a` is a vector of constants and `x` a vector
# of affine functions for instance.
function MA.mutable_operate!(
    op::MA.AddSubMul,
    f::TypedScalarLike,
    args::Vararg{Any,N},
) where {N}
    return MA.mutable_operate!(MA.add_sub_op(op), f, *(args...))
end


text:
"""
    Models.TestUtils

Provides test fakes, [`FakeTemplate`](@ref) and [`FakeModel`](@ref), that are useful for
testing downstream dependencies, and [`test_interface`](@ref) for testing the Models API has
been correctly implemented.
"""
module TestUtils
using Distributions
using Models
using NamedDims
using StatsBase
using Test

export FakeModel, FakeTemplate
export test_interface

"""
    FakeTemplate{E <: EstimateTrait, O <: OutputTrait} <: Template

This template is a [test double](https://en.wikipedia.org/wiki/Test_double) for testing
purposes. It should be defined (before fitting) with a `predictor`, which can be changed by
mutating the field.

## Fields
- `predictor::Function`: predicts the outputs of the [`FakeModel`](@ref).
   It is `(num_variates, inputs) -> outputs`, where the `num_variates` will be memorized
   during `fit`.

## Methods
- [`fit`](@ref) does not learn anything, it just creates an instance of the corresponding [`Model`](@ref).
- [`predict`](@ref) applies the `predictor` to the inputs.
"""
mutable struct FakeTemplate{E<:EstimateTrait, O<:OutputTrait} <: Template
    predictor::Function
end

"""
    FakeTemplate{PointEstimate, SingleOutput}()

A [`Template`](@ref) whose [`Model`](@ref) will accept real value variables to predict 0
for each observation.
"""
function FakeTemplate{PointEstimate, SingleOutput}()
    FakeTemplate{PointEstimate, SingleOutput}() do num_variates, inputs
        @assert(num_variates == 1, "$num_variates != 1")
        inputs = _handle_inputs(inputs)
        return NamedDimsArray{(:variates, :observations)}(
            zeros(1, size(inputs, :observations))
        )
    end
end

"""
    FakeTemplate{PointEstimate, MultiOutput}()

A [`Template`](@ref) whose [`Model`](@ref) will accept real value variables to predict a
vector of 0s for each observation. The input and output will have the same dimension.
"""
function FakeTemplate{PointEstimate, MultiOutput}()
    FakeTemplate{PointEstimate, MultiOutput}() do num_variates, inputs
        inputs = _handle_inputs(inputs)
        return NamedDimsArray{(:variates, :observations)}(
            zeros(num_variates, size(inputs, :observations))
        )
    end
end

"""
    FakeTemplate{DistributionEstimate, SingleOutput}()

A [`Template`](@ref) whose [`Model`](@ref) will accept real value variables to predict a
univariate normal distribution (with zero mean and unit standard deviation) for each
observation.
"""
function FakeTemplate{DistributionEstimate, SingleOutput}()
    FakeTemplate{DistributionEstimate, SingleOutput}() do num_variates, inputs
        @assert(num_variates == 1, "$num_variates != 1")
        inputs = _handle_inputs(inputs)
        return NoncentralT.(3.0, zeros(size(inputs, :observations)))
    end
end

"""
    FakeTemplate{DistributionEstimate, MultiOutput}()

A [`Template`](@ref) whose [`Model`](@ref) will accept real value variables to predict a
multivariate normal distribution (with zero-vector mean and identity covariance matrix) for
each observation.
"""
function FakeTemplate{DistributionEstimate, MultiOutput}()
    FakeTemplate{DistributionEstimate, MultiOutput}() do num_variates, inputs
        std_dev = ones(num_variates)
        inputs = _handle_inputs(inputs)
        return [Product(Normal.(0, std_dev)) for _ in 1:size(inputs, 2)]
    end
end

"""
    _handle_inputs(inputs::AbstractMatrix)
    _handle_inputs(inputs::AbstractVector{<:Sampleable})

Process the inputs to `predict` appropriately depending on whether they are real valued or
distributions over input variables.
"""
function _handle_inputs(inputs::AbstractVector{<:Sampleable})
    return NamedDimsArray{(:features, :observations)}(hcat([mean(inputs[i]) for i in 1:size(inputs, 1)]...))
end

_handle_inputs(inputs::AbstractMatrix) = NamedDimsArray{(:features, :observations)}(inputs)

"""
    FakeModel

A fake Model for testing purposes. See [`FakeTemplate`](@ref) for details.
"""
mutable struct FakeModel{E<:EstimateTrait, O<:OutputTrait} <: Model
    predictor::Function
    num_variates::Int
end

Models.estimate_type(::Type{<:FakeModel{E, O}}) where {E, O} = E
Models.output_type(::Type{<:FakeModel{E, O}}) where {E, O} = O

Models.estimate_type(::Type{<:FakeTemplate{E, O}}) where {E, O} = E
Models.output_type(::Type{<:FakeTemplate{E, O}}) where {E, O} = O

function StatsBase.fit(
    template::FakeTemplate{E, O},
    outputs,
    inputs,
    weights=uweights(Float32, size(outputs, 2))
) where {E, O}
    outputs = NamedDimsArray{(:variates, :observations)}(outputs)
    num_variates = size(outputs, :variates)
    return FakeModel{E, O}(template.predictor, num_variates)
end

StatsBase.predict(m::FakeModel, inputs::AbstractMatrix) = m.predictor(m.num_variates, inputs)
StatsBase.predict(m::FakeModel, inputs::AbstractVector{<:Sampleable}) = m.predictor(m.num_variates, inputs)

"""
    test_interface(template::Template; inputs=rand(5, 5), outputs=rand(5, 5))

Test that subtypes of [`Template`](@ref) and [`Model`](@ref) implement the expected API.
Can be used as an initial test to verify the API has been correctly implemented.
"""
function test_interface(
    template::Template;
    inputs=rand(5,5),
    outputs=_default_outputs(template),
    distribution_inputs=[MvNormal(5, m) for m in 1:5],
    kwargs...
)
    @testset "Models API Interface Test: $(nameof(typeof(template)))" begin
        predictions = test_common(template, inputs, outputs)
        test_estimate_type(estimate_type(template), predictions, outputs)
        test_output_type(output_type(template), predictions, outputs)
        test_predict_input_type(predict_input_type(template), template, outputs, inputs, distribution_inputs)
    end
end

_default_outputs(template) = _default_outputs(output_type(template))
_default_outputs(::Type{SingleOutput}) = rand(1, 5)
_default_outputs(::Type{MultiOutput}) = rand(2, 5)

function test_estimate_type(::Type{PointEstimate}, predictions, outputs)
    @test predictions isa NamedDimsArray{(:variates, :observations), <:Real, 2}
    @test size(predictions) == size(outputs)
end

function test_estimate_type(::Type{DistributionEstimate}, predictions, outputs)
    @test predictions isa AbstractVector{<:ContinuousDistribution}
    @test length(predictions) == size(outputs, 2)
end

function test_output_type(::Type{SingleOutput}, predictions, outputs)
    @test all(length.(predictions) .== size(outputs, 1))
    @test all(length.(predictions) .== 1)
end

function test_output_type(::Type{MultiOutput}, predictions, outputs)
    if  eltype(predictions) <: Distribution
        @test all(length.(predictions) .== size(outputs, 1))
        @test all(length.(predictions) .> 1)
    else
        @test size(predictions, 1) == size(outputs, 1)
        @test size(predictions, 1) > 1
    end
end

function test_predict_input_type(::Type{PointPredictInput}, template, outputs, inputs, distribution_inputs)
    model = fit(template, outputs, inputs)
    @test hasmethod(Models.predict, (typeof(model), typeof(inputs)))
end

function test_predict_input_type(::Type{PointOrDistributionPredictInput}, template, outputs, inputs, distribution_inputs)
    model = fit(template, outputs, inputs)
    @test hasmethod(Models.predict, (typeof(model), typeof(distribution_inputs)))
    predictions = predict(model, distribution_inputs)
    test_estimate_type(estimate_type(template), predictions, outputs)
    test_output_type(output_type(template), predictions, outputs)
end

function test_names(template, model)
    template_type_name = string(nameof(typeof(template)))
    template_base_name_match = match(r"(.*)Template", template_type_name)
    @test template_base_name_match !== nothing  # must have Template suffix

    model_type_name = string(nameof(typeof(model)))
    model_base_name_match = match(r"(.*)Model", model_type_name)
    @test model_base_name_match !== nothing  # must have Model suffix

    # base_name must agreee
    @test model_base_name_match[1] == template_base_name_match[1]
end

function test_common(template, inputs, outputs)

    model = fit(template, outputs, inputs)

    @test template isa Template
    @test model isa Model

    @testset "type names" begin
        test_names(template, model)
    end

    @testset "test fit/predict errors" begin
        @test_throws MethodError predict(template, inputs)  # can only predict on a model
        @test_throws MethodError fit(model, outputs, inputs)  # can only fit a template
    end

    @testset "test weights can also be passed" begin
        weights = uweights(Float32, size(outputs, 2))
        model_weights = fit(template, outputs, inputs, weights)
    end

    @testset "traits" begin
        @test estimate_type(template) == estimate_type(model)
        @test output_type(template) == output_type(model)
        @test predict_input_type(template) == predict_input_type(model)
    end

    @testset "submodels" begin
        @test length(submodels(template)) == length(submodels(model))
        foreach(test_names, submodels(template), submodels(model))
    end

    predictions = predict(model, inputs)
    return predictions
end

end


text:
# Run this in another process:
#using WinControllerToFile;WinControllerToFile.subscribe()
using FileToRigidTransform
FileToRigidTransform.subscribe(logger= FileToRigidTransform.log_by_channel, timeout = 150)

#FileToRigidTransform.subscribe(timeout = 5)

#import FileToRigidTransform.Bytevector
#import FileToRigidTransform.DevState
#v = UInt8[0,127, 255]
#bv = Bytevector(v)
#sv = String(v)
#cu = codeunits(sv)

#drophighbits(1, UInt8(255))


text:
module aoc_22
    using DataStructures

    export read_input, play!, score, part1

    player_r = r"Player (\d):"

    const Decks = Tuple{Deque{Int},Deque{Int}}

    function read_input(filename = "input.txt")
        decks = (Deque{Int}(), Deque{Int}())
        lines = readlines(filename)
        player = 0
        for line in lines
            if isempty(line)
                player = 0
                continue
            end
            if player == 0
                m = match(player_r ,line)
                player = parse(Int, m.captures[1])
                continue
            end
            push!( decks[ player ], parse( Int, line ) )
        end
        decks
    end

    function play!( decks::Decks = read_input() )
        while !any( isempty.(decks) )
            cards = popfirst!.(decks)
            if(>( cards ... ))
                push!.( (decks[1],) , cards )
            else
                push!.( (decks[2],) , reverse(cards) )
            end
            # @info "Turn"
            # println( decks[1] )
            # println( decks[2] )
        end
        decks
    end

    function score( decks::Decks )
        winner = decks[ [ .~isempty.(decks)... ] ][1]
        sum(winner .* (length(winner):-1:1))
    end

    part1( input = read_input() ) = score( play!() )

    # Part 2
    # const previous_decks = Set{ Decks }()
    # game_counter = 1
    # recurrence = false

    function play2!( decks::Decks = read_input(), game_counter = Ref(0) )
        # previous_decks = Vector{ Decks }()
        previous_decks = Vector{ UInt64 }()

        game_counter[] += 1
        this_game = game_counter[]
        #if game_counter[] % 1000 == 0
        #    @info "Game $(game_counter[])"
        #end

        #    @info "Start"
        #    println( decks[1] )
        #    println( decks[2] )

        round_counter = 1

        while !any( isempty.(decks) )
            # @info "Round $round_counter of Game $this_game"
            # println( decks[1] )
            # println( decks[2] )
            decks_hash = hash(decks)
            if decks_hash in previous_decks
                # @info "Player 1 wins game $this_game due to recurrence!"
                # global recurrence = true
                empty!( decks[2] )
                return decks
            else
                push!( previous_decks, decks_hash )
            end
            cards = popfirst!.(decks)

            decks_longer_than_cards = length.(decks) .>= cards

            if( all( decks_longer_than_cards) )
                # @info "Subgame!"
                sub_deck = (Deque{Int}(), Deque{Int}())
                push!.( (sub_deck[1],), collect(decks[1])[1:cards[1]] )
                push!.( (sub_deck[2],), collect(decks[2])[1:cards[2]] )
                subgame_decks = play2!( sub_deck, game_counter )
                # if recurrence
                #    empty!( decks[2] )
                #    return decks
                # end
                winner = findfirst( .~isempty.(subgame_decks) )
                if winner == 2
                    cards = reverse(cards)
                end
                push!.( (decks[ winner ],), cards )
                # @info "Player $winner wins round $round_counter in game $this_game"
                continue
            end

            if(>( cards ... ))
                push!.( (decks[1],) , cards )
                # @info "Player 1 wins round in game $this_game"
            else
                push!.( (decks[2],) , reverse(cards) )
                # @info "Player 2 wins round game $this_game"
            end
            # println( decks[1] )
            # println( decks[2] )
            # sleep( 1 )
            round_counter += 1
        end
        return decks
    end

    function part2( input = read_input() )
        # global recurrence = false
        # empty!(previous_decks)
        game_counter = Ref(0)
        s = score( play2!( input, game_counter ) )
        # @info "Game counter: $(game_counter[])"
        return s
        # 30498
    end

    function __init__()
        input = read_input()
        println( part1( input ) )
        println( part2( input ) )
    end
end


text:
export EventChannel

@doc (@doc EventChannel)
mutable struct EventChannel <: GphysChannel
  id    ::String                # id
  name  ::String                # name
  loc   ::InstrumentPosition    # loc
  fs    ::Float64               # fs
  gain  ::Float64               # gain
  resp  ::InstrumentResponse    # resp
  units ::String                # units
  az    ::Float64               # source azimuth
  baz   ::Float64               # backazimuth
  dist  ::Float64               # distance
  pha   ::PhaseCat              # phase catalog
  src   ::String                # src
  misc  ::Dict{String,Any}      # misc
  notes ::Array{String,1}       # notes
  t     ::Array{Int64,2}        # time
  x     ::FloatArray            # data

  function EventChannel(
                        id    ::String,               # id
                        name  ::String,               # name
                        loc   ::InstrumentPosition,   # loc
                        fs    ::Float64,              # fs
                        gain  ::Float64,              # gain
                        resp  ::InstrumentResponse,   # resp
                        units ::String,               # units
                        az    ::Float64,              # source azimuth
                        baz   ::Float64,              # backazimuth
                        dist  ::Float64,              # distance
                        pha   ::PhaseCat,             # phase catalog
                        src   ::String,               # src
                        misc  ::Dict{String,Any},     # misc
                        notes ::Array{String,1},      # notes
                        t     ::Array{Int64,2},       # time
                        x     ::FloatArray            # data
                        )
      return new(id, name, loc, fs, gain, resp, units, az, baz, dist, pha, src, misc, notes, t, x)
    end
end

EventChannel(;
  id    ::String                = "",
  name  ::String                = "",
  loc   ::InstrumentPosition    = GeoLoc(),
  fs    ::Float64               = zero(Float64),
  gain  ::Float64               = one(Float64),
  resp  ::InstrumentResponse    = PZResp(),
  units ::String                = "",
  az    ::Float64               = zero(Float64),    # source azimuth
  baz   ::Float64               = zero(Float64),    # backazimuth
  dist  ::Float64               = zero(Float64),    # distance
  pha   ::PhaseCat              = PhaseCat(),
  src   ::String                = "",
  misc  ::Dict{String,Any}      = Dict{String,Any}(),
  notes ::Array{String,1}       = Array{String,1}(undef, 0),
  t     ::Array{Int64,2}        = Array{Int64,2}(undef, 0, 2),
  x     ::FloatArray            = Array{Float32,1}(undef, 0)
  ) = EventChannel(id, name, loc, fs, gain, resp, units, az, baz, dist, pha, src, misc, notes, t, x)

function getindex(S::EventTraceData, j::Int)
  C = EventChannel()
  [setfield!(C, f, getfield(S,f)[j]) for f in tracefields]
  return C
end
setindex!(S::EventTraceData, C::EventChannel, j::Int) = (
  [(getfield(S, f))[j] = getfield(C, f) for f in tracefields];
  return S)

function isempty(Ch::EventChannel)
  q::Bool = Ch.gain == 1.0
  for f in (:az, :baz, :dist, :fs)
    q = min(q, getfield(Ch, f) == 0.0)
    (q == false) && return q
  end
  for f in (:id, :loc, :misc, :name, :notes, :pha, :resp, :src, :t, :units, :x)
    q = min(q, isempty(getfield(Ch, f)))
    (q == false) && return q
  end
  return q
end

# ============================================================================
# Conversion and push to EventTraceData
function EventTraceData(C::EventChannel)
 S = EventTraceData(1)
 for f in tracefields
   setindex!(getfield(S, f), getfield(C, f), 1)
 end
 return S
end

function push!(S::EventTraceData, C::EventChannel)
 for i in tracefields
   push!(getfield(S,i), getfield(C,i))
 end
 S.n += 1
 return nothing
end

function sizeof(C::EventChannel)
  s = 136
  for f in tracefields
    v = getfield(C,f)
    s += sizeof(v)
    if f == :notes
      if !isempty(v)
        s += sum([sizeof(j) for j in v])
      end
    elseif f == :misc || f == :pha
      k = collect(keys(v))
      s += sizeof(k) + 64 + sum([sizeof(j) for j in k])
      for p in values(v)
        s += sizeof(p)
        if typeof(p) == Array{String,1}
          s += sum([sizeof(j) for j in p])
        end
      end
    end
  end
  return s
end

function write(io::IO, S::EventChannel)
  write(io, Int64(sizeof(S.id)))
  write(io, S.id)                                                     # id
  write(io, Int64(sizeof(S.name)))
  write(io, S.name)                                                   # name
  write(io, loctyp2code(S.loc))
  write(io, S.loc)                                                    # loc
  write(io, S.fs)                                                     # fs
  write(io, S.gain)                                                   # gain
  write(io, resptyp2code(S.resp))
  write(io, S.resp)                                                   # resp
  write(io, Int64(sizeof(S.units)))
  write(io, S.units)                                                  # units
  write(io, S.az)                                                     # az
  write(io, S.baz)                                                    # baz
  write(io, S.dist)                                                   # dist
  write(io, S.pha)                                                    # pha
  write(io, Int64(sizeof(S.src)))
  write(io, S.src)                                                    # src
  write_misc(io, S.misc)                                              # misc
  write_string_vec(io, S.notes)                                       # notes
  write(io, Int64(size(S.t,1)))
  write(io, S.t)                                                      # t
  write(io, typ2code(eltype(S.x)))
  write(io, Int64(length(S.x)))
  write(io, S.x)                                                      # x
  return nothing
end

read(io::IO, ::Type{EventChannel}) = EventChannel(
  String(read(io, read(io, Int64))),                                    # id
  String(read(io, read(io, Int64))),                                    # name
  read(io, code2loctyp(read(io, UInt8))),                              # loc
  read(io, Float64),                                                    # fs
  read(io, Float64),                                                    # gain
  read(io, code2resptyp(read(io, UInt8))),                              # resp
  String(read(io, read(io, Int64))),                                    # units
  read(io, Float64),                                                    # az
  read(io, Float64),                                                    # baz
  read(io, Float64),                                                    # dist
  read(io, PhaseCat),                                                   # pha
  String(read(io, read(io, Int64))),                                    # src
  read_misc(io, getfield(BUF, :buf)),                                   # misc
  read_string_vec(io, getfield(BUF, :buf)),                             # notes
  read!(io, Array{Int64, 2}(undef, read(io, Int64), 2)),                # t
  read!(io, Array{code2typ(read(io,UInt8)),1}(undef, read(io, Int64))), # x
  )


text:
# This file seems to need to exist for Julia package manager to be happy
module CORBITS

export prob_of_transits_approx

if VERSION >= v"0.7.0-"
   using Libdl
   using Pkg
   using CORBITS
   global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[joinpath(dirname(pathof(CORBITS)),".."),"/usr/local/lib"] ) # WARNING: Assumes can find libcorbits.so
elseif VERSION >= v"0.4.0-"
   global const LIB_CORBITS = Libdl.find_library(["libcorbits.so"],[".",joinpath(Pkg.dir(),"ExoplanetsSysSim/"),joinpath(Pkg.dir(),"CORBITS/"),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
else
   global const LIB_CORBITS = find_library(["libcorbits.so"],[".",joinpath(Pkg.devdir(),"ExoplanetsSysSim/"),joinpath(Pkg.devdir(),"CORBITS/"),"/usr/local/lib"])  # WARNING: Assumes can find libcorbits.so
end

# Call CORBITS's function prob_of_transits_approx_arrays(a, r_star, r, e, Omega, omega, inc, use)
# Returns a Cdouble (aka Float64)
# For documentation see https://github.com/jbrakensiek/CORBITS
function prob_of_transits_approx(a::Vector{Cdouble},r_star::Cdouble,r::Vector{Cdouble}, e::Vector{Cdouble},
                                 Omega::Vector{Cdouble}, omega::Vector{Cdouble}, inc::Vector{Cdouble}, use::Vector{Cint})
  @assert(length(a) == length(r) == length(e) == length(Omega) == length(omega)
== length(inc) >= length(use) )
  @assert(length(use) >=1 )
  num_pl = length(use)
  return ccall( (:prob_of_transits_approx_arrays, LIB_CORBITS), Cdouble,
      (Ptr{Cdouble}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Cint),
       a, r_star, r, e, Omega, omega, inc, use, num_pl)
end

end # module


text:
module TM1

using Dates
using Base64

import HTTP, JSON, MbedTLS, Sockets

# potentially set some high level constants Here

include("utils/requests.jl")

include("utils/TM1Type.jl")

export name

include("utils/auth.jl")

export authenticate

#export # requests.jl -

# simple utility, not needing authentication

include("ping/ping.jl")

export ping

include("dimensions/elements.jl")
include("dimensions/hierarchies.jl")
include("dimensions/dimensions.jl")
include("dimensions/subsets.jl")

include("cubes/cubes.jl")

include("processes/processes.jl")
include("processes/chores.jl")

include("cellsets/cellsets.jl")

export Cube,
  get_all_cubes, get_all_control_cubes, get_all_model_cubes, delete_cube, create_cube

export Dimension,
  get_all_dimensions,
  get_all_model_dimensions,
  get_all_control_dimensions,
  get_dimension,
  delete_dimension

export Hierarchy, get_hierarchy, delete_hierarchy

export Element, get_element, delete_element

export Subset, get_subset, delete_subset

export Process, get_all_processes, get_process, delete_process

export Chore, get_all_chores, get_chore, delete_chore

export Cellset, get_cellset, delete_cellset

end #Module


text:
## Code for the CAMERAS #################################################################################################################

abstract type Renderer end

# Camera is the abstract type which the two different cameras are generated from.
abstract type Camera
end

"""
    OrthogonalCamera( aspect_ratio=1.0, transformation=Transformation())

It creates a **Orthogonal Camera**.

## Arguments:
- *a* -> aspect ratio;
- *T* -> generic Transformation.

When not specified in the constructor, a = 1 and T = Transformation(),
where *a* is the aspect ratio and *T* a generic transformation.
"""
struct OrthogonalCamera <: Camera
    aspect_ratio::Float32
    transformation::Transformation

    OrthogonalCamera( aspect_ratio=1.0f0, transformation=Transformation()) = new(aspect_ratio, transformation)
end

"""
#### Usage 1:
    fire_ray(camera, u, v)

It fires a ray from a camera (Orthogonal or Perspective) directed to a pixel with coordinates *u* and *v* (on the screen).

## Usage 2:
    fire_ray(im, col, row, u_pixel, v_pixel)

It fires a ray from a camera (contained in *im*) directed to (col, row)-pixel. It hits the pixel in hte point with coordinates
*u_pixel*, *v_pixel*.

## Arguments:
- *im* -> object of type ImageTracer;
- *col* & *row* -> integers for the coordinates of the pixel in the raster image;
- *u_pixel* and *v_pixel* -> integers for the coordinates inside the pixel.

If not specified *u_pixel* = *v_pixel* = 0.5.

"""
function fire_ray( camera::OrthogonalCamera, u, v)
    Ray_StandardFrame = Ray( Point(-1.0f0, (1.0f0-2*u)*camera.aspect_ratio, 2*v-1), Vec(1.0f0, 0.0f0, 0.0f0))
    return (camera.transformation * Ray_StandardFrame)
end

"""
    PerspectiveCamera(a, T, d)

It creates a **Perspective Camera**.

## Arguments:
- *a* -> aspect ratio;
- *T* -> generic Transformation;
- *d* -> distance between the observer and the screen.

When not specified in the constructor, a = 1 and T = Transformation(),
where *a* is the aspect ratio and *T* a generic transformation.
"""
struct PerspectiveCamera <: Camera
    aspect_ratio::Float32
    transformation::Transformation
    distance::Float32

    PerspectiveCamera(aspect_ratio=1.0f0, transformation=Transformation(), distance=1.0f0) = new(aspect_ratio, transformation, distance )
#    PerspectiveCamera(transformation=Transformation(), distance=1.0) = new(1.0, transformation, distance )
end

function fire_ray(camera::PerspectiveCamera, u, v)
    Ray_StandardFrame = Ray( Point(-camera.distance, 0.0f0, 0.0), Vec(camera.distance, (1.0f0-2*u)*camera.aspect_ratio, 2*v-1))
    return (camera.transformation * Ray_StandardFrame)
end

## Defining IMAGETRACER and its methods: fire_ray and fire_all_rays: ##################################################################à

"""
    ImageTracer()

It creates a **ImageTracer**.

## Arguments:
- image -> object of type HdrImage;
- camera;
- ray_per_side^2 is the number of rays fired in each pixel to reduce aliasing;

"""
struct ImageTracer
    image::HdrImage
    camera::Camera
    ray_per_side::Int
    pcg::PCG

    ImageTracer(image::HdrImage,
                camera::Camera,
                ray_per_side::Int=1,
                pcg::PCG=PCG()) =
                new(image, camera, ray_per_side, pcg)
end

function fire_ray(im::ImageTracer, col, row, u_pixel=0.5f0, v_pixel=0.5f0)
    ray_vector = []
    for i in 1:im.ray_per_side
        for j in 1:im.ray_per_side
            if im.ray_per_side > 1 # ANTIALIASING!!
                u_pixel = (j - 1 + pcg_randf(im.pcg)) / im.ray_per_side
                v_pixel = (i - 1 + pcg_randf(im.pcg)) / im.ray_per_side
            end
            u = Float32((col -1 + u_pixel)/(im.image.width))
            v = Float32(1.0f0 - (row -1 + v_pixel)/(im.image.height))
            ray = fire_ray(im.camera, u, v)
            push!(ray_vector, ray)
        end
    end

    return ray_vector
end

"""
    fire_all_rays(im, func)

It fires all rays, requiring a ImageTracer and a generic function (to assign colors to the pixels).
"""
function fire_all_rays(im::ImageTracer, func, renderer::Renderer)
    temp = -2
    for row ∈ 1:im.image.height
        for col ∈ 1:im.image.width
            cum_color = RGB(0.f0,0.f0,0.f0)
            rays_in_pixel= fire_ray(im, col, row)

            for ray in rays_in_pixel
                cum_color += func(ray, renderer)
            end
            im.image.pixels[get_pixel(im.image, col, row)] = cum_color * (1.f0 / im.ray_per_side^2)
        end

        # mini-script for printing the percentage of rays computed
        percentage= convert(Int32 ,floor(100*(row-1)/im.image.height)) ##print progress
        if percentage == temp+2
            i = convert(Int32 ,floor(percentage/2))
            print("\rComputed $(percentage)% of pixels [$("#"^i)$("."^(49-i))]")
            temp = percentage
        end

    end
    print("\rComputed 100% of pixels [$("#"^50)]")
end

function fire_all_rays(im::ImageTracer, func)
    for row ∈ 1:im.image.height
        for col ∈ 1:im.image.width
            ray = fire_ray(im, col, row)
            color = func(ray)
            im.image.pixels[get_pixel(im.image, col, row)] = color
        end
    end
end


text:
include("../../testcases.jl")
using PyCall
using DFTK

atoms = [ElementPsp(:Al, psp=load_psp("hgh/lda/al-q3")) => aluminium_primitive.positions]
model = model_LDA(aluminium_primitive.lattice, atoms, temperature=0.01,
                  smearing=Smearing.MethfesselPaxton2(), spin_polarization=:collinear)

abinitpseudos = [joinpath(pyimport("abipy.data").__path__[end], "hgh_pseudos/13al.3.hgh")]
DFTK.run_abinit_scf(model, @__DIR__; abinitpseudos=abinitpseudos,
                    Ecut=15, kgrid=[4, 4, 4], n_bands=8, tol=1e-10,
                    iscf=3)  # Use Anderson mixing instead of minimization


text:
tiny(x::Float64) = 0x1p-1023
huge(x::Float64) = 0x1p1023

# scalbn(x,n) = x * 2^n
function scalbn(x::Float64, n::Int)
    xu = reinterpret(Uint64,x)
    k = int(xu >> 52) & 0x07ff
    if k == 0 # x is subnormal
        x == zero(x) && return x
        x *= 0x1p54 # normalise significand
        xu = reinterpret(Uint64,x)
        k = int(xu >> 52) & 0x07ff
        k -= 54
        n < -50000 && return tiny(x)*x # so exponent doesn't underflow
    end
    k == 0x07ff && return x + x # NaN or Inf
    k += n
    if k > 0x07fe || n > 50000  # overflow result
        return huge(x)*copysign(huge(x),x)
    elseif k > 0 # normal result
        xu = (xu & 0x800f_ffff_ffff_ffff) | (int64(k) << 52)
        return reinterpret(Float64,xu)
    elseif k > -54 # subnormal result
        xu = (xu & 0x800f_ffff_ffff_ffff) | (int64(k+54) << 52)
        return reinterpret(Float64,xu)*0x1p-54
    else # underlow result
        return tiny(x)*copysign(tiny(x),x)
    end
end


text:
Random.seed!(42)
a = rand(SymmetricTensor{Float64, 2}, 3)
b = rand(SymmetricTensor{Float64, 3}, 3)
@everywhere using LinearAlgebra
@everywhere testf(a,b,bool)= det(a[bool,bool])
println("nworkers()")
println(nworkers())

@testset "greedesearch parallel implementation" begin
  g = greedesearchdata(a,b, testf, 3)

  if VERSION <= v"1.7"
    @test g[1][1] == [true, false, true]
    @test g[2][1] == [false, false, true]
    @test g[3][1] == [false, false, false]
    @test g[1][2] ≈ 0.48918301293211774
    @test g[2][2] ≈ 0.9735659798036858
    @test g[3][2] == 1.0
    @test g[1][3] == 2
    @test g[2][3] == 1
    @test g[3][3] == 3
  else
    @test g[1][1] == [true, true, false]
    @test g[2][1] == [true, false, false]
    @test g[3][1] == [false, false, false]
    @test g[1][2] ≈ 0.09764869605558585
    @test g[2][2] ≈ 0.6293451231426089
    @test g[3][2] == 1.0
    @test g[1][3] == 3
    @test g[2][3] == 2
    @test g[3][3] == 1
  end
end

Random.seed!(42)

@testset "cum2mat parallel tests on cumulants" begin
    @testset "small data sie" begin
        x = rand(200,12)
        c = cumulants(x,5,2)
        s = cum2mat(c[3])
        X = unfold(Array(c[3]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

        s = cum2mat(c[4])
        X = unfold(Array(c[4]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

        s = cum2mat(c[5])
        X = unfold(Array(c[5]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end

    @testset "larger data size" begin
        x = rand(200,27)
        c = cumulants(x,5,5)
        s = cum2mat(c[3])
        X = unfold(Array(c[3]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

        s = cum2mat(c[4])
        X = unfold(Array(c[4]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)

        s = cum2mat(c[5])
        X = unfold(Array(c[5]), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-10)
    end
end

@testset "cum2mat parallel tests on random tensors" begin
    @testset "m = 5, n = 24, d = 4" begin
        Random.seed!(42)
        t = rand(SymmetricTensor{Float64, 5}, 24, 4)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-8)
    end

    @testset "m = 5, n = 37, d = 6" begin
        t = rand(SymmetricTensor{Float64, 5}, 37, 6)
        s = cum2mat(t)
        X = unfold(Array(t), 1)
        M = X*transpose(X)
        @test maximum(abs.(M - Array(s))) ≈ 0 atol = 10^(-8)
    end

end


text:
function find_and_append_deps!(idxs, nodes, i)
    node = nodes[i]
    if node.type isa AbstractStructNodeType
        push!(idxs, i)
    elseif node.type isa AbstractTypedefNodeType
        if isempty(node.adj)
            return nothing
        else
            @assert length(node.adj) == 1
            find_and_append_deps!(idxs, nodes, first(node.adj))
        end
    end
    return nothing
end

function is_argtype_dep(struct_node, nodes, i)
    node = nodes[i]
    if node.type isa AbstractStructNodeType
        return struct_node == node
    elseif node.type isa AbstractTypedefNodeType
        if isempty(node.adj)
            return false
        else
            @assert length(node.adj) == 1
            return is_argtype_dep(struct_node, nodes, first(node.adj))
        end
    end
    return false
end

function should_tweak(nodes, idx)
    struct_node = nodes[idx]
    is_orphan = true
    is_all_non_pointer = true
    is_all_no_integer_pointer = true
    for (i, node) in enumerate(nodes)
        # only consider FunctionProto for now
        node.type isa FunctionProto || continue

        # do a quick-check
        is_dep = false
        for j in node.adj
            if is_argtype_dep(struct_node, nodes, j)
                is_dep = true
                break
            end
        end

        if is_dep
            is_orphan && (is_orphan = false;)
            cursor = node.cursor
            types = [getArgType(getCursorType(cursor), i - 1) for i in 1:getNumArguments(cursor)]
            has_integer_args = false
            has_a_pointer_dep = false

            for ty in types
                jlty = tojulia(ty)
                leaf_ty = get_jl_leaf_type(jlty)

                if is_jl_basic(leaf_ty)
                    is_jl_integer(jlty) && (has_integer_args = true;)
                    continue
                end

                if leaf_ty isa JuliaCtypedef && haskey(EXTRA_DEFINITIONS, leaf_ty.sym)
                    if is_jl_integer(EXTRA_DEFINITIONS[leaf_ty.sym])
                        has_integer_args = true
                    end
                end

                expr = first(filter(x->Meta.isexpr(x, :struct), struct_node.exprs))
                struct_name = string(expr.args[2])
                if string(leaf_ty.sym) == struct_name
                    is_jl_pointer(jlty) && (has_a_pointer_dep = true;)
                end
            end

            if is_all_non_pointer && has_a_pointer_dep
                is_all_non_pointer = false
            end

            if is_all_no_integer_pointer && has_a_pointer_dep && has_integer_args
                is_all_no_integer_pointer = false
            end
        end
    end

    return is_orphan || is_all_non_pointer || is_all_no_integer_pointer
end