text:
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    ident(p)

Description:
Create dense identity matrix of dimension p × p.

Inputs:
- p = number of columns and rows of identity matrix.

Outputs:
- ident_mat = p × p dense identity matrix with Float64.
"""
function ident(p)
    ident_mat = 1.0 * Matrix(I, p, p)
    return ident_mat
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    transp_dbl(A)

Description:
Transpose matrix, so that the output type is Array{Float64,2}.

Inputs:
- A = matrix.

Outputs:
- transp_A = dense transposed matrix A.
"""
function transp_dbl(A)
    transp_A = 1.0 * Matrix(A')
    return transp_A
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    transp_int(A)

Description:
Transpose matrix, so that the output type is Array{Int64,2}.

Inputs:
- A = matrix.

Outputs:
- transp_A = dense transposed matrix A.
"""
function transp_int(A)
    transp_A = 1 * Matrix(A')
    return transp_A
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    seqa(a, b, c)

Description:
Create a sequence with `c` # of inputs, starting at `a` with `b` increments.

Inputs:
- a = starting value of the sequence.
- b = increment size.
- c = number of increments.

Outputs:
- seq = sequence with `c` # of inputs, starting at `a` with `b` increments.

"""
function seqa(a, b, c)

    seq = zeros(c)

    for i = 1:c
        seq[i] = a + b * (i - 1)
    end

    return seq
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    sigbig(phi, p, capt)

Description:
Create T×T matrix S^(-1) (pg. 1003 of Otrok and Whiteman (1998)).

Inputs:
- phi   = Lag coefficients associated with idiosyncratic error autoregression.
- p     = Number of lags.
- capt  = total number of time periods in the sample.

Outputs:
- Si    = T×T matrix S^(-1).
"""
function sigbig(phi, p, capt)

    rotate = trunc.(Int, seqa(0, 1, capt - p))
    Siinv_upper = Hermitian(sigmat(phi, p))
    Siinv_upper = [inv(cholesky(Siinv_upper, Val(true), check = false).U)' zeros(p, capt - p)]
    Siinv_lower = [kron((-reverse(phi, dims = 1)'), ones(capt - p, 1)) ones(capt - p, 1) zeros(capt - p, capt - p - 1)]

    for s = 1:(capt-p)
        Siinv_lower[s, :] = circshift(Siinv_lower, (0, rotate[s]))[s, :]
    end

    Si = [Siinv_upper; Siinv_lower]

    return Si
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    rnchisq(m)

Description:
Draw sample from a chi-square distribution.

Inputs:
- m = degrees of freedom.

Outputs:
- g = draw from a chi-square distribution with m degrees of freedom.
"""
function rnchisq(m)

    g = rand(Chisq(m))

    return g
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    sigmat(phi, p)

Description:
Yield the stationary covariance matrix of the first p errors (pg 1001 of Otrok and Whiteman (1998)).

Inputs:
- phi  = Lag coefficients associated with idiosyncratic error autoregression.
- p    = Number of lags.

Outputs:
- Stationary covariance matrix of the first p errors.

"""
function sigmat(phi, p)

    r2 = p^2
    i = [ident(p - 1) zeros(p - 1, 1)]
    Pcap = [phi'; i]
    pp = ident(r2) - kron(Pcap, Pcap)
    e1 = [1; zeros(p - 1, 1)]
    sig = inv(pp) * vec(e1 * e1')
    sigma = reshape(sig', p, p)
    sigma = Matrix(sigma)

    return sigma
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    gendiff(z, phi)

Description:
Difference a series `z` using the lag coefficients provided in `phi`.

Inputs:
- z = time series.
- phi = lag coefficient vector.

Outputs:
- zgdiff = quasi-differenced version of the inputted series `z`.

"""
function gendiff(z, phi)

    phi = vec(phi)

    p = size(phi)[1]
    ztrim = z[(p+1):end, :]
    zgdiff = ztrim

    for i = 1:p
        zgdiff = zgdiff - phi[i, 1] * z[(p-i+1):(end-i), :]
    end

    return zgdiff
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
"""
function arfac(y, p, r0_, R0__, phi0, sig2, capt)

    # Generation of phi
    yp = y[1:p, 1]          # the first p observations
    e = y
    e1 = e[p+1:capt, 1]
    ecap = zeros(capt, p)

    for j = 1:p
        ecap[:, j] = lag(e, j, default = 0.0)
    end

    ecap = ecap[p+1:capt, :]

    V = invpd(R0__ + sig2^(-1) * ecap' * ecap)
    phihat = V * (R0__ * r0_ + sig2^(-1) * ecap' * e1)

    phi1 = sim_MvNormal(phihat, V)

    coef = [-reverse(phi1, dims = 1); 1]                # check stationarity
    root = roots(Polynomial(reverse(vec(coef), dims = 1)))                      # Find lag polynomial roots

    rootmod = abs.(root)
    accept = min(rootmod...) >= 1.01                     # all the roots bigger than 1

    if accept == 0                                      # doesn't pass stationarity
        phi1 = phi0
    else

        #=
        sigma1 = Hermitian(sigmat(vec(phi1), p))               # numerator of acceptance prob
        d = det(sigma1)
        psi1 = (d^(-1/2)) * exp((-0.5 / sig2) * (transp_dbl(yp)*invpd(sigma1)*(yp))[1])

        sigma1 = Hermitian(sigmat(vec(phi0), p))               # numerator of acceptance prob
        d = det(sigma1)
        psi0 = (d^(-1/2)) * exp((-0.5 / sig2) * (transp_dbl(yp)*invpd(sigma1)*(yp))[1])
        =#
        ##=
        sigma1 = sigmat(phi1, p)       # numerator of acceptance prob
        sigma1 = Hermitian(sigma1)
        sigroot = cholesky(sigma1, Val(true)).U
        p1 = transp_dbl(inv(sigroot))
        ypst = p1 * yp
        d = det(p1' * p1)
        psi1 = (d^(1 / 2)) * exp(-0.5 * (ypst)' * (ypst) / sig2)

        sigma1 = sigmat(phi0, p)       # numerator of acceptance prob
        sigma1 = Hermitian(sigma1)
        sigroot = cholesky(sigma1, Val(true)).U
        p1 = transp_dbl(inv(sigroot))
        ypst = p1 * yp
        d = det(p1' * p1)
        psi0 = (d^(1 / 2)) * exp(-0.5 * (ypst)' * (ypst) / sig2)
        ##=#

        if psi0 == 0
            accept = 1
        else
            u = rand(1)
            accept = u[1] <= psi1 / psi0
        end
        phi1 = phi1 * accept + phi0 * (1 - accept)
    end

    return phi1
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    seqm(a,b,c)

Description:
Prouce a sequence of values.

Inputs:
- a = initial value of sequence
- b = increment scaling size
- c = number of values in the sequence

Outputs:
- seq = geometric sequence
"""
function seqm(a, b, c)

    seq = zeros(c, 1)

    seq[1] = a

    if c > 1

        seq[2] = a * b

        for i = 3:c
            seq[i] = seq[i-1] * b
        end
    end

    return seq
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    invpd(X)

Description:
Invert matrix X using augmented eigenvalues if X is not positive-definite.

Inputs:
- X = matrix of interest

Output:
- X_inv = inverse of input matrix X
"""
function invpd(X)

    X_inv = similar(X)

    if isposdef(X) == true
        X_inv = inv(X)
    else
        n, m = size(X)
        U, dd, V = svd(X)
        dd = vec(dd)
        xchk = U * Diagonal(dd) * V'
        dd = dd .+ 1000 * 2.2204e-16
        di = ones(n, 1) ./ dd
        di = vec(di)
        X_inv = U * Diagonal(di) * V
    end

    return X_inv
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    diagrv(X,V)

Description:
Replace the diagonal of matrix X with vector V.

Inputs:
- X = matrix of interest
- V = new diagonal of X

Outputs:
- X_new = X with new diagonal V
"""
function diagrv(X, V)

    X_new = X
    X_new[diagind(X_new)] = V

    X_new
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
    ar(y, x, p, b0_, B0__, r0_, R0__, v0_, d0_, b0, s20, phi0, xvar, nfc, facts, capt, nreg, Size)

Description:
Generate β_i , ϕ_i, and σ^2_i estimates based on the posterior distributions (eqs. 6-8) on pg. 1002 of Otrok and Whiteman (1998).

Inputs:
- y = observed dependent variable / data.
- x = measurement equation regressors (intercept + factor).
- p = number of state eq. lags.
- b0_ = prior mean of factor loading (factor coefficient in measurement equation).
- B0__ = prior variance of factor loading (factor coefficient in measurement equation).
- r0_ = prior mean of phi (idiosyncratic AR polynomial).
- R0__ = prior precision of phi (idiosyncratic AR polynomial).
- v0_ = inverted gamma parameter of innovation variances.
- d0_ = inverted gamma parameter of innovation variances.
- b0 = old draw of factor loading.
- s20 = old draw of innovation variances.
- phi0 = old draw of phi (idiosyncratic AR polynomial).
- xvar = measurement eq. regressor index.
- nfc = number of factors.
- facts = factor draw.
- capt = total number of time periods / observations in the sample.
- nreg = number of regressors in each observable equation, constant plus K factors
- Size = number of variables each factor loads on

Outputs:
- b0 = β coefficient hyperparameter estimates.
- s20 = σ^2 coefficient hyperparameter estimate.
- phi1 = ϕ coefficient hyperparameter estimates.
- facts = factor estimate.
"""
function ar(y, x, p, b0_, B0__, r0_, R0__, v0_, d0_, b0, s20, phi0, xvar, nfc, facts, capt, nreg, Size)

    local xst, yst, b1, phi1

    n = capt
    signmax1 = 1
    signbeta1 = 1

    while signbeta1 >= 1

        # generation of phi1
        yp = y[1:p, 1]          # the first p observations
        xp = x[1:p, :]

        e = y - x * b0
        e1 = e[(p+1):n, 1]
        ecap = zeros(n, p)

        for j = 1:p
            ecap[:, j] = lag(e, j, default = 0.0)
        end

        ecap = ecap[(p+1):n, :]

        V = invpd(R0__ .+ inv(s20) * (transp_dbl(ecap) * ecap))
        phihat = V * (R0__ * r0_ + inv(s20) * (transp_dbl(ecap) * e1))

        phi1 = sim_MvNormal(vec(phihat), V)

        coef = [-reverse(vec(phi1), dims = 1); 1]                      # check stationarity
        root = roots(Polynomial(reverse(coef)))
        rootmod = abs.(root)
        accept = min(rootmod...) >= 1.01             # all the roots bigger than 1

        if accept == 0
            phi1 = phi0
        else
            sigma1 = Hermitian(sigmat(vec(phi1), p))               # numerator of acceptance prob
            d = det(sigma1)
            psi1 = (d^(0.5)) * exp((-0.5 / s20) * (transp_dbl(yp - xp * b0)*invpd(sigma1)*(yp-xp*b0))[1])

            sigma1 = Hermitian(sigmat(vec(phi0), p))               # numerator of acceptance prob
            d = det(sigma1)
            psi0 = (d^(1 / 2)) * exp((-0.5 / s20) * (transp_dbl(yp - xp * b0)*invpd(sigma1)*(yp-xp*b0))[1])

            if psi0 == 0
                accept = 1
            else
                u = rand(Uniform(0, 1))
                accept = u <= psi1 / psi0
            end
            phi1 = phi1 * accept + phi0 * (1 - accept)
        end

        # generation of beta
        sigma = Hermitian(sigmat(phi1, p))             # sigma = sigroot' * sigroot
        sigroot = cholesky(sigma, Val(true)).U                 # signber2v = p1' * p1
        p1 = transp_dbl(inv(sigroot))
        ypst = p1 * yp
        xpst = p1 * xp
        yst = [ypst; gendiff(y, phi1)]
        xst = [xpst; gendiff(x, phi1)]

        V = invpd(B0__ + s20^(-1) * xst' * xst)
        bhat = V * (B0__ * b0 + s20^(-1) * xst' * yst)
        b1 = sim_MvNormal(vec(bhat), Matrix(V))

        signbeta1 = (b1[2, 1] <= 0.0) * (xvar == 1)
        signmax1 = signmax1 + (1 * signbeta1)

        if signmax1 >= 100
            facts[:, 1] = (-1) * facts[:, 1]
            x[:, 2] = facts[:, 1]
            signmax1 = 1
        end
    end

    # generation of s2
    nn = n + v0_
    d = d0_ + transp_dbl(yst - xst * b1) * (yst - xst * b1)
    c = rnchisq(nn)
    t2 = c ./ d
    s21 = 1 ./ t2
    b0 = b1
    s20 = s21[1, 1]

    return b0, s20, phi1, facts
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
@doc """
"""
function ar_LJ(y, x, p, b0_, B0__, r0_, R0__, v0_, d0_, b0, s20, phi0, xvar, nfc, facts, capt, nreg, Size, varassign)

    # Make sure relevant objects are accessible
    # outside of the upcoming for-loops
    local xst, yst, b1, phi1

    # Save number of observations in the sample
    n = capt

    # Initialize parameter signs
    signmax1 = 1.0      #
    signbeta1 = 1.0     #
    signmax2 = 1.0      #
    signbeta2 = 0.0     #

    testind = 0

    while signbeta1 + signbeta2 >= 1.0 ##&& testind < 10000

        testind = testind + 1

        # generation of phi1
        yp = y[1:p, 1]          # the first p observations
        xp = x[1:p, :]

        e = y - x * vec(b0)
        e1 = e[(p+1):n, 1]
        ecap = zeros(n, p)

        for j = 1:p
            ecap[:, j] = lag(e, j, default = 0.0)
        end

        ecap = ecap[(p+1):n, :]

        V = invpd(R0__ .+ inv(s20) * (transp_dbl(ecap) * ecap))
        phihat = V * (R0__ * r0_ + inv(s20) * (transp_dbl(ecap) * e1))

        phi1 = sim_MvNormal(vec(phihat), Matrix(V))

        coef = [-reverse(vec(phi1), dims = 1); 1]                      # check stationarity
        root = roots(Polynomial(reverse(coef)))
        rootmod = abs.(root)
        accept = min(rootmod...) >= 1.01             # all the roots bigger than 1

        if accept == 0
            phi1 = phi0
        else

            #=
            sigma1 = Hermitian(sigmat(vec(phi1), p))               # numerator of acceptance prob
            d = det(sigma1)
            psi1 = (d^(-1/2)) * exp((-0.5 / s20) * (transp_dbl(yp - xp * b0')*invpd(sigma1)*(yp-xp*b0'))[1])

            sigma1 = Hermitian(sigmat(vec(phi0), p))               # numerator of acceptance prob
            d = det(sigma1)
            psi0 = (d^(-1/2)) * exp((-0.5 / s20) * (transp_dbl(yp - xp * b0')*invpd(sigma1)*(yp-xp*b0'))[1])
            =#

            sigma1 = sigmat(phi1, p)       # numerator of acceptance prob
            sigma1 = Hermitian(sigma1)
            sigroot = cholesky(sigma1, Val(true)).U
            p1 = transp_dbl(inv(sigroot))
            ypst = p1 * yp
            xpst = p1 * xp
            d = det(p1' * p1)
            psi1 = (d^(1 / 2)) * exp(-0.5 * (ypst - xpst * b0')' * (ypst - xpst * b0') / s20)

            sigma1 = sigmat(phi0, p)       # numerator of acceptance prob
            sigma1 = Hermitian(sigma1)
            sigroot = cholesky(sigma1, Val(true)).U
            p1 = transp_dbl(inv(sigroot))
            ypst = p1 * yp
            xpst = p1 * xp
            d = det(p1' * p1)
            psi0 = (d^(1 / 2)) * exp(-0.5 * (ypst - xpst * b0')' * (ypst - xpst * b0') / s20)

            if psi0 == 0
                accept = 1
            else
                u = rand(Uniform(0, 1))
                accept = u <= psi1[1, 1] / psi0[1, 1]
            end
            phi1 = phi1 * accept + phi0 * (1 - accept)
        end

        # generation of beta
        sigma = Hermitian(sigmat(phi1, p))              # sigma = sigroot' * sigroot
        sigroot = cholesky(sigma, Val(true)).U                       # signber2v = p1' * p1
        p1 = transp_dbl(inv(sigroot))
        ypst = p1 * yp
        xpst = p1 * xp
        yst = [ypst; gendiff(y, phi1)]
        xst = [xpst; gendiff(x, phi1)]

        V = invpd(B0__ + s20^(-1) * xst' * xst)
        bhat = V * (B0__ * vec(b0) + s20^(-1) * xst' * yst)

        b1 = sim_MvNormal(vec(bhat), Matrix(V))

        signbeta1 = (b1[2, 1] <= 0.0) * (xvar == 1)
        signmax1 = signmax1 + (1 * signbeta1)
        signbeta2 = (b1[3, 1] <= 0.0) * (varassign[nfc][1] == xvar)
        signmax2 = signmax2 + (1 * signbeta2)

        if signmax1 >= 100
            facts[:, 1] = (-1) * facts[:, 1]
            x[:, 2] = facts[:, 1]
            signmax1 = 1
        end
        if signmax2 >= 100
            facts[:, 1+nfc] = (-1) * facts[:, 1+nfc]
            x[:, 3] = facts[:, 1+nfc]
            signmax2 = 1
        end
    end

    # generation of s2
    nn = n + v0_
    d = d0_ + (transp_dbl(yst - xst * b1)*(yst-xst*b1))[1, 1]
    c = rnchisq(nn)
    t2 = c ./ d
    s21 = 1 ./ t2
    b0 = b1
    s20 = s21[1, 1]

    return b0, s20, phi1, facts
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################


text:
using ApproxFunBase, Test
    import ApproxFunBase: Block

@testset "RaggedMatrix" begin
    cols=Int[rand(1:k+2) for k=1:5]
    B=ApproxFunBase.rrand(Float64,maximum(cols),cols)
    cols=Int[rand(1:k+2) for k=1:size(B,1)]
    A=ApproxFunBase.rrand(Float64,maximum(cols),cols)
    @test Matrix(A)*Matrix(B) ≈ Matrix(A*B)

    @test ApproxFunBase.RaggedMatrix(B) === B
    @test ApproxFunBase.RaggedMatrix{Float64}(B) === B
    @test Matrix(ApproxFunBase.RaggedMatrix{ComplexF64}(B)) == Matrix{ComplexF64}(Matrix(B))

    B = ApproxFunBase.brand(10,10,2,3)
    @test Matrix(B) == Matrix(ApproxFunBase.RaggedMatrix(B))
    @test ApproxFunBase.RaggedMatrix(B) == ApproxFunBase.RaggedMatrix{Float64}(B)
    @test ApproxFunBase.RaggedMatrix(ApproxFunBase.BandedMatrix{ComplexF64}(B)) == ApproxFunBase.RaggedMatrix{ComplexF64}(B)
end

ext: jl


text:
using JLD2
using PyPlot
using SegyIO

run(`wget -r ftp://slim.gatech.edu/data//synth/Compass/final_velocity_model_ieee_6m.sgy`) # this might take a while

# get velocity
block = segy_read("slim.gatech.edu/data/synth/Compass/final_velocity_model_ieee_6m.sgy")

# original compass model is in 25m*25m*6m
n1 = (1911,2730,341)
d = (25f0,25f0,6f0)

n = (637,910,341)

sx = get_header(block, "SourceX")
sy = get_header(block, "SourceY")

v_nogas = zeros(Float32,n)

for i = 1:n[1]
    x = d[1].*(i-1)
    inds = findall(sx.==x)
    slice = block.data[:,inds[sortperm(sy[inds])]]

    v_nogas[i,:,:] = transpose(slice[:,1:Int(end/3)])
end

v = v_nogas[1:256,1:256,end-255:end]
JLD2.@save "256cube_v.jld2" v

max_stars_repo_name: slimgroup/TwoPhaseSpeedTest

max_issues_repo_name: slimgroup/TwoPhaseSpeedTest

max_forks_repo_name: slimgroup/TwoPhaseSpeedTest


text:
import Base: unsafe_string, show, print, isempty, length, *, ==, unsafe_string
export fixed_string
using StructIO

## Here we define various string-related utilities for MachO and COFF files
"""
    fixed_string{T}

A fixed-byte string, stored as an integer type (e.g. `T` = `UInt128`, or `T` =
`UInt64`) but displayed and treated as a string.
"""
@io struct fixed_string{T <: Integer}
    data::T
end

"""
    unsafe_string(x::fixed_string)

Convert a `fixed_string` object to a native-Julia `String`
"""
function unsafe_string(x::fixed_string{T}) where {T <: Integer}
    data_array = reinterpret(UInt8, [x.data])
    zero_idx = findfirst(isequal(0x00), data_array)
    if zero_idx === nothing
        zero_idx = sizeof(T) + 1
    end

    str = String(data_array[1:zero_idx-1])
    if !isvalid(str)
        return "$(zero_idx)-byte String of invalid UTF-8 data"
    end
    return str
end

## Various overrides for `fixed_string` that make it act like a `String`
show(io::IO, x::fixed_string) = show(io, unsafe_string(x))
print(io::IO, x::fixed_string) = print(io, unsafe_string(x))
isempty(x::fixed_string) = (x.data & 0xff) == 0
function length(x::fixed_string{T}) where {T <: Integer}
    for idx in 0:(sizeof(T)-1)
        if (x.data & (0xff << (idx*8))) == 0
            return idx
        end
    end
    return sizeof(T)
end

==(x::fixed_string, y::AbstractString) = unsafe_string(x) == y
==(x::AbstractString, y::fixed_string) = y==x
*(a::String, b::fixed_string) = a*unsafe_string(b)
*(a::fixed_string, b::String) = unsafe_string(a)*b

"""
    unsafe_string(io, max_len::Integer)

Read in a null-terminated string, stopping with a maximum length of `max_len`.
"""
function unsafe_string(io, max_len::T) where {T <: Integer}
    str = UInt8[]
    idx = 0
    c = read(io, UInt8)
    while c != 0x00 && idx < max_len
        push!(str, c)
        c = read(io, UInt8)
        idx += 1
    end
    return String(str)
end


text:
export MixPredictor

"""
    MixPredictor()

This uses a Padé-approximation of type (2,1) for prediction.
"""
struct MixPredictor{P1<:AbstractPredictor,P2<:AbstractPredictor} <: AbstractPredictor
    predictor1::P1
    predictor2::P2
end
struct MixPredictorCache{
    P1<:AbstractPredictorCache,
    P2<:AbstractPredictorCache,
} <: AbstractPredictorCache
    predictor1::P1
    predictor2::P2
    predictor1_active::Base.RefValue{Bool}
end

function cache(MP::MixPredictor, H, x, ẋ, t)
    p1 = cache(MP.predictor1, H, x, ẋ, t)
    p2 = cache(MP.predictor2, H, x, ẋ, t)
    MixPredictorCache(p1, p2, Ref(true))
end

function update!(
    cache::MixPredictorCache,
    H,
    x,
    ẋ,
    t,
    Jac::JacobianMonitor,
    eval_err::Float64,
)
    if cache.predictor1_active[]
        update!(cache.predictor1, H, x, ẋ, t, Jac, eval_err)
    else
        update!(cache.predictor2, H, x, ẋ, t, Jac, eval_err)
    end
    cache
end

function predict!(x̂, cache::MixPredictorCache, H::HomotopyWithCache, x, t, Δt, ẋ, J)
    if cache.predictor1_active[]
        predict!(x̂, cache.predictor1, H, x, t, Δt, ẋ, J)
    else
        predict!(x̂, cache.predictor2, H, x, t, Δt, ẋ, J)
    end
end

for f in [:order, :highest_derivative, :second_derivative, :third_derivative]
    @eval begin
        function $f(cache::MixPredictorCache)
            cache.predictor1_active[] ? $f(cache.predictor1) : $f(cache.predictor2)
        end
    end
end


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

JLLWrappers.@generate_main_file_header("Pangolin")
JLLWrappers.@generate_main_file("Pangolin", UUID("883017f5-2c5b-5b91-9ffe-927a00b93392"))
end  # module Pangolin_jll


text:
export FactorEquality

include(joinpath(@__DIR__, "../prob_operations.jl"))
include(joinpath(@__DIR__, "../graph_operations.jl"))
include(joinpath(@__DIR__, "../util.jl"))

mutable struct FactorEquality
    """
    Equality node.

    This node has the function:
    f(x,y,z) = δ(z - x) δ(z - y)
           _
    x --->[=]---> z
           ^
           |
           y
    """

    # Identifiers of edges/nodes in factor graph
    id::String
    time::Int64
    connected_edges::Dict{String, String}

    # Reaction parameters
    incoming::Queue{Tuple}
    beliefs::Dict{String, Any}
    energies::Dict{String, Float64}

    function FactorEquality(id; time=0, edges=[])

        # Keep track of recognition distributions
        connected_edges = Dict{String, String}()
        for edge in edges
            # Edge name is also edge key
            connected_edges[edge] = edge
        end

        # Initialize beliefs for all edges
        beliefs = Dict{String, Any}()

        # Initialize queue for incoming messages
        incoming = Queue{Tuple}()

        # Initialize energy placeholder
        energies = Dict{String,Float64}()
        for edge in edges
            energies[edge] = Inf
        end

        # Create instance
        self = new(id, time, connected_edges, incoming, beliefs, energies)
        return self
    end
end

function energy(node::FactorEquality)
    "Pass on energy from incoming messages."
    # Return the sum of stored energies
    return sum([node.energies[key] for key in keys(node.energies)])
end

function grad_energy(node::FactorEquality, edge_id::String)
    "Pass on energy gradients from other connected nodes."
    #TODO (probably via Automatic Differentiation)
    # Find set of other edges
    other_edges = setdiff(keys(node.energies), [edge_id])

    return 1.0, 1.0
end

function message(node::FactorEquality, edge_id::String)
    """
    Pass on messages from other edges.

    The message to z is computed by ∫ qx(x) qy(y) log[f(x,y,z)] dx dy. The
    sifting property of the delta distributions turns this integral into:
        qx(z)⋅qy(z).
    """

    # Find set of other edges
    other_edges = setdiff(keys(node.beliefs), [edge_id])

    # Construct array of incoming messages from other edges
    incoming_messages = [node.beliefs[edge] for edge in other_edges]

    # Outgoing message is the product of the incoming messages
    outgoing_message = prod(incoming_messages)

    # Pass message
    return outgoing_message
end

function act(node::FactorEquality, edge_id::String, graph::MetaGraph)
    "Send out message for one of the connecting edges"

    # Compute message for a particular edge
    outgoing_message = message(node, edge_id)

    # Extract edge from graph
    edge = eval(Symbol(edge_id))

    # Check if edge is blocked
    if !edge.block

        # Pass message to edge
        edge.messages[node.id] = outgoing_message
    end

    return Nothing
end

function react(node::FactorEquality, graph::MetaGraph)
    "React to incoming messages from edges"

    # Find edges attached to node
    N = neighbors(graph, graph[node.id, :id])
    edge_ids = [graph[n, :id] for n in N]

    # Loop over all edges that have produced incoming messages
    for n = 1:length(node.incoming)

        # Extract the id of the edge and the communicated delta free energy
        edge_id, delta_free_energy = dequeue!(node.incoming)

        # Store incoming energy
        node.energies[edge_id] = delta_free_energy

        # Act on other edges than incoming one
        for edge_out in setdiff(Set(edge_ids), Set([edge_id]))
            act(node, edge_out, graph)
        end
    end

    return Nothing
end

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
# --------------- Parameter Writing ---------------
function writemparams(m::ModelParams, trialname)
	open(trialname * "_mparams.txt", "w") do f write(f, "fs\t\t\t$(m.fs)\n")
		write(f, "num_neurons\t\t$(Float64(m.num_neurons))\n")
		write(f, "p_contact\t\t$(m.p_contact)\n")
		write(f, "p_exc\t\t\t$(m.p_exc)\n")
		write(f, "maxdelay\t\t$(m.maxdelay)\n")
		write(f, "synmax\t\t\t$(m.synmax)\n")
		write(f, "tau_pre\t\t\t$(m.tau_pre)\n")
		write(f, "tau_post\t\t$(m.tau_post)\n")
		write(f, "a_pre\t\t\t$(m.a_pre)\n")
		write(f, "a_post\t\t\t$(m.a_post)\n")
	end
end

function writetparams(t::TrialParams, trialname)
	open(trialname * "_tparams.txt", "w") do f
		write(f, "dur\t\t\t$(t.dur)\n")
		write(f, "lambda\t\t\t$(t.lambda)\n")
		write(f, "randspikesize\t\t$(t.randspikesize)\n")
		write(f, "randinput\t\t$(Float64(t.randinput))\n")
		write(f, "inhibition\t\t$(Float64(t.inhibition))\n")
        write(f, "stdp\t\t$(Float64(t.stdp))\n")
		write(f, "inputmode\t\t$(t.inputmode)\n")
		write(f, "multiinputmode\t\t$(t.multiinputmode)\n")
		write(f, "inputweight\t\t$(t.inputweight)\n")
		write(f, "recordstart\t\t$(t.recordstart)\n")
		write(f, "recordstop\t\t$(t.recordstop)\n")
		write(f, "lambdainput\t\t$(t.lambdainput)\n")
		write(f, "inputrefractorytime\t$(t.inputrefractorytime)\n")
	end
end

# --------------- IO for delnet Communication ---------------
function read(s::IO, ::Type{Spike})
    n = read(s,Int64)
    t = read(s,Float64)
	Spike(n,t)
end

function readneurons(trialname)
    neurons = Nothing
	open(trialname*"_neurons.bin", "r") do f
		n = Base.read(f, Int64)
        neurons = Array{SimpleNeuron,1}(undef,n)
		for k ∈ 1:n
            a = Base.read(f, Float64)
            b = Base.read(f, Float64)
            c = Base.read(f, Float64)
            d = Base.read(f, Float64)
            neurons[k] = SimpleNeuron(a,b,c,d)
		end
	end
    neurons
end

function writeneurons(neurons::Array{SimpleNeuron,1}, trialname)
	open(trialname*"_neurons.bin", "w") do f
		write(f, length(neurons))
		for neuron ∈ neurons
			write(f, neuron.a)
			write(f, neuron.b)
			write(f, neuron.c)
			write(f, neuron.d)
		end
	end
end

function writeinput(input::Array{Spike,1}, trialname; weight=18.0)
	open(trialname * "_input.bin", "w") do f
		write(f,Int64(1)) 			# number of inputs (only 1 here)
		write(f,length(input)) 		# write an Int64 (long int)
        for spike ∈ input
			write(f,spike.n-1)
			write(f,spike.t)
		end
        for _ ∈ 1:length(input)
            write(f,weight)
        end
	end
end

function writeinput(input::Array{Spike,1}, weights::Array{Float64,1}, trialname)
	open(trialname * "_input.bin", "w") do f
		write(f,Int64(1)) 			# number of inputs (only 1 here)
		write(f,length(input)) 		# write an Int64 (long int)
        @assert length(weights) == length(input)
        for spike ∈ input
			write(f,spike.n-1)
			write(f,spike.t)
		end
        for weight ∈ weights
            write(f,weight)
        end
	end
end

function writeinput(inputs::Array{Array{Spike,1},1},
                    weights::Array{Array{Float64,1},1}, trialname)
	open(trialname * "_input.bin", "w") do f
		write(f,length(inputs)) 			# number of inputs
        for (i, input) ∈ enumerate(inputs)
			write(f,length(input)) 		# write an Int64 (long int)
            for spike ∈ input
				write(f,spike.n-1)
				write(f,spike.t)
			end
            for weight ∈ weights[i]
                write(f,weight)
            end
		end
	end
end

function writeinput(inputs::Array{Array{Spike,1},1}, trialname; weight=18.0)
	open(trialname * "_input.bin", "w") do f
		write(f,length(inputs)) 			# number of inputs
        for input ∈ inputs
			write(f,length(input)) 		# write an Int64 (long int)
            for spike ∈ input
				write(f,spike.n-1)
				write(f,spike.t)
			end
            for _ ∈ 1:length(input)
                # println(weight)
                write(f,weight)
            end
		end
	end
end

function loadinput(trialname)
	open(trialname * "_input.bin", "r") do f
		numinputs = Base.read(f,Int64)
		inputs = Array{Array{Spike,1},1}(undef,0)
		for _ ∈ 1:numinputs
			len = Base.read(f,Int64)
			input = Array{Spike, 1}(undef,len)
			#input = read(f, input)
			for k ∈ 1:len
				input[k] = Spike(Base.read(f, Int64) + 1, Base.read(f, Float64))
			end
			push!(inputs, input)
		end
		return length(inputs) == 1 ? inputs[1] : inputs 	# type stability here...
	end
end

function writedelgraph(graph::Array{Int64,2}, trialname)
	open(trialname * "_graph.bin", "w") do f
		n = size(graph)[1]
		graphraw = reshape(graph', (n*n,))
		write(f, Int64(n)) 	# write an Int64 (long int)
		write(f, graphraw)
	end
end

function loaddelgraph(trialname)
	open(trialname * "_graph.bin", "r") do f
		n = Base.read(f, Int64)
		graphraw = Array{UInt64, 1}(undef, n*n)
		read!(f, graphraw);
		graph = reshape(graphraw, (n, n))'
		return graph
	end
end

function writesyngraph(graph::Array{Float64,2}, trialname)
	open(trialname * "_syngraph.bin", "w") do f
		n = size(graph)[1]
		graphraw = reshape(graph', (n*n,))
		write(f, Int64(n*n)) 	# write an Int64 (long int)
		write(f, graphraw)
	end
end

function loadbinvector(filename; datatype=Float64)
	open(filename, "r") do f
		len = Base.read(f, UInt64)
		input = Array{datatype, 1}(undef, Int64(len))
		read!(f, input)
		return input
	end
end

function loadsynapses(trialname)
	filename = trialname*"_synapses.bin"
	open(filename, "r") do f
		len = Int64(Base.read(f, UInt64))
		strengths = Array{Float64,1}(undef, len)
		sources = Array{UInt64,1}(undef, len)
		dests = Array{UInt64,1}(undef, len)
		delays = Array{UInt16,1}(undef, len)
		read!(f, strengths)
		read!(f, sources)
		read!(f, dests)
		read!(f, delays)
		return [Synapse(strengths[i], sources[i], dests[i], delays[i]) for i ∈ 1:len]
	end

end

function loadspikesmat(trialname)
	open(trialname * "_spikes.txt", "r") do f
		lines = readlines(f)
		ts = [parse(Float64, split(line, "  ")[1]) for line ∈ lines]
		ns = [parse(Int64, split(line, "  ")[2]) for line ∈ lines]
		return [ts ns]
	end
end

function loadinputtimes(trialname)
	open(trialname * "_instarttimes.txt", "r") do f
		lines = readlines(f)
		ts = [parse(Float64, split(line, "  ")[1]) for line ∈ lines]
		ids = [parse(Int64, split(line, "  ")[2])+1 for line ∈ lines]
		return ts, ids
	end
end

function loadspikes(trialname::String)
	open(trialname * "_spikes.txt", "r") do f
		lines = readlines(f)
		n = length(lines)
		spikes = Array{Spike,1}(undef, n)
		for k ∈ 1:n
			# Note "+1" for 1-based indexing in Julia
			spikes[k] = Spike( parse(Int64, split(lines[k], "  ")[2])+1,
							   parse(Float64, split(lines[k], "  ")[1]))
		end
		return spikes
	end
end

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
module SequenceJacobians

using ForwardDiff
using Graphs: AbstractGraph, Edge, SimpleDiGraphFromIterator, topological_sort_by_dfs
using LinearAlgebra: I, UniformScaling, LU, lu!, ldiv!, norm
using LinearMaps

import Base: ==, eltype, zero, show, convert
import Graphs: SimpleDiGraph, edgetype, nv, ne, vertices, edges, is_directed,
    has_vertex, has_edge, inneighbors, outneighbors, neighborhood

# Reexport objects from Graphs
export SimpleDiGraph, edgetype, nv, ne, vertices, edges, is_directed, has_vertex, has_edge,
    inneighbors, outneighbors

export Shift,
       Lag,
       Lead,
       ShiftMap,

       VarInput,
       var,
       lag,
       lead,
       name,
       shift,
       AbstractBlock,
       inputs,
       invars,
       ssinputs,
       outputs,
       SimpleBlock,
       block,
       steadystate!,
       jacobian,
       transition!,

       BlockOrVar,
       SequenceSpaceModel,
       model,
       srcs,
       sssrcs,
       dests,
       isblock,
       SteadyState,
       residuals!,
       criterion!,

       TotalJacobian,
       GEJacobian,
       getG!,

       Transition,
       update!,
       solve!,

       linirf,
       nlirf

include("shift.jl")
include("blocks.jl")
include("model.jl")
include("jacobian.jl")
include("transition.jl")
include("irf.jl")

end # module


text:
using GRIB, DataFrames, Query, DataStructures, H3, H3.API, H3Geometry, GZip, Dates
using Arrow, Parquet, TimeZones, Printf, Distributed, Dagger

"""
Reads downloaded MRMS weather data and arranges it in a dataframe with
columns for latitude, longitude, time, h3 key, tile id and precipitation.

Time is represented as milliseconds since 1970 instead of a proper ZonedDateTime
in order to make storing the result using Parquet easy. Likewise, the H3
fields (h3, tile) are converted to signed integers.

Due to memory constraints, it isn't practical to ingest all of the data
over a long time period into a DataFrame (even one day can be taxing).
As a result, the normal process is to read each hour of data separately,
group by tile and write the data one file per tile retaining only time, h3 and
precipitation in order to save space. Subsequently, the data for a particular
tile can be concatenated and sorted by h3 and time in order to build a useful
historical dataset.

The DataFrame that we return here uses RepeatVectors for latitude, longitude,

All of these steps come down to two functions. The first splits a
single GRIB file into files, one for each hex at a specifed level
(levels 3,4 and 5 result in 2025, 13709, and 94718 files respectively). The
second function merges corresponding files across multiple directories
into a composite output. This arrangement means that we only need to
store at most one GRIB file or data for one tile of output.

"""
function readData(inputFile::String, tileLevel::Int)
    t = @elapsed begin
        if endswith(inputFile, ".gz")
            mktemp() do tmp, tmp_io
                GZip.open(inputFile, "r") do fh
                    data = read(fh)
                    write(tmp_io, data)
                end
                dataset = GribFile(tmp)
            end
        else
            dataset = GribFile(inputFile)
        end
        m = first(dataset)
        lons, lats, v = data(m)
        t0 = ZonedDateTime((m[i] for i in ("year", "month", "day", "hour", "minute", "second"))..., tz"UTC")
        t1 = Int((t0 - ZonedDateTime(1970, 1, 1, tz"UTC")) / Millisecond(1))
    end

    @info "read data" size=length(lons) time=t file=inputFile

    t = @elapsed begin
        # small cleanups on the data
        lons = collapse(Vector{Float32}(reshape(lons, :)))
        lats = collapse(Vector{Float32}(reshape(lats, :)))
        v = Vector{Float32}(reshape(v, :))
        lons.base[lons.base .> 180.0] .-= 360.0
        time = RepeatVector([t1], length(lons), length(lons))
    end
    @info "cleaned up data" elapsed=t time="$(length(time.base)) $(time.each)" lons="$(length(lons.base)) $(lons.each)" lats="$(length(lats.base)) $(lats.each)"

    # some geo data added
    t = @elapsed begin
        tile = zeros(Int64, length(lons))
        pt = zeros(Int64, length(lons))
        Threads.@threads for i in 1:length(lons)
            tile[i] = Int64(geoToH3(lats[i], lons[i], tileLevel))
            pt[i] = Int64(geoToH3(lats[i], lons[i], 15))
        end
    end
    @info "geo time" t=t
    base = (longitude=lons, latitude=lats, tile=tile, h3=pt, t=time)
    base = (base..., precipitation=v)
    return DataFrame(Tables.columntable(Tables.Schema(keys(base), typeof.(values(base))), base))
end

"""
Output a data frame into files, each containing data for a single geo-tile.
The data frame is assumed to contain latitude, longitude, h3-id, tile-id,
time plus other data columns. Only h3-id, time and the other data are preserved
in order to minimize the size of the saved data. The data saved in each tile
file is ordered by h3 and then by time before being written out.
"""
function splitData(data, outputDir::String, hashDir::String, tileLevel)
    outCount = Threads.Atomic{Int}(0);
    step = 1
    scale = 1
    mkpath(outputDir)
    mkpath(hashDir)
    t0 = time()
    n0 = 0
    groups = [g for g in data |> @groupby(_.tile)]
    for td in groups
        tileData = DataFrame(td)
        @assert size(tileData)[1] > 0
        tile = tileData[1, :tile]
        fname = joinpath(outputDir, "x-$tileLevel-$tile")
        outputData = tileData |>
            @orderby(_.h3) |> @thenby(_.t) |>
            @select(-:longitude, -:latitude, -:tile)
        Parquet.write_parquet(fname, outputData)
        open(fname) do input
            h = hash(read(input))
            open(joinpath(hashDir, "x-$tileLevel-$tile"), "w") do output
                println(output, h)
            end
            try
                df0 = DataFrame(Parquet.Table(fname))
                @assert collapse(tileData[:,:t])==collapse(df0[:,:t])
                @assert length(unique(tileData[:,:h3]))==length(unique(df0[:,:h3]))
            catch e
                @warn "failed" fname
                throw(e)
            end
        end

        k = Threads.atomic_add!(outCount, 1)
        if k%(step*scale) == 0
            t1 = time()
            t = t1 - t0
            step = floor(step * 2.5)
            if step >= 10
                scale *= 10
                step = 1
            end
            @info "file" files=k time=t rate=(k-n0)/t s=step*scale
            n0 = k
        end
    end
end

"""
With roughly 25 million points per day and 100 day history, we need to
store about a 60 billion data points. Splitting these up into 1000
files means that each file only has a 25,000 histories containing
about 60 million data points which might be small enough for fast
access, but probably not. On the other hand, splitting into 100,000
files (which is still quite practical) results in each file having (on
average) 100 histories and about 240,000 data points. This is almost
certainly small enough to be fast to scan at query time.

Creating these files is done by splitting up each hourly file into
tile files and then merging all of the files from each hour into the
daily files when each hourly file becomes available. Due to space
and file count issues, it is a good thing to remove hourly splits
as soon as practical (they can always be recreated).

This merge function works by scanning all of the directories to
get a complete list of all of the files that exist. Then, for
each file name, we concatenate all of the data in the file by
that name in any directory and write that concatenation out under
the same name in the output directory. When that's all done, the
output directory can be atomically renamed to be the final output
directory if we want to avoid exposing half done work.

"""
function mergeData(output::String, dirs::String...)
    mkpath(output)
    dirs = sort([x for x in dirs])
    # compute list of all files
    files = union(readdir.(dirs)...) |> @filter(x -> startswith(x,"x-"))

    k = 0
    for f in files
        df = Vector()
        for fx in joinpath.(dirs, f)
            tbl = Parquet.read_parquet(fx)
            push!(df, DataFrame(tbl))
        end
        all = vcat(df...)
        Parquet.write_parquet(joinpath(output, f), all)
        k += 1
        if k % 100 == 0 || (k < 100 && k % 10 == 0)
            @info "progress" processed=k
        end
    end
end

"""
This scans input directories to find MRMS hourly data files. For
each file found, that file is split into geo-tiles (unless it has
already been split).

This is done across as many workers as are available because
it can be kind of slow. Once this is done, the hours are
available for merging.
"""
function processDays(source, dest, first, last, tileLevel)
    @info "before"
    files = []
    date = first
    while date <= last
        day = Dates.format(date, "yyyy/mm/dd")
        dailyDir = joinpath(source, day)
        @info "Day" dailyDir
        hour = 0
        if isdir(dailyDir)
            for f in sort(readdir(dailyDir, join=true))
                data = joinpath(dest, "hourly/data", day, @sprintf("%02d", hour))
                if !isdir(data)
                    push!(files, (f,day,hour))
                end
                hour += 1
            end
        end
        date += Dates.Day(1)
    end
    @info "Starting" length(files)
    @distributed for f in 1:length(files)
        fname,day,hour = files[f]
        data = joinpath(dest, "hourly/data", day, @sprintf("%02d", hour))
        hash = joinpath(dest, "hourly/hash", day, @sprintf("%02d", hour))
        mkpath(data)
        mkpath(hash)

        m = readData(fname, tileLevel)
        @info "Read file" size(m)
        @info "Split data" data hash
        splitData(m, data, hash, tileLevel)
    end
end

"""
A RepeatVector looks like a vector with a fixed cycle of repeated elements
but only requires a small amount of storage
"""
struct RepeatVector{T} <: AbstractVector{T}
    base::Vector{T}
    each::Int64
    n::Int64
end

"""
Assuming that a vector contains lots of very regular repetition,
collapse will extract a RepeatVector vector which requires less
space.
"""
function collapse(input::Vector{T})::RepeatVector{T} where (T)
    each = 0
    x1 = first(input)
    for x in input
        if x == x1
            each += 1
        else
            break
        end
    end
    base = Vector{T}()
    push!(base, x1)
    x2 = x1
    for x in input[(each+1):end]
        if x == x1
            break
        end
        if x != x2
            push!(base, x)
            x2 = x
        end
    end
    return RepeatVector(base, each, length(input))
end

function Base.getindex(v::RepeatVector, i::Int)
    i -= 1
    i = i ÷ v.each
    return v.base[i % length(v.base) + 1]
end

function Base.size(v::RepeatVector)
    return (v.n,)
end

max_stars_repo_name: agstack/weather-server

max_issues_repo_name: agstack/weather-server

max_forks_repo_name: agstack/weather-server


text:
using StringLiterals
using Format

@static VERSION < v"0.7.0-DEV" ? (using Base.Test) : (using Test)
eval_parse(s) = eval((@static VERSION < v"0.7.0-DEV.2995" ? parse : Meta.parse)(s))
const ErrorType = @static VERSION < v"0.7.0-DEV" ? ArgumentError : LoadError

ts(io) = String(take!(io))

@testset "LaTeX Entities" begin
    @test f"\<dagger>" == "†"
    #@test f"\<mscrl>" == "𝓁" # \U1f4c1
    @test f"\<c_l>" == "𝓁" # \U1f4c1
    @test f"\<nleqslant>" == "⩽̸" # \u2a7d\u338
end
@testset "Emoji Entities" begin
    @test f"\:sailboat:" == "\u26f5"
    @test f"\:ring:"     == "\U1f48d"
    @test f"\:flag-us:"  == "\U1f1fa\U1f1f8"
end
@testset "Unicode Entities" begin
    @test f"\N{end of text}" == "\x03" # \3
    @test f"\N{TIBETAN LETTER -A}" == "\u0f60"
    @test f"\N{LESS-THAN OR SLANTED EQUAL TO}" == "\u2a7d"
    @test f"\N{REVERSED HAND WITH MIDDLE FINGER EXTENDED}" == "\U1f595"
end
@testset "HTML Entities" begin
    @test f"\≰"    == "\u2270"
    @test f"\𝒫"   == "\U1d4ab"
    @test f"\≨"   == "\u2268\ufe00"
end
@testset "Unicode constants" begin
    @test f"\u{3}"     == "\x03"
    @test f"\u{f60}"   == "\u0f60"
    @test f"\u{2a7d}"  == "\u2a7d"
    @test f"\u{1f595}" == "\U0001f595"
end
@testset "Interpolation" begin
    scott = 123
    @test f"\(scott)" == "123"
end
@testset "\$ not interpolation" begin
    @test f"I have $10, $spj$" == "I have \$10, \$spj\$"
end
@testset "Valid quoted characters" begin
    @test f"\$" == "\$"
    @test f"\"" == "\""
    @test f"\'" == "\'"
    @test f"\\" == "\\"
    @test f"\0" == "\0"
    @test f"\a" == "\a"
    @test f"\b" == "\b"
    @test f"\e" == "\e"
    @test f"\f" == "\f"
    @test f"\n" == "\n"
    @test f"\r" == "\r"
    @test f"\t" == "\t"
    @test f"\v" == "\v"
end
@testset "Invalid quoted characters" begin
    for ch in "cdghijklmopqsuwxy"
        @test_throws ErrorType eval_parse(string("f\"\\", ch, '"'))
    end
    for ch in 'A':'Z'
        @test_throws ErrorType eval_parse(string("f\"\\", ch, '"'))
    end
end

@testset "Legacy mode only sequences" begin
    # Check for ones allowed in legacy mode
    for s in ("f\"\\x\"", "f\"\\x7f\"", "f\"\\u\"", "f\"\\U\"", "f\"\\U{123}\"")
        @test_throws ErrorType eval_parse(s)
    end
end

@testset "Legacy mode Hex constants" begin
    @test F"\x3"     == "\x03"
    @test F"\x7f"    == "\x7f"
    @test_throws ErrorType eval_parse("F\"\\x\"")
    @test_throws ErrorType eval_parse("F\"\\x!\"")
end

@testset "Legacy mode Unicode constants" begin
    @test F"\u3"     == "\x03"
    @test F"\uf60"   == "\u0f60"
    @test F"\u2a7d"  == "\u2a7d"
    @test F"\U1f595" == "\U0001f595"
    @test F"\u{1f595}" == "\U0001f595"
    @test_throws ErrorType eval_parse("F\"\\U\"")
    @test_throws ErrorType eval_parse("F\"\\U!\"")
    @test_throws ErrorType eval_parse("F\"\\U{123}\"")
end

@testset "Legacy mode valid quoted characters" begin
    @test f"\$" == "\$"
    @test f"\"" == "\""
    @test f"\'" == "\'"
    @test f"\\" == "\\"
    @test f"\0" == "\0"
    @test f"\a" == "\a"
    @test f"\b" == "\b"
    @test f"\e" == "\e"
    @test f"\f" == "\f"
    @test f"\n" == "\n"
    @test f"\r" == "\r"
    @test f"\t" == "\t"
    @test f"\v" == "\v"
end

@testset "Legacy mode \$ interpolation" begin
    scott = 123
    @test F"$scott" == "123"
    @test F"$(scott+1)" == "124"
end

@testset "Quoted \$ in Legacy mode" begin
    @test F"I have \$10, \$spj\$" == "I have \$10, \$spj\$"
end

@testset "C Formatting" begin
    @testset "int" begin
        @test f"\%d(typemax(Int64))" == "9223372036854775807"
        @test f"\%a(typemax(Int64))" == "0x7.fffffffffffffffp+60"
        @test f"\%A(typemax(Int64))" == "0X7.FFFFFFFFFFFFFFFP+60"
    end
    @testset "printing an int value" begin
        for num in (UInt16(42), UInt32(42), UInt64(42), UInt128(42),
                    Int16(42), Int32(42), Int64(42), Int128(42))
            @test f"\%i(num)"       == "42"
            @test f"\%u(num)"       == "42"
            @test f"Test: \%i(num)" == "Test: 42"
            @test f"\%#x(num)"      == "0x2a"
            @test f"\%#o(num)"      == "052"
            @test f"\%x(num)"       == "2a"
            @test f"\%X(num)"       == "2A"
            @test f"\% i(num)"      == " 42"
            @test f"\%+i(num)"      == "+42"
            @test f"\%4i(num)"      == "  42"
            @test f"\%-4i(num)"     == "42  "
            @test f"\%a(num)"       == "0x2.ap+4"
            @test f"\%A(num)"       == "0X2.AP+4"
            @test f"\%20a(num)"     == "            0x2.ap+4"
            @test f"\%-20a(num)"    == "0x2.ap+4            "
            @test f"\%f(num)"       == "42.000000"
            @test f"\%g(num)"       == "42"
        end
    end
    @testset "pointers" begin
        @static if Sys.WORD_SIZE == 64
                    @test f"\%20p(0)"  == "  0x0000000000000000"
                    @test f"\%-20p(0)" == "0x0000000000000000  "
                elseif Sys.WORD_SIZE == 32
                    @test f"\%20p(0)"  == "          0x00000000"
                    @test f"\%-20p(0)" == "0x00000000          "
                else
                    @test false
                end
    end
    @testset "float / BigFloat" begin
        for num in (1.2345, big"1.2345")
            @test f"\%7.2f(num)"  == "   1.23"
            @test f"\%-7.2f(num)" == "1.23   "
            @test f"\%07.2f(num)" == "0001.23"
            @test f"\%.0f(num)"   == "1"
            @test f"\%#.0f(num)"  == "1."
            @test f"\%.4e(num)"   == "1.2345e+00"
            @test f"\%.4E(num)"   == "1.2345E+00"
            @test f"\%.2a(num)"   == "0x1.3cp+0"
            @test f"\%.2A(num)"   == "0X1.3CP+0"
        end
    end
    @testset "Inf / NaN handling" begin
        bf = big"Inf"
        bn = big"NaN"
        @test f"\%f(Inf)" == "Inf"
        @test f"\%f(NaN)" == "NaN"
        @test f"\%f(bf)"  == "Inf"
        @test f"\%f(bn)"  == "NaN"
    end
    @testset "scientific notation" begin
        fv = 3e142
        bv = big"3e142"
        @test f"\%.0e(fv)"  == "3e+142"
        @test f"\%#.0e(fv)" == "3.e+142"
        @test f"\%.0e(bv)"  == "3e+142"
        @test f"\%#.0e(bv)" == "3.e+142"
        bv = big"3e1042"
        @test f"\%.0e(bv)"   == "3e+1042"
        @test f"\%e(3e42)"   == "3.000000e+42"
        @test f"\%E(3e42)"   == "3.000000E+42"
        @test f"\%e(3e-42)"  == "3.000000e-42"
        @test f"\%E(3e-42)"  == "3.000000E-42"
        @test f"\%a(3e4)"    == "0x1.d4cp+14"
        @test f"\%A(3e4)"    == "0X1.D4CP+14"
        @test f"\%.4a(3e-4)" == "0x1.3a93p-12"
        @test f"\%.4A(3e-4)" == "0X1.3A93P-12"
    end
    @testset "%g" begin
        for (val, res) in (
            (12345678., "1.23457e+07"),
            (1234567.8, "1.23457e+06"),
            (123456.78, "123457"),
            (12345.678, "12345.7"),
            (12340000.0, "1.234e+07"),
            (big"12345678.", "1.23457e+07"),
            (big"1234567.8", "1.23457e+06"),
            (big"123456.78", "123457"),
            (big"12345.678", "12345.7"))
            @test f"\%.6g(val)" == res
        end
        for num in (123.4, big"123.4")
            @test f"\%10.5g(num)"   == "     123.4"
            @test f"\%+10.5g(num)"  == "    +123.4"
            @test f"\% 10.5g(num)"  == "     123.4"
            @test f"\%#10.5g(num)"  == "    123.40"
            @test f"\%-10.5g(num)"  == "123.4     "
            @test f"\%-+10.5g(num)" == "+123.4    "
            @test f"\%010.5g(num)"  == "00000123.4"
        end
        bv = big"-123.4"
        b2 = big"12340000.0"
        @test f"\%10.5g(-123.4)"    == "    -123.4"
        @test f"\%010.5g(-123.4)"   == "-0000123.4"
        @test f"\%.6g(12340000.0)"  == "1.234e+07"
        @test f"\%#.6g(12340000.0)" == "1.23400e+07"
        @test f"\%10.5g(bv)"        == "    -123.4"
        @test f"\%010.5g(bv)"       == "-0000123.4"
        @test f"\%.6g(b2)"          == "1.234e+07"
        @test f"\%#.6g(b2)"         == "1.23400e+07"
        @test f"\%.5g(42)"          == "42"
        @test f"\%#.2g(42)"         == "42."
        @test f"\%#.5g(42)"         == "42.000"
    end

    @testset "hex float" begin
        bv = big"1.5"
        @test f"\%a(1.5)"    == "0x1.8p+0"
        @test f"\%a(1.5f0)"  == "0x1.8p+0"
        @test f"\%a(bv)"     == "0x1.8p+0"
        @test f"\%#.0a(1.5)" == "0x2.p+0"
        @test f"\%+30a(1/3)" == "         +0x1.5555555555555p-2"
    end

    @testset "chars" begin
        @test f"\%c(65)"    == "A"
        @test f"\%c('A')"   == "A"
        @test f"\%3c('A')"  == "  A"
        @test f"\%-3c('A')" == "A  "
        @test f"\%c(248)"   == "ø"
        @test f"\%c('ø')"   == "ø"
    end

    @testset "type width specifier parsing (ignored)" begin
        @test f"\%llf(1.2)" == "1.200000"
        @test f"\%Lf(1.2)"  == "1.200000"
        @test f"\%hhu(1)"   == "1"
        @test f"\%hu(1)"    == "1"
        @test f"\%lu(1)"    == "1"
        @test f"\%llu(1)"   == "1"
        @test f"\%Lu(1)"    == "1"
        @test f"\%zu(1)"    == "1"
        @test f"\%ju(1)"    == "1"
        @test f"\%tu(1)"    == "1"
    end

    @testset "strings" begin
        s1 = "test"
        s2 = "tést"

        @test f"\%s(s1)"   == s1
        @test f"\%s(s2)"   == s2

        @test f"\%8s(s1)"  == "    test"
        @test f"\%-8s(s1)" == "test    "

        # This was broken in v0.5
        @static if VERSION >= v"0.6-"
            @test f"\%8.3s(s1)"    == "     tes"
            @test f"\%#8.3s(s1)"   == "     \"te"
            @test f"\%-8.3s(s1)"   == "tes     "
            @test f"\%#-8.3s(s1)"  == "\"te     "
            @test f"\%.3s(s1)"     == "tes"
            @test f"\%#.3s(s1)"    == "\"te"
            @test f"\%-.3s(s1)"    == "tes"
            @test f"\%#-.3s(s1)"   == "\"te"
        end
    end

    @testset "type based formatting" begin
        fmt_default!()        # resets all defaults

        # some basic functionality testing
        x = 1234.56789

        @test f"\%(x)"             == "1234.567890"
        @test f"\%(x;prec=2)"      == "1234.57"
        @test f"\%(x,10,3)"        == "  1234.568"
        @test f"\%(x,10,3,:left)"  == "1234.568  "
        @test f"\%(x,10,3,:ljust)" == "1234.568  "
        @test f"\%(x,:commas)"     == "1,234.567890"

        i = 1234567

        @test f"\%(i)" == "1234567"
        @test f"\%(i,:commas)" == "1,234,567"

        fmt_default!(Int, :commas, width = 12)
        @test f"\%(i)" == "   1,234,567"
        @test f"\%(x)" == "1234.567890"  # default hasn't changed

        fmt_default!(:commas)
        @test f"\%(i)" == "   1,234,567"
        @test f"\%(x)" == "1,234.567890"  # width hasn't changed, but added commas

        fmt_default!(Int) # resets Integer defaults
        @test f"\%(i)" == "1234567"
        @test f"\%(i,:commas)" == "1,234,567"
    end
end

@testset "Python Formatting" begin
    @testset "Format string" begin
        s = "abc"
        @test f"\{}(s)"     == "abc"
        @test f"\{s}(s)"    == "abc"
        @test f"\{2s}(s)"   == "abc"
        @test f"\{5s}(s)"   == "abc  "
        @test f"\{>5s}(s)"  == "  abc"
        @test f"\{*>5s}(s)" == "**abc"
        @test f"\{*<5s}(s)" == "abc**"
    end

    @testset "Format Char" begin
        @test f"\{}('c')"     == "c"
        @test f"\{c}('c')"    == "c"
        @test f"\{3c}('c')"   == "c  "
        @test f"\{>3c}('c')"  == "  c"
        @test f"\{*>3c}('c')" == "**c"
        @test f"\{*<3c}('c')" == "c**"
    end

    @testset "Format integer" begin
        @test f"\{}(1234)" == "1234"
        @test f"\{d}(1234)" == "1234"
        @test f"\{n}(1234)" == "1234"
        @test f"\{x}(0x2ab)" == "2ab"
        @test f"\{X}(0x2ab)" == "2AB"
        @test f"\{o}(0o123)" == "123"
        @test f"\{b}(0b1101)" == "1101"

        @test f"\{d}(0)" == "0"
        @test f"\{d}(9)" == "9"
        @test f"\{d}(10)" == "10"
        @test f"\{d}(99)" == "99"
        @test f"\{d}(100)" == "100"
        @test f"\{d}(1000)" == "1000"

        @test f"\{06d}(123)" == "000123"
        @test f"\{+6d}(123)" == "  +123"
        @test f"\{+06d}(123)" == "+00123"
        @test f"\{ d}(123)" == " 123"
        @test f"\{ 6d}(123)" == "   123"
        @test f"\{<6d}(123)" == "123   "
        @test f"\{>6d}(123)" == "   123"
        @test f"\{*<6d}(123)" == "123***"
        @test f"\{*>6d}(123)" == "***123"
        @test f"\{< 6d}(123)" == " 123  "
        @test f"\{<+6d}(123)" == "+123  "
        @test f"\{> 6d}(123)" == "   123"
        @test f"\{>+6d}(123)" == "  +123"

        @test f"\{+d}(-123)" == "-123"
        @test f"\{-d}(-123)" == "-123"
        @test f"\{ d}(-123)" == "-123"
        @test f"\{06d}(-123)" == "-00123"
        @test f"\{<6d}(-123)" == "-123  "
        @test f"\{>6d}(-123)" == "  -123"
    end

    @testset "Format floating point (f)" begin

        @test f"\{}(0.125)" == "0.125"
        @test f"\{f}(0.0)" == "0.000000"
        @test f"\{f}(0.001)" == "0.001000"
        @test f"\{f}(0.125)" == "0.125000"
        @test f"\{f}(1.0/3)" == "0.333333"
        @test f"\{f}(1.0/6)" == "0.166667"
        @test f"\{f}(-0.125)" == "-0.125000"
        @test f"\{f}(-1.0/3)" == "-0.333333"
        @test f"\{f}(-1.0/6)" == "-0.166667"
        @test f"\{f}(1234.5678)" == "1234.567800"
        @test f"\{8f}(1234.5678)" == "1234.567800"

        @test f"\{8.2f}(8.376)" == "    8.38"
        @test f"\{<8.2f}(8.376)" == "8.38    "
        @test f"\{>8.2f}(8.376)" == "    8.38"
        @test f"\{8.2f}(-8.376)" == "   -8.38"
        @test f"\{<8.2f}(-8.376)" == "-8.38   "
        @test f"\{>8.2f}(-8.376)" == "   -8.38"

        @test f"\{<08.2f}(8.376)" == "00008.38"
        @test f"\{>08.2f}(8.376)" == "00008.38"
        @test f"\{<08.2f}(-8.376)" == "-0008.38"
        @test f"\{>08.2f}(-8.376)" == "-0008.38"
        @test f"\{*<8.2f}(8.376)" == "8.38****"
        @test f"\{*>8.2f}(8.376)" == "****8.38"
        @test f"\{*<8.2f}(-8.376)" == "-8.38***"
        @test f"\{*>8.2f}(-8.376)" == "***-8.38"

        @test f"\{.2f}(0.999)" == "1.00"
        @test f"\{.2f}(0.996)" == "1.00"
        # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
        @test (f"\{.2f}(0.995)" == "1.00" || f"\{.2f}(0.995)" == "0.99")
        @test f"\{.2f}(0.994)" == "0.99"
    end

    @testset "Format floating point (e)" begin

        @test f"\{E}(0.0)" == "0.000000E+00"
        @test f"\{e}(0.0)" == "0.000000e+00"
        @test f"\{e}(0.001)" == "1.000000e-03"
        @test f"\{e}(0.125)" == "1.250000e-01"
        @test f"\{e}(100/3)" == "3.333333e+01"
        @test f"\{e}(-0.125)" == "-1.250000e-01"
        @test f"\{e}(-100/6)" == "-1.666667e+01"
        @test f"\{e}(1234.5678)" == "1.234568e+03"
        @test f"\{8e}(1234.5678)" == "1.234568e+03"

        @test f"\{<12.2e}(13.89)" == "1.39e+01    "
        @test f"\{>12.2e}(13.89)" == "    1.39e+01"
        @test f"\{*<12.2e}(13.89)" == "1.39e+01****"
        @test f"\{*>12.2e}(13.89)" == "****1.39e+01"
        @test f"\{012.2e}(13.89)" == "00001.39e+01"
        @test f"\{012.2e}(-13.89)" == "-0001.39e+01"
        @test f"\{+012.2e}(13.89)" == "+0001.39e+01"

        @test f"\{.1e}(0.999)" == "1.0e+00"
        @test f"\{.1e}(0.996)" == "1.0e+00"
        # Floating point error can upset this one (i.e. 0.99500000 or 0.994999999)
        @test (f"\{.1e}(0.995)" == "1.0e+00" || f"\{.1e}(0.995)" == "9.9e-01")
        @test f"\{.1e}(0.994)" == "9.9e-01"
        @test f"\{.1e}(0.6)" == "6.0e-01"
        @test f"\{.1e}(0.9)" == "9.0e-01"
    end

    @testset "Format special floating point value" begin

        @test f"\{f}(NaN)" == "NaN"
        @test f"\{e}(NaN)" == "NaN"
        @test f"\{f}(NaN32)" == "NaN"
        @test f"\{e}(NaN32)" == "NaN"

        @test f"\{f}(Inf)" == "Inf"
        @test f"\{e}(Inf)" == "Inf"
        @test f"\{f}(Inf32)" == "Inf"
        @test f"\{e}(Inf32)" == "Inf"

        @test f"\{f}(-Inf)" == "-Inf"
        @test f"\{e}(-Inf)" == "-Inf"
        @test f"\{f}(-Inf32)" == "-Inf"
        @test f"\{e}(-Inf32)" == "-Inf"

        @test f"\{<5f}(Inf)" == "Inf  "
        @test f"\{>5f}(Inf)" == "  Inf"
        @test f"\{*<5f}(Inf)" == "Inf**"
        @test f"\{*>5f}(Inf)" == "**Inf"
    end

    @testset "test commas..." begin
        @test f"\%'d(1000)" == "1,000"
        @test f"\%'d(-1000)" == "-1,000"
        @test f"\%'d(100)" == "100"
        @test f"\%'d(-100)" == "-100"
        @test f"\%'f(Inf)" == "Inf"
        @test f"\%'f(-Inf)" == "-Inf"
        @test f"\%'s(1000.0)" == "1,000.0"
        @test f"\%'s(1234567.0)" == "1.234567e6"
    end # testset test commas
end

@testset "Print macro support" begin
    io = IOBuffer()
    scott = 123
    pr"\(io)This is a test with \(scott)"
    @test ts(io) == "This is a test with 123"
end

@testset "escape, unescape" begin
    @test s_escape_string(f"' \" \\ \u{7f} \u{20ac} \u{1f596} \u{e0000}") ==
        "' \\\" \\\\ \\u{7f} € 🖖 \\u{e0000}"
    @test s_unescape_string(f"' \\\" \\\\ \\u{7f} € 🖖 \\u{e0000}") ==
        "' \" \\ \x7f € 🖖 \Ue0000"
    io = IOBuffer()
    s_print_escaped(io, f"' \" \\ \u{7f} \u{20ac} \u{1f596} \u{e0000}", "")
    @test ts(io) == f"' \" \\\\ \\u{7f} € 🖖 \\u{e0000}"
    io = IOBuffer()
    s_print_unescaped(io, f"' \" \\\\ \\u{7f} € 🖖 \\u{e0000}")
    @test ts(io) == "' \" \\ \x7f € 🖖 \Ue0000"
end


text:
import ..Embeddings:
    AbstractEmbedding

"""
    minima_and_stepsizes(points, ϵ) -> (Vector{Float}, Vector{Float})

Find the minima along each axis of the embedding, and computes appropriate
`stepsizes` given `ϵ`, which provide instructions on how to grid the space.
Assumes each point is a column vector.

Specifically, the binning procedure is controlled by the type of `ϵ`:

1. `ϵ::Int` divides each axis into `ϵ` intervals of the same size.
2. `ϵ::Float` divides each axis into intervals of size `ϵ`.
3. `ϵ::Vector{Int}` divides the i-th axis into `ϵᵢ` intervals of the same size.
4. `ϵ::Vector{Float64}` divides the i-th axis into intervals of size `ϵᵢ`.
"""
function minima_and_stepsizes(points, ϵ)
    if size(points, 1) > size(points, 2)
        error("The dimension of the dataset exceeds the number of points.")
    end
    D = size(points, 1)
    n_pts = size(points, 2)

    axisminima = Vector{Float64}(undef, D)
    top = Vector{Float64}(undef, D)

    for i = 1:D
        axisminima[i] = minimum(points[i, :])
        top[i] = maximum(points[i, :])
    end
    axisminima = axisminima - (top - axisminima) / 100
    top = top + (top - axisminima) / 100

    stepsizes = Vector{Float64}(undef, D)
    if ϵ isa Float64
        stepsizes = [ϵ for i in 1:D]
    elseif ϵ isa Vector{Float64}
        stepsizes = ϵ
    elseif ϵ isa Int
        stepsizes = (top - axisminima) / ϵ
    elseif ϵ isa Vector{Int}
        stepsizes = (top - axisminima) ./ ϵ
    elseif ϵ isa Tuple{Vector{Tuple{Float64, Float64}}, Int}
        # We have predefined axis minima and axis maxima.
        n_bins = ϵ[2]
        dim = length(ϵ[1])
        stepsizes = zeros(Float64, dim)
        axisminima = zeros(Float64, dim)

        for i = 1:dim
            stepsizes[i] = (maximum(ϵ[1][i]) - minimum(ϵ[1][i]))/n_bins
            axisminima[i] = minimum(ϵ[1][i])
        end
    end

    axisminima, stepsizes
end

"""
    assign_integer_bin_label_to_eachpoint!(
        A::Array{Int, 2},
        embedding_points::Array{Float64, 2},
        axisminima::Vector{Float64},
        stepsizes::Vector{Float64},
        npts::Int) -> Array{Int, 2}

Assign integer bin labels to the points of an embedding, storing the
labels in the preallocated array `A`.
"""
function assign_integer_bin_label_to_eachpoint!(
            A::AbstractArray{Int, 2},
            pts::AbstractArray{T, 2},
            axisminima::Vector{Float64},
            δs::Vector{Float64},
            npts::Int) where T
    @inbounds for i = 1:npts
        A[:, i] .= floor.(Int, abs.(axisminima .- pts[:, i]) ./ δs)
    end
end

"""
    assign_bin_labels(points, ϵ)

Consider a rectangular grid specified by ϵ. Assign bin labels to the provided
`points` by checking which bins each point fall into. Each points is given a
label which is an integer encoding of the origin of the bin it falls into.

The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.

The points are assumed to be provided as an array where each point is a column.
"""
function assign_bin_labels(points, ϵ)
    # Make sure that the array contains points as columns.
    if size(points, 1) > size(points, 2)
        error("The dimension of the dataset exceeds the number of points.")
    end

    # How many points are there in the embedding, and what is its dimension?
    D = size(points, 1)
    npts = size(points, 2)

    if npts <= D
        error("Number of pts is less than dim=$D. Is every column a point?")
    end
    # Find the minima of the embedding and generate step
    # sizes along each axis
    bottom, stepsizes = minima_and_stepsizes(points, ϵ)
    if typeof(stepsizes) <: Union{Float64, Int}
        stepsizes = [stepsizes]
    end

    # Each points of the embedding gets assigned to one bin.
    # The coordinates of each point of the original embedding are
    # assigned an integer number indicating which bin along the respective
    # dimension it falls into.
    visited_bins_inds = zeros(Int, D, npts)
    assign_integer_bin_label_to_eachpoint!(
        visited_bins_inds,
        points,
        bottom,
        stepsizes,
        npts)

    return visited_bins_inds
end

"""
    assign_bin_labels(E::AbstractEmbedding, ϵ) -> Array{Int, 2}

Find bins visited by the reconstructed orbit and assign
to them integer tuples uniquely identifying the bins.

Each point of the embedding gets assigned one integer tuple label,
stored as column vectors in the returned array. If some bin
is visited more than once, there will be repeated columns.

# Details
Given a bin size ϵ, whics is `Union{Int, Float}` for cubic bins,
and `Union{Vector{Int}, Vector{Float64}}` for rectangular bins,
assign a tuple of integers uniquely identifying that bin.

Take, for example, the pair `(pᵢ, (n₁, n₂, n₃)`. Having a
3-element tuple associated with pᵢ means that the reconstruction
is 3D, so we have three coordinate axes - x₁, x₂ and x₃ - to
consider.

This generic example should be read as:
for point `pᵢ`, one must jump n₁ steps along x₁, n₂ steps
along x₂ and n₃ steps along x₃ (in units of ϵ₁, ϵ₂ and ϵ₃)
from the minima of each axis to reach the coordinates of `pᵢ`.

Each visited bin is thus uniquely identified by a tuple of
integers. `visited_bins_indices` assigns one such tuple
to each of the points in the embedding. These are gathered
in an array where column vectors represent the tuples
associated to each point.
"""
function assign_bin_labels(E::AbstractEmbedding, ϵ)
    assign_bin_labels(E.points, ϵ)
end

"""
    assign_coordinate_labels_to_eachpoint!(
            A::Array{Float64, 2},
            pts::Array{Float64, 2},
            axisminima::Vector{Float64},
            δs::Vector{Float64},
            npts::Int)

Assign bin labels to the points of an embedding as coordinate
relative to the minima along each coordinate axis. Store the
coordinate labels in the preallocated array `bin_origins`.
"""
function assign_coordinate_labels_to_eachpoint!(
        bin_origins::AbstractArray{Float64, 2},
        visited_bin_inds::AbstractArray{Int, 2},
        pts::AbstractArray{Float64, 2},
        axisminima::Vector{Float64},
        δs::Vector{Float64},
        n_orbit_pts::Int)
    for i = 1:n_orbit_pts
        bin_origins[:, i] = axisminima .+ (δs .* visited_bin_inds[:, i])
    end
end

"""
    assign_coordinate_labels(points, ϵ) -> Array{Float64, 2}

Consider a rectangular grid specified by ϵ. Assume the bin labels to the provided
`points` by checking which bins each point fall into. Each points is given a
label which is the coordinates of the origin of the bin it falls into.

Together with `ϵ`, the coordinates of the origins of the bins provide
complete information about a coarse graining that covers all visited
states of the system (when `points` represent the orbit).

The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.

The points are assumed to be provided as an array where each point is a column.
"""
function assign_coordinate_labels(points, ϵ)
    # Make sure that the array contains points as columns.
    if size(points, 1) > size(points, 2)
        error("The dimension of the dataset exceeds the number of points.")
    end

    # How many points are there in the embedding, and what is its dimension?
    D = size(points, 1)
    npts = size(points, 2)

    if npts <= D
        error("Number of pts is less than dim=$D. Is every column a point?")
    end

    # Find the minima of the embedding and generate step
    # sizes along each axis.
    axisminima, stepsizes = minima_and_stepsizes(points, ϵ)
    visited_bin_inds = assign_bin_labels(points, ϵ)
    bin_origins = zeros(Float64, D, npts)

    n_orbit_points = size(points, 2)
    assign_coordinate_labels_to_eachpoint!(
        bin_origins,
        assign_bin_labels(points, ϵ),
        points,
        axisminima,
        stepsizes,
        n_orbit_points
    )

    return bin_origins
end

"""
    assign_coordinate_labels(E::AbstractEmbedding, ϵ) -> Array{Float64, 2}

Consider a rectangular grid specified by ϵ. Assume the bin labels to the
provided embedding by checking which bins each point fall into. Each points is
given a label which is the coordinates of the origin of the bin it falls into.

Together with `ϵ`, the coordinates of the origins of the bins provide
complete information about a coarse graining that covers all visited
states in the reconstruction/embedding.

The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.
"""
function assign_coordinate_labels(E::AbstractEmbedding, ϵ)
    assign_coordinate_labels(E.points, ϵ)
end

"""
    assign_coordinate_labels(visited_bin_inds, points, ϵ)  -> Array{Float64, 2}

Consider a rectangular grid specified by ϵ. Assume that, given `ϵ`,
integer bin  labels have been assigned to each point and are stored in
`visited_bin_inds`. This array contains one column of integer bin
labels each of the points in `points`.

Here, we return an array of the same size as `visited_bin_inds`. However,
instead of the labels being integers, they are converted into the
coordinates of the origin of the corresponding bins.

Together with `ϵ`, the coordinates of the origins of the bins provide
complete information about a coarse graining that covers all visited
states of the system (when `points` represent the orbit).

The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.

The points are assumed to be provided as an array where each point is a column.
"""
function assign_coordinate_labels(visited_bin_inds, points, ϵ)
    # Make sure that the array contains points as columns.
    if size(points, 1) > size(points, 2)
        error("The dimension of the dataset exceeds the number of points.")
    end

    # How many points are there in the embedding, and what is its dimension?
    D = size(points, 1)
    npts = size(points, 2)

    if npts <= D
        error("Number of pts is less than dim=$D. Is every column a point?")
    end

    # Find the minima of the embedding and generate step
    # sizes along each axis.
    axisminima, stepsizes = minima_and_stepsizes(points, ϵ)

    bin_origins = zeros(Float64, size(visited_bin_inds))

    n_orbit_points = size(points, 2)
    assign_coordinate_labels_to_eachpoint!(
        bin_origins,
        visited_bin_inds,
        points,
        axisminima,
        stepsizes,
        n_orbit_points
    )

    return bin_origins
end

"""
    assign_coordinate_labels(visited_bin_inds,
        E::AbstractEmbedding, ϵ)  -> Array{Float64, 2}

Consider a rectangular grid specified by ϵ. Assume that, given `ϵ`,
integer bin  labels have been assigned to each point and are stored in
`visited_bin_inds`. This array contains one column of integer bin
labels each of the points of the embedding.

Here, we return an array of the same size as `visited_bin_inds`. However,
instead of the labels being integers, they are converted into the
coordinates of the origin of the corresponding bins.

Together with `ϵ`, the coordinates of the origins of the bins provide
complete information about a coarse graining that covers all visited
states of the reconstructed state space.

The following `ϵ` will work:

* `ϵ::Int` divide each axis into `ϵ` intervals of the same size.
* `ϵ::Float` divide each axis into intervals of size `ϵ`.
* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size.
* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`.
"""
function assign_coordinate_labels(visited_bin_inds, E::AbstractEmbedding, ϵ)
    assign_coordinate_labels(visited_bin_inds, E.points, ϵ)
end

export
minima_and_stepsizes,
assign_integer_bin_label_to_eachpoint!,
assign_coordinate_labels_to_eachpoint!,
assign_bin_labels,
assign_coordinate_labels


text:
# radial plot

using VegaLite

@vlplot(
    data={values=[12, 23, 47, 6, 52, 19]},
    theta={:data, stack=true},
    radius={:data, scale={type=:sqrt, zero=true, range=[20,100]}},
    color={"data:n", legend=nothing},
    view={stroke=nothing}
) +
@vlplot(mark={:arc, innerRadius=20, stroke="#fff"}) +
@vlplot(mark={:text, radiusOffset=10}, text=:data)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
export MaxFlow
module MaxFlow

export Graph, PathWeight, Path, WeightedPath

# solve Kattis problem "https://open.kattis.com/problems/maxflow" using weighted graphs techniques

struct Graph{T}
    n::Int
    node::Vector{Vector{Int}}
    weight::Vector{Vector{T}}
end

struct PathWeight{T}
    weight::T
    length::Int
    PathWeight(w::T, d::Integer) where T = new{T}(w, d)
end

struct Path
    id::Int
    next::Union{Path,Nothing}
    Path(id::Integer) = new(id, nothing)
    Path(id::Integer, p::Path) = new(id, p)
end

struct WeightedPath{T}
    path::Path
    cweight::PathWeight{T}
    WeightedPath(id::Integer, w::T) where T = new{T}(Path(id), PathWeight(w, 0))
    WeightedPath(id::Integer, w::T, wp::WeightedPath{T}) where T = new{T}(Path(id, wp.path), combine_weights(w, wp.cweight))
end

struct Stack{T}
    a::Vector{WeightedPath{T}}
    Stack(wp::WeightedPath{T}) where T = new{T}([wp])
end

function Graph(n, ii, jj, aa::AbstractVector{T}) where T
    m = length(ii)
    @assert m == length(jj) == length(aa)
    node = Vector{Vector{Int}}(undef, n)
    weight = Vector{Vector{T}}(undef, n)
    for i = 1:n
        node[i] = zeros(Int, 0)
        weight[i] = Vector{T}(undef, 0)
    end
    kk = 0
    for k = 1:m
        i, j, a = ii[k], jj[k], aa[k]
        @assert 0 <= i < n
        @assert 0 <= j < n
        nodei = node[i+1]
        weighti = weight[i+1]
        ix = findfirst(isequal(j), nodei)
        if ix === nothing
            push!(nodei, j)
            push!(weighti, a)
        else
            weighti[ix] += a
        end
    end
    Graph(n, node, weight)
end

function main(ioin::IO=stdin, ioout::IO=stdout)
    n, s, t, graph0 = read_data(ioin)
    # print_results(ioout, n, 0, graph0)
    f, graph = run(n, s, t, graph0)
    print_results(ioout, n, f, graph)
end

function run(n, s, t, graph::Graph{T}) where T
    @assert 0 <= s < n
    @assert 0 <= t < n
    graph = deepcopy(graph)
    akkug = Graph(n, Int[], Int[], T[])
    sum = 0
    wp = bestpath(graph, s, t)
    while is_adding(wp, s, t)
        sum += wp.cweight.weight
        add!(akkug, wp)
        sub!(graph, wp)
        wp = bestpath(graph, s, t)
    end
    sum, akkug
end

"""
    bestpath(g::Graph, s, t)

best path in `g` from node `s` to node `t`
"""
function bestpath(graph::Graph{T}, s, t) where T
    w = typemax(T)
    res = WeightedPath(s, w)
    stack = Stack(res)

    while !isempty(stack)
        top = pop!(stack)
        if top.path.id == t
            res = top
            break
        end
        for c in children(top, graph)
            if is_eligible(c)
                insert_sorted!(stack, c)
            end
        end
    end
    res
end

function is_adding(wp::WeightedPath, s, t)
    wp.path.id == t && !iszero(wp.cweight.weight)
end

is_eligible(wp::WeightedPath) = !is_cyclictop(wp.path)

function is_cyclictop(p::Path)
    ( p.next === nothing || p.next.next == nothing ) && return false
    is_cyclictop(p.next, p.id)
end
is_cyclictop(::Nothing, id) = false
is_cyclictop(p::Path, id) = p.id == id || is_cyclictop(p.next, id)

sub!(g::Graph{T}, wp::WeightedPath{T}) where T = _add(g, -1, wp)
add!(g::Graph{T}, wp::WeightedPath{T}) where T = _add(g, 1, wp)
function _add(g::Graph{T}, sig, wp::WeightedPath{T}) where T
    w = wp.cweight.weight * sig
    p = wp.path
    t = p.id
    while p.next !== nothing
        p = p.next
        s = p.id
        _add(g, s, t, w)
        t = s
    end
    nothing
end

function _add(g, s, t, w)
    node, weight = g.node[s+1], g.weight[s+1]
    k = findfirst(isequal(t), node)
    if k === nothing
        if !iszero(w)
            push!(node, t)
            push!(weight, w)
        end
    else
        nw = weight[k] + w
        if iszero(nw)
            deleteat!(node, k)
            deleteat!(weight, k)
        else
            weight[k] = nw
        end
    end
    nothing
end

function Base.pop!(s::Stack)
    popat!(s.a, 1)
end

function children(wp::WeightedPath{T}, g::Graph{T}) where T
    id = wp.path.id
    node = g.node[id+1]
    weight = g.weight[id+1]
    [WeightedPath(node[c], weight[c], wp) for c in 1:length(node)]
end

function insert_sorted!(st::Stack{T}, wp::WeightedPath{T}) where T
    insert_sorted!(st.a, 1, wp)
end
function insert_sorted!(awp::AbstractArray{WeightedPath{T}}, ix::Integer, wp::WeightedPath{T}) where T
    if ix > length(awp)
        insert!(awp, ix, wp)
    else
        top = awp[ix]
        if wp.path.id == top.path.id
            if wp.cweight < top.cweight
                awp[ix] = wp
            end
        elseif wp.cweight < top.cweight
            insert!(awp, ix, wp)
            delete_sorted!(awp, ix+1, wp)
        else
            insert_sorted!(awp, ix+1, wp)
        end
    end
    nothing
end

function delete_sorted!(awp::AbstractArray{WeightedPath{T}}, ix::Integer, wp::WeightedPath{T}) where T
    ix > length(awp) && return
    top = awp[ix]
    if wp.path.id == top.path.id
        deleteat!(awp, ix)
    else
        delete_sorted!(awp, ix+1, wp)
    end
end

function combine_weights(w::T, pw::PathWeight{T}) where T
    PathWeight(min(w, pw.weight), pw.length + 1)
end

function Base.isless(a::PathWeight{T}, b::PathWeight{T}) where T
    a.weight > b.weight || isequal(a.weight, b.weight) && a.length < b.length
end

Base.isempty(s::Stack) = isempty(s.a)
weight(g::Graph) = sum(sum.(g.weight))

# utilities
function read_data(io::IO=stdin)
    line = split(readline(io))
    n, m, s, t = parse.(Int, line)
    ii, jj, aa = Vector{Int}.(undef, (m, m, m))
    for k = 1:m
        i, j, a = parse.(Int, split(readline(io)))
        ii[k] = i
        jj[k] = j
        aa[k] = a
    end
    n, s, t, Graph(n, ii, jj, aa)
end

function print_results(io::IO, n, f, g::Graph)
    n = length(g.node)
    m = sum(length.(g.weight))
    println(io, "$n $f $m")
    for s = 0:n-1
        node = g.node[s+1]
        weight = g.weight[s+1]
        perm = sortperm(node)
        for k = keys(perm)
            t = node[perm[k]]
            w = weight[perm[k]]
            println(io, "$s $t $w")
        end
    end
end

function Base.show(io::IO, p::Path)
    print(io, "Path(", p.id)
    _showrest(io, p)
end
function _showrest(io::IO, p::Path)
    if p.next === nothing
        print(io, ")")
    else
        print(io, " ← ", p.next.id)
        _showrest(io, p.next)
    end
end

Base.show(io::IO, pw::PathWeight) = print(io, "Weight", (pw.weight, pw.length))
Base.show(io::IO, wp::WeightedPath) = print(io, (wp.path, wp.cweight))

end #module


text:
function step_RK(
    particle_info,
    gr,
    particle_weights,
    Ucartesian,
    xp0,
    zp0,
    Δt,
    θThermal,
    rThermal,
    coordinates,
    IntC,
    multiplier,
)

    e2n_p1 = gr.e2n_p1

    # (a) velocity @ particles
    particle_info = velocities2particle(particle_info,
                                        e2n_p1,
                                        particle_weights,
                                        Ucartesian)

    # velocities2particle_cubic!(gr, particle_info, Ucartesian)

    # (b) advect particles
    particle_info = particlesadvection(xp0, zp0, particle_info, Δt, multiplier=multiplier)

    particle_info, particle_weights, found = tsearch_parallel(particle_info,
                                    particle_weights,
                                    θThermal,
                                    rThermal,
                                    coordinates,
                                    gr,
                                    IntC)

    # println(length(particle_info) - sum(found), " lost particles")

    # (c) store particles velocity
    Uxp, Uzp = particle_velocity(particle_info)

    return Uxp, Uzp, particle_info, particle_weights, found
end

function advection_RK2!(Particles, xp0, zp0, Uxp1, Uzp1, Uxp2, Uzp2, Δt)
    cte = Δt/2
    Threads.@threads for i in eachindex(Particles)
        @inbounds Particles[i].CCart.x = @muladd xp0[i] + (Uxp1[i] + Uxp2[i])*cte
        @inbounds Particles[i].CCart.z = @muladd zp0[i] + (Uzp1[i] + Uzp2[i])*cte

        @inbounds Particles[i].CPolar.x = atan(Particles[i].CCart.x, Particles[i].CCart.z)
        @inbounds Particles[i].CPolar.z = sqrt(Particles[i].CCart.x^2 + Particles[i].CCart.z^2)

        # Fix θ
        @inbounds if Particles[i].CPolar.x < 0
            Particles[i].CPolar.x += 2π
        elseif Particles[i].CPolar.x > 2π
            Particles[i].CPolar.x -= 2π
        end

        # Fix radius
        @inbounds if Particles[i].CPolar.z < 1.22
            Particles[i].CPolar.z = 1.22001
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)

        elseif Particles[i].CPolar.z > 2.22
            Particles[i].CPolar.z = 2.21999
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)
        end

    end
end

function final_advection_RK2!(Particles, xp0, zp0, Uxp1, Uzp1, Uxp2, Uzp2, Δt)
    cte = Δt/2
    Threads.@threads for i in eachindex(Particles)
        @inbounds Particles[i].CCart.x = @muladd xp0[i] + (Uxp1[i] + Uxp2[i])*cte
        @inbounds Particles[i].CCart.z = @muladd zp0[i] + (Uzp1[i] + Uzp2[i])*cte

        @inbounds Particles[i].CPolar.x = atan(Particles[i].CCart.x, Particles[i].CCart.z)
        @inbounds Particles[i].CPolar.z = sqrt(Particles[i].CCart.x^2 + Particles[i].CCart.z^2)

        # Fix θ
        @inbounds if Particles[i].CPolar.x < 0
            Particles[i].CPolar.x += 2π
        elseif Particles[i].CPolar.x > 2π
            Particles[i].CPolar.x -= 2π
        end

        # Fix radius
        @inbounds if Particles[i].CPolar.z < 1.22
            Particles[i].CPolar.z = 1.22001
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)

        elseif Particles[i].CPolar.z > 2.22
            Particles[i].CPolar.z = 2.21999
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)
        end

    end
end

function advection_RK2(particle_info::Vector{PINFO}, gr, particle_weights,
    Ucartesian, Δt, θThermal, rThermal, coordinates, IntC, to)

    # Allocate a couple of buffer arrays
    xp0, zp0 = particle_coordinates(particle_info) # initial particle position

    # STEP 1
    Uxp1, Uzp1, particle_info, particle_weights, found = step_RK(
        particle_info,
        gr,
        particle_weights,
        Ucartesian,
        xp0,
        zp0,
        Δt,
        θThermal,
        rThermal,
        coordinates,
        IntC,
        1.0,
    );

    # STEP 2
    # _, _, _, _, particle_info, particle_weights, found = step_RK(
    Uxp2, Uzp2, particle_info, particle_weights, found = step_RK(
        particle_info,
        gr,
        particle_weights,
        Ucartesian,
        xp0,
        zp0,
        Δt,
        θThermal,
        rThermal,
        coordinates,
        IntC,
        1.0,
    );

    # (c) advect particles
    advection_RK2!(particle_info, xp0, zp0, Uxp1, Uzp1, Uxp2, Uzp2, Δt)

    return particle_info, particle_weights, to
end

function advection_RK4(particle_info::Vector{PINFO}, gr, particle_weights,
    Ucartesian, Δt, θThermal, rThermal, IntC, to)

    # Allocate a couple of buffer arrays
    xp0, zp0 = particle_coordinates(particle_info) # initial particle position

    @timeit to "Advection" begin
        # STEP 1
        xp1, zp1, Uxp1, Uzp1, particle_info, particle_weights = step_RK(
            particle_info,
            gr,
            particle_weights,
            Ucartesian,
            xp0,
            zp0,
            Δt,
            θThermal,
            rThermal,
            IntC,
            0.5,
        )

        # STEP 2
        xp2, zp2, Uxp2, Uzp2, particle_info, particle_weights = step_RK(
            particle_info,
            gr,
            particle_weights,
            Ucartesian,
            xp0,
            zp0,
            Δt,
            θThermal,
            rThermal,
            IntC,
            0.5,
        )

        # STEP 3
        xp3, zp3, Uxp3, Uzp3, particle_info, particle_weights = step_RK(
            particle_info,
            gr,
            particle_weights,
            Ucartesian,
            xp0,
            zp0,
            Δt,
            θThermal,
            rThermal,
            IntC,
            1.0,
        )

        # STEP 4
        _, _, Uxp4, Uzp4, particle_info, particle_weights = step_RK(
            particle_info,
            gr,
            particle_weights,
            Ucartesian,
            xp0,
            zp0,
            Δt,
            θThermal,
            rThermal,
            IntC,
           1.0,
        )

        # (c) advect particles
        advection_RK_step4!(particle_info, xp0, zp0, Uxp1, Uzp1, Uxp2, Uzp2, Uxp3, Uzp3, Uxp4, Uzp4, Δt)

    end

    return particle_info, to
end

function advection_RK_step4!(Particles, xp0, zp0, Uxp1, Uzp1, Uxp2, Uzp2, Uxp3, Uzp3, Uxp4, Uzp4, Δt)
    cte = Δt/6
    Threads.@threads for i in eachindex(Particles)
        @inbounds Particles[i].CCart.x = xp0[i] + (Uxp1[i] + 2*(Uxp2[i] + Uxp3[i]) + Uxp4[i])*cte
        @inbounds Particles[i].CCart.z = zp0[i] + (Uzp1[i] + 2*(Uzp2[i] + Uzp3[i]) + Uzp4[i])*cte

        @inbounds Particles[i].CPolar.x = atan(Particles[i].CCart.x, Particles[i].CCart.z)
        @inbounds Particles[i].CPolar.z = sqrt(Particles[i].CCart.x^2 + Particles[i].CCart.z^2)

        # Fix θ
        @inbounds if Particles[i].CPolar.x < 0
            Particles[i].CPolar.x += 2π

        elseif Particles[i].CPolar.x > 2π
            Particles[i].CPolar.x -= 2π
        end

        # Fix radius
        @inbounds if Particles[i].CPolar.z < 1.22
            Particles[i].CPolar.z = 1.22001
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)

        elseif Particles[i].CPolar.z > 2.22
            Particles[i].CPolar.z = 2.21999
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)
        end

    end
end

@inline function velocities2particle(particle_info, e2n_p1, particle_weights, Ucartesian)
    Threads.@threads for i in axes(particle_info,1)
        # -- Reshape velocity arrays
        # Ux_p1, Uz_p1 = getvelocity(Ucartesian, view(e2n_p1,:,particle_info[i].t) )
        # -- Velocity interpolation (directly using linear SF)
        @inbounds particle_info[i].UCart.x = node2particle(getvelocityX(Ucartesian, view(e2n_p1,:,particle_info[i].t)),
                                                 particle_weights[i].barycentric)
        @inbounds particle_info[i].UCart.z = node2particle(getvelocityZ(Ucartesian, view(e2n_p1,:,particle_info[i].t)),
                                                 particle_weights[i].barycentric)
    end
    return particle_info
end

@inline function velocities2particle!(particle_info, e2n_p1, particle_weights, Ucartesian, i)
    # -- Reshape velocity arrays
    Ux_p1, Uz_p1 = getvelocity(Ucartesian, view(e2n_p1,:,particle_info[i].t) )
    # -- Velocity interpolation (directly using linear SF)
    particle_info[i].UCart.x = node2particle(Ux_p1,
                                             particle_weights[i].barycentric)
    particle_info[i].UCart.z = node2particle(Uz_p1,
                                             particle_weights[i].barycentric)
end

function particlesadvection(xp0, zp0, Particles, Δt; multiplier=1)
    # np = length(Particles)
    # xp = Vector{Float64}(undef, np)
    # zp = similar(xp)
    Δtmultiplier = Δt * multiplier
    @fastmath Threads.@threads for i in eachindex(Particles)
        @inbounds Particles[i].CCart.x = @muladd xp0[i] + Particles[i].UCart.x * Δtmultiplier
        @inbounds Particles[i].CCart.z = @muladd zp0[i] + Particles[i].UCart.z * Δtmultiplier
        @inbounds Particles[i].CPolar.x = atan(Particles[i].CCart.x, Particles[i].CCart.z)
        @inbounds Particles[i].CPolar.z =
            @muladd √(Particles[i].CCart.x*Particles[i].CCart.x + Particles[i].CCart.z*Particles[i].CCart.z)

        # Fix θ
        @inbounds if Particles[i].CPolar.x < 0
            Particles[i].CPolar.x += 2π
        elseif Particles[i].CPolar.x > 2π
            Particles[i].CPolar.x -= 2π
        end

        # Fix radius
        @inbounds if Particles[i].CPolar.z < 1.22
            Particles[i].CPolar.z = 1.22001
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)
        elseif Particles[i].CPolar.z > 2.22
            Particles[i].CPolar.z = 2.21999
            Particles[i].CCart.x = polar2x(Particles[i].CPolar)
            Particles[i].CCart.z = polar2z(Particles[i].CPolar)
        end
    end
    return Particles
end

function particlesadvection!(xp, zp, xp0, zp0, particle_info, Δt; multiplier=1)

    Threads.@threads for i in eachindex(particle_info)
        @inbounds particle_info[i].CCart.x = xp[i] = xp0[i] + particle_info[i].UCart.x * Δt * multiplier
        @inbounds particle_info[i].CCart.z = zp[i] = zp0[i] + particle_info[i].UCart.z * Δt * multiplier
        @inbounds particle_info[i].CPolar.x = atan(particle_info[i].CCart.x, particle_info[i].CCart.z)
        @inbounds particle_info[i].CPolar.z = sqrt(particle_info[i].CCart.x^2 + particle_info[i].CCart.z^2)

        # Fix θ
        @inbounds if particle_info[i].CPolar.x < 0
            particle_info[i].CPolar.x += 2π

        elseif particle_info[i].CPolar.x > 2π
            particle_info[i].CPolar.x -= 2π
        end

        # Fix radius
        @inbounds if particle_info[i].CPolar.z < 1.22
            particle_info[i].CPolar.z = 1.22001
            particle_info[i].CCart.x = xp[i] = polar2x(particle_info[i].CPolar)
            particle_info[i].CCart.z = zp[i] = polar2z(particle_info[i].CPolar)

        elseif particle_info[i].CPolar.z > 2.22
            particle_info[i].CPolar.z = 2.21999
            particle_info[i].CCart.x = xp[i]= polar2x(particle_info[i].CPolar)
            particle_info[i].CCart.z = zp[i]= polar2z(particle_info[i].CPolar)
        end

    end

end

function advection!(particle_info,Δt)
    particlesadvection!(particle_info, Δt)  # advect cartesian coordinates
    cartesian2polar!(particle_info)
    return particle_info
end

function particle_velocity(P::Vector{PINFO})
    np = length(P)
    uxp = Vector{Float64}(undef, np)
    uzp = similar(uxp)
    Threads.@threads for i in 1:np
        @inbounds uxp[i] = P[i].UCart.x
        @inbounds uzp[i] = P[i].UCart.z
    end

    return uxp, uzp
end

function particle_velocity(P::PINFO)
    uxp = P.UCart.x
    uzp = P.UCart.z
    return uxp, uzp
end

function particle_velocity!(uxp, uzp, P::Vector{PINFO})
    Threads.@threads for i in eachindex(P)
        @inbounds uxp[i] = P[i].UCart.x
        @inbounds uzp[i] = P[i].UCart.z
    end
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone