text:
import StateSpaceReconstruction: embed

"""
    te(driver, response;
        estimator = :tetogrid,
        E = nothing, v = nothing,
        ϵ = nothing, n_ϵ = 5, ν = 1, τ = 1, dim = 3,
        k1 = 2, k2 = 3, distance_metric = Chebyshev(),
        which_is_surr = :none,
        surr_func = aaft,
        min_numbins = nothing,
        max_numbins = nothing)

Compute transfer entropy from a `driver` time series to a `response` time
series.

## Method
`method` sets the transfer entropy estimator to use. There are two types of
estimators: grid-based approaches, and nearest-neighbor based approaches.

## Embedding instructions
Specify embedding and instructions for how to compute marginal
entropies are with `E::AbstractEmbedding` and `v::TEVars`. If either
is missing, the data is embedded using the provided forward prediction lag `ν`,
with embedding dimension `dim` and embedding lag `η`.

## Grid-based estimators

- `:transferentropy_transferoperator_grid`, or `:tetogrid`. Grid-based transfer entropy
    estimator based on the transfer operator.
- `::transferentropy_visitfreq`, or `:tefreq`. Simple visitation frequency based estimator.

For the grid based TE estimators, `ϵ` sets the bin size. If `ϵ == nothing`, then
the algorithm uses a bin size corresponding to a number of subdivisions `N`
along each axis so that `N ≦ npoints^(1/(dim+1))`. Transfer entropy is then
computed as an average over `n_ϵ` different bin sizes corresponding to
`N-1` to `N` subdivisions along each axis. If `ϵ` is specified and consists
of multiple bin sizes, an average
TE estimate over the bin sizes is returened. If `ϵ` is a single value,
then TE is estimated for that bin size only.

## Nearest neighbor based estimator.
- `:transferentropy_kraskov`, `:tekraskov`, or `:tekNN`. A nearest-neighbor
    based estimator that computes transfer entropy as the sum of two mutual
    information terms. `k1` and `k2` sets the
    number of nearest neighbors used to estimate the mutual information terms.
    `distance_metric` sets the distance metric (must be an instance
    of a valid metric from `Distances.jl`).

## Surrogate analysis
A surrogate analysis will be run if `which_is_surr` is set to
    `:both`, `:driver`, `:response`. Default is `:none`.
- `which_is_surr = :both` will replace both time series with a surrogate.
- `which_is_surr = :driver` will replace the driver time series with a surrogate.
- `which_is_surr = :response` will replace the response time series with a surrogate.

The type of surrogate must be either `randomshuffle`, `randomphases`,
`randomamplitudes`, `aaft` or `iaaft`.

"""
function transferentropy(driver, response;
        dim = 3,
        τ = 1,
        ν = 1,
        estimator = :tetogrid,
        E = nothing, v = nothing,
        ϵ = nothing, n_ϵ = 5,
        k1 = 2, k2 = 3, distance_metric = Chebyshev(),
        which_is_surr = :none,
        surr_func = aaft,
        min_numbins = nothing,
        max_numbins = nothing)

    estims_te_grid = [:transferentropy_transferoperator_grid, :tetogrid]
    estims_te_kNN = [:transferentropy_kraskov, :tekraskov, :tekNN]
    estims_te_freq = [:transferentropy_visitfreq, :tefreq]
    valid_estimators = estims_te_grid ∪ estims_te_kNN ∪ estims_te_freq

    if dim < 3
        error("Dimension $dim is too low. Must be at least 3.")
    end

    if !(estimator ∈ valid_estimators)
        error("Transfer entropy method $estimator not valid.")
    end

    valid_surr_funcs = [TimeseriesSurrogates.randomshuffle,
                        TimeseriesSurrogates.randomphases,
                        TimeseriesSurrogates.randomamplitudes,
                        TimeseriesSurrogates.aaft,
                        TimeseriesSurrogates.iaaft]

    if !(surr_func ∈ valid_surr_funcs)
        error("Surrogate type $surr_func is not valid.")
    end
    if which_is_surr == :driver
        driver[:] = surr_func(driver)
    elseif which_is_surr == :response
        response[:] = surr_func(response)
    elseif which_is_surr == :both
        driver[:] = surr_func(driver)
        response[:] = surr_func(response)
    elseif !(which_is_surr == :none)
        error("$which_is_surr is not a valid surrogate choice." )
    end

    # Compute transfer entropy analogously to how Krakovskà et al. (2018)
    # computes conditional mutual information.
    if !(typeof(E) == Embeddings.AbstractEmbedding) || !(typeof(v) == TEVars)
        if dim == 3
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 1], [ν, 0, 0])
            v = TEVars([1], [2], [3])
        elseif dim == 4
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 2, 1], [ν, 0, -τ, 0])
            v = TEVars([1], [2, 3], [4])
        elseif dim == 5
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 2, 2, 1], [ν, 0, -τ, -2*τ, 0])
            v = TEVars([1], [2, 3, 4], [5])
        elseif dim == 6
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 2, 2, 2, 1], [ν, 0, -τ, -2*τ, -3*τ, 0])
            v = TEVars([1], [2, 3, 4, 5], [6])
        elseif dim == 7
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 2, 2, 2, 2, 1], [ν, 0, -τ, -2*τ, -3*τ, -4*τ, 0])
            v = TEVars([1], [2, 3, 4, 5, 6], [7])
        elseif dim == 8
            E = StateSpaceReconstruction.embed([driver, response], [2, 2, 2, 2, 2, 2, 2, 1], [ν, 0, -τ, -2*τ, -3*τ, -4*τ, -5*τ, 0])
            v = TEVars([1], [2, 3, 4, 5, 6, 7], [8])
        end
    end

    if ϵ == nothing
        ϵ = zeros(Float64, n_ϵ, dim)
        # We a range of bin sizes along each axis.
        for d in 1:dim
            dmin = Base.minimum(E.points[d, :])
            dmax = Base.maximum(E.points[d, :])
            maxrange = dmax - dmin

            # If minimum and maximum number of bins are not provided,
            # use the Palus recommendation of N^(1/dim) as a max.
            if !(typeof(min_numbins) == Int && typeof(max_numbins) == Int)
                max_numbins = ceil(Int, npoints(E)^(1/(dim+1)))
                min_numbins = max_numbins - 1
            end

            max_bin_size = maxrange/max_numbins
            min_bin_size = maxrange/min_numbins

            ϵ[:, d] = LinRange(log(10, min_bin_size), log(10, max_bin_size), n_ϵ)
        end
        binsizes = [ϵ[i, :] for i = 1:n_ϵ]
        ϵs = binsizes
    else
        if length(ϵ) == 1
            ϵs = [ϵ]
        else
            ϵs = ϵ
        end
    end

    tes = zeros(Float64, length(ϵs))

    # Loop over the bin sizes
    for i = 1:length(ϵs)
        if estimator ∈ (:transferentropy_transferoperator_grid,
                        :tetogrid)
            tes[i] = transferentropy_transferoperator_grid(E, ϵs[i], v)

        elseif estimator ∈ (:transferentropy_kraskov,
                            :tekraskov,
                            :tekNN)
            tes[i] = transferentropy_kraskov(E, k1, k2, v, metric = distance_metric)

        elseif estimator ∈ (:transferentropy_visitfreq,
                            :tefreq)
            tes[i] = transferentropy_visitfreq(E, ϵs[i], v)

        elseif estimator == :transferoperator_triang
                warn("$estimator is not implemented yet.")
                #warn("Invariantizing embedding before calculating transfer matrix...")
                #t = triangulate(invariantize(E))
                #TO = transferoperator_approx(t, n_pts = n_randpts, sample_randomly = sample_randomly)
                #invdist = left_eigenvector(TO)
                #tes[i] = transferentropy_transferoperator_triang(t, invdist, n_ϵ, n_reps)
        end

    end

    # If we only use one binsize, we can't take the integral. This is always
    # the case for the kNN based estimators, because they don't use a bin size.
    if estimator ∈ estims_te_kNN
        return tes[1]
    else
        return mean(tes)
    end
end
te = transferentropy
export transferentropy, te


text:
"""
    bessel_moment_G(λ, ν, a)

Computes the standard Bessel moment ``G_{\\lambda, \\nu}(a)`` defined as
```math
G_{\\lambda, \\nu}(a) = \\int_0^a t^{\\lambda} J_{\\nu}(t)\\,dt.
```
"""
function bessel_moment_G(λ, ν, a::T) where T
    return HypergeometricFunctions.pFq(
            SA[(1 + λ + ν)/2],SA[(3+λ+ν)/2, 1+ν], (-a^2/4)) * exp(
        log(a) * (1 + λ + ν) - log(2) * (1 + ν) +
        lgamma((1+λ+ν)/2) - lgamma((3+λ+ν)/2) - lgamma((1+ν)))
end

# modified moment boundary conditions (Piessens & Branders 1983 eq. 15-18)
M₀(a, ν) = bessel_moment_G(0, ν, a) / a
M₁(a, ν) = ((2/a) * bessel_moment_G(1, ν, a) - bessel_moment_G(0, ν, a)) / a
M₂(a, ν) = ((8/a^2)*bessel_moment_G(2, ν, a) -
    (8/a)*bessel_moment_G(1, ν, a) + bessel_moment_G(0, ν, a)) / a
M₃(a, ν) = ((32/a^3) * bessel_moment_G(3, ν, a) -
    (48/a^2) * bessel_moment_G(2, ν, a) +
    (18/a) * bessel_moment_G(1, ν, a) - bessel_moment_G(0, ν, a)) / a

@doc raw"""
    Mk_asymptotic(a::T, ν, k) where T

Computes a modified moment using asymptotics as ``k \rightarrow \infty``.
"""
function Mk_asymptotic(a::T, ν, k) where T

    ϕ0 = SA[
        J(ν, a),
        zero(T),
        ((3ν-2) * J(ν, a) - 3a * J(ν-1, a)) / 2,
        zero(T),
        ( 15a * J(ν-1, a) + (4 - 15a^2 + 15ν * (ν-1)) * J(ν, a) ) / 4,
        0,
        (21a * (-2 + 5a^2 - 5ν^2) * J(ν-1, a) +
            (-8 - 105a^2 * (ν-3) + 21ν * (2 + 5 * (ν-1) * ν)) * J(ν, a)) / 8
    ]

    ϕπ = SA[
        -J(ν, zero(T)),
        zero(T),
        J(ν, zero(T)) - (3a/4) * J(ν-1, zero(T)),
        zero(T),
        (-15a^2 * J(ν-2, zero(T)) + 60a * J(ν-1, zero(T)) +
            2 * (15a^2 - 8) * J(ν, zero(T))) / 16,
        zero(T),
        (-105a^3 * J(ν-3, zero(T)) +
            840a^2 * J(ν-2, zero(T)) -
            (315a^3 - 1008a) * J(ν-1, zero(T)) +
            (64 - 1680a^2) * J(ν, zero(T))
        ) / 64
    ]

    Σ = zero(T)
    for j in 0:3
        Σ = Σ + (-1)^j * ((-1)^k * ϕπ[2j+1] - ϕ0[2j+1]) * T(k)^(-2j-2)
    end

    return Σ / 2
end

"""
    BrandersPiessensProblem(a, ν)

Construct a boundary value problem from the modified moment recurrence relation.

`N` and `M` correspond to the total number of specified initial boundary conditions,
and the number specified on the end of the recurrence, respectively. T is the type
we're computing on, `a` is the parameter in the Bessel function Jᵥ(ax), and `ν` is
the order.
"""
struct BrandersPiessensProblem{T,N,M} <: OliverProblem{T,N,M}
    a::T
    ν::T
end

@doc raw"""
    ClenshawCurtisBessel.OliverP(OP::BrandersPiessensProblem{T}, s, k) where T

Computes the modifed moments from Branders & Piessens 1987.

```math
\begin{aligned}
\frac{a^2}{16} &M_{k}(a, \nu) + \left[ (k+1)^2 - \nu^2 - \frac{a^2}{4} \right]
M_{k + 2}(a, \nu) + \left[ 4 \nu^2 - 2(k+4) + 4 \right] M_{k+3}(a,\nu) \\
&- \left[ 2(k+4)^2 - 6 + 6\nu^2 - \frac{3a^2}{8} \right] M_{k+4}(a,\nu) +
(4\nu^2 + 2(k+4) + 4)M_{k+5}(a,\nu) \\
&+ \left[ (k+7)^2 - \nu^2 - \frac{a^2}{4} \right]M_{k+6} + \frac{a^2}{16} M_{k+8}(a,\nu) \quad = 0
\end{aligned}
```

Note that the coefficient indices decrease ``P_8, P_7, \cdots, P_0`` as the ``M_k``
index increases.
"""
function OliverP(
        OP::BrandersPiessensProblem{T}, s, k) where T
    a, ν = OP.a, OP.ν
    if s == 8
        return a^2/16
    elseif s == 7
        return zero(T)
    elseif s == 6
        return (k+1)^2 - ν^2 - a^2/4
    elseif s == 5
        return 4ν^2 - 2(k+4) + 4
    elseif s == 4
        return -(2 * (k+4)^2 - 6 + 6ν^2 - 3 * a^2 / 8)
    elseif s == 3
        return 4ν^2 + 2(k+4) + 4
    elseif s == 2
        return (k+7)^2 - ν^2 - a^2/4
    elseif s == 1
        return zero(T)
    elseif s == 0
        return a^2/16
    end
end

"""
    OliverR(OP::BrandersPiessensProblem{T}, i)

Computes ``R(i)``, the right side of the recurrence relation.
"""
OliverR(OP::BrandersPiessensProblem{T}, i) where T = zero(T)

"""
    generate_BC(OP::BrandersPiessensProblem{T,8,2}, kmax::Int) where T

Generates the boundary conditions required for Oliver's method to solve
the boundary value problem posed Branders & Piessens.

# Arguments:
- `OP::BrandersPiessensProblem{T,8,2}`: problem parameters
- `kmax`: maximum order ``k`` to compute.

# Returns:
- `Tuple{Int,OffsetVector}`: the maximum index, and a sparse offset array
     with the boundaries.
"""
function generate_BC(OP::BrandersPiessensProblem{T,8,2}, kmax::Int) where T
    a, ν = OP.a, OP.ν
    ℓ = max(kmax+1, Int(ceil(a + 10)))  # added +1, never use asymptotic for kmax
    Mbc = OffsetVector(spzeros(T,ℓ+5), -3:(ℓ+1))

    Mbc[-3] = M₃(a, ν)
    Mbc[-2] = M₂(a, ν)
    Mbc[-1] = M₁(a, ν)
    Mbc[ 0] = M₀(a, ν)
    Mbc[ 1] = M₁(a, ν)
    Mbc[ 2] = M₂(a, ν)
    Mbc[ 3] = M₃(a, ν)

    Mbc[ℓ] = Mk_asymptotic(a, ν, ℓ)
    Mbc[ℓ+1] = Mk_asymptotic(a, ν, ℓ+1)

    return ℓ+1, Mbc
end

function momentM(Tout::Type{<:Number}, a::T, nu::T, kmax::Int, kbuffer=100) where T
    bpp = BrandersPiessensProblem{T,8,2}(a, nu)
    upperindex, Mbc = generate_BC(bpp, kmax+kbuffer)
    ρ = assembleρ(bpp, -3, upperindex, Mbc)
    𝐏 = assembleP(bpp, -3, upperindex)

    sol = OffsetVector(Array{Tout}(undef,kmax+1), 0:kmax)
    for i in 0:3
        sol[i] = Mbc[i]
    end
    res = qr(𝐏) \ ρ
    for i in 4:kmax
        sol[i] = Tout(res[i-3])
    end
    return sol
end

ext: jl


text:
"""
    nonlinear objective
"""
struct NonlinearObjective <: Objective
    f # objectivce
    g # objective gradient
    gradient_type
end

function nonlinear_objective(f)
    NonlinearObjective(f, x -> zero(x), :auto)
end

function nonlinear_objective(f, g)
    NonlinearObjective(f, g, :gradient)
end

function objective(Z, obj::NonlinearObjective, idx, T)
    return obj.f(Z)
end

function objective_gradient!(∇J, Z, obj::NonlinearObjective, idx, T)
    if obj.gradient_type == :auto
        ∇J .+= ForwardDiff.gradient(obj.f, Z)
    else
        ∇J .+= obj.g(Z)
    end
    return nothing
end


text:
#list of added packages

#routines for adding the packages
add Atom Juno
add DataFrames Revise Distributions StatsBase GLM Dates NLopt ForwardDiff Formatting
add DataStructures LinearAlgebra Gadfly Measures Formatting CodecZlib
add CSV GZip NLopt CategoricalArrays Serialization Compose Distributed Random Cairo
add Fontconfig Compose StaticArrays CodecLz4
add https://github.com/clintonTE/Finometrics

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Tables

const lookup = Dict(
        ColumnTypeByte => Int8,
        ColumnTypeUByte => UInt8,
        ColumnTypeBool => Bool,
        ColumnTypeShort => Int16,
        ColumnTypeUShort => UInt16,
        ColumnTypeInt => Int32,
        ColumnTypeUInt => UInt32,
        ColumnTypeLong => Int64,
        ColumnTypeULong => UInt64,
        ColumnTypeFloat => Float32,
        ColumnTypeDouble => Float64,
        ColumnTypeString => String,
        ColumnTypeJson => String,
        ColumnTypeDateTime => String,
        ColumnTypeBinary => String
    )

Tables.istable(::Type{<:FlatGeobuffer}) = true
Tables.rowaccess(::Type{<:FlatGeobuffer}) = true
Tables.rows(fgb::FlatGeobuffer) = fgb
Tables.schema(fgb::FlatGeobuffer) = Tables.Schema((map(c -> Symbol(c.name), fgb.header.columns)..., :geom), (map(x -> x.nullable ? Union{Missing,lookup[x.type]} : lookup[x.type], fgb.header.columns)..., Geometry))


text:
module Omniscape
using ArchGDAL
using Base.Threads
using Circuitscape
using DelimitedFiles
using LinearAlgebra.BLAS
using ProgressMeter
using Random
using StatsBase
using Statistics

include("structs.jl")
include("config.jl")
include("io.jl")
include("consts.jl")
include("utils.jl")
include("main.jl")
include("errors_warnings.jl")
export run_omniscape

end


text:
# Executable

using LinearAlgebra, Random, Plots
using AutomatonEdification.NeuralNetworks, AutomatonEdification.ActivationFunctions, AutomatonEdification.LossFunctions

VECTORIZED = false

# Parameters
batch_size = 32
learning_rate = 0.2
epochs = 96000

nn = NeuralNetworks.create_nn([2,1], ActivationFunctions.sigmoid, LossFunctions.mse, 2, 1)

# Learn XOR
# [a, b] -> [a XOR b]
data = shuffle(repeat([[0,0],[0,1],[1,0],[1,1]], epochs))
# Each label is a vector (even if a singleton)
labels = map(point -> [float(point[1] ⊻ point[2])], data)

@time begin
    stats = VECTORIZED ?
        NeuralNetworks.train_vectorized!(nn, data, labels, learning_rate, batch_size) :
        NeuralNetworks.train!(nn, data, labels, learning_rate, batch_size)
end

for input in [[0,0],[0,1],[1,0],[1,1]]
    print("Input: $input / " *
          "Prediction: $(NeuralNetworks.predict(nn, input)) / " *
          "Actual: $(input[1] ⊻ input[2])\n")
end

plt = plot(1:100:length(stats.losses),
           stats.losses[1:100:end],
           seriestype = :scatter)
display(plt)

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
    maps::As
    function LinearCombination{T, As}(maps::As) where {T, As}
        N = length(maps)
        sz = size(maps[1])
        for n = 1:N
            size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination"))
            promote_type(T, eltype(maps[n])) == T || throw(InexactError())
        end
        new{T, As}(maps)
    end
end

LinearCombination{T}(maps::As) where {T, As} = LinearCombination{T, As}(maps)

# basic methods
Base.size(A::LinearCombination) = size(A.maps[1])
LinearAlgebra.issymmetric(A::LinearCombination) = all(issymmetric, A.maps) # sufficient but not necessary
LinearAlgebra.ishermitian(A::LinearCombination) = all(ishermitian, A.maps) # sufficient but not necessary
LinearAlgebra.isposdef(A::LinearCombination) = all(isposdef, A.maps) # sufficient but not necessary

# adding linear maps
function Base.:(+)(A1::LinearCombination, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1.maps..., A2.maps...))
end
function Base.:(+)(A1::LinearMap, A2::LinearCombination)
    size(A1) == size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2.maps...))
end
Base.:(+)(A1::LinearCombination, A2::LinearMap) = +(A2, A1)
function Base.:(+)(A1::LinearMap, A2::LinearMap)
    size(A1)==size(A2) || throw(DimensionMismatch("+"))
    T = promote_type(eltype(A1), eltype(A2))
    return LinearCombination{T}(tuple(A1, A2))
end
Base.:(-)(A1::LinearMap, A2::LinearMap) = +(A1, -A2)

# comparison of LinearCombination objects, sufficient but not necessary
Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A) == eltype(B) && A.maps == B.maps)

# special transposition behavior
LinearAlgebra.transpose(A::LinearCombination) = LinearCombination{eltype(A)}(map(transpose, A.maps))
LinearAlgebra.adjoint(A::LinearCombination)   = LinearCombination{eltype(A)}(map(adjoint, A.maps))

# multiplication with vectors
function A_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    A_mul_B!(y, A.maps[1], x)
    l = length(A.maps)
    if l>1
        z = similar(y)
        for n=2:l
            A_mul_B!(z, A.maps[n], x)
            y .+= z
        end
    end
    return y
end
function At_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    At_mul_B!(y, A.maps[1], x)
    l = length(A.maps)
    if l>1
        z = similar(y)
        for n = 2:l
            At_mul_B!(z, A.maps[n], x)
            y .+= z
        end
    end
    return y
end
function Ac_mul_B!(y::AbstractVector, A::LinearCombination, x::AbstractVector)
    # no size checking, will be done by individual maps
    Ac_mul_B!(y, A.maps[1], x)
    l = length(A.maps)
    if l>1
        z = similar(y)
        for n=2:l
            Ac_mul_B!(z, A.maps[n], x)
            y .+= z
        end
    end
    return y
end


text:
@testset "300.longest-increasing-subsequence.jl" begin
    @test length_of_LIS([7, 7, 7, 7, 7, 7, 7]) == 1
    @test length_of_LIS([0, 1, 0, 3, 2, 3]) == 4
    @test length_of_LIS([10, 9, 2, 5, 3, 7, 101, 18]) == 4
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    unix,
    linux,
    LIBSWSCALE_VERSION_MAJOR,
    LIBSWSCALE_VERSION_MINOR,
    LIBSWSCALE_VERSION_MICRO,
    # LIBSWSCALE_BUILD,
    SWS_FAST_BILINEAR,
    SWS_BILINEAR,
    SWS_BICUBIC,
    SWS_X,
    SWS_POINT,
    SWS_AREA,
    SWS_BICUBLIN,
    SWS_GAUSS,
    SWS_SINC,
    SWS_LANCZOS,
    SWS_SPLINE,
    SWS_SRC_V_CHR_DROP_MASK,
    SWS_SRC_V_CHR_DROP_SHIFT,
    SWS_PARAM_DEFAULT,
    SWS_PRINT_INFO,
    SWS_FULL_CHR_H_INT,
    SWS_FULL_CHR_H_INP,
    SWS_DIRECT_BGR,
    SWS_ACCURATE_RND,
    SWS_BITEXACT,
    SWS_ERROR_DIFFUSION,
    SWS_MAX_REDUCE_CUTOFF,
    SWS_CS_ITU709,
    SWS_CS_FCC,
    SWS_CS_ITU601,
    SWS_CS_ITU624,
    SWS_CS_SMPTE170M,
    SWS_CS_SMPTE240M,
    SWS_CS_DEFAULT,
    SWS_CS_BT2020,
    SwsVector,
    SwsFilter,
    SwsContext

const unix = 1
const linux = 1
const LIBSWSCALE_VERSION_MAJOR = 5
const LIBSWSCALE_VERSION_MINOR = 3
const LIBSWSCALE_VERSION_MICRO = 100

# Skipping MacroDefinition: LIBSWSCALE_VERSION_INT AV_VERSION_INT ( LIBSWSCALE_VERSION_MAJOR , LIBSWSCALE_VERSION_MINOR , LIBSWSCALE_VERSION_MICRO )
# Skipping MacroDefinition: LIBSWSCALE_VERSION AV_VERSION ( LIBSWSCALE_VERSION_MAJOR , LIBSWSCALE_VERSION_MINOR , LIBSWSCALE_VERSION_MICRO )

# const LIBSWSCALE_BUILD = LIBSWSCALE_VERSION_INT

# Skipping MacroDefinition: LIBSWSCALE_IDENT "SwS" AV_STRINGIFY ( LIBSWSCALE_VERSION )
# Skipping MacroDefinition: FF_API_SWS_VECTOR ( LIBSWSCALE_VERSION_MAJOR < 6 )

const SWS_FAST_BILINEAR = 1
const SWS_BILINEAR = 2
const SWS_BICUBIC = 4
const SWS_X = 8
const SWS_POINT = 0x10
const SWS_AREA = 0x20
const SWS_BICUBLIN = 0x40
const SWS_GAUSS = 0x80
const SWS_SINC = 0x0100
const SWS_LANCZOS = 0x0200
const SWS_SPLINE = 0x0400
const SWS_SRC_V_CHR_DROP_MASK = 0x00030000
const SWS_SRC_V_CHR_DROP_SHIFT = 16
const SWS_PARAM_DEFAULT = 123456
const SWS_PRINT_INFO = 0x1000
const SWS_FULL_CHR_H_INT = 0x2000
const SWS_FULL_CHR_H_INP = 0x4000
const SWS_DIRECT_BGR = 0x8000
const SWS_ACCURATE_RND = 0x00040000
const SWS_BITEXACT = 0x00080000
const SWS_ERROR_DIFFUSION = 0x00800000
const SWS_MAX_REDUCE_CUTOFF = 0.002
const SWS_CS_ITU709 = 1
const SWS_CS_FCC = 4
const SWS_CS_ITU601 = 5
const SWS_CS_ITU624 = 5
const SWS_CS_SMPTE170M = 5
const SWS_CS_SMPTE240M = 7
const SWS_CS_DEFAULT = 5
const SWS_CS_BT2020 = 9

struct SwsVector
    coeff::Ptr{Cdouble}
    length::Cint
end

struct SwsFilter
    lumH::Ptr{SwsVector}
    lumV::Ptr{SwsVector}
    chrH::Ptr{SwsVector}
    chrV::Ptr{SwsVector}
end

struct SwsContext
end


text:
module HessMain
using Hecke
using ..HessQRModule, ..GenericRound2
import ..HessQRModule: HessQR

function GenericRound2.integral_closure(S::HessQR, F::Generic.FunctionField{T}) where {T}
  return GenericRound2._integral_closure(S, F)
end

function _gcdx(a::fmpq, b::fmpq)
  l = lcm(denominator(a), denominator(b))
  g, e, f = gcdx(numerator(a*l), numerator(b*l))
  return g//l, e, f
end

#=
  base case:
  given
    a 0
    b c
  where a, b, c are polynomials, deg b < deg c
  do Q[x]-col transforms and Z<x>-row transforms (HessQR) to get diagonal.

  Several steps.
  Input
  a*alpha 0
  b*beta  c*gamma
  where a, b, c in Q, alpha, beta, gamma in Z[x], primitive

  Step 1:
  alpha is a Z<x> unit, so
  ->
    a       0
    b*beta  c*gamma
  is valid row-transform

  Step 2:
    g = gcd(a, b) = ea + fb (via common denominator of a, b, g), so in particular
    a/g, b/g, e, f in Z

    more row-tranforms:
    e*beta       f          a      0           g*beta   c*f*gamma
    -b/g*beta    a/g    *   b*beta c gamma  =  0        a/g*c*gamma

    det(trans) = (ea+fb)/g * beta = beta is a Z<x> unit

  Step 3:
    Q[x] col. operations: since deg beta < deg gamma we get
    g*beta (c*f*gamma mod g*beta)
    0      a/g*c*gamma

  Step 4: row and col swap
    a/g*c*gamma 0
    d*delta     g*beta  (d*delta :=  (c*f*gamma mod g*beta))

    and deg delta < deg beta

  This is iterated until delta == 0
=#

function two_by_two(Q::MatElem{<:Generic.Rat{_T}}, R::PolyRing{_T}, S::HessQR) where {_T}
  @assert size(Q) == (2,2)
  @assert iszero(Q[1,2])

  Qt = base_ring(Q)
  T1 = identity_matrix(Qt, 2)
  T2 = identity_matrix(Qt, 2)

  d = reduce(lcm, [denominator(x, R) for x = Q]; init = R(1))
  c = content(d)
  # c is a unit in R
  # d//c is a unit in S
  @assert isunit(S(d*inv(c)))
  T1 = T1 * d*inv(c)
  T2 = divexact(T2, Qt(c))
  Q = T1*Q*T2

  @assert all(x->isone(denominator(x, R)), Q)

  if degree(numerator(Q[2,1])) >= degree(numerator(Q[2,2]))
    q = div(numerator(Q[2,1]), numerator(Q[2,2]))
    T = matrix(Qt, 2, 2, [Qt(1), Qt(0), Qt(-q), Qt(1)])
    T2 = T2*T
    Q = Q*T
  end

  while !iszero(Q[2,1])
    @assert all(x->isone(denominator(x, R)), Q)
    @assert iszero(Q[1,2])
    @assert degree(numerator(Q[2,1])) < degree(numerator(Q[2,2]))
    c = reduce(gcd, map(x->content(numerator(x)), Q), init = fmpq(1,1))
    Q = Q*inv(c)
    T2 = T2*inv(c)
#    @show degree(numerator(Q[2,1])), degree(numerator(Q[2,2]))
#    @show map(x->(degree(numerator(x)), degree(denominator(x))), Array(T1))
#    @show map(x->(degree(numerator(x)), degree(denominator(x))), Array(T2))
#    @show map(x->(degree(numerator(x)), degree(denominator(x))), Array(Q))
#    @show map(x->max(nbits(maximum(map(height, coefficients(numerator(x))), init = fmpz(1))), nbits(maximum(map(height, coefficients(denominator(x))), init = fmpz(1)))), Array(Q))

    n, d = integral_split(Q[1,1], S)
    @assert isconstant(d)
    a = n.c//d.c
    Q[1,1] = a

    c = Qt(n)//Qt(d)*inv(Q[1,1])
    @assert isunit(c)
    T1[1,:] *= inv(c)
    n, d = integral_split(Q[2,1], S)
    b = n.c//d.c
    beta = Q[2,1]//b

    g, e, f = _gcdx(a, b)
    T = matrix(Qt, 2, 2, [e*beta, Qt(f), -Q[2,1]//g, Q[1,1]//g])
    Q = T*Q

    @assert iszero(Q[2,1])
    T1 = T*T1

    if iszero(Q[1,2])
      return Q, T1, T2
    end

    n, d = integral_split(Q[1,1], R)
    @assert isone(d)
    nn, d = integral_split(Q[1,2], R)
    @assert isone(d)
    @assert degree(nn) > degree(n)
    q, r = divrem(nn, n)

    T = matrix(Qt, 2, 2, [1, -q, Qt(0), Qt(1)])
    Q = Q*T
    T2 = T2 * T

    swap_cols!(Q, 1, 2)
    swap_cols!(T2, 1, 2)
    swap_rows!(Q, 1, 2)
    swap_rows!(T1, 1, 2)
  end

  @assert isone(integral_split(T1, S)[2])
  @assert isone(integral_split(T2, R)[2])
  @assert isunit(integral_split(det(T1), S)[1])
  @assert isunit(integral_split(det(T2), R)[1])

  return Q, T1, T2
end

function clear_triangle!(T, TT1, TT2, R, S, l::Int)
  n = nrows(T)
  Qt = base_ring(T)
  @show :clear_triang, l
  for i=n-1:-1:l
    for j=i+1:n
      @assert i<j
      @show i, j
      q, t1, t2 = two_by_two(T[ [i,j], [i,j]], R, S)
      @assert t1*T[ [i,j], [i,j]]*t2 == q

      TT = identity_matrix(Qt, n)
      TT[[i,j], [i,j]] = t1
      TT1 = TT*TT1
      T = TT*T
      TT[[i,j], [i,j]] = t2
      TT2 = TT2 * TT
      T = T*TT
#      display(map(Int, map(iszero, T)))
#      @show :herew
    end
    if i<n-2
      clear_triangle!(T, TT1, TT2, R, S, i+2)
    end
  end
  display(map(Int, map(iszero, T)))
  display(map(degree, map(numerator, T)))
  display(map(degree, map(denominator, T)))
  @show :herew, l, n
  @assert all(i->all(j -> iszero(T[i,j]), i+1:n), l:n)
end

function GenericRound2.integral_closure(Zx::FmpzPolyRing, F::Generic.FunctionField)
  Qt = base_ring(F)
  t = gen(Qt)
  S = HessQR(Zx, Qt)
  R = parent(numerator(t))
  o1 = integral_closure(S, F)
  o2 = integral_closure(R, F)
  if isdefined(o1, :trans)
    T = o1.trans
  else
    T = identity_matrix(Qt, degree(F))
  end
  if isdefined(o2, :itrans)
    T = T * o2.itrans
  end
  q, w = integral_split(T, R)
  h, T2 = Hecke._hnf_with_transform(q', :upperright)
#    display(map(Int, map(iszero, h)))
#    @show :hallo
  T = map_entries(Qt, h')
    display(map(Int, map(iszero, T)))
  T = divexact(T, Qt(w))
#TODO: we don't need TT2 other than to debug assertions
# make it optional? tricky to also do this in two_by_two...
  TT2 = map_entries(Qt, T2')
  TT1 = identity_matrix(Qt, degree(F))
  cnt = 0

  for i=degree(F)-2:-1:1
    clear_triangle!(T, TT1, TT2, R, S, i)
  end

  @assert isdiagonal(T)

  # the diagonal in Q(t) is split into a/b * alpha/beta where
  #  a/b in Q (hence a unit there)
  # and alpha, beta in Z[x] primitive, so alpha/beta is a unit in Z<x>
  for i=1:degree(F)
    n, d = integral_split(T[i,i], S)
    @assert isconstant(d)
    u = Qt(n.f)//Qt(n.g)
#    @assert n.c//d.c*u == T[i,i]
    TT2[:, i] *= Qt(d.c)*inv(Qt(n.c))
    TT1[i, :] *= inv(u)
    T[i,i] = 1
#  @assert TT1*o1.trans*o2.itrans*TT2 == T
  end

  TT1 = TT1
  n, d = integral_split(TT1, Zx)
  @assert map_entries(Qt, n) == TT1 * Qt(d)
  o3 = GenericRound2.Order(Zx, F)
  if isdefined(o1, :trans)
    return GenericRound2.Order(o3, integral_split(map_entries(Qt, TT1)*o1.trans, Zx)..., check = false)
  else
    return GenericRound2.Order(o3, integral_split(map_entries(Qt, TT1), Zx)..., check = false)
  end
  return GenericRound2.Order(o1, TT1, one(S)), GenericRound2.Order(o2, inv(TT2'), one(base_ring(TT2)))
end

function Base.denominator(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  return integral_split(a, S)[2]
end

function Base.numerator(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  return integral_split(a, S)[1]
end

function Hecke.integral_split(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  #TODO: feels too complicated....
  if iszero(a)
    return zero(S), one(S)
  end
  n = numerator(a)
  d = denominator(a)
  dn = reduce(lcm, map(denominator, coefficients(n)), init = fmpz(1))
  dd = reduce(lcm, map(denominator, coefficients(d)), init = fmpz(1))
  zn = S(n*dn)
  zd = S(d*dd)
  cn = content(zn)
  cd = content(zd)
  zn = divexact(zn, cn)
  zd = divexact(zd, cd)
  cn *= dd
  cd *= dn
  g = gcd(cn, cd)
  cn = divexact(cn, g)
  cd = divexact(cd, g)
  return cn*zn, cd*zd
end

function (S::FmpzPolyRing)(a::Generic.Rat{fmpq})
  n, d = integral_split(a, S)
  @assert isone(d)
  return n
end

end

using .HessMain

#=
  this should work:

Hecke.example("Round2.jl")

?GenericRound2

Qt, t = RationalFunctionField(QQ, "t")
Qtx, x = PolynomialRing(Qt, "x")
F, a = FunctionField(x^6+27*t^2+108*t+108, "a")
integral_closure(parent(denominator(t)), F)
integral_closure(Localization(Qt, degree), F)
integral_closure(Hecke.Globals.Zx, F)
basis(ans, F)
derivative(F.pol)(gen(F)) .* ans #should be integral

k, a = wildanger_field(3, 8*13)
integral_closure(ZZ, k)
integral_closure(Localization(ZZ, 2), k)

more interesting and MUCH harder:

G, b = FunctionField(x^6 + (140*t - 70)*x^3 + 8788*t^2 - 8788*t + 2197, "b")

=#


text:
using Revise
using HybridFunnelGraphs

# graph = create_funnel_graph("test/hpd/domain.hpd", "test/hpd/problem.hpd"; max_levels=170)
graph = create_funnel_graph("test/hpd/dom2.hpd", "test/hpd/prob2.hpd"; max_levels=300, has_placement_constraint=false)
graph.num_levels


text:
function initialize_fields(particles::Particles{N}, grid::Grid{N}) where {N}
    padded_grid = padded(grid)

    potential_source = zeros(size(padded_grid))
    sample_charge_density!(potential_source, padded_grid, particles)
    potential_source .*= VACUUM_PERMITTIVITY⁻¹

    potential = zeros(size(padded_grid))
    solve_poisson!(
        potential,
        potential_source,
        padded_grid,
        Free,
        GaussSeidel(),
    )
    potential .*= -1

    fields = Fields(grid)
    grad_padded!(fields.E, potential, padded_grid, Up)
    fields
end

function advance!(particles::Particles{N}, fields::Fields{N}, Δt) where {N}
    update_accelerations!(particles, fields)
    advance_positions!(particles, 0.5 * Δt)
    advance_velocities!(particles, 0.5 * Δt)
    update_j!(fields, particles)
    update_curl_E!(fields)
    advance_B!(fields, 0.5 * Δt)
    update_curl_B!(fields)
    advance_E!(fields, Δt)
    advance_B!(fields, 0.5 * Δt)
    advance_positions!(particles, 0.5 * Δt)
    advance_velocities!(particles, 0.5 * Δt)
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


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

JLLWrappers.@generate_main_file_header("libmseed")
JLLWrappers.@generate_main_file("libmseed", UUID("ccc2c699-b150-5417-88f2-a95a0d1581d9"))
end  # module libmseed_jll


text:
abstract MatchSteering
immutable NotGreedy <: MatchSteering
	x
end

isgreedy(x::NotGreedy) = false
isgreedy(x) = true

trymatch(ms::MatchSteering, val) = trymatch(ms.x, val)
trymatch(f::Function, val) = f(val)
trymatch(val1, val2) = val1 == val2

function ismatch{N}(x, atoms::Vararg{Any, N}; matchall = false)
	matches = NTuple{N, Vector{Any}}[]
	current_atom = 1
	states, tmp_matches = [], ntuple(x->Any[], Val{N})
	elem, state = next(x, start(x))
	matched_state = state
	last_state = state
	while true
		atom = atoms[current_atom]
		if !trymatch(atom, elem) # reset if not matching
			current_atom = 1
			atom = atoms[current_atom]
			tmp_matches = ntuple(x->Any[], Val{N})
		end
		while !trymatch(atom, elem) # continue until done or found match
			done(x, state) && return matches, states
			elem, state = next(x, state)
		end
		# if at first atom, we remember the first match, else don't modify it!
		matched_state = ifelse(current_atom == 1, state, matched_state)
		ismatched = false
		while trymatch(atom, elem)
			ismatched = true
			push!(tmp_matches[current_atom], elem)
			done(x, state) && break
			elem, state = next(x, state)
			next_atom = atoms[min(current_atom+1, length(atoms))]
			isgreedy(atom) || trymatch(next_atom, elem) && break
		end
		current_atom = ifelse(ismatched, current_atom+1, 1)
		if current_atom > length(atoms)
			push!(states, matched_state)
			push!(matches, tmp_matches)
			matchall || break
			tmp_matches = ntuple(x->Any[], Val{N})
			current_atom = 1
		end
	end
	matches, states
end

function forward(x, state, elem, n)
	for i=1:n
		done(x, state) && break
		elem, state = next(x, state)
	end
	elem, state
end

function _replace(f, x, atoms...; matchall = false)
	matches, states = ismatch(x, atoms..., matchall=matchall)
	if isempty(matches)
		return copy(x)
	end
	result = similar(x, 0)
	(elem, state), i = next(x, start(x)), 1
	while !done(x, state)
		isreplace = i <= length(states) && state == states[i]
		replacements, n = if isreplace
			i += 1
			n = sum(map(length, matches[i-1]))
			tmp = f(matches[i-1]...)
			r = isa(tmp, Tuple) ? tmp : (tmp,)
			r, n
		else
			(elem,), 1
		end
		elem, state = forward(x, state, elem, n)
		for r in replacements
			push!(result, r)
		end
	end
	result

end
[1,2]
[(1,2,3),4,5]

isunless(x::Expr) = x.head == :gotoifnot
isunless(x) = false
islabelnode(x::LabelNode) = true
islabelnode(x) = false
function remove_goto(ast)
	_replace(ast, isunless, (x->true)*(*), islabelnode) do unless, body, label
		condition = unless[].args[1]
		ifbody = Expr(:block, remove_goto(body)...)
		Expr(:if, condition, ifbody)
	end
end


text:
__precompile__()

module SRTM15plus

    using NetCDF
    include("elevations.jl")

end # module


text:
import Base: ==

function decLULu(iscirc, band, lastcols, lastrows)
    wd = size(band, 1)
    kl, ku = get_kl_ku(wd)
    szb = size(band, 2)
    n = iscirc ? size(lastrows, 2) : szb
    begrow = n - ku
    begcol = n - kl
    for k = 1:(iscirc ? begrow : n)
        pivot = band[ku+1, k]
        for i = ku+2:wd
            band[i, k] /= pivot
        end
        if iscirc
            for i = 1:ku
                lastrows[i, k] /= pivot
            end
        end
        #        println("trace1 k=$k")
        for i = 1:kl, j = 1:ku
            if k + j <= szb
                #                println("i=$i j=$j avant band=$(band[ku+1+i-j, k+j]) b1=$(band[ku+1+i, k]) b2=$(band[ku+1+j, k+j])")
                # A[k+i,k+j] -= A[k+i, k]*A[k, k+j]
                band[ku+1+i-j, k+j] -= band[ku+1+i, k] * band[ku+1-j, k+j]
                #                println("i=$i j=$j apres band=$(band[ku+1+i-j, k+j])")
            end
        end
        if iscirc
            bkl = min(kl, begrow - k)
            for i = 1:kl, j = 1:bkl
                # A[k+j, begcol+i] = A[k+j, k]*A[k, begcol+i]
                lastcols[k+j, i] -= band[ku+1+j, k] * lastcols[k, i]
            end
            bku = min(ku, szb - k)
            for i = 1:bku, j = 1:ku
                # A[begrow+j, k+i] = A[begrow+j, k]*A[k, k+i]
                lastrows[j, k+i] -= lastrows[j, k] * band[ku+1-i, k+i]
            end
            for i = 1:ku, j = 1:kl
                # A[begrow+i, begcol+j] = A[begrow+i, k]*A[k, begcol+j]
                lastrows[i, begcol+j] -= lastrows[i, k] * lastcols[k, j]
            end
        end
    end
    if iscirc
        for k = begrow+1:n
            i_k = k - begrow
            pivot = lastrows[i_k, k]
            for i = i_k+1:ku
                lastrows[i, k] /= pivot
            end
            for i = i_k+1:ku, j = k+1:n
                # A[i+begrow, j] -= A[i+begrow, k]*A[k, j]
                lastrows[i, j] -= lastrows[i, k] * lastrows[k-begrow, j]
            end
        end
    end
end

"""
    struct LuSpline{T}
    LuSpline(n, t::Vector{T}; iscirc=true, isLU=true) where{T}

Structure of a LU decomposition of circular banded matrix,
a LU decomposition can be stored in a Matrix which is equal to L + U - I.
For a circular Banded matrix all non zero coefficients are in the band and in the last columns and lines

# Implementation
- band::Matrix{T} : matrix of size (kl+ku+1, n-kl)
- ku : size of the band upper the diagonal
- kl : size of the band lower the diagonal
- iscirc : true if and only if original matrix is circular
- isLU : true if LU decomposition has been perform
- lastcols : only in circular case, Matrix of size (n-ku, kl) that represents teh last columns of matrix
- lastrows : only in circular case, Matrix of size (n, ku) that represents the last rows of matrix

# Arguments
- 'n' : size of the matrix
- 't::Vector{T}` : vector of all values, the size is order+1, where order is the order of the spline.

"""
struct LuSpline{T}
    band::Matrix{T}
    ku::Int64
    kl::Int64
    iscirc::Bool
    isLU::Bool
    lastcols::Union{Matrix{T},Missing} # missing when iscirc=false
    lastrows::Union{Matrix{T},Missing} # missing when iscirc=false
    function LuSpline(n, t::Vector{T}; iscirc = true, isLU = true) where {T}
        wd = size(t, 1) # band width
        kl, ku = get_kl_ku(wd)
        szb = iscirc ? n - kl : n
        band = zeros(T, wd, szb)
        for i = 1:wd
            jbeg = i <= ku + 1 ? ku - i + 2 : 1
            if iscirc
                jend = i >= kl + 2 ? szb - i + kl + 1 : szb
            else
                jend = i >= ku + 2 ? n - i + ku + 1 : n
            end
            for j = jbeg:jend
                band[i, j] = t[i]
            end
        end
        if iscirc
            lastrows = zeros(T, ku, n)
            lastcols = zeros(T, n - ku, kl)
            for i = 1:ku
                for ind = 1:wd
                    j = n - wd + i + ind
                    lastrows[i, (j-1)%n+1] = t[ind]
                end
            end
            for i = 1:kl
                for j = 1:kl+1-i
                    lastcols[j, j+i-1] = t[i]
                    lastcols[n-kl-ku+i+j-1, j] = t[i]
                end
            end
        else
            lastcols = lastrows = missing
        end
        if isLU
            decLULu(iscirc, band, lastcols, lastrows)
        end
        return new{T}(band, ku, kl, iscirc, isLU, lastcols, lastrows)
    end
    """
        LuSpline( A::Matrix{T}, ku, kl; iscirc=true, isLU=false) where {T}

    Contructor from a matrix, for test only

    """
    function LuSpline(A::Matrix{T}, ku, kl; iscirc = true, isLU = false) where {T}
        n = size(A, 1)
        if iscirc
            lastrows = copy(A[end-ku+1:end, :])
            lastcols = copy(A[1:end-ku, end-kl+1:end])
        else
            lastcols = lastrows = missing
        end
        szb = iscirc ? n - kl : n
        wd = kl + ku + 1
        band = zeros(T, wd, szb)
        for j = 1:szb
            for i = j-ku:j+kl
                if 1 <= i <= n - (iscirc ? ku : 0)
                    band[ku+i+1-j, j] = A[i, j]
                end
            end
        end
        return new{T}(band, ku, kl, iscirc, isLU, lastcols, lastrows)
    end
end
function ==(la::LuSpline{T}, lb::LuSpline{T}) where {T}
    return (
        la.ku == lb.ku &&
        la.kl == la.kl &&
        la.iscirc == lb.iscirc &&
        la.isLU == lb.isLU &&
        la.band == lb.band &&
        (!la.iscirc || (la.lastrows == lb.lastrows && la.lastcols == lb.lastcols))
    )
end
function sol!(X::AbstractVector{I}, spA::LuSpline{T}, Y::AbstractVector{I}) where {T,I}
    szb = size(spA.band, 2)
    n = spA.iscirc ? size(spA.lastrows, 2) : szb
    begrow = n - spA.ku
    begcol = n - spA.kl
#    copyto!(Y, b)
    endmat = spA.iscirc ? begrow : n
    endmat2 = spA.iscirc ? begcol : n
    for i = 2:endmat
        fin = i - 1
        deb = max(1, i - spA.kl)
        Y[i] -= sum([Y[j] * spA.band[spA.ku+1+i-j, j] for j = deb:fin])
    end
    if spA.iscirc
        for i = begrow+1:n
            Y[i] -= sum(Y[1:i-1] .* spA.lastrows[i-begrow, 1:i-1])
        end
    end
#    X = zeros(T, n)
    fill!(X, zero(eltype(X)))
    if spA.iscirc
        for i = n:-1:begrow+1
            X[i] =
                (Y[i] - sum(X[i+1:n] .* spA.lastrows[i-begrow, i+1:n])) /
                spA.lastrows[i-begrow, i]
        end
    end
    for i = endmat:-1:1
        deb = i + 1
        fin = min(i + spA.ku, endmat2)
        if deb <= fin
            s = sum(spA.band[spA.ku+1+i-j, j] * X[j] for j = deb:fin)
        else
            s = zero(eltype(Y))
        end
        if spA.iscirc
            s += sum(spA.lastcols[i, 1:spA.kl] .* X[end-spA.kl+1:end])
        end
        X[i] = (Y[i] - s) / spA.band[spA.ku+1, i]
    end
    return X, Y
end
sol(spA::LuSpline{T}, b::AbstractVector) where {T} = sol!(zeros(eltype(b), size(b, 1)), spA, copy(b))
# get_n(sp::LuSpline)=sp.iscirc ? size(sp.lastrows, 2) : size(sp.band, 2)
get_order(sp::LuSpline) = sp.ku + sp.kl + 1

"""
    B_SplineLU{T, edge, order} <: AbstractInterpolation{T, edge, order}

Type containing spline coefficients for b-spline interpolation

# Type parameters
- `T` : the type of data that is interpolate
- `edge::TypeEdge=CircEdge` : true if function is circular
- `order::Int` : order of lagrange interpolation

# Implementation :
- `ls::LuSpline{T}` : the LU matrix
- `tabfct::Vector{Polynomial{T}}` : function table for interpolation

# Arguments :
- `n` : size of the matrix
- `order` : the order of interpolation
- `[T::DataType=Float64]` : The type values to interpolate

"""
struct B_SplineLU{T,edge,order} <: B_Spline{T,edge,order}
    ls::LuSpline{T}
    tabfct::Vector{Polynomial{T}}
    function B_SplineLU(order::Int, n::Int, T::DataType = Float64)
        (order % 2 == 0) && throw(
            ArgumentError(
                "order=$order B_SplineLU for even  order is not implemented n=$n",
            ),
        )
        bspline = getbspline(order, 0)
        tabfct_rat = map(x -> bspline[order-x](Polynomial([order - x, 1])), 0:order)
        # bspline = SplineInt(order)
        # N = typeof(bspline.fact_order)
        # tabpol = map(x -> bspline[order-x](Polynomial([order-x,1])), 0:order)
        ls = LuSpline(n, convert.(T, bspline.(1:order)), iscirc = true, isLU = true)
        return new{T,CircEdge,order}(ls, convert.(Polynomial{T}, tabfct_rat))
    end
    B_SplineLU(o::Int, n::Int, elt::T; kwargs...) where {T<:Number} =
        B_SplineLU(o, n, T; kwargs...)
end

sol!(X::AbstractVector{I}, bsp::B_SplineLU{T}, Y::AbstractVector{I}) where {T,I} = sol!(X, bsp.ls, Y)[1]

sol(bsp::B_SplineLU{T}, b::AbstractVector{I}) where {T<:Number,I} = sol!(zeros(I, length(b)), bsp, copy(b))


text:
main() = let part1 = 0, part2 = 0
    groups = split.(split(read("input.txt", String), "\n\n"))
    part1 = sum(length(union(g...)) for g in groups)
    part2 = sum(length(intersect(g...)) for g in groups)
    part1, part2
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
# Autogenerated wrapper script for Opus_jll for i686-w64-mingw32
export libopus

JLLWrappers.@generate_wrapper_header("Opus")
JLLWrappers.@declare_library_product(libopus, "libopus-0.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libopus,
        "bin\\libopus-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PodPresetSpec is a description of a pod preset.

    IoK8sApiSettingsV1alpha1PodPresetSpec(;
        env=nothing,
        envFrom=nothing,
        selector=nothing,
        volumeMounts=nothing,
        volumes=nothing,
    )

    - env::Vector{IoK8sApiCoreV1EnvVar} : Env defines the collection of EnvVar to inject into containers.
    - envFrom::Vector{IoK8sApiCoreV1EnvFromSource} : EnvFrom defines the collection of EnvFromSource to inject into containers.
    - selector::IoK8sApimachineryPkgApisMetaV1LabelSelector : Selector is a label query over a set of resources, in this case pods. Required.
    - volumeMounts::Vector{IoK8sApiCoreV1VolumeMount} : VolumeMounts defines the collection of VolumeMount to inject into containers.
    - volumes::Vector{IoK8sApiCoreV1Volume} : Volumes defines the collection of Volume to inject into the pod.
"""
mutable struct IoK8sApiSettingsV1alpha1PodPresetSpec <: SwaggerModel
    env::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvVar} } # spec name: env
    envFrom::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvFromSource} } # spec name: envFrom
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    volumeMounts::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeMount} } # spec name: volumeMounts
    volumes::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Volume} } # spec name: volumes

    function IoK8sApiSettingsV1alpha1PodPresetSpec(;env=nothing, envFrom=nothing, selector=nothing, volumeMounts=nothing, volumes=nothing)
        o = new()
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("env"), env)
        setfield!(o, Symbol("env"), env)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("envFrom"), envFrom)
        setfield!(o, Symbol("envFrom"), envFrom)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("volumeMounts"), volumeMounts)
        setfield!(o, Symbol("volumeMounts"), volumeMounts)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("volumes"), volumes)
        setfield!(o, Symbol("volumes"), volumes)
        o
    end
end # type IoK8sApiSettingsV1alpha1PodPresetSpec

const _property_map_IoK8sApiSettingsV1alpha1PodPresetSpec = Dict{Symbol,Symbol}(Symbol("env")=>Symbol("env"), Symbol("envFrom")=>Symbol("envFrom"), Symbol("selector")=>Symbol("selector"), Symbol("volumeMounts")=>Symbol("volumeMounts"), Symbol("volumes")=>Symbol("volumes"))
const _property_types_IoK8sApiSettingsV1alpha1PodPresetSpec = Dict{Symbol,String}(Symbol("env")=>"Vector{IoK8sApiCoreV1EnvVar}", Symbol("envFrom")=>"Vector{IoK8sApiCoreV1EnvFromSource}", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("volumeMounts")=>"Vector{IoK8sApiCoreV1VolumeMount}", Symbol("volumes")=>"Vector{IoK8sApiCoreV1Volume}")
Base.propertynames(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }) = collect(keys(_property_map_IoK8sApiSettingsV1alpha1PodPresetSpec))
Swagger.property_type(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiSettingsV1alpha1PodPresetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }, property_name::Symbol) =  _property_map_IoK8sApiSettingsV1alpha1PodPresetSpec[property_name]

function check_required(o::IoK8sApiSettingsV1alpha1PodPresetSpec)
    true
end

function validate_property(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }, name::Symbol, val)
end


text:
include("./make.jl")

deploydocs(
    repo = "github.com/$(ENV["GITHUB_REPOSITORY"]).git",
)


text:
function conAIC(m::LMMScalar1)
  q=length(m.L); theta=m.theta; p,n = size(m.Xt);
  ## note
  ## L*Rzx =lambdat *ZtX  => Rzx =theta*L^-1 *ZtX
  ## Hl = [ L ,0 ; Rzx , Rx]
  Hlinv=inv(tril(hvcat((2,2),diagm(m.L),(theta ./ m.L).*m.XtZ',(theta ./ m.L)'.*m.XtZ, m.RX.UL')))
  H= Hlinv' * Hlinv
  #making Zt
  Zt=zeros(Float64,q,n)
  for i in 1:n
    j = m.Ztrv[i]; z = m.Ztnz[i]
    Zt[j,i] = z
  end
  # rho= tr ( [Z*lambda; X] H [lambda] [lambdat*Zt, Xt] )
  rho = trace( hcat(theta*Zt',m.Xt') * H * vcat(theta*Zt, m.Xt))
  #REML
  #rho=(n-p-1.0)*(rho+1.0)/(n-p-2.0) + (p+1.0)/(n-p-2.0)
  #MLE
  rho=n/(n-p-2.0) * ((rho+1) -(rho-p)/(n-p))
  (objectivecon(m)+2.0*rho, rho)
end

function conAICb(m::LinearMixedModel)
  n,p,q,k = size(m)
  spλtZt=sparse(m.s.λtZt)
  k12= spλtZt * sparse(m.X.m)
  #Hlinv = inv(hvcat((2,2),m.s.λtZt *m.s.λtZt' +eye(q) ,k12,k12',sparse(m.X.m)'*sparse(m.X.m)))
  Hlinv = inv(tril(full([ spλtZt *spλtZt'+speye(q) k12 ; k12' sparse(m.X.m)'*sparse(m.X.m)])))
  rho = trace( vcat(spλtZt,sparse(m.X.m)')*vcat(spλtZt,sparse(m.X.m)')' *Hlinv)
  #REML
  #rho=(n-p-1.0)*(rho+1.0)/(n-p-2.0) + (p+1.0)/(n-p-2.0)
  #MLE
  rho=n/(n-p-2.0) * ((rho+1) -(rho-p)/(n-p))
  (-2.0*condll(m)+2.0*rho, rho)
end

function conAICa(m::LinearMixedModel)
  n,p,q,k = size(m)
  #k12= m.s.λtZt * sparse(m.X.m)
  #Hlinv = inv(hvcat((2,2),m.s.λtZt *m.s.λtZt' +eye(q) ,k12,k12',sparse(m.X.m)'*sparse(m.X.m)))
  #Hlinv = inv([ m.s.λtZt *m.s.λtZt'+eye(q) k12 ; k12' sparse(m.X.m)'*sparse(m.X.m)])
  spλtZt=sparse(m.s.λtZt)
  RZXa = spλtZt * m.X.m
  for j in 1:size(RZXa,2)
            permute!(view(RZXa,:,j),m.s.perm)
  end
  RZXa = Base.LinAlg.CHOLMOD.solve(m.s.L, RZXa, CHOLMOD_L)
  Hlinv2=inv(tril(hvcat((2,2),sparse(m.s.L),RZXa,RZXa',m.s.RX.UL')))
  Hlinv= Hlinv2' * Hlinv2
  rho = trace( vcat(spλtZt,sparse(m.X.m)')*vcat(spλtZt,sparse(m.X.m)')' *Hlinv)
  #REML
  #rho=(n-p-1.0)*(rho+1.0)/(n-p-2.0) + (p+1.0)/(n-p-2.0)
  #MLE
  rho=n/(n-p-2.0) * ((rho+1) -(rho-p)/(n-p))
  (-2.0*condll(m)+2.0*rho, rho)
end

function objectivecon(m::LinearMixedModel)
    n,p,q,k = size(m); fn = float64(n - (m.REML ? p : 0))
    sig2 = std(m)[end][1]^2
    n*log(2*pi*sig2)+sumfdiff(Abs2Fun(), m.y, m.mu)/sig2
end

function TRCderHat(f::LinearMixedModel, fo::Formula, fr::AbstractDataFrame)
    rhohat=0.0
    n,p,q,k = size(f)

    for i in 1:n
      x=f.y[i]
      h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
      xph = x + h
      dx = xph - x

      newdf = deepcopy(fr)
      newdf[fo.lhs]=convert(DataArray{Float64,1},newdf[fo.lhs])
      newdf[fo.lhs][i]=xph
      f1= fit(lmm(fo,newdf,true)).mu[i]
      f0 = f.mu[i]
      rhohat += (f1 - f0) / dx
    end

    rhohat
end

function HAIC(m::LMMGeneral)
  n,p,q,k = size(m)
  k12= m.LambdatZt * sparse(m.X.m)
  #Hlinv = inv(hvcat((2,2),m.LambdatZt *m.LambdatZt' +eye(q) ,k12,k12',sparse(m.X.m)'*sparse(m.X.m)))
  Hlinv = inv([ m.LambdatZt *m.LambdatZt'+eye(q) k12 ; k12' sparse(m.X.m)'*sparse(m.X.m)])
  vcat(m.LambdatZt,sparse(m.X.m)')' *Hlinv*vcat(m.LambdatZt,sparse(m.X.m)')
end

function TRCderHat2(f::LinearMixedModel, G::Array{Float64},
                    fo::Formula, fr::AbstractDataFrame)

    n,p,q,k = size(f)
    Hout = eye(n)

    for i in 1:n
      for j in 1:n
      x=f.y[j]
      # pick a small value for h
      h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
      # floating point arithmetic gymnastics
      xph = x + h
      dx = xph - x

      # evaluate f at x + h
      newdf = deepcopy(fr)
      #newdf[fo.lhs]=convert(DataArray{Float64,1},newdf[fo.lhs])
      newdf[fo.lhs][j]=xph
      f1= HAIC(fit(lmm(fo,newdf,true)))[i,j]
      # evaluate f at x
      f0 = G[i,j]
      # divide the difference by h
      Hout[i,j]= (f1 - f0) / dx
      #Hout[j,i]=Hout[i,j]
    end
    end
    ones(n)'*Hout*f.y
    #Hout
end

#????? ==> newdf = deepcopy(fr)???
function derHat(f::LinearMixedModel, idx::Int64, fo::Formula, fr::AbstractDataFrame)

    return function(x)

      # pick a small value for h
      h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
      # floating point arithmetic gymnastics
      xph = x + h
      dx = xph - x

      # evaluate f at x + h
      newdf = fr
      #newdf[fo.lhs]=convert(DataArray{Float64,1},newdf[fo.lhs])
      newdf[fo.lhs][idx]=xph
      f1= fit(lmm(fo,newdf,true)).mu[idx]
      # evaluate f at x
      f0 = f.mu[idx]
      # divide the difference by h
      (f1 - f0) / dx
    end

end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
module Lorentz

using Vector2d
using Docile

#using ImmutableArrays
#import Base.convert

export run_Lorentz_gas
export run_many_particles

export Particle, Disc, collision_time, Plane

type Particle
    x::Vector2D
    v::Vector2D
end

immutable Disc
    centre::Vector2D
    radius::Float64
end

@doc """Compute *time of collision* of a particle and a disc,
        assuming the particle starts outside the disc and the
        particle speed is one (using the quadratic formula).

        Returns -1 if no collision.
        """ ->
function collision_time(p::Particle, disc::Disc)

    disp = p.x-disc.centre

    B = dot(p.v, disp)
    C = dot(disp, disp) - disc.radius*disc.radius

    discriminant = B*B - C

    if discriminant < 0
        return -1.
    end

    return -B - sqrt(discriminant)  # NB: this supposes that the particle is *outside* the disc

end

@doc """Compute normal vector to disc boundary, at a point x that must lie *on the boundary*""" ->
function normal(disc::Disc, x)
    return (x - disc.centre) / disc.radius
end

immutable Plane
    c::Vector2D  # point on plane
    normal::Vector2D
end

#Plane(c::Vector{Float64}, normal::Vector{Float64}) =
#    Plane(Vector2D(c), Vector2D(normal))

collision_time(p::Particle, plane::Plane) =
    dot(plane.c - p.x, plane.normal) / dot(p.v, plane.normal)

normal(plane::Plane, x) = plane.normal

function collision(p::Particle, current_cell, boundaries, jump_directions, previous_hit)

    min_collision_time = 100.0
    which_hit = -1

    for i in 1:5
        if i == previous_hit
            continue
        end

        @inbounds collision_t = collision_time(p, boundaries[i])

        if collision_t > 0.0 && collision_t < min_collision_time
            which_hit, min_collision_time = i, collision_t
        end
    end

    x_new = p.x + p.v*min_collision_time
    new_cell = current_cell

    @assert which_hit > 0

    if which_hit == 1  # reflect off disc
        @inbounds n = normal(boundaries[1], x_new)

        v_new = p.v - 2.0*n*dot(n, p.v)  # reflejar solo si es disco; si es plano, pasar a traves

        speed = sqrt(dot(v_new, v_new))
        v_new /= speed

    else
        # hit plane, so periodise
        @inbounds jump = jump_directions[which_hit]
        # print "jump = ", jump

        new_cell += jump
        x_new -= jump  # this is supposing that we are in exactly unit cell

        v_new = p.v

        which_hit += 2
        if which_hit > 5
            which_hit -= 4
        end
    end

    return x_new, v_new, new_cell, min_collision_time, which_hit
end

function initial_condition(radius)

    x, y = 0.0, 0.0

    while true
        x, y = rand(2) .- 0.5

        if (x*x + y*y) >= radius^2
            break
        end
    end

    # velocity in 2D:
    theta = rand() * 2*pi

    xx = Vector2D(x, y)
    vv = Vector2D(cos(theta), sin(theta))

    return xx, vv
end

# fast version:
function collision_time_with_disc(x, v, r)
    B = dot(v, x);
    C = dot(x, x) - r*r;

    discriminant = B*B - C;

    if (discriminant < 0)
        return -1.0;

    else
        return -B - sqrt(discriminant);
    end
end

intersection_horiz_boundary(xx, vv, y) =
    @inbounds return (y - xx[2]) / vv[2]

intersection_vert_boundary(xx, vv, x) =
    @inbounds return (x - xx[1]) / vv[1]

function new_collision(xx, vv, r, cell, previous_hit)

    min_collision_time = 1000.0

    collision_times =
        ( collision_time_with_disc(xx, vv, r),
            intersection_vert_boundary(xx, vv, -0.5),
            intersection_horiz_boundary(xx, vv, 0.5),
            intersection_vert_boundary(xx, vv, 0.5),
            intersection_horiz_boundary(xx, vv, -0.5)
        )

    object_hit = -1

    for i in 1:5  # num of collision times

        if i==previous_hit
            continue
        end

        @inbounds coll_time = collision_times[i]

        if coll_time > 0.0 && coll_time < min_collision_time

            min_collision_time = coll_time
            object_hit = i
        end
    end

    #@assert object_hit > 0

    #println(min_collision_time)

    xx += vv*min_collision_time

    @inbounds begin
        x, y = xx[1], xx[2]
        x_cell, y_cell = cell[1], cell[2]
    end

    if object_hit == 1
        n = xx / r  # vector normal

        vv -= 2. * dot(vv, n) * n
        speed = sqrt(dot(vv, vv))
        vv /= speed

    elseif object_hit == 2
        x = 0.5
        x_cell -= 1
        object_hit = 4

    elseif object_hit == 3
        y = -0.5
        y_cell += 1
        object_hit = 5

    elseif object_hit == 4
        x = -0.5
        x_cell += 1
        object_hit = 2

    elseif object_hit == 5
        y = 0.5
        y_cell -= 1
        object_hit = 3
    end

    cell = Vector2D(x_cell, y_cell)
    xx = Vector2D(x, y)

    return xx, vv, cell, min_collision_time, object_hit

end

### Main program

function make_boundaries(radius)

    boundaries = {}

    push!(boundaries, Disc([0., 0.], radius) )

    push!(boundaries, Plane([-0.5, -0.5], [0., -1.]) )
    push!(boundaries, Plane([0.5, -0.5], [1., 0.]) )
    push!(boundaries, Plane([0.5, 0.5], [0., 1.]) )
    push!(boundaries, Plane([-0.5, 0.5], [-1., 0.]) )

    jump_directions = [Vector2D(0.,0.)]  # dummy

    for boundary in boundaries[2:end]
        push!(jump_directions, normal(boundary, 0))
    end

    println("# Boundaries: ", boundaries)
    println("# Jump directions: ", jump_directions)

    return boundaries, jump_directions
end

function Lorentz_gas(radius, boundaries, jump_directions, t_final)

    current_cell = Vector2D(0., 0.)

    #collision_times = -ones(5)

    x0, v0 = initial_condition(radius)
    p = Particle(x0, v0)

    total_time = 0.0

    xs = [x0[1]]
    ys = [x0[2]]

    which_hit = -1
    time_since_disc_collision = 0.0
    disc_collision_times = Float64[]

    output_time = 10.
    next_output_time = output_time
    num_data = int(round(t_final / output_time))

    println("num_data = ", num_data)

    step = 0

    positions = [x0]
    times = [0.0]

    #while total_time <= t_final
    while step <= num_data+1
        x_new, v_new, new_cell, collision_t, which_hit =
            #collision(p, current_cell, boundaries,
            #        jump_directions, which_hit)
        new_collision(p.x, p.v, radius, current_cell, which_hit)

        next_time = total_time + collision_t

        #println("next_time = ", next_time, "; next_output_time = ", next_output_time)

        #while (next_time > next_output_time)
        if (next_time >= next_output_time)
            # not necessary to use a while

            new_position = p.x + (next_output_time - total_time) * p.v
            new_position += new_cell

            #println(next_output_time, "\t", new_position[1], "\t", new_position[2])

            push!(positions, new_position)
            push!(times, next_output_time)

            next_output_time += output_time
            step += 1

            if step == num_data
                break
            end

        end

        total_time = next_time

        if total_time + collision_t < t_final
            total_time += collision_t
            current_cell = new_cell

            p.x, p.v = x_new, v_new

            #println(x_new[1], "\t", x_new[2], "\t", v_new[1], "\t", v_new[2])

            time_since_disc_collision += collision_t

            if which_hit == 1  # disc

                x_new += new_cell

                push!(xs, x_new[1])
                push!(ys, x_new[2])

                push!(disc_collision_times, time_since_disc_collision)
                time_since_disc_collision = 0.0
            end

        else
            collision_t = t_final - total_time

            x_new = p.x + p.v*collision_t
            x_new += current_cell

            println("# Position at time ", t_final, " = ", x_new)
            println("# Cell: ", current_cell)
            println("# Distance from origin = ", (dot(x_new, x_new))^0.5)
#            println("# Distance squared from origin = ", dot(x_new, x_new))

            break
        end

    end

    return xs, ys, disc_collision_times, times, positions
end

function many_particles(N, radius, boundaries, jump_directions, t_final)

    xs, ys, disc_collision_times, times, positions = Lorentz_gas(radius, boundaries, jump_directions, t_final)

 println("times = ", times)
    println(length(positions))

    all_positions = { { positions[i] } for i in 1:length(positions) }

    for n in 1:N-1
        xs, ys, disc_collision_times, times, positions = Lorentz_gas(radius, boundaries, jump_directions, t_final)

         println("times = ", times)
    println("all_positions = ", all_positions)
    println(length(positions))
    println(length(all_positions))

        for i in 1:length(positions)
            push!(all_positions[i], positions[i])
        end
    end

    all_positions
end

# function continuous_time_trajectory()

#     p = Particle(initial_condition()...)

#     t = 0.0;
#     next_output_time = output_time;

#     step = 1;  /// location in data array

#     while (t < max_time) {

#         next_t = t + min_collision_time;

#         while (next_t > next_output_time) {
#             /// use a while in case have very long path

#             new_position = p.position + (next_output_time - t)*p.velocity;

#             new_position += Vec(x_cell, y_cell);

#             displacement = new_position - p.initial_position;
#             disp_norm = dot(displacement, displacement);

#             /// |Delta x|
#             basic_moment = sqrt(sqrt(disp_norm)); /// |Delta x|^{1/2}
#             /// the basic unit is |Delta x|^{1/2}

#             moment = basic_moment;
#             /// moment is multiplied over and over by basic_moment to give the powers

#             for (int m=1; m<=num_moments; m++) {
#                 moments_data[m][step] += moment;
#                 moment *= basic_moment;
#             }

#             vel_vel_correlation_data[step] += dot(p.velocity, p.initial_velocity);
#             vel_pos_correlation_data[step] += dot(displacement, p.initial_velocity);
#             current_vel_pos_correlation_data[step] += dot(displacement, p.velocity);

#             step++;
#             next_output_time += output_time;

#             if (step > num_data) return;

#         }

#         implementCollision();

#         t = next_t;

#     //cout << "After collision:\n";
#     //cout << "Pos: ";

#         findCollision();

#     }

# }

function run_Lorentz_gas(radius, max_time)

    println("# Using radius: ", radius)

    boundaries, jump_directions = make_boundaries(radius)
    @time xs, ys = Lorentz_gas(radius, boundaries, jump_directions, max_time)
end

function run_many_particles(N, radius, max_time)

    println("# Using radius: ", radius)
    println("# Using $N particles")

    boundaries, jump_directions = make_boundaries(radius)
    @time all_positions = many_particles(N, radius, boundaries, jump_directions, max_time)
end

if length(ARGS) > 1
    radius = float(ARGS[1])
    run_Lorentz_gas(radius)
end

# using PyPlot
# plot(xs, ys)

end

function elastic_collision(obstacle, particle)
end

function periodic_boundary(obstacle, particle)  # needs *more* information and returns more information
end

# There's a fundaental difference betweeen "nromal"elastic objects and periodic boundaries
# There could also be other types of boundary, e.g. ones where crossing it randomizes the position of the disc in the new cell or just updates which cell the particle
# lieves in.

# E.g. think of a sequence of cells with pre-randomised (quenched) disc positions.

max_stars_repo_name: dpsanders/BilliardModel

max_issues_repo_name: dpsanders/BilliardModel

max_forks_repo_name: dpsanders/BilliardModel


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

JLLWrappers.@generate_main_file_header("SentencePiece")
JLLWrappers.@generate_main_file("SentencePiece", UUID("0cc836cc-86cb-5b0e-b290-caf6c31794d8"))
end  # module SentencePiece_jll


text:
import Base: bits, show

bits{N, ES, mode}(x::Sigmoid{N, ES, mode}) = bits(reinterpret(@UInt, x))[1:N]
function bits{N, ES, mode}(x::Sigmoid{N, ES, mode}, separator::AbstractString)
  #we're going to create this as a string array, then join() it at the end.
  stringarray = Vector{String}()

  bitstring = bits(x)

  push!(stringarray, bitstring[1:1])
  seek_idx = 2
  term_idx = N - (mode == :ubit) #store the index of termination.

  #calculate how many regime bits there are
  r_length = regimebits(x)

  push!(stringarray, bitstring[seek_idx:(seek_idx + r_length - 1)])
  seek_idx += r_length

  (seek_idx > term_idx) && @goto finish

  if (ES > 0)
    if ((seek_idx + ES) < term_idx)
      push!(stringarray, bitstring[seek_idx:seek_idx + ES - 1])
      seek_idx += ES
    else
      push!(stringarray, bitstring[seek_idx:term_idx])
      @goto finish
    end
  end

  push!(stringarray, bitstring[seek_idx:term_idx])

  @label finish

  if mode == :ubit
    push!(stringarray, bitstring[end:end])
  end

  return join(stringarray, separator)
end

function show{N, ES, mode}(io::IO, x::Sigmoid{N, ES, mode})

  if mode == :ubit
    #cast to Exact or ULP, depending on what the number looks like.
    if (reinterpret((@UInt), x) & increment(Sigmoid{N, ES, mode})) != 0
      print(io, "ULP{$N,$ES}")
    else
      print(io, "Exact{$N,$ES}")
    end
  else
    show(io, typeof(x))
  end

  innerval = reinterpret((@UInt), x)

  if N == 16
    print(io, "(0x", hex(innerval,16)[1:4], ")")
  elseif N == 32
    print(io, "(0x", hex(innerval,16)[1:8], ")")
  elseif N == 64
    print(io, "(0x", hex(innerval,16)[1:16], ")")
  else
    print(io, "(0b", bits(innerval)[1:N], ")")
  end
end

@generated function show{N, ES, mode}(io::IO, ::Type{Sigmoid{N, ES, mode}})
  if (mode == :guess)
    sig = "{$N,$ES}"
    :(print(io, "Posit", $sig))
  else
    sig = "{$N,$ES,$mode}"
    :(print(io, "Sigmoid", $sig))
  end
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
using SequentialMonteCarlo
import SMCExamples.LinearGaussian.defaultLGModel
import SMCExamples.MVLinearGaussian.defaultMVLGModel
import SMCExamples.SMCSampler.defaultSMCSampler
import SMCExamples.Lorenz96.defaultLorenzModel

using Dates
using BenchmarkTools

function benchMachineInfo()
  s::String = "$(Sys.MACHINE)\n" *
  "$(Sys.cpu_info()[1].model) ($(Sys.CPU_NAME)) ; " *
    "$(Sys.CPU_THREADS) Logical cores\n" *
    "Julia $VERSION using $(Threads.nthreads()) threads"
  return s
end

function dateandtime()
  datetimenow = now()
  s::String = Dates.format(datetimenow, "e, dd u yyyy HH:MM:SS")
  return s
end

function closestN(N::Int64, nt::Int64)
  mod(N, nt) == 0 && return N
  return round(Int64, N / nt) * nt
end

function bench(model, range)
  nt = Threads.nthreads()
  println("log₂N  Benchmark")
  for k in range
    print("$(lpad(string(k),2,string(0)))   ")
    N::Int64 = closestN(2^k, nt)
    smcio = SMCIO{model.particle, model.pScratch}(N, model.maxn, nt, false)
    @btime smc!($model, $smcio)
  end
end

println("Running Benchmarks: $(dateandtime())")
println(benchMachineInfo())

range = [10, 15, 20]

println("Linear Gaussian Model, n = 10")
model, theta, ys, ko = defaultLGModel(10)
bench(model, range)

println("Multivariate Linear Gaussian Model, d = 10, n = 10")
model, theta, ys, ko = defaultMVLGModel(10, 10)
bench(model, range)

println("SMC Sampler Example, n = 12")
model, ltarget = defaultSMCSampler()
bench(model, range)

println("Lorenz96 Example, d = 8, n = 10")
model, theta, ys = defaultLorenzModel(8, 10)
bench(model, range)

println("Finished: $(dateandtime())")


text:
struct QuasiNewtonSolver{T, FT, TJ, TL} <: AbstractNewtonSolver{T}
    @newton_solver_variables
    tx::Vector{T}
    ty::Vector{T}

    refactorize::Int

    function QuasiNewtonSolver{T,FT,TJ,TL}(x, Fparams, Jparams, linear_solver) where {T,FT,TJ,TL}
        J  = zero(linear_solver.A)
        x₀ = zero(x)
        x₁ = zero(x)
        y₁ = zero(x)
        y₀ = zero(x)
        δx = zero(x)
        δy = zero(x)
        tx = zero(x)
        ty = zero(x)

        nls_params = NonlinearSolverParameters(T)
        nls_status = NonlinearSolverStatus{T}(length(x))

        new(x, J, x₀, x₁, y₀, y₁, δx, δy, Fparams, Jparams, linear_solver,
            nls_params, nls_status, tx, ty, get_config(:quasi_newton_refactorize))
    end
end

const DEFAULT_LINESEARCH_nmax=50
const DEFAULT_ARMIJO_λ₀ = 1.0
const DEFAULT_ARMIJO_σ₀ = 0.1
const DEFAULT_ARMIJO_σ₁ = 0.5
const DEFAULT_ARMIJO_ϵ  = 0.5

function QuasiNewtonSolver(x::AbstractVector{T}, F!::Function; J!::Union{Function,Nothing}=nothing) where {T}
    n = length(x)
    Jparams = getJacobianParameters(J!, F!, T, n)
    linear_solver = getLinearSolver(x)
    QuasiNewtonSolver{T, typeof(F!), typeof(Jparams), typeof(linear_solver)}(x, F!, Jparams, linear_solver)
end

function solve!(s::QuasiNewtonSolver{T}; n::Int=0) where {T}
    local λ::T
    local λₜ::T
    local y₀norm::T
    local y₁norm::T
    local lsiter::Int
    local p₀::T
    local p₁::T
    local p₂::T
    local nmax::Int = n > 0 ? nmax = n : s.params.nmax
    local refactorize::Int = s.refactorize

    s.F!(s.x, s.y₀)
    residual_initial!(s.status, s.x, s.y₀)
    s.status.i  = 0

    if s.status.rₐ ≥ s.params.atol² || n > 0 || s.params.nmin > 0
        computeJacobian(s)
        s.linear.A .= s.J
        factorize!(s.linear)

        for s.status.i = 1:nmax
            s.x₀ .= s.x
            y₀norm = l2norm(s.y₀)

            # b = - y₀
            s.linear.b .= -one(T) .* s.y₀

            # solve J δx = -f(x)
            solve!(s.linear)

            # TODO Separate line search algorithms into independent modules.

            # δx = b
            s.δx .= s.linear.b

            # δy = Jδx
            simd_mult!(s.δy, s.J, s.δx)

            # set λ to default initial value
            λ = DEFAULT_ARMIJO_λ₀

            # use simple line search to determine a λ for which there is not Domain Error
            for lsiter in 1:DEFAULT_LINESEARCH_nmax
                # x₁ = x₀ + λ δx
                s.x₁ .=  s.x₀ .+ λ .* s.δx

                try
                    # y₁ = f(x₁)
                    s.F!(s.x₁, s.y₁)

                    break
                catch DomainError
                    # in case the new function value results in some DomainError
                    # (e.g., for functions f(x) containing sqrt's or log's),
                    # decrease λ and retry

                    println("WARNING: Quasi-Newton Solver encountered Domain Error (lsiter=", lsiter, ", λ=", λ, ").")

                    λ *= DEFAULT_ARMIJO_σ₁
                end
            end

            # x₁ = x₀ + λ δx
            s.x₁ .= s.x₀ .+ λ .* s.δx

            # y₁ = f(x₁)
            s.F!(s.x₁, s.y₁)

            # compute norms of solutions
            y₀norm = l2norm(s.y₀)
            y₁norm = l2norm(s.y₁)

            if y₁norm < (one(T)-DEFAULT_ARMIJO_σ₀*λ)*y₀norm
                nothing
            else
                # determine coefficients of polynomial p(λ) = p₀ + p₁λ + p₂λ²
                p₀ = y₀norm^2
                p₁ = 2(⋅(s.y₀, s.δy))
                p₂ = (y₁norm^2 - y₀norm^2 - p₁*λ)/(λ^2)

                # compute minimum λₜ of p(λ)
                λₜ = - p₁/(2p₂)

                if λₜ < DEFAULT_ARMIJO_σ₀ * λ
                    λ = DEFAULT_ARMIJO_σ₀ * λ
                elseif λₜ > DEFAULT_ARMIJO_σ₁ * λ
                    λ = DEFAULT_ARMIJO_σ₁ * λ
                else
                    λ = λₜ
                end
            end

            s.δx .*= λ

            # simple Armijo line search
            #
            # λ = DEFAULT_ARMIJO_λ₀
            #
            # # δx = λ b
            # simd_copy_scale!(λ, s.linear.b, s.δx)
            #
            # for lsiter in 1:DEFAULT_LINESEARCH_nmax
            #     # x₁ = x₀ + λ δx
            #     simd_wxpy!(s.x₁, s.δx, s.x₀)
            #
            #     try
            #         # y₁ = f(x₁)
            #         function_stages!(s.x₁, s.y₁, s.Fparams)
            #
            #         # y₁norm = ||y₁||₂
            #         y₁norm = vecnorm(s.y₁, 2)
            #
            #         # exit loop if ||f(x₁)||₂ < (1-λϵ) ||f(x₀)||₂
            #         if y₁norm < (one(T)-λ*DEFAULT_ARMIJO_ϵ)*y₀norm
            #             break
            #         end
            #     catch DomainError
            #         # in case the new function value results in some DomainError
            #         # (e.g., for functions f(x) containing sqrt's or log's),
            #         # decrease λ and retry
            #         println("WARNING: Quasi-Newton Solver encountered Domain Error.")
            #         # if lsiter == DEFAULT_LINESEARCH_nmax
            #         #     λ *= 0
            #         #     simd_scale!(s.δx, 0)
            #         # end
            #     end
            #
            #     λ *= DEFAULT_ARMIJO_σ₁
            #     simd_scale!(s.δx, DEFAULT_ARMIJO_σ₁)
            # end

            s.x .+= s.δx
            s.F!(s.x, s.y₀)
            residual!(s.status, s.δx, s.x, s.y₀)

            if check_solver_converged(s.status, s.params) && s.status.i ≥ s.params.nmin && !(n > 0)
                if s.params.nwarn > 0 && s.status.i ≥ s.params.nwarn
                    println("WARNING: Quasi-Newton Solver took ", s.status.i, " iterations.")
                end
                break
            end

            if mod(s.status.i, refactorize) == 0
                computeJacobian(s)
                s.linear.A .= s.J
                factorize!(s.linear)
            end
        end
    end
end


text:
export add_update!

tr2blas(X::BlasArray) = 'N', X
tr2blas(X::BlasTranspose) = 'T', parent(X)
tr2blas(X::BlasAdjoint) = 'C', parent(X)
tr2blas(X::BlasAdjoint{T}) where T <: BlasReal = 'T', parent(X)

#  GEMM and SYMM =======================

do_gemm!(α, TA, A, TB, B, β, C) = BLAS.gemm!(TA, TB, α, A, B, β, C)      # C ← αAB + βC  (with transpositions as specified)
do_gemm!(α, TA, A::Symmetric, TB, B, β, C) = BLAS.symm!('L', A.uplo, α, blasnode(A), B, β, C)
do_gemm!(α, TA, A, TB, B::Symmetric, β, C) = BLAS.symm!('R', B.uplo, α, blasnode(B), A, β, C)

function gemm_αABβC!(α, A, B, β, C::BlasMatrix{T}) where T
    (C===blasnode(A) || C===blasnode(B)) && ip_error("multiplicative and additive updates cannot be combined.")
    try
        α, β = convert.(T, (α, β))
        do_gemm!(α, tr2blas(A)..., tr2blas(B)..., β, C)
    catch err
       ip_error(err)
   end
end

add_update0!(C::BlasMatrix{T}, β::Number,  α::Number, A::BlasMatrix{T}, B::BlasMatrix{T}) where T =
    gemm_αABβC!(α, A, B, β, C)
#

# AXPY and AXPBY and copyto!  =====================================================================
for (pm,sα) in ( (:+, :α), (:-, :(-α)) )
    @eval begin
        function add_update!(y::AxpyVec{T}, β::Number, ::typeof($pm), α::Number, x::AxpyVec{T}) where T
            try
                β==0 && α==1 && return copyto!(y,x)
                α, β = convert.(T, ($sα, β))
                β==1 ? BLAS.axpy!(α,x,y) : BLAS.axpby!(α,x,β,y)
            catch err
                ip_error(err)
            end
        end
    end
end
# add_update!(y::AxpyVec{T}, β::Number, pm::Function, x::AxpyVec{T}, α::Number) where T =
#     add_update!(y::AxpyVec{T}, β::Number, pm::Function, α::Number, x::AxpyVec{T})
# #
add_update!(y::AxpyVec{T}, β::Number, pm::Function, x::AxpyVec{T}) where T =
    add_update!(y::AxpyVec{T}, β::Number, pm::Function, 1, x::AxpyVec{T})
#

# Updating triangular and Symmetric matrices  =======================
function add_update!(y::BlasTriangular{T}, β::Number, ::typeof(+), α::Number, x::BlasTriangular{T}) where T
    uplo_eff(y) == uplo_eff(x) || ip_error("uplo must agree for additive triangular updates.")
    try
        if β==0
           y .= α .* x
        elseif β==1
            y .+= α .* x
        else
            y .= β .* y .+ α .* x
        end
    catch err
        ip_error(err)
    end
    return y
end

function add_update!(y::BlasTriangular{T}, β::Number, ::typeof(+), x::BlasTriangular{T}) where T
    uplo_eff(y) == uplo_eff(x) || ip_error("uplo must agree for additive triangular updates.")
    try
        if β==0
            y .= x
        elseif β==1
            y .+= x
        else
            y .= β .* y .+ x
        end
    catch err
        ip_error(err)
    end
    return y
end

function add_update!(y::BlasTriangular{T}, β::Number, ::typeof(-), α::Number, x::BlasTriangular{T}) where T
    uplo_eff(y) == uplo_eff(x) || ip_error("uplo must agree for additive triangular updates.")
    try
        if β==0
           y .= α .* x
        elseif β==1
            y .-= α .* x
        else
            y .= β .* y .- α .* x
        end
    catch err
        ip_error(err)
    end
    return y
end

function add_update!(y::BlasTriangular{T}, β::Number, ::typeof(-), x::BlasTriangular{T}) where T
    uplo_eff(y) == uplo_eff(x) || ip_error("uplo must agree for additive triangular updates.")
    try
        if β==0
            y .= x
        elseif β==1
            y .-= x
        else
            y .= β .* y .- x
        end
    catch err
        ip_error(err)
    end
    return y
end

function add_update!(y::Symmetric{T}, β::Number, pm::Function, α::Number, x::Symmetric{T}) where T <: BlasFloat
    if y.uplo=='L' && x.uplo=='L'
        add_update!(LowerTriangular(parent(y)), β, pm, α, LowerTriangular(parent(x)))
    elseif y.uplo=='U' && x.uplo=='U'
        add_update!(UpperTriangular(parent(y)), β, pm, α, UpperTriangular(parent(x)))
    else
        ip_error("uplo must agree for additive Symmetric updates.")
    end
    return y
end

function add_update!(y::Symmetric{T}, β::Number, pm::Function, x::Symmetric{T}) where T <: BlasFloat
    if y.uplo=='L' && x.uplo=='L'
        add_update!(LowerTriangular(parent(y)), β, pm, LowerTriangular(parent(x)))
    elseif y.uplo=='U' && x.uplo=='U'
        add_update!(UpperTriangular(parent(y)), β, pm, UpperTriangular(parent(x)))
    else
        ip_error("uplo must agree for additive Symmetric updates.")
    end
    return y
end

# GEMV and SYMV =====================================================================================
import LinearAlgebra: lmul!,rmul!

do_gemv!(α, TA, A, B, β, C) = BLAS.gemv!(TA, α, A, B, β, C)      # C ← αAB + βC  (with transposition as specified)
do_gemv!(α, TA, A::Symmetric, B, β, C) = BLAS.symv!(A.uplo, α, blasnode(A), B, β, C)

function gemv_αABβC!(α, A, B, β, C::BlasVector{T}) where T
    (C===blasnode(A) || C===blasnode(B)) && ip_error("multiplicative and additive updates cannot be combined.")
    try
        α, β = convert.(T, (α, β))
        do_gemv!(α, tr2blas(A)..., B, β, C)
    catch err
        ip_error(err)
    end
end

add_update0!(C::BlasVector{T}, β::Number, α::Number, A::BlasMatrix{T}, B::BlasVector{T}) where T =
    gemv_αABβC!(α, A, B, β, C)
#

# SYRK =============================================
function do_syrk!(α, A, B, β, C::Symmetric{T}) where T
    TA, A = tr2blas(A)
    TB, B = tr2blas(B)
    cond1 = A===B
    cond2 = (TA=='N' && TB=='T') || (TA=='T' && TB=='N')
    ( cond1 && cond2) || ip_error("conditions violated for update to Symmetric LHS.")
    try
        α, β = convert.(T, (α, β))
        BLAS.syrk!(C.uplo, TA, α, A, β, blasnode(C))
        return C  # syrk! returns parent(C), which is usually not symmetric
    catch err
        ip_error(err)
    end
end

add_update0!(C::Symmetric{T}, β::Number, α::Number, A::BlasMatrix{T}, B::BlasMatrix{T}) where T =
    do_syrk!(α, A, B, β, C)
#

# HERK =============================================
function do_herk!(α, A, B, β, C::Hermitian{T}) where T  <: Complex{F} where F
    TA, A = tr2blas(A)
    TB, B = tr2blas(B)
    cond1 = A===B
    cond2 = (TA=='N' && TB=='C') || (TA=='C' && TB=='N')
    ( cond1 && cond2) || ip_error("conditions violated for update to Hermitian LHS.")
    try
        α, β = convert.(F, (α, β))
        BLAS.herk!(C.uplo, TA, α, A, β, blasnode(C))
        return C  # herk! returns parent(C)
    catch err
        ip_error(err)
    end
end

add_update0!(C::Hermitian{T}, β::Number, α::Number, A::BlasMatrix{T}, B::BlasMatrix{T}) where T <: Complex =
    do_herk!(α, A, B, β, C)
#

# SYR =============================================
function do_syr!(α, A, B, β, C::Symmetric{T}) where T
    β==1 || ip_error("symmetric rank 1 update cannot prescale LHS (β==1 required).")
    TB, B = tr2blas(B)
    (A===B && TB=='T') || ip_error("conditions violated for update to Symmetric LHS.")
    try
        α = convert(T, α)
        BLAS.syr!(C.uplo, α, A, blasnode(C))
        return C  # syr! returns parent(C), which is usually not symmetric
    catch err
        ip_error(err)
    end
end

add_update0!(C::Symmetric{T}, β::Number, α::Number, A::BlasVector{T}, B::BlasRow{T}) where T =
    do_syr!(α, A, B, β, C)
#

# HER =============================================
function do_her!(α, A, B, β, C::Hermitian{T}) where T  <: Complex{F} where F
    β==1 || ip_error("Hermitian rank 1 update cannot prescale LHS (β==1 required).")
    A===parent(B) || ip_error("conditions violated for update to Hermitian LHS.")
    try
        α = convert(F, α)
        BLAS.her!(C.uplo, α, A, blasnode(C))
        return C  # her! returns parent(C)
    catch err
        ip_error(err)
    end
end

add_update0!(C::Hermitian{T}, β::Number, α::Number, A::BlasVector{T}, B::BlasAdjRow{T}) where T <: Complex =
    do_her!(α, A, B, β, C)
#

function add_update!(C::BlasMatrix{T}, β::Number, ::typeof(+), α::Number, U::UniformScaling) where T
    try
        αλ,β = convert.(T, (α*U.λ,β))
        n = checksquare(C)
        if β==0
            for i=1:n
                @inbounds C[i,i] = αλ
            end
        elseif β==1
            for i=1:n
                @inbounds C[i,i] += αλ
            end
        else
            for i=1:n
                @inbounds C[i,i] = β*C[i,i] + αλ
            end
        end
        return C
    catch err
        ip_error(err)
    end
end
add_update!(C::BlasMatrix, β::Number, ::typeof(-), α::Number, U::UniformScaling) = add_update!(C, β, +, -α, U)
add_update!(C::BlasMatrix, β::Number, ::typeof(+), U::UniformScaling) = add_update!(C, β, +, 1, U)
add_update!(C::BlasMatrix, β::Number, ::typeof(-), U::UniformScaling) = add_update!(C, β, +, -1, U)

add_update!(C, β, pm::Function, α, A, B) = ip_error("inplace assignment not available for this combination of types.")
add_update0!(C, β::Number, α::Number, A, B) = ip_error("inplace assignment not available for this combination of types.")
for (pm, sα, s) in ( (:+, :α, 1), (:-, :(-α), -1) )
    @eval begin
        add_update!(C, β::Number, ::typeof($pm), α::Number, A, B) =
            add_update0!(C, β, $sα, A, B)
        #
        # add_update!(C, β::Number, ::typeof($pm), A, α::Number, B) =
        #     add_update0!(C, β, $sα, A, B)
        # #
        # add_update!(C, β::Number, ::typeof($pm), A, B, α::Number) =
        #     add_update0!(C, β, $sα, A, B)
        # #

        add_update!(C::BlasMatrix{T}, β::Number, ::typeof($pm), A::BlasArray{T}, B::BlasArray{T}) where T =
            add_update0!(C, β, $s, A, B)
        #

    end
end

#Do we want the following error behaviour, or do we just go ahead and do the broadcast?
add_update!(C, ::Number, ::Function, ::Number) = ip_error("additive update with scalar RHS not available. Use broadcasting instead.")


text:
module RandomCircuits

using QuantumOpticsBase, QuantumInformation
using LinearAlgebra, SparseArrays, ElasticArrays, Kronecker
using OrdinaryDiffEq, DiffEqBase, DiffEqCallbacks

include("types.jl")
export AbstractGate, HaarGate, LazyHaarGate, hamiltonian, RowGate, LazyRowGate, generate_random_circuit, generate_lazy_random_circuit

include("timeevolution.jl")
export compute_tstops, time_evolution

end