text:
module PolyhedralRelaxations

using DataStructures
using JuMP

import ForwardDiff
import Memento

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via
# getlogger(PolyhedralRelaxations)`.
# NOTE: If this line is not included then the precompiled `PolyhedralRelaxations.LOGGER` won't
# be registered at runtime.
__init__() = Memento.register(_LOGGER)

"""
Suppresses information and warning messages output by PolyhedralRelaxations. For fine grained
control use the Memento package.
"""
function silence()
    Memento.info(
        _LOGGER,
        string(
            "Suppressing information and warning messages for the rest ",
            "of this session.  Use the Memento package for more fine-grained control of logging.",
        ),
    )
    Memento.setlevel!(Memento.getlogger(PolyhedralRelaxations), "error")
end

"Allows users to set the logging level without adding Memento."
function logger_config!(level)
    Memento.config!(Memento.getlogger("PolyhedralRelaxations"), level)
end

const EPS = 1e-6
const INF = 1e12

include("common.jl")
include("univariate_lp_relaxation.jl")
include("univariate_milp_relaxation.jl")
include("bilinear_relaxation.jl")
include("api.jl")

end


text:
using FixedEffects, Random, Statistics
Random.seed!(1234)
N = 10000000
K = 100
id1 = rand(1:div(N, K), N)
id2 = rand(1:K, N)
fes = [FixedEffect(id1), FixedEffect(id2)]
x = rand(N)

# simple problem
@time solve_residuals!(deepcopy(x), fes)

ext: jl

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects


text:
#
# Device functionality
#

## execution

struct ROCArrayBackend <: AbstractGPUBackend end

struct ROCKernelContext <: AbstractKernelContext end

function GPUArrays.gpu_call(::ROCArrayBackend, f, args, threads::Int, blocks::Int;
                            name::Union{String,Nothing})
    groupsize, gridsize = threads, blocks*threads
    wait(@roc groupsize=groupsize gridsize=gridsize f(ROCKernelContext(), args...))
end
function GPUArrays.gpu_call(::ROCArrayBackend, f, args; elements::Int, name::Union{String,Nothing}=nothing)
    wait(@roc groupsize=min(elements, 64) gridsize=elements f(ROCKernelContext(), args...))
end

## on-device

# indexing

for (f, froc) in (
        (:blockidx, :blockIdx),
        (:blockdim, :blockDim),
        (:threadidx, :threadIdx),
        (:griddim, :gridDimWG)
    )
    @eval @inline GPUArrays.$f(::ROCKernelContext) = AMDGPU.$froc().x
end

# math

@inline GPUArrays.cos(ctx::ROCKernelContext, x) = cos(x)
@inline GPUArrays.sin(ctx::ROCKernelContext, x) = sin(x)
@inline GPUArrays.sqrt(ctx::ROCKernelContext, x) = sqrt(x)
@inline GPUArrays.log(ctx::ROCKernelContext, x) = log(x)

# memory

@inline function GPUArrays.LocalMemory(ctx::ROCKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}) where {T,dims,id}
    ptr = AMDGPU.alloc_special(Val(id), T, AMDGPU.AS.Local, Val(prod(dims)))
    ROCDeviceArray(dims, ptr)
end

# synchronization

@inline function GPUArrays.synchronize_threads(::ROCKernelContext)
    AMDGPU.sync_workgroup()
    return
end

#
# Host abstractions
#

mutable struct ROCArray{T,N} <: AbstractGPUArray{T,N}
    buf::Mem.Buffer
    own::Bool

    dims::Dims{N}
    offset::Int

    syncstate::SyncState

    function ROCArray{T,N}(buf::Mem.Buffer, dims::Dims{N}; offset::Integer=0, own::Bool=true) where {T,N}
        @assert isbitstype(T) "ROCArray only supports bits types"
        xs = new{T,N}(buf, own, dims, offset, SyncState())
        if own
            hsaref!()
            Mem.retain(buf)
            finalizer(unsafe_free!, xs)
        end
        return xs
    end
end

function unsafe_free!(xs::ROCArray)
    Mem.release(xs.buf) && Mem.free(xs.buf)
    hsaunref!()
    return
end

wait!(x::ROCArray) = wait!(x.syncstate)
mark!(x::ROCArray, s) = mark!(x.syncstate, s)
wait!(xs::Vector{<:ROCArray}) = foreach(wait!, xs)
mark!(xs::Vector{<:ROCArray}, s) = foreach(x->mark!(x,s), xs)
wait!(xs::NTuple{N,<:ROCArray} where N) = foreach(wait!, xs)
mark!(xs::NTuple{N,<:ROCArray} where N, s) = foreach(x->mark!(x,s), xs)
function Adapt.adapt_storage(::WaitAdaptor, x::ROCArray)
    wait!(x.syncstate)
    x
end
function Adapt.adapt_storage(ma::MarkAdaptor, x::ROCArray)
    mark!(x.syncstate, ma.s)
    x
end

## aliases

const ROCVector{T} = ROCArray{T,1}
const ROCMatrix{T} = ROCArray{T,2}
const ROCVecOrMat{T} = Union{ROCVector{T},ROCMatrix{T}}

# anything that's (secretly) backed by a ROCArray
AnyROCArray{T,N} = Union{ROCArray{T,N}, WrappedArray{T,N,ROCArray,ROCArray{T,N}}}
AnyROCVector{T} = AnyROCArray{T,1}
AnyROCMatrix{T} = AnyROCArray{T,2}
AnyROCVecOrMat{T} = Union{AnyROCVector{T}, AnyROCMatrix{T}}

## constructors

# type and dimensionality specified, accepting dims as tuples of Ints
function ROCArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
    buf = Mem.alloc(prod(dims)*sizeof(T))
    ROCArray{T,N}(buf, dims)
end

# type and dimensionality specified, accepting dims as series of Ints
ROCArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = ROCArray{T,N}(undef, dims)

# type but not dimensionality specified
ROCArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
ROCArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    ROCArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# from Base arrays
function ROCArray{T,N}(x::Array{T,N}, dims::Dims{N}) where {T,N}
    r = ROCArray{T,N}(undef, size(x))
    Mem.upload!(r.buf, pointer(x), sizeof(x))
    return r
end

# type as first argument
# FIXME: Remove me!
#ROCArray(::Type{T}, dims::Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)

# empty vector constructor
ROCArray{T,1}() where {T} = ROCArray{T,1}(undef, 0)

Base.similar(a::ROCArray{T,N}) where {T,N} = ROCArray{T,N}(undef, size(a))
Base.similar(a::ROCArray{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
Base.similar(a::ROCArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)

## array interface

Base.elsize(::Type{<:ROCArray{T}}) where {T} = sizeof(T)

Base.size(x::ROCArray) = x.dims
Base.sizeof(x::ROCArray) = Base.elsize(x) * length(x)

## interop with Julia arrays

ROCArray{T,N}(x::AbstractArray{S,N}) where {T,N,S} =
    ROCArray{T,N}(convert(Array{T}, x), size(x))

# underspecified constructors
ROCArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = ROCArray{T,N}(xs)
(::Type{ROCArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = ROCArray{S,N}(x)
ROCArray(A::AbstractArray{T,N}) where {T,N} = ROCArray{T,N}(A)

# idempotency
ROCArray{T,N}(xs::ROCArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: ROCArray = x

## memory operations

function Base.copyto!(dest::Array{T}, d_offset::Integer,
                      source::ROCArray{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    wait!(source)
    Mem.download!(pointer(dest, d_offset),
                  Mem.view(source.buf, source.offset + (s_offset-1)*sizeof(T)),
                  amount*sizeof(T))
    dest
end
function Base.copyto!(dest::ROCArray{T}, d_offset::Integer,
                      source::Array{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    wait!(dest)
    Mem.upload!(Mem.view(dest.buf, dest.offset + (d_offset-1)*sizeof(T)),
                pointer(source, s_offset),
                amount*sizeof(T))
    dest
end
function Base.copyto!(dest::ROCArray{T}, d_offset::Integer,
                      source::ROCArray{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    wait!(dest)
    wait!(source)
    Mem.transfer!(Mem.view(dest.buf, dest.offset + (d_offset-1)*sizeof(T)),
                  Mem.view(source.buf, source.offset + (s_offset-1)*sizeof(T)),
                  amount*sizeof(T))
    dest
end

# TODO: Workaround for hanging copy() broadcast kernel
function Base.copy(X::ROCArray{T}) where T
    Xnew = ROCArray{T}(undef, size(X))
    copyto!(Xnew, 1, X, 1, length(X))
    Xnew
end

## views

# optimize view to return a ROCArray when contiguous

struct Contiguous end
struct NonContiguous end

# NOTE: this covers more cases than the I<:... in Base.FastContiguousSubArray
ROCIndexStyle() = Contiguous()
ROCIndexStyle(I...) = NonContiguous()
ROCIndexStyle(::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::Colon, ::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::AbstractUnitRange, ::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::Colon, I...) = ROCIndexStyle(I...)

rocviewlength() = ()
@inline rocviewlength(::Real, I...) = rocviewlength(I...) # skip scalar
@inline rocviewlength(i1::AbstractUnitRange, I...) = (Base.unsafe_length(i1), rocviewlength(I...)...)
@inline rocviewlength(i1::AbstractUnitRange, ::Base.ScalarIndex...) = (Base.unsafe_length(i1),)

@inline function Base.view(A::ROCArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(Array, j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(ROCArray, j), J)
    unsafe_view(A, J_gpu, ROCIndexStyle(I...))
end

@inline function unsafe_view(A, I, ::Contiguous)
    unsafe_contiguous_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I, rocviewlength(I...))
end
@inline function unsafe_contiguous_view(a::ROCArray{T}, I::NTuple{N,Base.ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M}
    offset = Base.compute_offset1(a, 1, I) * sizeof(T)

    Mem.retain(a.buf)
    b = ROCArray{T,M}(a.buf, dims, offset=a.offset+offset, own=false)
    finalizer(unsafe_free!, b)
    return b
end

@inline function unsafe_view(A, I, ::NonContiguous)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I...)
end

#FIXME: add pointer conversions

## reshape

# optimize reshape to return a ROCArray

function Base.reshape(a::ROCArray{T,M}, dims::NTuple{N,Int}) where {T,N,M}
    if prod(dims) != length(a)
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $len"))
    end

    if N == M && dims == size(a)
        return a
    end

    Mem.retain(a.buf)
    b = ROCArray{T,N}(a.buf, dims, offset=a.offset, own=false)
    finalizer(unsafe_free!, b)
    return b
end

# allow missing dimensions with Colon()
if VERSION < v"1.6.0-DEV.1358"
    Base.reshape(parent::ROCArray, dims::Tuple{Vararg{Union{Int,Colon}}}) =
    Base.reshape(parent, Base._reshape_uncolon(parent, dims))
end

## fft

#=
using AbstractFFTs

# defining our own plan type is the easiest way to pass around the plans in FFTW interface
# without ambiguities

struct FFTPlan{T}
    p::T
end

AbstractFFTs.plan_fft(A::ROCArray; kw_args...) = FFTPlan(plan_fft(A.data; kw_args...))
AbstractFFTs.plan_fft!(A::ROCArray; kw_args...) = FFTPlan(plan_fft!(A.data; kw_args...))
AbstractFFTs.plan_bfft!(A::ROCArray; kw_args...) = FFTPlan(plan_bfft!(A.data; kw_args...))
AbstractFFTs.plan_bfft(A::ROCArray; kw_args...) = FFTPlan(plan_bfft(A.data; kw_args...))
AbstractFFTs.plan_ifft!(A::ROCArray; kw_args...) = FFTPlan(plan_ifft!(A.data; kw_args...))
AbstractFFTs.plan_ifft(A::ROCArray; kw_args...) = FFTPlan(plan_ifft(A.data; kw_args...))

function Base.:(*)(plan::FFTPlan, A::ROCArray)
    x = plan.p * A.data
    ROCArray(x)
end
=#

## GPUArrays interfaces

GPUArrays.device(x::ROCArray) = x.buf.agent

GPUArrays.backend(::Type{<:ROCArray}) = ROCArrayBackend()

function Base.convert(::Type{ROCDeviceArray{T,N,AS.Global}}, a::ROCArray{T,N}) where {T,N}
    ptr = Base.unsafe_convert(Ptr{T}, a.buf)
    ROCDeviceArray{T,N,AS.Global}(a.dims, AMDGPU.LLVMPtr{T,AS.Global}(ptr+a.offset))
end
Adapt.adapt_storage(::AMDGPU.Adaptor, x::ROCArray{T,N}) where {T,N} =
    convert(ROCDeviceArray{T,N,AS.Global}, x)

## interop with CPU arrays

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{ROCArray}, xs::AT) where {AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(ROCArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:ROCArray{T}}, xs::AT) where {T, AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(ROCArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::ROCArray) = convert(Array, xs)

## Float32-preferring conversion

struct Float32Adaptor end

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray) =
  isbits(xs) ? xs : convert(ROCArray, xs)

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{<:AbstractFloat}) =
  isbits(xs) ? xs : convert(ROCArray{Float32}, xs)

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{<:Complex{<:AbstractFloat}}) =
  isbits(xs) ? xs : convert(ROCArray{ComplexF32}, xs)

# not for Float16
Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{Float16}) =
  isbits(xs) ? xs : convert(ROCArray, xs)

roc(xs) = adapt(Float32Adaptor(), xs)

Base.unsafe_convert(::Type{Ptr{T}}, x::ROCArray{T}) where T =
    Base.unsafe_convert(Ptr{T}, x.buf)

# some nice utilities

ones(dims...) = ones(Float32, dims...)
ones(T::Type, dims...) = fill!(ROCArray{T}(undef, dims...), one(T))
zeros(dims...) = zeros(Float32, dims...)
zeros(T::Type, dims...) = fill!(ROCArray{T}(undef, dims...), zero(T))

# create a derived array (reinterpreted or reshaped) that's still a ROCArray
# TODO: Move this to GPUArrays?
@inline function _derived_array(::Type{T}, N::Int, a::ROCArray, osize::Dims) where {T}
  Mem.retain(a.buf)
  b = ROCArray{T,N}(a.buf, osize, offset=a.offset, own=false)
  finalizer(unsafe_free!, b)
  return b
end

## reinterpret

function Base.reinterpret(::Type{T}, a::ROCArray{S,N}) where {T,S,N}
  err = _reinterpret_exception(T, a)
  err === nothing || throw(err)

  if sizeof(T) == sizeof(S) # for N == 0
    osize = size(a)
  else
    isize = size(a)
    size1 = div(isize[1]*sizeof(S), sizeof(T))
    osize = tuple(size1, Base.tail(isize)...)
  end

  return _derived_array(T, N, a, osize)
end

function _reinterpret_exception(::Type{T}, a::AbstractArray{S,N}) where {T,S,N}
  if !isbitstype(T) || !isbitstype(S)
    return _ROCReinterpretBitsTypeError{T,typeof(a)}()
  end
  if N == 0 && sizeof(T) != sizeof(S)
    return _ROCReinterpretZeroDimError{T,typeof(a)}()
  end
  if N != 0 && sizeof(S) != sizeof(T)
      ax1 = axes(a)[1]
      dim = length(ax1)
      if Base.rem(dim*sizeof(S),sizeof(T)) != 0
        return _ROCReinterpretDivisibilityError{T,typeof(a)}(dim)
      end
      if first(ax1) != 1
        return _ROCReinterpretFirstIndexError{T,typeof(a),typeof(ax1)}(ax1)
      end
  end
  return nothing
end

struct _ROCReinterpretBitsTypeError{T,A} <: Exception end
function Base.showerror(io::IO, ::_ROCReinterpretBitsTypeError{T, <:AbstractArray{S}}) where {T, S}
  print(io, "cannot reinterpret an `$(S)` array to `$(T)`, because not all types are bitstypes")
end

struct _ROCReinterpretZeroDimError{T,A} <: Exception end
function Base.showerror(io::IO, ::_ROCReinterpretZeroDimError{T, <:AbstractArray{S,N}}) where {T, S, N}
  print(io, "cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a different size")
end

struct _ROCReinterpretDivisibilityError{T,A} <: Exception
  dim::Int
end
function Base.showerror(io::IO, err::_ROCReinterpretDivisibilityError{T, <:AbstractArray{S,N}}) where {T, S, N}
  dim = err.dim
  print(io, """
      cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
      The resulting array would have non-integral first dimension.
      """)
end

struct _ROCReinterpretFirstIndexError{T,A,Ax1} <: Exception
  ax1::Ax1
end
function Base.showerror(io::IO, err::_ROCReinterpretFirstIndexError{T, <:AbstractArray{S,N}}) where {T, S, N}
  ax1 = err.ax1
  print(io, "cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first.")
end

## reinterpret(reshape)

function Base.reinterpret(::typeof(reshape), ::Type{T}, a::ROCArray) where {T}
  N, osize = _base_check_reshape_reinterpret(T, a)
  return _derived_array(T, N, a, osize)
end

# taken from reinterpretarray.jl
# TODO: move these Base definitions out of the ReinterpretArray struct for reuse
function _base_check_reshape_reinterpret(::Type{T}, a::ROCArray{S}) where {T,S}
  isbitstype(T) || throwbits(S, T, T)
  isbitstype(S) || throwbits(S, T, S)
  if sizeof(S) == sizeof(T)
      N = ndims(a)
      osize = size(a)
  elseif sizeof(S) > sizeof(T)
      d, r = divrem(sizeof(S), sizeof(T))
      r == 0 || throwintmult(S, T)
      N = ndims(a) + 1
      osize = (d, size(a)...)
  else
      d, r = divrem(sizeof(T), sizeof(S))
      r == 0 || throwintmult(S, T)
      N = ndims(a) - 1
      N > -1 || throwsize0(S, T, "larger")
      axes(a, 1) == Base.OneTo(sizeof(T) ÷ sizeof(S)) || throwsize1(a, T)
      osize = size(a)[2:end]
  end
  return N, osize
end

@noinline function throwbits(S::Type, T::Type, U::Type)
  throw(ArgumentError("cannot reinterpret `$(S)` as `$(T)`, type `$(U)` is not a bits type"))
end

@noinline function throwintmult(S::Type, T::Type)
  throw(ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got $(sizeof(T))) and `sizeof(eltype(a))` (got $(sizeof(S))) be an integer multiple of the other"))
end

@noinline function throwsize0(S::Type, T::Type, msg)
  throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a $msg size"))
end

@noinline function throwsize1(a::AbstractArray, T::Type)
    throw(ArgumentError("`reinterpret(reshape, $T, a)` where `eltype(a)` is $(eltype(a)) requires that `axes(a, 1)` (got $(axes(a, 1))) be equal to 1:$(sizeof(T) ÷ sizeof(eltype(a))) (from the ratio of element sizes)"))
end


text:
###################################################
## network/iterative/limitedpaths.jl
## Iterative state with fixed subse
###################################################

"""
    `LimitedPaths`, Iterative state optimizing on static trip list (subset of data)
    - each trip has a fixed max number of paths
    - Default number of paths per trip is infinite
"""
mutable struct LimitedPaths <: IterativeState
    # inherited attributes
    data::NetworkData
    timings::NetworkTimings
    trips::Vector{NetworkTrip}
    paths::Vector{Vector{Dict{Edge, Float64}}}  # for each trip, a vector of paths
    pathDiff::Float64

    "max number of paths per trip"
    pathsPerTrip::Int64
    "only trips longer than tripLength get extra paths"
    tripLength::Float64
    "distances of network, used to estimate trip length"
    roadDistances::NetworkTimings
end

"""
    Creates initial fixed number of paths per trip state
    takes in NetworkData object, initial timings (as link times or full timings),
    and max number of paths per trip as integer
"""
function LimitedPaths(data::NetworkData, startSolution::NetworkTimings; pathsPerTrip::Int = typemax(Int), tripLength::Float64 = 0., maxTrip::Int=length(data.trips), yenksp::Bool=false)
    if pathsPerTrip < 1
        error("Must have at least one path per trip")
    end
    # randomly select the trips
    Random.seed!(1991)
    trips = shuffle(data.trips)[1:min(maxTrip,length(data.trips))]
    if yenksp
        paths = yenKSP(startSolution, trips, pathsPerTrip)
    else
        # One path per trip: the initial shortest path
        paths = [Dict{Edge,Float64}[getFullPathEdges(t, startSolution)] for t in trips]
    end
    roadDistances = NetworkTimings(data.network)
    return LimitedPaths(data,startSolution,trips,paths,Inf,pathsPerTrip,tripLength,roadDistances)
end

LimitedPaths(data::NetworkData, initTimes::AbstractArray{Float64, 2}; args...) =
    LimitedPaths(data, NetworkTimings(data.network, initTimes); args...)

heuristicPaths(data::NetworkData; maxTrip::Int=1000) = LimitedPaths(data, uniformTimes(data.network, 15.0), pathsPerTrip=1, maxTrip=maxTrip)

"""
    Update paths of LimitedPaths object given new times
"""
function updateState!(s::LimitedPaths, times::AbstractArray{Float64, 2}; updatePaths::Bool=true)
    # update the timings and compute shortest paths
    s.timings = NetworkTimings(s.data.network, times)

    if updatePaths
        s.pathDiff = 0.
        for (d,t) in enumerate(s.trips)
            # get new shortest path with extra edges if necessary
            sp = getFullPathEdges(t, s.timings)
            s.pathDiff += length(symdiff(collect(keys(sp)), collect(keys(s.paths[d][1]))))/2
            if s.pathsPerTrip == 1 || traveltime(s.roadDistances, t.orig[2], t.dest[1]) < s.tripLength # short trips get one path
                s.paths[d][1] = sp
            else
                index = findfirst(isequal(sp), s.paths[d])
                if index != nothing   # if so, put it in first position
                    s.paths[d][index], s.paths[d][1] = s.paths[d][1], s.paths[d][index] #swap
                elseif length(s.paths[d]) < s.pathsPerTrip # if not, and if enough room to add, add it in first position
                    pushfirst!(s.paths[d], sp)
                else        # replace least useful path
                    worstIndex = argmax([pathEdgesTime(s.timings, collect(keys(p))) for p in s.paths[d]])
                    s.paths[d][worstIndex] = s.paths[d][1]
                    s.paths[d][1] = sp
                end
            end
        end
        s.pathDiff /= length(s.trips)
    end
end


text:
# * Term symbols

struct Term
    L::HalfInteger
    S::HalfInteger
    parity::Parity
end
Term(L::Real, S::Real, parity::Integer) =
    Term(convert(HalfInteger, L), convert(HalfInteger, S), convert(Parity, parity))

function Base.parse(::Type{Term}, s::AbstractString)
    m = match(r"([0-9]+)([A-Z]|\[[0-9/]+\])([oe ]{0,1})", s)
    isnothing(m) && throw(ArgumentError("Invalid term string $s"))
    L = lowercase(m[2])
    L = if L[1] == '['
        L = strip(L, ['[',']'])
        if occursin("/", L)
            Ls = split(L, "/")
            length(Ls) == 2 && length(Ls[1]) > 0 && length(Ls[2]) > 0 ||
                throw(ArgumentError("Invalid term string $(s)"))
            Rational(parse(Int, Ls[1]),parse(Int, Ls[2]))
        else
            parse(Int, L)
        end
    else
        findfirst(L, spectroscopic)[1]-1
    end
    denominator(L) ∈ [1,2] || throw(ArgumentError("L must be integer or half-integer"))
    S = (parse(Int, m[1]) - 1)//2
    Term(L, S, m[3] == "o" ? p"odd" : p"even")
end

macro T_str(s::AbstractString)
    parse(Term, s)
end

Base.zero(::Type{Term}) = T"1S"

multiplicity(t::Term) = convert(Int, 2t.S + 1)
weight(t::Term) = (2t.L + 1) * multiplicity(t)

import Base.==
==(t1::Term, t2::Term) = ((t1.L == t2.L) && (t1.S == t2.S) && (t1.parity == t2.parity))

import Base.<
<(t1::Term, t2::Term) = ((t1.S < t2.S) || (t1.S == t2.S) && (t1.L < t2.L)
                         || (t1.S == t2.S) && (t1.L == t2.L) && (t1.parity < t2.parity))
import Base.isless
isless(t1::Term, t2::Term) = (t1 < t2)

Base.hash(t::Term) = hash((t.L,t.S,t.parity))

include("xu2006.jl")

# This function calculates the term symbol for a given orbital ℓʷ
function xu_terms(ℓ::Int, w::Int, p::Parity)
    ts = map(((w//2 - floor(Int, w//2)):w//2)) do S
        S′ = 2S |> Int
        map(L -> repeat([Term(L,S,p)], Xu.X(w, ℓ, S′, L)), 0:w*ℓ)
    end
    vcat(vcat(ts...)...)
end

function terms(orb::Orbital, occ::Int)
    ℓ = orb.ℓ
    g = degeneracy(orb)
    occ > g && throw(ArgumentError("Invalid occupancy $occ for $orb with degeneracy $g"))
    (occ > g/2 && occ != g) && (occ = g - occ)

    p = parity(orb)^occ
    if occ == 1
        [Term(ℓ,1//2,p)] # Single electron
    elseif ℓ == 0 && occ == 2 || occ == degeneracy(orb) || occ == 0
        [Term(0,0,p"even")] # Filled ℓ shell
    else
        xu_terms(ℓ, occ, p) # All other cases
    end
end

function terms(config::Configuration{O}) where {O<:AbstractOrbital}
    ts = map(config) do (orb,occ,state)
        terms(orb,occ)
    end

    final_terms(ts)
end

"""
    count_terms(orb, occ, term)

Count how many times `term` occurs among the valid terms of
`orb`^`occ`. Example:

```
count_terms(o"1s", 2, T"1S") # 1
```
"""
function count_terms(orb::Orbital, occ::Int, term::Term)
    ℓ = orb.ℓ
    g = degeneracy(orb)
    occ > g && throw(ArgumentError("Invalid occupancy $occ for $orb with degeneracy $g"))
    (occ > g/2 && occ != g) && (occ = g - occ)

    p = parity(orb)^occ
    if occ == 1
        term == Term(ℓ,1//2,p) ? 1 : 0
    elseif ℓ == 0 && occ == 2 || occ == degeneracy(orb) || occ == 0
        term == Term(0,0,p"even") ? 1 : 0
    else
        S′ = convert(Int, 2term.S)
        Xu.X(occ, orb.ℓ, S′, convert(Int, term.L))
    end
end

function write_L(io::IO, term::Term)
    if isinteger(term.L)
        write(io, uppercase(spectroscopic_label(convert(Int, term.L))))
    else
        write(io, "[$(numerator(term.L))/$(denominator(term.L))]")
    end
end

function Base.show(io::IO, term::Term)
    write(io, to_superscript(multiplicity(term)))
    write_L(io, term)
    write(io, to_superscript(term.parity))
end

# * Intermediate terms, seniority

struct IntermediateTerm
    term::Term
    seniority::Int
    function IntermediateTerm(term::Term, seniority::Int)
        iseven(multiplicity(term)) ⊻ iseven(seniority) ||
            throw(ArgumentError("Invalid seniority $(seniority) for term $(term)"))
        new(term, seniority)
    end
end

function Base.show(io::IO, iterm::IntermediateTerm)
    # This is the notation by Giulio Racah, p.377:
    # - Racah, G. (1943). Theory of complex spectra. iii. Physical Review,
    #   63(9-10), 367–382. http://dx.doi.org/10.1103/physrev.63.367
    write(io, to_subscript(iterm.seniority))
    show(io, iterm.term)
end

Base.isless(a::IntermediateTerm, b::IntermediateTerm) =
    a.seniority < b.seniority ||
    a.seniority == b.seniority && a.term < b.term

function intermediate_terms(orb::Orbital, w::Int=one(Int))
    ts = terms(orb, w)
    its = map(unique(ts)) do t
        its = IntermediateTerm[]
        previously_seen = 0
        # The seniority number is defined as the minimum occupancy
        # number ν ∈ n:-2:0 for which the term first appears, e.g. the
        # ²D term first occurs in the d¹ configuration, then twice in
        # the d³ configuration (which will then have the terms ₁²D and
        # ₃²D).
        #
        # We have to loop in reverse, since odd occupation numbers
        # should go from 1 and even from 0.
        for ν ∈ reverse(w:-2:0)
            nn = count_terms(orb, ν, t) - previously_seen
            previously_seen += nn
            append!(its, repeat([IntermediateTerm(t, ν)], nn))
        end
        its
    end
    sort(vcat(its...))
end

function intermediate_terms(config::Configuration)
    map(config) do (orb,occ,state)
        intermediate_terms(orb,occ)
    end
end

export Term, @T_str, multiplicity, weight, terms, count_terms, IntermediateTerm, intermediate_terms


text:
# wrappers for functionality provided by the Metal library toolkit

include("metal/arguments.jl")
include("metal/math.jl")
include("metal/synchronization.jl")


text:
using Test
using GradientBoost.GBDecisionTree
using GradientBoost.LossFunctions

using DecisionTree

instances = [
  1 1;
  1 8;
  1 10
]
labels = [
  1;
  9;
  9
]

# Create DT.
function create_tree(instances, labels)
  model = build_tree(labels, instances)
  model
end

@testset "GB Decision Tree" begin
  @testset "build_base_func works" begin
    gb = GBDT(;
      loss_function=LeastSquares(),
      sampling_rate=0.5,
      learning_rate=0.01,
      num_iterations=100
    )
    prev_func_pred =
      fill(minimizing_scalar(gb.loss_function, labels), size(instances, 1))
    psuedo = negative_gradient(
      gb.loss_function,
      labels,
      gb.learning_rate .* prev_func_pred
    )

    base_func = build_base_func(
      gb,
      instances,
      labels,
      prev_func_pred,
      psuedo
    )

    predictions = base_func(instances)
    expected = [6.2700000000000005; 6.2700000000000005; 6.2700000000000005 ]
    @test predictions ≈ expected
  end

  @testset "instance_to_node indexes" begin
    model = create_tree(instances, labels)
    inst_node_index = GBDecisionTree.InstanceNodeIndex(model, instances)
    num_instances = size(instances, 1)

    @test length(inst_node_index.n2i) == num_instances - 1

    n2i_values = reduce(union, collect(values(inst_node_index.n2i)))
    @test length(n2i_values) == num_instances

    for i = 1:num_instances
      @test in(i, inst_node_index.n2i[inst_node_index.i2n[i]])
    end
  end

  @testset "update_regions! updates terminal regions on tree" begin
    model = create_tree(instances, labels)
    function val_func(leaf::Leaf)
      leaf.majority / 2.0
    end
    n2v = Dict{Leaf, Any}()
    GBDecisionTree.update_regions!(n2v, model, val_func)

    num_instances = size(instances, 1)
    for i = 1:num_instances
      old_pred = apply_tree(model, instances[i,:])
      pred_node = GBDecisionTree.instance_to_node(model, instances[i,:])
      new_pred = n2v[pred_node]

      @test old_pred ./ 2.0 == new_pred
    end
  end

  @testset "LeastAbsoluteDeviation fit_best_constant works" begin
    lf = LeastAbsoluteDeviation()
    dummy_vec = [0.0,0.0,0.0,0.0]
    labels = [0.0,1.0,2.0,3.0]
    prev_func_pred = [3.0,2.0,1.0,0.0]
    expected = 0.0

    actual = GBDecisionTree.fit_best_constant(
      lf, labels, dummy_vec, dummy_vec, prev_func_pred
    )
    @test actual ≈ expected
  end
  @testset "BinomialDeviance fit_best_constant works" begin
    lf = BinomialDeviance()
    dummy_vec = [0.0,0.0,0.0,0.0]

    labels = [0.0,0.0,1.0,1.0]
    psuedo = [0.0,0.5,0.0,0.5]
    expected = -2.0
    actual = GBDecisionTree.fit_best_constant(
      lf, labels, psuedo, dummy_vec, dummy_vec
    )
    @test actual == expected

    labels = [1.0,1.0,1.0,1.0]
    psuedo = [1.0,1.0,1.0,1.0]
    expected = 0.0
    actual = GBDecisionTree.fit_best_constant(
      lf, labels, psuedo, dummy_vec, dummy_vec
    )
    @test actual == expected
  end
end


text:
module VectorizedStatistics

    using LoopVectorization, Static

    _dim(::Type{StaticInt{N}}) where {N} = N::Int

    # Implemented by reduction, recursively
    include("vreducibles.jl")

    # Implemented with @generated functions, single- and multithreaded
    include("vmean.jl")
    include("vsum.jl")
    include("vvar.jl")
    include("vstd.jl")
    include("vcov.jl")

    # Sorting-based statistics
    include("quicksort.jl")
    include("vsort.jl")
    include("vmedian.jl")
    include("vquantile.jl")

end


text:
# # Promises and Futures
#
# [Blocking containers](@ref ex-blockingcontainers) tutorial demonstrated how to
# wait on the arrival (`put!`) of a new element.  However, `put!` is not the
# only interesting event on a concurrent data structure.  By defining `Promise`
# and `Future`, the following code demonstrates a strategy for signaling and
# responding to additional events such as a call to `close`.

using Reagents
using Reagents: Block, CAS, Computed, Map, PostCommit, Read, Return

# ## Promise
#
# Let us implement
# [*Promise*](https://en.wikipedia.org/wiki/Futures_and_promises) with
# the following API:

function test_promise_fetches()
    #=
    We can create a promise, possibly with a specific element type:
    =#
    p = Promise{Int}()

    #=
    Calling `fetch(p::Promise)` will wait for `p` to be fullfiled:
    =#
    task = @task fetch(p)
    yield(task)
    @test !istaskdone(task)  # the task is suspended

    #=
    We can set the value of the promise with the `Ref`-like interface:
    =#
    p[] = 111

    #=
    Once the value is set, all the calls to `fetch` are unblocked:
    =#
    @test fetch(task) == 111  # Note: `task` is calling `fetch(p)`

    #=
    `fetch(p::Promise)` can be called multiple times and it does not block after
    the value is set:
    =#
    @test fetch(p) == 111
end

# We also implement `close(::Promise)`, which unblock `fetch` but with
# exception.

function test_promise_close_before_fetches()
    #=
    Suppose we created a promise and there is a task waiting for it:
    =#
    p = Promise{Int}()
    t = @task fetch(p)
    yield(t)

    #=
    ... but the promise is closed before setting the value
    =#
    close(p)

    #=
    Then, previously blocked `fetch(::Promise)` rasies an exception:
    =#
    err = try
        wait(t)
        nothing
    catch err
        err
    end
    @test err isa TaskFailedException
    @test occursin("promise is closed", sprint(showerror, err))

    #=
    Subsequent call to `fetch(::Promise)` also throws an exception:
    =#
    @test_throws ErrorException("promise is closed") fetch(p)
end

# ### Implementing `Promise`
#
# We store the state of `Promise` in a single `Ref` by using the `Union` type.

struct Closed end

const PromiseRef{T} = Reagents.Ref{Union{
    Nothing,  # indicates the value is not set and the promise is not closed
    Some{T},  # indicates that the value of type `T` is set
    Closed,   # indicates that the promise is closed
}}

# The `Promise` type also contains a channel for sending and receiving signals
# on the state change:

struct Promise{T,Ref<:PromiseRef{T}}
    value::Ref
    send::typeof(Reagents.channel(Nothing)[1])
    receive::typeof(Reagents.channel(Nothing)[2])
end

Promise() = Promise{Any}()
function Promise{T}() where {T}
    send, receive = Reagents.channel(Nothing)
    return Promise(PromiseRef{T}(nothing), send, receive)
end

# Since setting value and closing the channel are similar, we define an internal
# function that tries to set `p.value::Reagents.Ref` if it's not already set and
# then, upon success, notify all the waiters:

tryputting_internal(p::Promise) =
    Computed() do x
        CAS(p.value, nothing, x)
    end ⨟ PostCommit() do _
        while Reagents.trysync!(p.send) !== nothing
        end
    end

# Then, we can define a reagent for setting a value and a reagent for closing
# the promise as simple wrappers:

tryputting(p::Promise{T}) where {T} = Map(Some{T}) ⨟ tryputting_internal(p)
closing(p::Promise) = Return(Closed()) ⨟ tryputting_internal(p)

# The reagent for fetching the promise needs to first listen to the putting and
# closing events (to avoid missing the notification) and *then* check if the
# value is set:

fetching(p::Promise{T}) where {T} =
    (p.receive ⨟ Read(p.value) ⨟ Map(something)) |
    (Read(p.value) ⨟ Map(x -> x === nothing ? Block() : something(x)))

# We check the returned value of `fetching` outside reagent. If it is the
# `Closed` sentinel value, the exception is thrown:

function check_promise_closed(@nospecialize(value))
    if value isa Closed
        error("promise is closed")
    end
    return value
end

# It is now straightforward to define the API mentioned above:

Base.fetch(p::Promise) = check_promise_closed(fetching(p)())

Base.close(p::Promise) = closing(p)()
Base.isopen(p::Promise) = !(p.value[] isa Closed)

function Base.setindex!(p::Promise{T}, x) where {T}
    x = convert(T, x)
    if Reagents.trysync!(tryputting(p), x) === nothing
        check_promise_closed(p.value[])
        error("promise already has a value")
    end
end

# Since we defined underlying synchronization mechanisms as reagents, we can
# compose them. For example, to wait for two promises to be ready, we can use
# the combinator `&`:

function test_promise_fetch_all()
    p1 = Promise{Int}()
    p2 = Promise{Int}()
    t = @task (fetching(p1) & fetching(p2))()
    yield(t)
    p1[] = 222
    @test !istaskdone(t)
    p2[] = 333
    @test fetch(t) == (222, 333)
end

# Or to wait for the first available promise, use `|`

function test_promise_fetch_any()
    p1 = Promise{Int}()
    p2 = Promise{Int}()
    t = @task (fetching(p1) | fetching(p2))()
    yield(t)
    p1[] = 444
    @test fetch(t) == 444
end

# ## Future
#
# Let us define a `Future` as a `Promise` and a thunk that generates the value
# to be stored in the `Promise`. That is to say, we'd like to have the following
# API:

function test_future_fetch_calls_thunk()
    thunk() = 111 + 222
    f = Future{Int}(thunk)
    @test fetch(f) == 333
end

# Importantly, `Future(thunk)` calls `thunk` at most once.

function test_future_thunk_is_called_once()
    #=
    To define this behavior, consider that we have a `thunk` that has a
    side-effect (which is not an intended use-case but useful for describing the
    behavior):
    =#
    ncalled = Ref(0)
    function thunk()
        ncalled[] += 1
        return 111 + 222
    end
    f = Future{Int}(thunk)

    #=
    The first `fetch` will call the `thunk`:
    =#
    ncalled[] = 0
    @test fetch(f) == 333
    @test ncalled[] == 1

    #=
    The Subsequent call to `fetch` does not call the `thunk` and uses the value
    internally stored:
    =#
    ncalled[] = 0
    @test fetch(f) == 333
    @test ncalled[] == 0
end

# Like `Promise`, `fetch`ing `close`d `Future` throws an exception:

function test_future_close()
    ncalled = Ref(0)
    function thunk()
        ncalled[] += 1
        return 111 + 222
    end
    f = Future{Int}(thunk)
    close(f)
    @test_throws ErrorException("promise is closed") fetch(f)
    @test_throws ErrorException("promise is closed") fetch(f)
    #=
    Furthermore, the thunk is not called when the future is `close`d before the
    first `fetch` call:
    =#
    @test ncalled[] == 0
end

# ### Implementing `Future`
#
# A `Future{T}` wraps a `Promise{T}` and a thunk that produces a value of type
# `T`. We also have an auxiliary state `started` tracking if the call to thunk
# is already started or not.

struct Future{T,F,Value<:Promise{T}}
    thunk::F
    value::Value
    started::Threads.Atomic{Bool}
end

Future(f) = Future{Any}(f)
Future{T}(f) where {T} = Future(f, Promise{T}(), Threads.Atomic{Bool}(false))

# We also use this example to demonstrate that not all states have to be
# expressed through Reagents.jl API.  Here, we use a simple
# `Threads.Atomic{Bool}` flag for the `started` state.
#
# The core functionality of `Future` is the ability to run the `thunk` (at most)
# once. Let us define an internal function that assures this invariance. The
# following function `tryrun!(f::Future)` returns `nothing` when the thunk is
# already called. Otherwise, it calls the thunk and then set its value. However,
# if there is a call to `close` before this function returns (i.e., setting the
# value to the promise failed), it also returns `nothing`. If it successfully
# sets the `value` to the promise, it returns `Some(nothing)`.

function tryrun!(f::Future{T}) where {T}
    if Threads.atomic_cas!(f.started, false, true) === false
        y = f.thunk()
        y = convert(T, y)
        # Set the value, if it hasn't been set:
        if Reagents.trysync!(tryputting(f.value), y) === nothing
            return nothing  # already closed
        end
        # Successfully stored the value:
        return Some(y)
    else
        # Lost the race:
        return nothing
    end
end

# We can then wrap this in a reagent. If the call to thunk is successfully, the
# computed value is returned as-is (`Return(something(y))`).  Otherwise,

fetching(f::Future) =
    Computed() do _
        # Optimization: if already closed, not need to call the thunk:
        isopen(f.value) || return fetching(f.value)
        # If still open, try to compute the value:
        y = tryrun!(f)
        if y === nothing
            fetching(f.value)
        else
            Return(something(y))
        end
    end

# The future can be closed by simply closing the underlying promise:

closing(f::Future) = closing(f.value)

# Finally, we can wrap these reagents into the blocking API mentioned above:

Base.fetch(f::Future) = check_promise_closed(fetching(f)())
Base.close(f::Future) = closing(f)()


text:
# test numericRootGenericRandomized
using IncrementalInference
using Test

import IncrementalInference: getSample

##

mutable struct LineResidual <: AbstractRelativeRoots
  m::Float64
  c::Float64
end

getSample(cf::CalcFactor{<:LineResidual}, N::Int=1) = (reshape(rand(Normal(),N),1,N),)

# y = mx + c
# res = y      -      m*x  -  c
#       meas          variables and fixed values
function (cr::CalcFactor{<:LineResidual})(z,
                                          x,
                                          y  )
  #
  return z[1] - (y[1] - (cf.factor.m*x[1] + cf.factor.c))
end

##

@testset "test CommonConvWrapper{T}, solve of residual functions..." begin

##

function assembleConvType(dfgvrs::Vector{<:DFGVariable},
                          functor::T,
                          xDim::Int,
                          zDim::Int,
                          nvars::Int) where {T <: AbstractRelativeRoots}
  #
  # @info "assembleConvType -- development testing function only, not intended for production."
  N = 3

  vars = Array{Array{Float64,2},1}()
  for i in 1:nvars
    push!(vars, zeros(xDim, N))
  end

  fmd = FactorMetadata(dfgvrs, getLabel.(dfgvrs), vars, :null, nothing)

  CommonConvWrapper(functor,vars[1],zDim,vars, fmd, measurement=(zeros(zDim,N),))
end

##

# TODO -- expand testing to include subcomponent tests from _solveCCWNumeric!
lr1 = LineResidual(2.0, 3.0)

fg = initfg()
X0 = addVariable!(fg, :x0, ContinuousScalar)
X1 = addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0,:x1], lr1, graphinit=false)

ccw = assembleConvType([X0;X1], lr1, 1, 1, 2)

res = zeros(1)
ccw(res, zeros(1))

# gwp(x, res)
ccw.cpt[1].particleidx = 1
_solveCCWNumeric!( ccw )

@test abs(ccw.cpt[1].Y[1] + 1.50) < 1e-5

##

end

# println("Test shuffling function")
#
# function testshuffle!(res, x)
#   # println("testshuffle!(x,res) gets x=$(x), and length(res)=$(length(res))")
#   @show res
#   @show x
#   nval = sum(abs.(x-[11.0;12.0;13.0]))
#   # trick the solver so that f/f' = 0/1e10, and therefore < tol
#   # resulting in termination of solver
#   if nval < 1e-8
#     res[1:2] = [0.0;0.0] # fake at root evaluation
#   else
#     res[1:2] = 1e10*randn(2) # fake large gradient
#   end
#   nothing
# end
#
# # test shuffling function
# ALLTESTS = Bool[]
# @show x0 = collect(11:13)+0.0
# for i in 1:10
#   y = numericRootGenericRandomizedFnc(
#           testshuffle!,
#           2, 3, x0,
#           perturb=0.0    )
#   @show y
#   @show y.%10
#   @show yy = y.%10.0
#   push!(ALLTESTS, norm(yy-collect(1:3)) < 0.1)
# end
# @test sum(ALLTESTS) > 9
#
#
# println("Test if shuffling results in correct mapping for solving")
#
#
# function testshuffle2!(res, x)
#   # println("testshuffle2!(x,res) gets x=$(x), and length(res)=$(length(res))")
#   res[1:2] = x - [1.0;2.0] # fake at root evaluation
#   nothing
# end
#
# # test shuffling function
# x0 = collect(11:12)+0.0
# for i in 1:10
#   println("starting with x0=$(x0)")
#   y = numericRootGenericRandomizedFnc(
#           testshuffle2!,
#           2, 2, x0,
#           perturb=0.0,
#           testshuffle=true    )
#   println("result y=$(y)")
#   @test norm(y-collect(1:2)) < 0.1
# end
#
#
#
# # x is dimension 3, z dimension is 2
# function residualrandtest!(res, x)
#   val = norm(x[1:2])
#   res[1] = 10.0 - val
#   nothing
# end
#
# for i in 1:10
#   x0 = ones(2)
#   y = numericRootGenericRandomizedFnc(
#           residualrandtest!,
#           1, 2, x0   )
#   #
#   @test abs(norm(y) - 10.0) < 1e-4
# end


text:
module TestUtils

using Test, Random, DataFrames

@testset "make_unique" begin
    @test DataFrames.make_unique([:x, :x, :x_1, :x2], makeunique=true) == [:x, :x_2, :x_1, :x2]
    @test_throws ArgumentError DataFrames.make_unique([:x, :x, :x_1, :x2], makeunique=false)
    @test DataFrames.make_unique([:x, :x_1, :x2], makeunique=false) == [:x, :x_1, :x2]
end

@testset "repeat count" begin
    df = DataFrame(a=1:2, b=3:4)
    ref = DataFrame(a=repeat(1:2, 2),
                    b=repeat(3:4, 2))
    @test repeat(df, 2) == ref
    @test repeat(view(df, 1:2, :), 2) == ref

    @test size(repeat(df, 0)) == (0, 2)
    @test size(repeat(df, false)) == (0, 2)
    @test_throws ArgumentError repeat(df, -1)
end

@testset "repeat inner_outer" begin
    df = DataFrame(a=1:2, b=3:4)
    ref = DataFrame(a=repeat(1:2, inner=2, outer=3),
                    b=repeat(3:4, inner=2, outer=3))
    @test repeat(df, inner=2, outer=3) == ref
    @test repeat(view(df, 1:2, :), inner=2, outer=3) == ref

    @test size(repeat(df, inner=2, outer=0)) == (0, 2)
    @test size(repeat(df, inner=0, outer=3)) == (0, 2)
    @test size(repeat(df, inner=2, outer=false)) == (0, 2)
    @test size(repeat(df, inner=false, outer=3)) == (0, 2)
    @test_throws ArgumentError repeat(df, inner=2, outer=-1)
    @test_throws ArgumentError repeat(df, inner=-1, outer=3)
end

@testset "repeat! count" begin
    df = DataFrame(a=1:2, b=3:4)
    ref = DataFrame(a=repeat(1:2, 2),
                    b=repeat(3:4, 2))
    a = df.a
    b = df.b
    repeat!(df, 2)
    @test df == ref
    @test a == 1:2
    @test b == 3:4

    for v in (0, false)
        df = DataFrame(a=1:2, b=3:4)
        repeat!(df, v)
        @test size(df) == (0, 2)
    end

    df = DataFrame(a=1:2, b=3:4)
    @test_throws ArgumentError repeat(df, -1)
    @test df == DataFrame(a=1:2, b=3:4)

    @test_throws MethodError repeat!(view(df, 1:2, :), 2)
end

@testset "repeat! inner_outer" begin
    df = DataFrame(a=1:2, b=3:4)
    ref = DataFrame(a=repeat(1:2, inner=2, outer=3),
                    b=repeat(3:4, inner=2, outer=3))
    a = df.a
    b = df.b
    repeat!(df, inner = 2, outer = 3)
    @test df == ref
    @test a == 1:2
    @test b == 3:4

    for v in (0, false)
        df = DataFrame(a=1:2, b=3:4)
        repeat!(df, inner=2, outer=v)
        @test size(df) == (0, 2)

        df = DataFrame(a=1:2, b=3:4)
        repeat!(df, inner=v, outer=3)
        @test size(df) == (0, 2)
    end

    df = DataFrame(a=1:2, b=3:4)
    @test_throws ArgumentError repeat(df, inner = 2, outer = -1)
    @test_throws ArgumentError repeat(df, inner = -1, outer = 3)
    @test df == DataFrame(a=1:2, b=3:4)

    @test_throws MethodError repeat!(view(df, 1:2, :), inner = 2, outer = 3)
end

@testset "funname" begin
    @test DataFrames.funname(sum ∘ skipmissing ∘ Base.div12) ==
          :sum_skipmissing_div12
end

@testset "pre-Julia 1.3 @spawn replacement" begin
    t = @sync DataFrames.@spawn begin
        sleep(1)
        true
    end
    @test fetch(t) === true
end

@testset "split_indices" begin
    for len in 1:100, basesize in 1:10
        x = DataFrames.split_indices(len, basesize)

        @test length(x) == max(1, div(len, basesize))
        @test reduce(vcat, x) == 1:len
        vmin, vmax = extrema(length(v) for v in x)
        @test vmin + 1 == vmax || vmin == vmax
        @test len < basesize || vmin >= basesize
    end

    @test_throws AssertionError DataFrames.split_indices(0, 10)
    @test_throws AssertionError DataFrames.split_indices(10, 0)

    # Check overflow on 32-bit
    len = typemax(Int32)
    basesize = 100_000_000
    x = collect(DataFrames.split_indices(len, basesize))
    @test length(x) == div(len, basesize)
    @test x[1][1] === 1
    @test x[end][end] === Int(len)
    vmin, vmax = extrema(length(v) for v in x)
    @test vmin + 1 == vmax || vmin == vmax
    @test len < basesize || vmin >= basesize
end

@testset "split_to_chunks" begin
    for lg in 1:100, nt in 1:11
        if lg < nt
            @test_throws AssertionError DataFrames.split_to_chunks(lg, nt)
            continue
        end
        x = collect(DataFrames.split_to_chunks(lg, nt))
        @test reduce(vcat, x) == 1:lg
        @test sum(length, x) == lg
        @test first(x[1]) == 1
        @test last(x[end]) == lg
        @test length(x) == nt
        for i in 1:nt-1
            @test first(x[i+1])-last(x[i]) == 1
        end
    end

    @test_throws AssertionError DataFrames.split_to_chunks(0, 10)
    @test_throws AssertionError DataFrames.split_to_chunks(10, 0)
    @test_throws AssertionError DataFrames.split_to_chunks(10, 11)
end

@testset "_findall(B::BitVector)" begin
    Random.seed!(1234)
    BD = Dict(
        "Empty" => (BitVector([]), UnitRange{Int}),
        "T Big" => (trues(100000), UnitRange{Int}),
        "F Big" => (falses(100000), UnitRange{Int}),
        "T64 F64" => ([trues(64); falses(64)], UnitRange{Int}),
        "F64 T64" => ([falses(64); trues(64)], UnitRange{Int}),
        "F80 T100" => ([falses(85); trues(100)], UnitRange{Int}),
        "F256 T32" => ([falses(256); trues(32)], UnitRange{Int}),
        "F260 T32" => ([falses(260); trues(32)], UnitRange{Int}),
        "TF Big" => ([trues(100000); falses(100000)], UnitRange{Int}),
        "FT Big" => ([falses(100000); trues(100000)], UnitRange{Int}),

        # some edge cases
        "TFT small" => ([trues(85); falses(100); trues(85)], Vector{Int}),
        "FTFFT small" => ([falses(64 + 32); trues(32); falses(128); trues(32)], Vector{Int}),
        "TFTF small" => ([falses(64); trues(64); falses(64); trues(64)], Vector{Int}),
        "TFT small2" => ([trues(64); falses(10); trues(100)], Vector{Int}),

        "FTF Big" => ([falses(8500); trues(100000); falses(65000)], UnitRange{Int}),
        "TFT Big" => ([trues(8500); falses(100000); trues(65000)], Vector{Int}),
        "FTFTFTF Big" => ([falses(65000); trues(65000); falses(65000); trues(65000); falses(65000); trues(65000); falses(65000)], Vector{Int}),

        "FTFR small" => ([falses(85); trues(100); falses(65); rand([true, false], 2000)], Vector{Int}),
        "R Big" => (BitVector(rand([true, false], 2000000)), Vector{Int}),
        "RF Big" => ([BitVector(rand([true, false], 1000000)); falses(1000000)], Vector{Int}),
        "RT Big" => ([BitVector(rand([true, false], 1000000));  trues(1000000)], Vector{Int}),
        "FR Big" => ([falses(1000000); BitVector(rand([true, false], 1000000))], Vector{Int}),
        "TR Big" => ([trues(1000000);  BitVector(rand([true, false], 1000000))], Vector{Int}),
        "FRT Big" => ([falses(1000000); BitVector(rand([true, false], 1000000)); trues(1000000)], Vector{Int}),
        "TRF Big" => ([trues(1000000); BitVector(rand([true, false], 1000000)); falses(1000000)], Vector{Int}),
        "FRF Big" => ([falses(1000000); BitVector(rand([true, false], 1000000)); falses(1000000)], Vector{Int}),
        "TRT Big" => ([trues(1000000); BitVector(rand([true, false], 1000000)); trues(1000000)], Vector{Int}),
        "RFR Big" => ([BitVector(rand([true, false], 1000000)); falses(1000000); BitVector(rand([true, false], 1000000))], Vector{Int}),
        "FTFR Big" => ([falses(65000);  trues(65000);  falses(65000); rand([true, false], 20000)], Vector{Int}),
        "T256 R100" => ([trues(256);  rand([true, false], 2000)], Vector{Int}),
        "F256 R100" => ([falses(256); rand([true, false], 2000)], Vector{Int}),
    )
    for (_, (B, T)) in BD
        res = DataFrames._findall(B)
        @test Base.findall(B) == res
        @test res isa T
        res = DataFrames._findall(Vector{Bool}(B))
        @test Base.findall(B) == res
        @test res isa T
    end

    # 1:200 is to test all small cases
    # 1000 is to test skipping multiple 64-bit blocks of 0 and 1
    for n in [1:200; 1000], i in 1:n, j in i:n
        x = falses(n)
        x[i:j] .= true
        res = DataFrames._findall(x)
        @test res == i:j
        @test res isa UnitRange{Int}
        res = DataFrames._findall(Vector{Bool}(x))
        @test res == i:j
        @test res isa UnitRange{Int}
        if j + 1 < n
            x[j + 2] = true
            # add one false and then true
            @test  DataFrames._findall(x) == [i:j; j+2]
            @test  DataFrames._findall(Vector{Bool}(x)) == [i:j; j+2]
            # sprinkle trues randomly after one false
            rand!(view(x, j + 2:n), Bool)
            @test  DataFrames._findall(x) == findall(x)
            @test  DataFrames._findall(Vector{Bool}(x)) == findall(x)
        end
    end
end

end # module


text:
"""
    DriverModel{DriveAction}
A DriverModel represents a specific driving behavior.
It specifies the action taken by the agent at a given scene.
The ation will be of type `DriveAction`.
It can be interpreted as a distribution, the likelihood of taking a certain action
in a given scene.
The DriverModel type is an abstract type! Custom driver models should inherit from it.
"""
abstract type DriverModel{DriveAction} end

"""
    action_type(::DriverModel{A}) where {A}
returns the type of the actions that are sampled from the model
"""
action_type(::DriverModel{A}) where {A} = A

"""
    set_desired_speed!(model::DriverModel, v_des::Float64)
Sets a desired speed.
This method is relevant for models like IDM where the vehicle tracks a nominal speed.
"""
function set_desired_speed! end

"""
    reset_hidden_state!(model::DriverModel)
Resets the hidden states of the model.
"""
function reset_hidden_state! end

"""
    reset_hidden_states!(models::Dict{I,M}) where {M<:DriverModel}
reset hidden states of all driver models in `models`
"""
function reset_hidden_states!(models::Dict{I,M}) where {I, M<:DriverModel}
    for model in values(models)
        reset_hidden_state!(model)
    end
    return models
end

"""
    observe!(model::DriverModel, scene, roadway, egoid)
Observes the scene and updates the model states accordingly.
"""
function observe! end

"""
    rand(model::DriverModel)
    rand(rng::AbstractRNG, model::DriverModel)
Samples an action from the model.
"""
Base.rand(model::DriverModel) = rand(Random.GLOBAL_RNG, model)
Base.rand(rng::AbstractRNG, model::DriverModel) = error("AutomotiveDrivingModelsError: Base.rand(::AbstractRNG, ::$(typeof(model))) not implemented")

####

"""
    StaticDriver{A,P<:ContinuousMultivariateDistribution} <: DriverModel{A}

A driver model where actions are always sampled by the same distribution specified
by the field `distribution`.

# Fields
- `distribution::P`
"""
struct StaticDriver{A,P<:ContinuousMultivariateDistribution} <: DriverModel{A}
    distribution::P
end

function Base.rand(rng::AbstractRNG, model::StaticDriver{A,P}) where {A,P}
    a = rand(rng, model.distribution)
    return convert(A, a)
end

Distributions.pdf(model::StaticDriver{A}, a::A) where {A} = pdf(model.distribution, convert(Vector{Float64}, a))
Distributions.logpdf(model::StaticDriver{A}, a::A) where {A} = logpdf(model.distribution, convert(Vector{Float64}, a))
set_desired_speed!(model::StaticDriver, v::Float64) = model
observe!(model::StaticDriver, scene::S, roadway::Roadway, egoid::I) where {S, I} = model


text:
"""
constants : small set of variables related to translation
"""

aa_table = Dict{String, Array{String, 1}}("A" => ["GCA", "GCC", "GCG", "GCT"], "C" => ["TGT", "TGC"],
            "E" => ["GAG", "GAA"], "D" => ["GAT", "GAC"], "G" => ["GGT", "GGG", "GGA", "GGC"],
            "F" => ["TTT", "TTC"], "I" => ["ATC", "ATA", "ATT"], "H" => ["CAT", "CAC"],
            "K" => ["AAG", "AAA"], "*" => ["TAG", "TGA", "TAA"], "M" => ["ATG"],
            "L" => ["CTT", "CTG", "CTA", "CTC", "TTA", "TTG"], "N" => ["AAC", "AAT"],
            "Q" => ["CAA", "CAG"], "P" => ["CCT", "CCG", "CCA", "CCC"],
            "S" => ["AGC", "AGT", "TCT", "TCG", "TCC", "TCA"],
            "R" => ["AGG", "AGA", "CGA", "CGG", "CGT", "CGC"],
            "T" => ["ACA", "ACG", "ACT", "ACC"], "W" => ["TGG"],
            "V" => ["GTA", "GTC", "GTG", "GTT"], "Y" => ["TAT", "TAC"])

first_letter_codons = Dict{String, Array{String, 1}}("T" => ["TTT", "TTC", "TTA", "TTG", "TCT", "TCC", "TCA", "TCG", "TAT", "TAC", "TGT", "TGC", "TGG"],
                       "G" => ["GTT", "GTC", "GTA", "GTG", "GCT", "GCC", "GCA", "GCG", "GAT", "GAC", "GAA", "GAG", "GGT", "GGC", "GGA", "GGG"],
                       "A" => ["ATT", "ATC", "ATA", "ATG", "ACT", "ACC", "ACA", "ACG", "AAT", "AAC", "AAA", "AAG", "AGT", "AGC", "AGA", "AGG"],
                       "C" => ["CTT", "CTC", "CTA", "CTG", "CCT", "CCC", "CCA", "CCG", "CAT", "CAC", "CAA", "CAG", "CGT", "CGC", "CGA", "CGG"])

last_letter_codons = Dict{String, Array{String, 1}}("C" => ["TTC", "TCC", "TAC", "TGC", "CTC", "CCC", "CAC", "CGC", "ATC", "ACC", "AAC", "AGC", "GTC", "GCC", "GAC", "GGC"],
                      "T" => ["TTT", "TCT", "TAT", "TGT", "CTT", "CCT", "CAT", "CGT", "ATT", "ACT", "AAT", "AGT", "GTT", "GCT", "GAT", "GGT"],
                      "G" => ["TTG", "TCG", "TGG", "CTG", "CCG", "CAG", "CGG", "ATG", "ACG", "AAG", "AGG", "GTG", "GCG", "GAG", "GGG"],
                      "A" => ["TTA", "TCA", "CTA", "CCA", "CAA", "CGA", "ATA", "ACA", "AAA", "AGA", "GTA", "GCA", "GAA", "GGA"])

first_letter_possibilities = Dict{String, Array{String, 1}}("A" => ["I", "M", "T", "N", "K", "S", "R"],
                              "C" => ["L", "P", "H", "Q", "R"],
                              "T" => ["F", "L", "S", "Y", "C", "W"],
                              "G" => ["V", "A", "D", "E", "G"])

top_codon = Dict{Char, AbstractString}('A'=>"GCG", 'C'=>"TGC", 'D'=>"GAC", 'E'=>"GAA", 'F'=>"TTT",
											'G'=>"GGT", 'H'=>"CAT", 'I'=>"ATT", 'K'=>"AAA", 'L'=>"CTG",
											'M'=>"ATG", 'N'=>"AAC", 'P'=>"CCG", 'Q'=>"CAG", 'R'=>"CGT",
											'S'=>"AGC", 'T'=>"ACC", 'V'=>"GTG", 'W'=>"TGG", 'Y'=>"TAT")


text:
export StandardMass

mutable struct StandardMass{K} <: MassComponent{Vector{Float},K}
    "Smoothing kernel."
    W::K
    "Mass of the fluid particles."
    m::Float
    mass_densities::Vector{Float}

    """
        StandardMass(
            W::Kernel{N},
            mass_distribution::InitialMassDistribution{N};
            η::Number = 1.2,
        ) -> (StandardMass, Vector{SVectorF{N}})
    """
    function StandardMass(
        W::K,
        mass_distribution::InitialMassDistribution{N};
        η::Number = 1.2,
    ) where {N,K<:Kernel{N}}
        @check !isa(W, ConstantWidthKernel)
        @check η > 0
        data = mass_distribution()
        m = data.m
        mass_densities = data.mass_densities
        positions = data.positions
        @check data.m > 0
        @checkeq(length(mass_densities), length(positions))
        h = @show estimate_initial_kernel_width(
            mean(mass_densities),
            m,
            η,
            Dim{N}(),
        )
        @check(h > 0)
        constant_width_W = ConstantWidthKernel(W, h)
        new{ConstantWidthKernel{N,K}}(constant_width_W, m, mass_densities),
        positions
    end
end

get_kernel(standard_mass::StandardMass) = standard_mass.W
get_particle_mass(standard_mass::StandardMass) = standard_mass.m
get_mass_densities(standard_mass::StandardMass) = standard_mass.mass_densities
get_kernel_widths(standard_mass::StandardMass) =
    fill(get_kernel_width(standard_mass), length(standard_mass.mass_densities))

"""
    get_kernel_width(standard_mass::StandardMass) -> Float

Return the constant kernel width of the standard mass component.
"""
get_kernel_width(standard_mass::StandardMass) =
    get_kernel_width(standard_mass.W)

"""
    get_kernel_widths(standard_mass::StandardMass, ::Tuple{Unsigned,Unsigned})
        -> Float

Return the constant kernel width of the standard mass component.
"""
get_kernel_widths(standard_mass::StandardMass, ::Tuple{Unsigned,Unsigned}) =
    get_kernel_width(standard_mass)

"""
    get_kernel_width(standard_mass::StandardMass, ::Unsigned) -> Float

Return the constant kernel width of the standard mass component.
"""
get_kernel_width(standard_mass::StandardMass, ::Unsigned) =
    get_kernel_width(standard_mass)

"""
    normalize_distance(
        standard_mass::StandardMass,
        r::Number,
        ::Number,
    ) -> Float

Normalize the distance by the constant kernel width of the standard mass component.
"""
normalize_distance(standard_mass::StandardMass, r::Number, ::Number) =
    r * get_inverse_kernel_width(standard_mass.W)

initderivatives(standard_mass::StandardMass) =
    zeros(length(standard_mass.mass_densities))

initderivatives!(dρdt::AbstractVector{Float}, ::StandardMass) = fill!(dρdt, 0.0)

function updatederivatives!(
    ::AbstractVector{Float},
    ::Tuple{Unsigned,Unsigned},
    ::Number,
    ::Tuple{Number,Number},
    ::StandardMass,
)
    throw(InvalidStateException)
end

function updatederivatives!(
    dρdt::AbstractVector{Float},
    (i, j)::Tuple{Unsigned,Unsigned},
    vᵣ::Number,
    m_ddr_W_r_h::Number,
    ::StandardMass,
)
    Δdρdt = compute_Δdρdt(m_ddr_W_r_h, vᵣ)
    dρdt[i] += Δdρdt
    dρdt[j] += Δdρdt
end

function updatederivative!(
    dρdt::AbstractVector{Float},
    i::Unsigned,
    vᵣ::Number,
    m_ddr_W_r_h::Number,
    ::StandardMass,
)
    dρdt[i] += compute_Δdρdt(m_ddr_W_r_h, vᵣ)
end

compute_Δdρdt(m_ddr_W_r_h::Number, vᵣ::Number) = m_ddr_W_r_h * vᵣ

evolvevariables!(
    standard_mass::StandardMass,
    dρdt::AbstractVector{Float},
    Δt::Number,
) = standard_mass.mass_densities .+= Δt * dρdt

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl